**3. Discussion**

Knowledge about a protein–protein interaction interface at the atomic level is essential for the development or design of substances intended to treat diseases. Similarly, the understanding of epitope–paratope interactions in an antibody–antigen complex at the atomic level is crucial for rational development of effective therapeutics.

An important step for the characterization of a functional antibody is accurate delineation of an epitope [31]. Insights about an epitope from a therapeutic perspective can be applied to rational design of mAbs [32]. On the other hand, a detailed epitope analysis is a major restraint because individual amino acid residues constituting an epitope are difficult to pinpoint. Given that experimental methods are time and resource consuming, and success is not guaranteed, computational techniques represent a rapid alternative for validation and characterization of an epitope and affinity maturation of mAbs [17].

In the present study, the epitope reported by the Elson group [22] was characterized and validated by in silico alanine mutagenesis. Epitopic residues in TLR4 were mutated to alanine to investigate the impact of each residue on TLR4 conformation and binding to an antibody. Among the six mutated residues, three, i.e., Y328, N329, and K349, were found to cause a conformational change in the TLR4 structure. RIN analysis suggested that these mutations induce a shift in residue topology owing to a loss of contact with the nearest neighbor or the formation of a new network in the protein structure (Figure S1). Because the epitope of TLR4 is conformational, the loss of the contact with the surrounding residues owing to a mutation will make the epitope ambiguous. The compactness of these muteins, as denoted by *Rg*, showed a significant increase (Figure 2C). The higher *Rg* values sugges<sup>t</sup> that the mutations decreased structural stability. The MD simulation analysis of these mutations revealed the protein's unstable behavior because RMSDs of the muteins uncovered variable behavior during the simulation, and residues from 175 to 400 (Hu 15C1–binding epitope) [26] of the mutated structures oscillated with a higher amplitude, as seen in the RMSF plot (Figure 1C). Gibbs FEL analysis indicated the energetic importance of these mutated residues for the stability of the TLR4 structure. The FEL plot contains splitting of the energy landscape with higher energy barriers, and similar trends were observed for all three mutations. The different conformations assumed by TLR4 owing to the mutations led to conformational changes in the epitopes and resulted in abrogation of the antibody binding, as illustrated by the PCA (Figure S4). The experimental data on the site-directed mutagenesis showed that one of these mutations, i.e., the K349 mutation, alters Hu 15C1 binding [26]. This result may be due to the conformational change in TLR4, which makes the Hu 15C1–binding epitope ambiguous. The impact of the mutation on TLR4 conformation, flexibility, and stability, as evidenced by NMA, showed that the protein folding energy may be increased and stability decreased. The motion of the termini of the muteins, causing twisting and deformity in the protein structure, leads to an increase in internal energy (Figure 2). Steric hindrance or torsion depending upon the direction of the motion induces an increase in the tension on intramolecular bonding, because of which internal hydrogen bonds between the residues break, which results in instability of the TLR4 solenoid structure. From thermodynamic and biophysical points of view, a conformational change in the TLR4 structure is due to the mutation and can make the conformational epitope ambiguous for Hu 15C1 binding. These observed impacts including steric clashes, hydrogen bond network disruption or abrogation of antibody binding, can lead to gain or loss in function of the receptor [33,34]. Moreover, Tsukamoto et al., assessed the functional role of TLR4 activation by a dual-luciferase NF-κB reporter assay, and demonstrated that the overexpression of TLR4w<sup>t</sup> induced NF-κB activation was reduced upon transfection of the mutated plasmid in which LLR13 residues (TLR4 residing epitope) were mutated [27]. From these studies, we can assume that the mutation-driven structure change makes the antibody unable to recognize the epitope, and these residues play a functional role in TLR4 activation as they are essential for dimerization. We believe that the evidence is in support of our findings that mutations are responsible for structural change, which impacts the function of the protein [34,35].

To ascertain whether the mutations affected the binding affinity of Hu 15C1 and overall stability of the Hu 15C1–TLR4 complex, Hu 15C1 was docked with muteins TLR4Y328A, TLR4N329A and TLR4K349A, and with TLR4wt. The backbone RMSD of the simulated mutein mAb-docked complexes uncovered an increasing deviation trend and consequently revealed the unstable nature of these mutein complexes (Figure 3A). Moreover, the bindingenergy data showed a significant decline of binding affinity as compared to Ctl (Table 1). The hampering of mAb binding can be explained by the conformational changes in the TLR4 structure because of these mutations. The findings were next confirmed by PCA and FEL data, which suggested that the decrease in binding affinity may be due to the epitopemasking phenomenon. Several studies have shown that minimal structural changes in an antigen are sufficient to prevent antibody binding [36,37]. The effect of these mutations on binding affinity can influence the function of mAbs. Hu 15C1 did not manifest cross-reactivity in other species because of the ambiguity of the epitope owing to several mutations in TLR4 in these species [26]. Computationally assisted mutagenesis of an antigen or antibody, and a dynamic analysis of these muteins, are beneficial during antibody design for mutated epitopes and for maturation of mAbs. Recently, Takefumi et al. investigated a single-amino-acid substitution to improve antigen–antibody interaction by alanine scanning of interfacial residues. Furthermore, using MD, they found that the enthalpic improvement can be attributed to the stabilization of antigen–antibody interfaces [38]. The

dynamics at atomic-level resolution of antigen–antibody interactions provide accurate insights for antibody development. To further investigate the allosteric binding of mAbs Hu 15C1 and C2E3 to TLR4, computationally assisted re-epitoping was performed. TLR4 is an important protein for the immune response. It may possess multiple binding sites because immune-response proteins have multiple and overlapping binding sites [39]. The Ctl epitope is located in the LLR13 region near the TLR4 dimeric interface [27]. One of the predicted epitopes, EP1, overlaps with Ctl. In our analysis of the paratope–epitope interaction interface for therapeutic and specificity purposes, Hu 15C1 and C2E3 were docked to two predicted epitopes. RMSD data on MD simulation trajectories of these complexes revealed that Hu 15C1 and C2E3 are highly stable and buried at EP1 (Figure 3C–F). In contrast to Ctl, this site provides a maximal surface for the interaction with mAbs Hu 15C1 and C2E3. The interface interaction analysis of Hu 15C1 and C2E3 at sites Ctl, EP1 and EP2 indicated that the molecular association of TLR4-mAb residues ensured the highest-peak formation at a shorter distance at the EP1 site, which means that these contacts are energetically preferred for both mAbs. Moreover, ionic and hydrogen bonds persisted during the simulation at the EP1 site, suggesting that Hu 15C1 and C2E3 are strongly bonded at the EP1 site in contrast to the other epitopic sites (Figure 4). The bonding of contacting residues is energetically favored because the distance between the interacting residues remains <4 Å during the simulation (Figure S9). The C2E3–Ctl interface interaction data showed an extended distance between the interacting residues during the simulation and energetically unfavorable interactions because the distribution peaks were short and formed at a distance of >5 Å, which could be the reason for the increased RMSD values. The binding-energy data support this finding. We observed the highest binding energy of Hu 15C1 at Ctl in contrast to the other epitopic sites; similarly, for C2E3, binding energy was the lowest. The salt bridge formation at the EP1 site for Hu 15C1 and C2E3 was greater, which provided extra stability to the complex. Salt bridge formation with a favorable geometrical orientation gives stability to the overall structure [40]. The binding energy data (Table 1) uncovered higher binding energy of C2E3 in comparison with Hu 15C1. It has been experimentally proved that C2E3 has stronger affinity than Hu 15C1 [41]. In short, the epitope analyses provided a comprehensive view of the TLR4 conformational changes that occur due to the mutations and indicated that they may yield epitope ambiguity. These outcomes will make the development of new mAbs easier and faster and will elucidate epitope structure for designing highly specific mAbs. For designing CDRs and their rationalization, the analyses of antibody–antigen interface interaction are important. Investigating the dynamic effect at the atomic level in the paratope–epitope interaction at three different TLR4 epitopic sites, i.e., Ctl, EP1, and EP2, for Hu 15C1 and C2E3, provided an in-depth understanding of these interactions. We believe that these findings will be helpful for rational design of high-affinity mAbs and for affinity maturation of already known mAbs Hu 15C1 and C2E3. Moreover, this study gives insights into a novel epitope at the TLR4-MD2 interface, indicating that heterodimerization of TLR4 and MD2 may be disrupted and that TLR4 signaling can be inhibited. In 2011, a group reported a MD2 mimicking peptide (MD2-I), which binds to the same site and disrupts the TLR4/MD2 interactions. MD2-I specificity has been evaluated in macrophages and animal models [42]. In addition to this, a few studies have also shown the importance of the druggability of this site [43–45].

mAbs have become an essential therapeutic tool, but their rationalization is a challenging phase and a crucial step in the development of highly specific mAbs. Comparative epitope analysis and antigen–antibody interface analysis via MD simulation offers an atomic-level understanding of the dynamic behavior of the antibody–antigen interactions, which is essential for antibody development. Furthermore, in comparison with other methods, the MMPBSA technique has an edge for the evaluation of affinity maturation and for clarification of the changes in binding affinity. However, TLR4 is a multidomain protein, but its full-length crystal structure is not available; therefore, we used the crystal structure of an isolated external domain in our study. The isolated domains can be less stable than the full-length multidomain protein but are fully functional [46]. In this study, we used

computational methods for epitope verification and affinity maturation of TLR4-targeting antibodies. However, due to the limited facilities, we could not validate our findings experimentally, which can be performed in future studies. Therefore, for efficient mAb design, it is important to couple MD simulations with machine learning and deep learning methods. The structure from conformational space of the MD simulations can be utilized to find new plausible mAb conformations complementing pre-existing ones. We expect that the coupling of these computational methods with experimental methods will not only reduce the cost and time but also offer an efficient approach to the development of rationalized mAbs.

In this study, we demonstrated a computational-driven approach for the epitope verification and affinity maturation of TLR4 antibodies. By using network analysis and implicit MD simulation with FEL, we were able to determine the change in hydrogen bond network and van der Waals interactions in Y328, N329 and K349 among six mutated epitopic residues affecting structural integrity and the energy landscape of the TLR4 and, consequently, the antibody affinity. Further, we predicted the novel epitope located in the MD2 binding site region, which could be explored for new therapeutic antibodies. The uniqueness of our approach is that we successfully employed a wide array of computational-driven techniques to determine the dynamics of mAb-TLR4 interactions at the atomic level, and comparative epitope analysis by computing the binding affinity of the TLR4 antibodies. We believe that this computational-driven approach will accelerate the rational design of therapeutic antibodies. However, experimental validation of the present computational methods remains to be observed until more mAb-antigen are studied and experimentally verified. Moreover, coupling the present pipeline with deep learning and machine learning techniques can further enhance the feasibility of the method. Overall, this work shows that network analysis, MD simulation and MMPBSA techniques can be used as an exploratory methodology for the study of antibody rationalization.

#### **4. Materials and Methods**

#### *4.1. In Silico Structure Reconstruction and Mutation of TLR4*

The 3D coordinates of TLR4 (Protein Data Bank [PDB] ID: 3FXI) and TLR4-specific mAbs Hu 15C1 (PDB ID: 4R7D) and C2E3 (PDB ID: 4R7N) were retrieved from the RCSB PDB database. A monomeric TLR4 structure was first extracted from the TLR4–MD2 complex and optimized in Molecular Operating Environment (MOE) 2019.01. After the removal of water related to crystallization, hydrogen atoms and partial charges were added, incomplete and broken side chains were remodeled and hybridization events were adjusted at default structure preparation parameters in MOE. TLR4 muteins were created by the computational alanine-scanning method. The residues with low energy rotamers were selected and replaced with the amino acid alanine by means of a protein builder tool in the MOE suite. A similar protocol was used to generate TLR4–mAb complexes for the muteins followed by energy minimization to optimize the mutein structures.

#### *4.2. Construction of the Interaction Network*

Network description is widely used to analyze network topology and dynamics of complex systems. The RIN was constructed from a graph-based model of protein structure and topological residue connectivity, where residues are nodes and detected interactions are represented as edges. The RIN was constructed in RINalyzer with the help of the structure viz module implemented in Cytoscape 3.7.2 as per the protocol of Guilaume et al. [47]. The cutoff distance was set to <5 Å for residue interaction. The topological network parameters were calculated as an undirected network using a network analyzer [48].

#### *4.3. Epitope Prediction and CDR Annotations*

Conformational epitopes were predicted with the Epipred tool of the SAbPred web server [49,50]. Based on combined conformational matching of the antibody–antigen structures via geometric fitting and knowledge-based asymmetrical antibody–antigen scoring, epitopes of the antigen are listed rank-wise. The score of an epitope is given by:

$$Equitope\,\mathrm{Score} = \sum d(n)\Pr\left(T\_{ab}, T\_{\mathrm{ag}}\right) \tag{1}$$

where *Tab* and *Tag* are the amino acid types of the antibody and antigen residues, respectively, that belong to node *n*. The stereochemistry and geometry of selected epitopes was verified by Ramachandran plot via MOE's inbuilt utility.

CDRs of antibodies retrieved from the PDB database are annotated according to the numbering scheme of Chothia and Lesk [51,52]. Instead of whole Fab regions, CDRs of the mAbs were considered ligand sites in the docking and simulation.

#### *4.4. Molecular Docking and MD Simulation*

Molecular docking of Hu 15C1 and C2E3 Fabs to TLR4 was performed on the Haddock 2.2 web server on the basis of the information derived from experimental and bioinformatics data [53]. The two mAbs were docked with the two newly predicted epitopes (discussed above) and a previously reported epitope (Ctl). In all three cases, the docked solutions generated by the Haddock server were clustered by their RMSDs. The RMSD measures the average distance of backbone atoms of the superimposed docked solutions, signifying the proximity between solutions. The most populated cluster with a lower RMSD was selected, and solutions with the lowest RMSD within the selected cluster were chosen for the MD analysis. The docking results obtained via Haddock were next validated by means of an antigen–antibody docking package implemented in the MOE suite.

MD simulations capture the dynamic behavior of proteins and other biomolecules in full atomic detail and provide an accurate insight into antigen-antibody interactions at very fine temporal resolution. The MD simulations were performed in GROMACS 2019.3 [54] on Intel E5-2680 and Intel E3-1275 with an Nvidia GeForce GTX 1060 graphics processing unit. TLR4, its muteins, and mAb–TLR4 complexes were solvated using the TIP3P water model in a cubic box; the dimensions of the box boundaries were extended by 10 Å from protein atoms. Na+ and Cl− counter ions were added to neutralize the charge of the simulation system and energy minimization was performed using the CHARMM36 (Chemistry at Harvard macromolecular mechanics) and AMBER99SB-ILDN force-fields and a steepdescent algorithm [55]. V-rescale and Berendsen coupling schemes were employed for temperature and pressure equilibration procedures, respectively. After the temperature and pressure equilibration, MD simulations were carried out for 150 ns for each system in the CHARMM36 and AMBER99SB-ILDN force-fields.

## *4.5. Binding-Free-Energy Calculations*

The MMPBSA approach was used to analyze free-energy interactions between the mAbs and TLR4. The enthalpy of the system was computed via molecular mechanics of MMPBSA, whereas the effect of the polar and nonpolar part of the solvent effect on free energy was determined via the Poisson–Boltzmann equation and calculation of the surface area. The basic equation is

$$
\Delta G\_{bind} = \Delta E\_{MM} + \Delta \Delta G\_{sol} - T\Delta S \tag{2}
$$

where **Δ** *Gbind* is 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 entropy change. In GROMACS, the built-in g\_mmpbsa tool and APBS were called for the MMPBSA calculations, and the last 50 ns of an MD simulation trajectory of each complex with 10-frame intervals was extracted for the energy calculations. For the g\_mmpbsa run, the dielectric constant of the aqueous solvent was set to 80, and the interior dielectric constant was set to 4; the surface tension constant g was set to 0.022 kJ/mol.

#### *4.6. Principal Component Analysis and Free-Energy Landscape*

PCA reduces a multidimensional dataset with many variables to a few 'principal components' that still preserve most of the differences between the data. In MD, PCA helps to elucidate the essential dynamics of the system in a low-dimensional FEL. PCA ultimately provides a view of the atoms (in the MD simulations) that move anisotropically to maximize the variance. To assess conformational flexibility of TLR4 and its muteins, PCA was performed in GROMACS, and the results were analyzed by means of the bio3D R package [56]. Before the PCA, translational and rotational motions were eliminated from the MD trajectory of all systems toward the average geometric center of the protein by leastsquare fit superimposition onto a reference structure [57]. To generate a covariance matrix, configurational space is reconstructed using a simple linear transformation in Cartesian coordinate space. Covariance matrix diagonalization generates an eigenvector set, where its vector components describe the direction of the motion, and the corresponding eigenvalue represents the magnitude of its energetic contribution to the motion. The FEL characterizes the state of a protein by exploring its conformations via an MD sampling technique and finds where Gibbs free energy for a sample is minimal. The gmx\_sham tool implemented in GROMACS was used to predict the FEL through plotting of principal components PC1 and PC2 in the academic version of the Mathematica software (version 11.3; Wolfram Research, Inc., Champaign, IL, USA).

The DynaMut [58] machine-learning–based web server was employed to predict the impact of the mutations on TLR4 structure stability and flexibility. This tool predicts the harmonic motion of the wild-type and mutant TLR4 systems. By means of normal mode analysis (NMA), the nontrivial mode of a molecular motion of TLR4 and its mutants was studied and its vector representation was visualized in Pymol 2.3.

**Supplementary Materials:** Supplementary materials can be found at https://www.mdpi.com/ article/10.3390/ijms22115989/s1.

**Author Contributions:** B.A. and S.C.: Conceptualization, methodology, data curation, and writingoriginal draft preparation. B.A. and M.B.: Investigation and formal analysis. M.-S.K. and S.C.: Supervision. All the authors: writing, reviewing and editing. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Research Foundation of Korea (NRF-2020R1F1A1 071517, 2019M3D1A1078940, and 2019R1A6A1A11051471).

**Data Availability Statement:** The structure models and simulation trajectories are available upon request (sangdunchoi@ajou.ac.kr). All remaining data are contained within the manuscript.

**Acknowledgments:** We gratefully acknowledge Asma Achek for the grea<sup>t</sup> help she provided in manuscript editing, and her valuable suggestions that improved the quality of the paper.

**Conflicts of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
