**2. Results**

#### *2.1. Computational TLR4 Epitope Mutagenesis and Network Analysis*

Hu 15C1 is a potent anti-TLR4 mAb, which, owing to its efficacy, has reached clinical trials and is being evaluated in RA patients [25]. The epitope of mAb Hu 15C1 has been mapped by alanine mutagenesis and has been found to be located in the ectodomain of TLR4 near the dimerization interface [27]. This epitope has also been investigated in another study [26], in which a new derivative antibody, C2E3, was designed. By exploring the TLR4 structure, researchers hypothesize that the epitope-constituting residues suggested by these studies are crucial for preserving the structural integrity of the TLR4 backbone (Figure 1A). With the aim of gaining insights into TLR4 structure and mAb binding at the atomic level of specific residues, computational alanine scanning, and an MD simulation were performed here. Six hot spot residues, i.e., Y328, N329, K349, K351, E369 and D371, were individually mutated to alanine and vetted for changes in the interactions between TLR4 and a mAb.

Moreover, to comprehend the multiple interaction types involved in the residue network in wild-type TLR4 (TLR4wt) and its muteins, network analysis was performed. Substantial changes in the hydrogen bond network and in van der Waals interactions were observed at the mutation sites, and these changes led to a distinct conformational change in the epitope (Figure S1). Topological parameters of the residue interaction network (RIN) are given in Figure S1H. In the RIN of TLR4w<sup>t</sup> and its muteins, there was a considerable variation in the node degree distribution. Thus, mutating these residues could distort the structure of TLR4 and make the precise Hu 15C1 epitope ambiguous.

#### *2.2. MD of TLR4w<sup>t</sup> and Muteins*

Behaviors of TLR4 and its muteins were evaluated through conventional MD. Root mean square deviations (RMSDs) of Cα atoms of muteins and TLR4w<sup>t</sup> were calculated as a function of time. Comparative RMSD plots suggested that TLR4w<sup>t</sup> remained stable, while all the muteins underwent rigorous fluctuations (Figure 1B). The muteins manifested RMSDs in the following order: TLR4Y328A > TLR4N329A > TLR4D371A > TLR4K349A > TLR4E369A > TLR4K351A > TLR4wt. Additionally, to probe the effect of these mutations on TLR4 structure, root mean square fluctuation (RMSF) analysis was carried out, which confirmed that muteins TLR4Y328A, TLR4N329A, and TLR4K349A underwent more fluctuations as compared to the other muteins and TLR4w<sup>t</sup> structures (Figure 1D). Detailed analyses then suggested (discussed below) that these mutations affect the overall structure of TLR4, particularly the proposed epitope region. Owing to these mutation-driven substantial deviations in structural integrity, we concluded that these residues are more important for preserving TLR4 folding. Furthermore, the compactness of TLR4w<sup>t</sup> and its muteins was determined by measuring the radius of gyration (*Rg*). Differences in *Rg* between TLR4w<sup>t</sup> and muteins are presented in Figure 1C. TLR4Y328A, TLR4N329A and TLR4D371A showed significant increases in *Rg*, indicating a decrease in the compactness of the system. In the second run of MD simulation the RMSD and *Rg* (Figure S2) showed almost a similar trend except TLR4K349A which underwent a bit higher RMSD fluctuation during 125–150 ns time-point (Figure S2D). To investigate the robustness of these results, an amber AMBER99SB-ILDN force field was implemented in all cases. All the trajectories showed a similar trend but with a lower RMSD than in the CHARMM36 force-field except TLR4Y328A which showed fluctuation from the 65–85 ns time-point (Figure S2B). In compactness, a similar trend was followed with lower *Rg* value as compared to the CHARMM36 force-field

(Figure S2E,F). The two force-fields generally produce very similar overall RMSD and Rg. It has been found that CHARMM36 tends to produce more structure variation shifted towards disordered conformations, while as AMBER99SB-ILDN over stabilizes the native fold [28,29].

**Figure 1.** Characterization of an epitope using MD. (**A**) Epitope location in the TLR4–MD2 dimer interface and the residues of the epitope are labelled. (**B**) The mutational effect of these residues on TLR4 stability was analyzed by RMSD, where TLR4Y328A (red), TLR4N329A (green), and TLR4D371A (grey) deviated substantially from TLR4w<sup>t</sup> (black) while RMSD of TLR4K349A (blue) remained constant. (**C**) *Rg* determining the compactness of a protein was analyzed for TLR4w<sup>t</sup> and its muteins. For TLR4Y328A (red) and TLR4N329A (green), the compactness was lost because their R*g* was higher than that of TLR4w<sup>t</sup> (black). (**D**) RMSF of TLR4w<sup>t</sup> and its muteins was measured. TLR4N328A, TLR4Y329A, and TLR4K349A showed a fluctuation more vividly in the epitopic region as illustrated in the zoomed view.

#### *2.3. Mutation-Induced Conformational Changes in TLR4 Structure*

To assess harmonic motions, NMA of TLR4w<sup>t</sup> and muteins was performed. TLR4w<sup>t</sup> showed an outward-motion tendency at both N and C termini; by contrast, TLR4Y328A featured an inward movement at both termini (Figure 2A,B). The Y328A mutation seemingly changes the rigidity of the convex surface of TLR4 and allows the N and C termini to bend inward and come closer. A conformational change of such magnitude could affect the dimerization interface of TLR4 monomers in the (TLR4–MD2)<sup>2</sup> complex (where MD2 is myeloid differentiation factor 2). The terminal regions of TLR4N329A showed a vertical motion, inducing structural torsion, which resulted in a reduction in the number of internal hydrogen bonds and deformation of the solenoid structure of TLR4 (Figure 2C). TLR4K349A underwent motions similar to those of TLR4w<sup>t</sup> (Figure 2D). This finding was in line with the RMSF and RMSD plots, where TLR4K349A did not show significant deviation and featured TLR4wt-like behavior. Other muteins, TLR4K351A, TLR4E369A and TLR4D371A, showed a peculiar motion; however, the magnitude of these movements was not significant (Figure S3).

**Figure 2.** Motion and FEL of TLR4w<sup>t</sup> and its muteins. The motion of N and C termini in a porcupine plot of (**A**) TLR4w<sup>t</sup> is outward and in (**D**) it is inward with large magnitude, and for (**G**) TLR4Y329A and (**J**) TLR4K349A, it is irregular. The FEL was computed using PC1 and PC2 as reaction coordinates for (**B**) TLR4wt, where all conformations remained confined to one minimum. (**E**) TLR4Y328A conformations split into two minima with higher energy barrier. (**H**) TLR4Y329A and (**K**) TLR4K349A split into three minima. A representative lowest-energy structure from local minima indicated by red circles was superimposed (**C**,**F**,**I**,**L**) with the respective first or 0 ns frame structure in TLR4w<sup>t</sup> and mutants shown along with the RMSD values. The black circle exemplifies the variation in structure at the epitopic region or dimer interface and termini.

#### *2.4. Visualization and Identification of Principal Motion of TLR4w<sup>t</sup> and Muteins*

The dominant motions in TLR4w<sup>t</sup> and muteins were determined through principal component analysis (PCA) in which the first five eigenvectors captured most of MD motions. Most of the dynamic structural information for each system was captured by 5–6 eigenvectors with considerable fluctuations, while fluctuations of the remaining eigenvectors had low amplitude. The first three eigenvectors accounted for 63.9%, 67.4%, 75.2%, 63.6%, 75.0%, 65.2% and 74.7% of principal motions for TLR4wt, TLR4Y328A, TLR4N329A, TLR4K349A, TLR4K351A, TLR4E369A, and TLR4D371A, respectively (Figure S4). The first three eigenvectors of TLR4w<sup>t</sup> and muteins were projected to determine interconformer relations (Figure S4). In Figure 2 and Figure S2, continuous color representation from red to blue indicates the change in conformation. The distribution of conformers in TLR4w<sup>t</sup> remained mostly in two subspaces, i.e., red and blue. In the case of TLR4Y328A, the conformers were scattered into a subspace owing to periodic jumps between different conformations. Unlike TLR4wt, TLR4N329A conformers dispersed and covered a large subspace, as was the case for TLR4Y328A. In the RMSD graph, TLR4Y328A and TLR4N329A showed instability because greater fluctuations were observed in PCA. The graph revealed continuous periodic jumps from one state to another, representing the instability of the alanine muteins of TLR4Y328 and TLR4N329. The transition between different conformations in these two muteins makes the epitope less exposed due to conformational changes. The other muteins showed different behavior in the PCA plots (Figure S4).

#### *2.5. Energetics of Conformation Transitions*

PCA was carried out to extract the predominant motions of the TLR4w<sup>t</sup> and mutein structures; these data can provide a better picture of conformational changes and structural evolution along the trajectory of the MD simulation. The first two principal components were employed to calculate Gibbs free energy of each system via plotting in the free-energy landscape (FEL). A plateau in an FEL plot indicates the state of a protein corresponding to its energy and structural transition to the lowest-energy state along a simulation path. An overall view suggested that TLR4w<sup>t</sup> remains in a single large global minima basin in the dark blue area (Figure 2B). Structure coordinates from the energy minima were randomly sampled to check conformational changes. Three samples were superimposed over the first frame. There was no substantial change in these structures, and their average RMSD was recorded and found to be 1.939 Å, which meant that TLR4w<sup>t</sup> retains its conformation during the simulation.

In contrast, the energy plateau for TLR4Y328A, TLR4N329A and TLR4K349A was found to split into two, four and three energy minima, respectively, separated by high energy barriers (Figure 2E,H,K). Besides, the representative structures randomly sampled from these energy minima showed a significant change in their N and C termini. Moreover, conformational changes in leucine-rich repeat (LRR) regions were observed. LRR13 is in the proximity to the dimeric interface of TLR4 and manifested a significant difference in the muteins. The α-helix near the dimer interface also underwent substantial conformational changes. Nonetheless, TLR4K351A, TLR4E369A, and TLR4D371A structures extracted at different time points from corresponding deep blue regions of the FEL did not show any differenceinstructuraltopologyandRMSD(FigureS3).

#### *2.6. The Influence of Mutations on the TLR4-mAb Interaction and Binding Affinity*

We noticed that the mutations caused various structural and conformational changes in TLR4, thereby increasing internal energy, and decreasing stability. Moreover, we noted that mutations TLR4Y328A, TLR4N329, and TLR4K349A, which mask the epitopes, influenced the structural integrity of TLR4. Taking into consideration the effect of these three mutations on TLR4, the mutated structures were docked with Hu 15C1 Fab and were simulated for 150 ns to gain insights into the detailed structural dynamics and comparative changes in binding affinity. Unlike the TLR4wt–Hu 15C1 complex, mutant complexes yielded a continuous rise in the RMSD plot (Figure 3A). The RMSD rise in these mutant complexes showed

a similar trend when the AMBER99SB-ILDN force-field was implemented (Figure S5). This observation suggested that these site-specific point mutations led to structural changes in TLR4 that subsequently affect the binding and stability of the complex. Given that these mutations break the hydrogen bond network in the surroundings and partially mask epitope hotspots, we propose that these residues have an indirect impact on the binding to Hu 15C1.

**Figure 3.** Epitope prediction and docking mode and stability analysis. (**A**) RMSDs of four complexes in which mutated complexes showed an increasing trend. (**B**) Residues for the construction of three conformational epitopes. Modes of docking of Hu 15C1 with TLR4 at sites (**C**) Ctl, (**D**) EP1, and (**E**) EP2 with superimposition of the human TLR4–MD2 dimer. RMSDs of six complexes are illustrated in (**F**), where C2E3–Ctl (red) showed a maximum deviation, and all remaining complexes were stable. The number of hydrogen bonds during the simulation was determined in (**G**), where C2E3–EP1 (magenta) featured the highest, C2E3–Ctl (red) and Hu 15C1–EP1 (green) moderate, and Hu 15C1–Ctl (black), Hu 15C1–EP2 (blue), and C2E3–EP2 (brown) the lowest number of hydrogen bonds formed during the simulation.

To investigate the influence of these mutations on the binding affinity of TLR4 for Hu 15C1, binding free energies of TLR4–mAb complexes were calculated. As expected, we observed a decrease in the binding energy of mutein complexes as compared to the

TLR4wt–Hu 15C1 complex (Table 1). Total binding free energy in each mutant system diminished by ~7–9%, which confirmed the decrease in the affinity of the mAb. The decrease in binding energy of the mutein complexes occurred because the interaction was reduced when residues, i.e., Y328, N329, and K349, were mutated into alanine in mutein complexes. The mutated residues do not interact with the antibody as in TLR4wt–Hu 15C1, which entails that the energy contributed by the mutated residue is lost. Thus, the decrease in the total binding energy of the mutein complexes was observed, suggesting that the original amino acids in this position is energetically vital for antigen binding.


**Table 1.** Binding free energy of Hu 15C1 and C2E3 toward TLR4.

#### *2.7. Epitope Prediction and Interaction Analysis of the Epitope–Paratope Interface for Rational Design of Antibodies*

The hotspot residues, as mentioned above, affect the conformation of TLR4, thereby reducing mAb-binding affinity. Therefore, we further investigated TLR4 epitopes and remapped them using conformational CDR information on Hu 15C1 and its derivative antibody, C2E3. Conformational epitopes on TLR4 were predicted based on an antibody– antigen score. Among the identified epitopes, EP1 and EP2 were selected and compared with a control (referred as Ctl) site. The residues of these epitopes are listed in Figure 3B and the Ramachandran plot of Ctl, EP1, and EP2 is shown in Figure S6. Ctl residues overlap with a few residues of EP1. Ctl is located near the TLR4 dimeric interface, whereas EP2 is located at the TLR4-MD2 interface. Hu 15C1 and C2E3 Fabs were docked to EP1, EP2, and Ctl epitope patches. In each resultant complex, the docked pose of a representative member from the most populated cluster was selected, as depicted in Figure 3C–E and Figure S5A–F. Superposition of the docking complexes revealed overlapping of one of mAb's chain at sites Ctl and EP1 (Figure S7G,H).

To assess the stability of the TLR4–mAb complexes, MD simulations were performed and then analyzed. Hu 15C1–Ctl, Hu 15C1–EP1 and C2E3–EP1 had a similar pattern of the RMSD of backbone atoms, indicating that these complexes are highly stable. Nevertheless, C2E3–EP2 yielded a deviation in the last 75 ns, whereas C2E3–Ctl manifested a continuous increase in RMSD throughout the trajectory (Figure 3F). The marked RMSD increase could be explained by the finding that C2E3 induces a surge in the energy of epitopic or CDR residues upon TLR4 binding (data not shown). A similar tendency was observed in the RMSF plots of these complexes (Figure S7I,J). Residues 200–410 of TLR4 and 150–200 of the antibody in TLR4–C2E3 underwent more fluctuation compared to similar residues in the other complexes. The robustness of these finds was investigated under the effect of the AMBER99SB-ILDN force-field. The RMSD graph (Figure S8) of all the complexes showed a similar trend except C2E3–Ctl, although with a continuous increase in RMSD but much lower than when simulated with CHARMM36 (Figure S8B). Again, the lower RMSD in AMBER99SB-ILDN could be due to the overstabilization of native conformation. Nevertheless, C2E3–EP1 showed a similar trend in RMSD but slightly higher than in the CHARMM36 force field (Figure S8D).

To investigate the antigen–antibody interaction, we performed a dynamic analysis of the interaction of Hu 15C1 and C2E3 with TLR4; this approach provides the basis for rational design of mAbs. Primarily, the number of hydrogen bonds between each mAb and TLR4 was determined with a cutoff of <5 Å (Figure 3G). In Hu 15C1–Ctl, Hu 15C1–EP1, and C2E3–EP1, the total number of hydrogen bonds was 12–15 and persisted during the simulation. By contrast, in Hu 15C1–EP2 and C2E3–EP2, the hydrogen bond number during the simulation was within the range of 7–12.

To illustrate the distribution of interacting residues of CDRs in the antibodies toward TLR4, the radial distribution function and minimum distance were computed. In a comparison between the interactions of the two mAbs with the TLR4, the CDRs of the C2E3 antibody showed stronger affinity because the distribution of CDR residues interacting with epitopes was higher in comparison with Hu 15C1 (Figure 4). The molecular association of pairwise residue contacts within 5 Å toward TLR4 at three epitopic sites was greater for C2E3, indicating its stronger binding in comparison with Hu 15C1. The peak distribution of Hu 15C1 was higher at Ctl with a peak amplitude above 20 for the Ser368– Gly91 molecular association. The height of distribution peaks for Gly394–Lys310, Glu396– Asp310, and Thr373–Pro312 reached 20 and formed at shorter distances (Figure 4A). In comparison with Hu 15C1, the molecular association peaks formed for C2E3 toward TLR4 at Ctl had a lower amplitude. The pronounced peak of Lys376–Ser101 was above 15 and formed within the 5 Å distance. On the other hand, the rest of the peaks were smaller and formed at above 5 Å; therefore, the pairwise antibody–antigen residue contacts were weak (Figure 4B). The minimum distance for most of the interactions of the Hu 15C1–Ctl complex persisted during the simulation (Figure S9A and Table S1).

**Figure 4.** Distribution of minimum distances between the residues of Hu 15C1 and C2E3 interacting with TLR4 at sites (**A**,**B**) Ctl, ( **C**,**D**) EP1, and (**E**,**F**) EP2.

The minimum distance during the simulation in the C2E3–Ctl complex remained above 4 Å for some interacting residues (Figure S9B and Table S1). Both mAbs Hu 15C1 and C2E3 at the EP1 site were found to be buried inside, as evidenced by the docking pose (Figure 3C–E) where the maximum surface interacts with TLR4 at this site; therefore, a maximal interaction between the interacting residues was recorded. The highest peaks in the C2E3–EP1 complex formed at the shortest distance of <5 Å. The most pronounced peaks were rendered by Arg231–Asp100 and Arg201–Gly278 with an amplitude of 30, whereas the peaks formed by Asp353–His53, Glu260–Lys267, and Lys328–Gly278 were >20 in height (Figure 4D). In the Hu 15C1–EP1 complex, the molecular association peaks of Glu350–Lys278, Gln307–Gly91, and Lys328–Asp270 were above 10 and formed at a distance of <5 Å, while the amplitude of the rest of the contacts remained below 10 and formed at a distance of >5 Å (Figure 4C). The interacting residues in both complexes mostly remained intact except for Glu152–Thr241, Arg201–Tyr318, and Glu307–Gly91 in Hu 15C1–EP1 and Arg201–Ser278 in C2E3–EP1 during the simulation (Figure S9C,D and Table S2). The Hu 15C1 distribution of interactions remained below 10 except for the pronounced Arg296–Glu1 peak, which showed an amplitude up to 40 at the EP1 site. C2E3 featured a higher distribution of residue interactions at the EP2 site because Arg238–Tyr54, Glu295–Ser101 and Lys245–Asn311 peak amplitudes were approximately 30, and the peaks formed within the 5 Å distance (Figure 4E,F). In Hu 15 C1–EP2 and C2E3–EP2, the minimum distance between the interacting residues was registered above 5 Å for a few interactions (Figure S9E,F and Table S3). The radial distribution functions for Hu 15C1 and C2E3 at the EP1 site were more vigorous as the pronounced peaks formed at <5 Å distance, and the minimum distance of most of the interacting residues of Hu 15C1 and C2E3 at the EP2 site was >4 Å. In contrast to the EP1 site, the gap remained <4 Å during the simulation (Figure 4 and Figure S6). The interface interaction analysis revealed that at the EP1 site, most of the mAb–TLR4 interactions are energetically favorable, as suggested by the peak distribution at this epitope.

The interactions associated with the formation of paratope–epitope interface residues resulted in a free-energy change for the binding reaction. To compare the binding affinity for TLR4 between Hu 15C1 and C2E3 at epitopic sites Ctl, EP1 and EP2, total binding energies were computed by the molecular mechanics Poisson–Boltzmann surface area (MMPBSA) method. The total binding energies of complexes Hu 15C1–TLR4 and C2E3–TLR4 were −3358.261 and −5068.979 kJ/mol, respectively, at the EP1 site (Table 1), thus ensuring their stability in comparison to the Ctl site with a total binding energy of −2868.266 kJ/mol toward Hu 15C1 and −2797.202 kJ/mol toward C2E3. In Table 1, the binding free energies imply that the affinity of Hu 15C1 and C2E3 for the EP2 site is less than that for the EP1 site but greater than that for the control site. The binding free energy of both mAbs toward TLR4 at EP1 is much higher than that at EP2. The more negative the binding energy value, the stronger is the binding affinity of the antibody for its specific target antigen [30].
