Unraveling the Tomaralimab Epitope on the Toll-like Receptor 2 via Molecular Dynamics and Deep Learning.

Tomaralimab (OPN-305) is the first humanized immunoglobulin G4 monoclonal antibody against TLR2 and is designed to prevent inflammation that is driven by inappropriate or excessive activation of innate immune pathways. Here, we constructed a homology model of Tomaralimab and its complex with TLR2 at different mapped epitopes and unraveled their behavior at the atomistic level. Furthermore, we predicted a novel epitope (leucine-rich region 9-12) near the lipopeptide-binding site that can be targeted and studied for the utility of therapeutic antibodies. A geometric deep learning algorithm was used to envisage Tomaralimab binding affinity changes upon mutation. There was a significant difference in binding affinity for Tomaralimab following epitope-mutated alanine substitutions of Val266, Pro294, Arg295, Asn319, Pro326, and His372. Using deep learning-based ΔΔG prediction, we computationally contrasted human TLR2-TLR2, TLR2-TLR1, and TLR2-TLR6 dimerization. These results reveal the mechanism that underlies Tomaralimab binding to TLR2 and should help to design structure-based mimics or bispecific antibodies that can be used to inhibit both lipopeptide-binding and TLR2 dimerization.

Introduction

Toll-like receptor (TLR) dysregulation due to pathogens causes inflammation or hyperresponsiveness of the immune system and may affect an adaptive immune response.[1] The nature and magnitude of the adaptive immune response have a great impact owing to the defense against pathogens via a TLR because TLR activation guides T-cell differentiation into CD4+ T helper cells or CD8+ cytotoxic T lymphocytes and promotes dendritic-cell maturation into fully competent antigen-presenting cells.[2] TLR signaling is initiated by the dimerization of intracellular Toll/IL-1 receptor (TIR) domains. Except for TLR3, all TLRs recruit the myeloid differentiation primary response 88 (MyD88) protein to the TLR TIR domains, resulting in downstream signaling that culminates in the production of proinflammatory cytokines. Depending on a ligand, TLR2 interacts with TLR1 (triacylated lipopeptides) or TLR6 (diacylated lipopeptides) to form two distinct heterodimers, TLR2–TLR1 and TLR2–TLR6, which lead to MyD88-dependent activation of NF-κB.[3] The pivotal role in the immune responses that are at risk of dysregulation makes TLRs an attractive therapeutic target.

TLR2 dysregulation has been implicated in numerous diseases, for example, atherosclerosis, arthritis, asthma, sepsis, septic shock, tumor metastasis, and autoimmunity.[4−7] Modulation of the TLR2 signaling pathway and the development of TLR2 antagonists to inhibit cytokine production in inflammatory diseases and autoimmune diseases are therapeutically worthwhile. Nonetheless, TLR2 agonists are among the most effective vaccine adjuvants against human immunodeficiency viruses, hepatitis B virus, and human papillomavirus.[8−10]N-methyl-4-nitro-2-(4-(4-(trifluoromethyl) phenyl)-1H-imidazol-1-yl) aniline (CU-T12-9), a small-molecule agonist, stabilizes the TLR1–TLR2 heterodimer and activates downstream signaling by invoking tumor necrosis factor α (TNF-α), interleukin 10 (IL-10), and inducible nitric oxide synthase (iNOS) through NF-κB signaling.[11] Synthetic small-molecule agonists called diprovocims induce heterodimerization of TLR2 and TLR1 as well as the formation of a TLR2 homodimer in vitro.[12] C29 and a derivative of ortho-vanillin inhibit TLR2 signaling in vitro and in vivo.[13] Staphylococcal superantigen-like protein 3 (SSL3) interferes with the TLR2 ligand binding by blocking its binding site and interacts with the complex of TLR2 with lipopeptides, thus preventing the formation of heterodimers TLR2–TLR1 and TLR2–TLR6.[14] On the other hand, neither TLR2 agonists nor antagonists are in clinical or preclinical development as small-molecule drugs because of the high selectivity of the ligands and the stabilization of the protein–protein interaction (PPI).

Nevertheless, the design of monoclonal antibodies (mAbs) competing with PPI formation (necessitating humanization of the immunoglobulins to prevent an immune reaction) has been remarkably successful. Accordingly, mAbs have become a key class of therapeutic agents for the treatment of many human disorders, particularly cancers and immunological, infectious, neurological, and metabolic diseases, because of the active development of antibody medicines in recent decades.[15] The unique structure of antibodies makes them highly antigen-specific and popular biological tools for the precise probing of individual molecules. Their therapeutic application is also highly desirable owing to their high specificity to disease-associated molecules and good safety. More recently, next-generation mAb immunotherapies have boosted the creation of antibody therapies. Computational techniques for antibody discovery have the potential to advance this discipline by generating findings faster than current standard laborious experimental procedures.[16] By using in silico methods rather than in vivo maturation or experimental selection procedures, one may search a much larger space, perhaps discovering bigger and more beneficial evolutionary steps. Iterative computational methods can improve the binding affinity of mAbs beyond that achieved by in vivo maturation.[17] For rational antibody design, well-established structural bioinformatic approaches such as homology modeling,[18,19] protein–protein docking,[20] and protein interface prediction[21] are already in use.

In the present study, we first designed an antibody model and studied its dynamics in detail. Next, we used a TLR2 crystal structure and a Tomaralimab-modeled structure to map the likely TLR2 epitope. As there is no experimentally resolved structure of the TLR2–Tomaralimab complex, we constructed it computationally via flexible protein–protein coupling. Molecular dynamics (MD) simulation trajectories and binding-energy calculation revealed that it is the PE1 site that favors Tomaralimab binding. To evaluate the effect of mutation on Tomaralimab binding affinity, we compared the ΔΔG of the wide type and mutant complex by means of a deep learning algorithm. Furthermore, we evaluated molecular interactions within complexes TLR2–TLR1, TLR2–TLR6, and TLR2–TLR2. We provide a comparative analysis of the most important residues from the dimeric interface of TLR2 that drive the homodimerization and heterodimerization (with TLR1 or TLR6). The deep learning-based ΔΔG prediction enabled us to determine that six dimeric interface residues are responsible for dimerization.

Materials and Methods

Tomaralimab Sequence Identification and Model Building

The sequence was retrieved from the Therapeutic Antibody Database (TABS) a unique database of therapeutic antibodies, and the AntiBodies Chemically Defined (ABCD) database.[22] After verifying the antibody sequence, a three-dimensional (3D) model was built in the MOE2020.09[23] built-in antibody homology modeling application and the PyRosetta program.[24,25] Antibody sequences were separated into a light chain variable domain (VL) and a heavy chain variable domain (VH). After that, for homology modeling, a suitable framework template was selected based on scoring and a structural fit of the integrity of the backbone. After determining the VL–VH framework of Tomaralimab, complementarity-determining regions (CDRs) in loops were assigned through grafting onto the Fv framework. Finally, we generated 100 single-chain variable fragment (scFv) models for the most likely orientations of the CDRs to preserve their distinct loop conformations. The CDRs of the antibody were annotated according to Chothia and Lesk’s numbering scheme.[26,27] The top model based on the lowest energy was chosen for further analysis.

Epitope Prediction

Epitopes were predicted as described previously.[28] The predicted epitopes were ranked by their score, which is given bywhere, Tab and Tag are the amino acid types of antibody and antigen residues, respectively, that belong to the node n.

Molecular Docking and MD Simulation

The Tomaralimab 3D model was docked to TLR2 (Protein Data Bank ID: 2xi) via the reported and predicted epitopes by an antibody–antigen docking protocol described in our previous report.[28] After the clustering of the docked solutions, three solutions were selected from the most populous clusters on the basis of root mean square deviation (RMSD) and the docking score in each case.

The stability of the docked complex was validated by MD simulations performed in GROMACS 2020.2[29] as described in our previous study.[28]

Binding Energy Calculations

The binding energy between the mAb and TLR2 was calculated using molecular mechanics Poisson–Boltzmann surface area (MMPBSA). The effects of polar and nonpolar parts of the solvent on the free energy were determined by means of the Poisson–Boltzmann equation and the calculation of the surface area, while the enthalpy of the system was computed from the MMPBSA. The basic equations arewhere ΔGbind is the binding free energy, ΔEMM represents the intramolecular energy difference in a vacuum, ΔΔGsol is the solvation energy difference, T denotes absolute temperature, and ΔS is the change in entropy. The MMPBSA calculations were performed using the g_mmpbsa tool and the adaptive Poisson–Boltzmann Solver. Frames from the last one-third of the MD simulation were extracted for each complex at 10-frame intervals. For the g_mmpbsa run, the dielectric constant of the aqueous solvent was set to 80, the interior dielectric constant was set to 4, and the surface tension constant g was set to 0.022 kJ/mol. The average contribution of the residues to the binding energy was calculated for each complex. The entropy contribution was ignored because the cost of computations for these large protein–antibody systems was too high. Additionally, we performed a per-residue decomposition analysis to elucidate the individual energy contribution of mAb and TLR2 amino acids to the overall binding energy. A binding free energy decomposition was carried out using the gmx_MMPBSA tool.[30] We calculated the energetically important residues within 5 Å at the interface by decomposing each residue using the effective free energy decomposition method.[31]

Deep Learning Framework for Predicting Binding Affinity Changes Upon Mutations

For the prediction of the binding affinity change upon mutation in TLR2-Tomaralimab, the single-point mutations in the TLR2 epitope were enumerated by means of PyRosetta implemented in the jupyter notebook. To estimate the mutation(s) effect, the resulted complexes’ sidechains were repacked around the mutation, and energy was minimized. The geometric deep learning method[32] was utilized to predict and identify the effect of mutation on Tomaralimab binding to TLR2. The geometrical neural network encodes the residues in wide type and the mutant complex with local coordinates as input. Let fwt represent the feature of the i th residue in the wild-type complex, and fmut represent the feature of its counterpart in the mutant complex. The features shared across residues in the complex are used as an input in the multilayer perceptron to predict the difference in binding affinity between the two complexes

Where MLP1 stands for standard multilayer perceptron network, while W stands for trainable weight matrix.

Furthermore, TLR2 homo- and heterodimers were constructed in PyRosetta. An alanine scan was performed using PyRosetta scripts, where computational models of the alanine variants were first generated, with energy minimization. We performed binding-energy calculations using a deep learning approach[32] to determine ΔΔG of alanine mutants.

Results and Discussion

The anti-TLR2 humanized immunoglobulin (Ig) G4 mAb in question has potential anti-inflammatory and antineoplastic activities.[33] Upon intravenous administration, Tomaralimab binds to the ligand-binding site on the receptor (TLR2) and blocks the activation of TLR2-mediated innate-immunity signaling.[34]

Characterization of the Structure and Dynamics of Tomaralimab

Tomaralimab is an anti-TLR2 humanized immunoglobulin G4 mAb created by Opsona Therapeutics. Tomaralimab is being studied in two patient groups: a Phase I/II trial against lower-risk myelodysplastic syndrome (NCT02363491) and a Phase II trial testing the ability to prevent delayed renal graft function (NCT01794663), a condition that can occur after a kidney transplant.[35]

Antibody structure prediction has been widely used in many biological analyses. The 3D structure of Tomaralimab has not yet been reported. The 3D homological model of the Tomaralimab scFv region was constructed here in the MOE2020 antibody modeler and PyRosettaAb from an IgG template (Protein Data Bank ID: 2NY7) that consists of one VL and one VH domain (Figure a). The quality of the top 3D model was evaluated on the MolProbity web server (http://molprobity.biochem.duke.edu/),[36] which determines the clash score (the number of unfavorable steric overlaps at >0.4 Å per 1,000 atoms), the percentage of backbone conformations in the favored Ramachandran region, and the MolProbity score that combines the clash score, the percentage of side-chain conformations classified as rotamer outliers, and the percentage of backbone Ramachandran conformations outside the favored region (Table ). The overall steric hindrance of the structure owing to clashes was removed through the refining of the model by maintaining the appropriate orientation of the CDRs.

Figure 1

Annotation and dynamics of the Tomaralimab homology model. (a) Framework regions are highlighted in green (VL) and marine blue (VH), and hypervariable loops are red (H1 and L1), orange (H2 and L2), and indigo (H3, L3). Root mean square (b) deviation and (c) fluctuation of VL and VH. (d) Compactness of the mAb.

Table 1

Validation of Tomaralimab Models

		Ramachandran data
model no.	clash score	outliers (%)	favored (%)	MolProbity score
1	5.05	1.40	92.09	2.55
2	6.12	2.79	89.07	2.62
3	6.74	3.49	83.72	2.75
4	5.51	3.02	86.74	2.69
5	6.43	3.26	85.81	2.61
6	8.57	1.86	85.35	2.79
7	7.66	3.26	85.35	2.76
8	6.28	3.02	86.98	2.57

To assess the flexibility, mobility, and accuracy of the modeled structure, the PyRosetta 3D model with the least outliers on the Ramachandran plot was subjected to MD simulations. By simply examining the simulation of the mAb, we can quantify the extent to which various regions of the molecule move at equilibrium and can determine which types of structural dynamics it undergoes. The stability of the simulation system was evaluated by determining the RMSD of Cα atoms of the 3D model (Figure a); initially, it was assumed that the system achieves stability after 50 ns, but the graph of the 3D model showed fluctuation between nanoseconds 105 and 125. Finally, the conformational stability of the mAb was reached after 125 ns because the RMSD remained constant thereafter. On the other hand, both the VL and VH showed the same behavior (Figure b). The fluctuation in RMSD from nanoseconds 105 to 125 may be due to the flexible nature of the loops. To obtain sufficient sampling of the mAb conformation sampling, we increased the sampling time to 450 ns (Figure b). The RMSD remained constant for further sampling time. Further, to check the reproducibility and reliability of the data, we performed the 3 replicas of mAb. The RMSD showed a similar trend in 3 replicas (Figure S2a). The displacement of individual atoms at an instant of the simulation was measured by determining root mean square fluctuation (RMSF). An RMSF graph (Figure c) revealed marked oscillation of loop residues 130–140 and 187–200 in the VH chain. Compared to other regions in both VL and VH, the loop residues oscillated with a higher amplitude. Furthermore, we assessed the mAb’s compactness by determining its radius of gyration (Rg). This parameter (Figure d) underwent a significant decrease between nanoseconds 50 and 80 and then stayed constant during the remainder of the simulation run. The second and third replicas showed slightly higher Rg values than the first up to 250 ns and remained constant thereafter (Figure S2b). Rg spanned the range from ∼24 to ∼25 Å, revealing that the mAb retained compactness during the simulation.

Mapping of the Tomaralimab-Binding Site on TLR2

To construct the TLR2–mAb complex, the epitope of the mAb should be known. A probable epitope of the anti-TLR2 mAb in question was predicted using Tomaralimab and TLR2 structure coordinate files, as demonstrated elsewhere.[29] Among the ranked epitopes, two were selected for docking analyses and compared with the reported epitope[36] represented by R hereafter in this study. A Ramachandran plot was analyzed to validate the epitopes (Figure S1). The mapping of these epitopes onto TLR2 is shown in Figure . PE1 overlaps with the reported epitope (R) because many residues are common between them and because both are at the TLR2 dimerization interface. These data support the published evidence that Tomaralimab interferes with TLR2 dimerization.[37] PE2 is located on the convex surface of TLR2 in the leucine-rich region at positions 9–12 (LRR9-12; Figure a).

Figure 2

Dynamics of the mapped epitopes. (a) Mapping of the epitopic residues of reported epitope R and of predicted epitopes PE1 and PE2. Tomaralimab docked in (b) R, (c) PE1, or (d) PE2. Structural fluctuations were measured as (e) RMSD and (f) root mean square fluctuation (RMSF).

Understanding the mode of binding of an immunoglobulin to its antigen has immense medical, industrial, and biological implications. To illustrate the mode of binding of the mAb at hand, molecular coupling of Fab (Tomaralimab) to TLR2 was performed in three different epitopes. The docked solutions were clustered, and representative complexes from the most populated clusters were chosen. The preferred solution was superimposed with the structures of TLR1 and TLR6 bound to TLR2 in the same orientation as seen in the respective complex (Figures b–d and S3a–f). The superimposition in the case of the reported epitope and PE1 indicated that OPN-305 disrupts dimerization because it binds to TLR2 in the same region as TLR1 and TLR6 do (Figures b,c and S3a,b and S3d,e). As revealed by the superimposition of the TLR2 dimer and the TLR–mAb complexes, the orientation of these two complexes favors blockage of the dimerization site. Therefore, Tomaralimab should block the heterodimerization of TLR2 with TLR1 and TLR6. According to surface plasma resonance analysis, TLR2 and the immunostimulatory lipopeptide had a direct and specific interaction that was blocked in a dose-dependent manner by OPN-301 (a murine analogue of Tomaralimab).[37] They used OPN-301 to treat HEK293 cells that had been overexpressed with mutant human TLR2 construct that lacked the respective portion of the wild-type extracellular domain.[38] The absence of NF-κB-dependent reporter gene activation in response to lipopeptide exposure after OPN-301 administration demonstrated that the epitope identified by OPN-301 is located within the LLR9-12 region of TLR2.[37] The entrance to the lipopeptide binding pocket in TLR2 is located between LRRs 11 and 12,[39] which is crucial to both the TLR–TLR (TLR2–TLR1 or TLR2–TLR6) dimerization surface and the antibody epitope (Figure b,c). The multiple sequence alignment of human, mouse, and monkey TLR2 (Figure S4) showed 72 and 96% sequence identity, respectively. Typically, the amino acids of TLR2 which are interacting with the antibody are highly conserved. However, OPN-301 specifically inhibits mammalian TLR2 activation and cross-reacts with human, pig, and monkey TLR2, indicating that this antibody is specific for a critical epitope.[40] Furthermore, in vitro studies show that Tomaralimab inhibits TLR2 signaling in mice, pigs,[41] cynomolgus monkeys, and human cells.[34] The data suggest that Tomaralimab is specific to a critical epitope yet conserved. Consequently, these data implicate that Tomaralimab binds to a specific epitope within the TLR2 extracellular domain to abrogate the TLR2 heterodimerization, resulting in the silencing of downstream signaling cascades. Furthermore, in the third complex, the mAb binds to the convex face of LRR10-13, partially covering the entrance of the lipopeptide-binding pocket (Figures d and S3c,f). We superposed this complex with the crystal structure of the TLR2–SSL3 complex (Figure S3g). We found that the binding of the mAb overlaps with the binding of SSL3. These data suggested that the antibody binding at the PE2 site accompanies a conformational change in TLR2 and prevents dimerization, which is crucial for the activation of downstream signaling. These results implied that the mAb designed based on the PE2 epitope interferes with lipopeptide binding and disrupts an already-formed TLR2–lipopeptide complex, thereby, preventing TLR heterodimerization and downstream signaling. These findings are backed up by a study on SSL3 binding to TLR2 on the convex face.[14]

Dynamics of Tomaralimab at the Mapped Sites on TLR2

To determine the correct mapping for locating more precisely the true binding site of Tomaralimab, we performed an MD simulation of the docked complexes. Structural fluctuation of the simulation systems was measured by means of RMSD and RMSF of Cα atoms. The RMSD of the PE1 and PE2 complexes proved to be slightly higher than that of the R complex up to 110 ns (Figure e). Meanwhile, to achieve sufficient sampling of conformations, we enhanced the sampling time up to 300 ns of each complex. There was an increase in RMSD of the R complex after 150 ns but remained constant (Figure e). However, the RMSD of the PE1 and PE2 complexes did not change. To make the data conclusive, we performed the 3 replicas of each complex (Figure S2c–e). All complexes showed similar behaviors in 3 replicas; however, R and PE1 underwent an increase in RMSD in the third replica with a similar trend as in their first and second run. To assess the protein–mAb interactions during the simulations, a number of hydrogen bonds were assessed. As the time trajectory progressed, the mAb formed 10 to 15 hydrogen bonds at the PE1 site, 7 to 12 at R and 10 to 12 at the PE2 site (Figure S2f). In addition to the RMSD of the complexes, the mAb and TLR2 in the complex were compared with their apo forms. The bound forms of both showed higher RMSDs than their apo forms did (Figure S3h,i). We determined the displacement of the residues—or the extent to which the residues of the mAb and TLR2 fluctuated during the simulation during the complex formation—by measuring their RMSF. The TLR2 residues in the PE2 complex were found to oscillate with a higher amplitude, especially the residues in the central region, where the antibody binds in this complex. On the contrary, in the other two complexes, only the terminal residues showed fluctuation, which was also seen in the apo form (Figure f).

Electrostatic Interactions within the TLR2–mAb Complex

To understand the structural origin of TLR2–mAb binding and affinity differences between different mapping sites, we focused on amino acid interactions of the binding contact regions (Figure ). The dimerization interface of TLR2–TLR1 has a small hydrophobic core in the center surrounded by ionic and hydrogen-bonding interactions.[39,42] In the case of the first and second complexes, the antibody binds to the dimeric region by forming hydrophobic, hydrogen, and ionic bonds with TLR2. Tomaralimab has several interlocking hydrogen bonds in R (Figure a,b), including Glu29–His292 (H), Ser67–Asp205 (H), Ser80–Asp159 (H), and Arg18–Met133 (H), while additional hydrophobic contacts with the antibody can be found in the Val58–Thr262-type ionic bond; the average distance is 2 Å.

Figure 3

Fluctuation of interatomic distances during interfacial interactions. (a) Computed averages and fluctuations of interatomic distances for selected residues. Three epitopic patches (b) R (c) PE1 and (d) PE2 are shown with structural details of key interfacial interactions between TLR2 and the antibody.

In the PE1 patch, both light and heavy chains of the antibody are involved, in sharp contrast to the TLR1 or TLR6 binding (Figure S3b,e). The hydrophobic core at the dimeric interface is overhauled by CDRs at this location along with Tyr32 and Tyr31, which form hydrogen bonds with Asp393 and Ser418 of the receptor (Figure c). There is also the consistency of such hydrogen bonds as Glu59–Cys321 (H), Asn57–Asp260 (H), and Thr103–Gln370 (H), with an average distance of >5 Å (Figure a,c). An ionic pair was found between Glu29 and His292 with an average distance of 2.5 Å.

Figure 4

Analysis of the binding affinity of Tomaralimab using deep learning. (a) Effect of epitope mutation on the Tomaralimab binding affinity measured by ΔΔG using deep learning. Tomaralimab binding affinity affected by TLR2 residues (b–h) is visualized.

The contact area between the antibody and TLR2 at the PE2 site is smaller (Figures d and S3c,f). The convex nature of TLR2 makes this patch less exposed to antibody binding because we see that fewer areas encounter antibodies. Furthermore, the binding of the antibody to the PE2 patch (LRR9-12) (Figure S3c) of the receptor does not interfere directly with the dimerization site. The hydrogen and ionic interactions of TLR2 and the mAb are presented in Figure d. Hydrogen bonds Tyr32–Le324, Ser60–Arg270, Arg98–Ser272, and Thr28–Asp279 turned out to be transient, with an average distance of >5 Å (Figure a,d). Notably, antibody binding at this site overlaps with SSL3 binding (Figure S3g). The structure of the TLR2–SSL3 complex indicates that SSL3 binds to LRR11-13 on the convex surface of TLR2.[14] SSL3 covers the entrance intended for lipopeptides, and therefore, it should be covered by antibody binding. Finding specific residue–residue interactions between an antibody and antigen and integrating them into a machine learning approach is used to design biologically active antibodies.[43]

Energetic Analysis of the TLR2–mAb Complex

Free energy describes the strength of antibody affinity for an antigen. The paratope–epitope binding reaction results in a free-energy change. The binding-free-energy change accompanying the interaction of the antibody at different epitopes of TLR2 was computed by the MMPBSA technique. Table summarizes the binding free energies of Tomaralimab in three epitopes. The ΔGbinding for the reported epitope proved to be −1158.740 ± 197.572 kJ/mol; notably, for the PE1 epitope, this parameter was found to be −1324.291 ± 191.206 kJ/mol, that is, 166 kJ/mol higher than that of R. In contrast, in PE2, ΔGbinding is 1080.718 ± 197.593 kJ/mol. Regarding the binding-energy data, we concluded that the antibody has a stronger affinity for the PE1 epitope than other epitopes.

Table 2

Binding Energy of Tomaralimab Toward TLR2

epitope position	Vdw energy	electrostatic energy	polar solvation	SASA	binding energy	P value
R	–325.962 ± 29.409	–2006.454 ± 185.514	1218.181 ± 172.642	–44.506 ± 8.320	–1158.740 ± 197.572	0.0001
PE1	–405.136 ± 33.046	–2000.633 ± 249.873	1135.067 ± 213.475	–53.589 ± 8.403	–1324.291 ± 191.206	0.0001
PE2	–452.370 ± 49.961	–1557.42 ± 193.853	988.298 ± 181.932	–59.217 ± 9.354	–1080.718 ± 197.593	0.0003

The lower free energy of the antibody in the PE2 patch means that Tomaralimab binds specifically to the dimerization site. As shown in Figures d and S3c, only a small area comes into contact during the formation of the antigen–antibody complex, and this area is also the loop region of the receptor, which is highly flexible. On the contrary, we observed notable binding energy for the antibody in the PE1 patch, indicating its stronger affinity for this epitope. The antibody was geometrically aligned with TLR2 in a stable configuration, and the area of contact was larger than that of other epitopic patches, as displayed in Figures c and S3b. MMPBSA analysis calculates the approximate free energy of binding within an antibody–antigen complex.[44] Furthermore, we performed per-residue decomposition analysis to elucidate the individual energy contributions of amino acids of mAb and TLR2 toward the overall binding energy. Figure S5 illustrates the energy contribution of the top ten residues in Tomaralimab and TLR2 toward total binding energy at R, PE1, and PE2. His292, Tyr297, and Met133 of TLR2 contribute −3.5 to −2.5 kcal/mol, and Arg18, Ser80, and Ser67 of the mAb contribute −4.6 to −3.6 kcal/mol toward total energy at R (Figure S5a). A total of −4.8 to −3.2 kcal/mol is added to the total energy at the PE1 epitope by amino acids Leu345, Tyr297, Val347, and His372 of TLR2, and −8.7 to −4.4 kcal/mol is added by Asn57, Thr31, Tyr32, and Lys32 of mAb (Figure S5b). In terms of energy content, Ser272, Val277, Tyr306, and Arg270 of TLR2 contribute −3.6 to −3.0 kcal/mol, whereas Tyr32, Thr34, and Arg54 of mAb contribute −4.9 to −4.0 kcal/mol to PE2 (Figure S5c).

Analysis of the Binding Affinity of Tomaralimab Using the Deep Learning Framework

To estimate the effect of mutation on Tomaralimab, we constructed the single-point alanine substitution in the Tomaralimab binding epitope. We computed the difference in binding affinity between the wide type and the mutant complex by measuring ΔΔG using the deep learning framework. The deep learning methods have been validated on split-by-complex fivefold cross-validation over the Structural Kinetic and Energetic database of the Mutant Protein Interactions (SKEMPI) V2.0 dataset. A subset of 1,131 single-point mutations and 1707 multipoint mutations was used to benchmark the model and other baselines. Pearson’s correlations between the predicted ΔΔG and the real ΔΔG values are 0.65. The model makes predictions with moderate to high correlation with experimental binding data.[32] The ΔΔG of the alanine mutants are shown in Figure .

Among epitope-mutated residues, the binding affinity of Tomaralimab was affected by Val266Ala, Pro294Ala, Arg295Ala, Asn319Ala, Pro326Ala, and His372Ala. Alanine substitution of Arg295, Asn319, and His372 showed a significant ΔΔG value of −0.6 to −4.32 kcal/mol. Figure b–g shows the structure of these mutations. Mutations in these residues showed negative free-energy (−ΔΔG) changes, which means that these mutations can reduce the efficacy of Tomaralimab. A potential CDR mutation can be evaluated to improve the efficacy of the antibody. An understanding of protein–protein binding affinity values is vital to understanding biological phenomena in a cell, such as how missense mutations alter the protein–protein binding. With the deep learning algorithm, the changes in binding affinity upon changes in amino acids can be modeled quickly and accurately.[45]

Energy Contribution of TLR2 Dimeric-Interface Residues

To assess the importance of the dimeric-interface residues of TLR2, we analyzed the dimeric interfaces within complexes TLR2–TLR1, TLR2–TLR6, and TLR2–TLR2. The surface corresponding to the PPI within these complexes contains a hydrophobic interaction reinforced by surrounding hydrophilic residues that form ionic and hydrogen bonds (Figure ). To find out the energy contribution of the dimeric-interface residues, we computationally mutated each to alanine, and the resultant complexes were energy-minimized and ΔΔG was determined by using the deep learning approach. To investigate how essential the native residues are, we substituted them with alanine and tested whether this alteration significantly affects the dimerization. The results of the alanine scan (Figure d) revealed that 63% (12 out of 18) of dimeric-interface residues of TLR2 are important for dimerization. We noticed that residues Lys347, Phe349, Leu371, Glu375, Tyr376, and Asn379 make a major contribution to the interactions at the dimerization interface. Each residue contributes −2.5 to −25 kcal/mol energy to the interaction. This is because they help establish strong bonding during TLR2 heterodimerization with TLR6 or TLR1 as well as in homodimerization. Glu375 engages in a strong ionic bond with Lys313 of TLR6. Lys347 and Tyr376 are involved in hydrogen bonding with Thr366 and a hydrophobic contact with TLR6 Pro342. These findings support crystallographic analysis data, which revealed that these residues participate in hydrogen bonding, hydrophobic contacts, and ionic interactions during TLR2–TLR6 and TLR1–TLR2 heterodimerization.[39,42]

Figure 5

Energy contribution of interfacial interactions of TLR2 during homo- and heterodimerization. Interfacial interactions within (a) TLR2–TLR2, (b) TLR2–TLR1, and (c) TLR2–TLR6. (d) Bars represent the energy contribution of TLR2’s interfacial residues to homodimerization and heterodimerization after alanine substitution. TLR2 appears to be highly optimized for heterodimerization, as evidenced by the binding-energy values after a single substitution.

Aberrant TLR2 activation is associated with atherosclerosis, arthritis, asthma, sepsis, septic shock, tumor metastasis, and autoimmunity. Unraveling the biophysical characterization of the Tomaralimab interaction with TLR2 points to a binding mechanism of mAb and provides new insights into the development of TLR2-based therapeutics. Our mapping data on the novel epitope offer a new approach to TLR2 inhibition via blockage of both ligand binding and dimerization by a bispecific antibody to disrupt heterodimerization, which is crucial for TLR activation and downstream signaling. A deep learning framework can identify and predict changes to the CDR that enhance the antibody’s efficacy by estimating the effect of mutations on Tomaralimab binding affinity. In a recent study, the CDR of the P36-5D2 antibody was optimized, and its potency increased ∼10 to 600-fold against SARS-CoV-2 variants, including Delta.[32] Altogether, deep learning and MD simulation methods can efficiently optimize the antibody and potentially develop a new antibody candidate with broader and more potent neutralization.