A stepwise docking molecular dynamics approach for simulating antibody recognition with substantial conformational changes.

Conformational changes or rearrangements are common events during inter-biomolecular recognition. Tracking these changes are essential for exploring the allosteric mechanism and it is usually achieved by molecular dynamics simulation in silico. We previously identified a broad-neutralizing antibody against H5 influenza virus, 13D4, and solved the crystal structures of the free 13D4 Fab and its complex with hemagglutinin (HA). Structural comparison of the unbound and bound 13D4 Fabs showed that the heavy chain complementarity-determining region 3 (HCDR3) undergoes a substantial conformational rearrangement when it recognizes the receptor-binding site (RBS). Here, we used molecular dynamics (MD) to simulate the conformational changes that occur during antibody recognition. We showed that neither conventional MD nor steered MD could recapitulate the loop fitting of the RBS structure contour. Consequently, to simulate these challenging conformational changes, we engaged a stepwise docking MD method that allowed for the gradual docking of the ligand to receptor. This new method recapitulates the bound shape of the HCDR3 and provides the best approximation of the shape rendered by the co-crystal structure, with an RMSD of 0.926 Å. This strategy affords a flexible MD approach for simulating complicated conformational changes that occur during molecular recognition, and helps to provide an understanding of the involved allosteric mechanism.

Introduction

Protein conformational changes are common in nature, and often endow the protein with special structural-based functions, such as allosteric enzyme regulation or control over structural modulations required for the dynamicity of ligand-gated ion channel proteins [1]. Conformational rearrangements are also necessary for protein–protein binding, particularly at the antibody–antigen interface. Such changes can be as simple as the rotation of side chains or may involve more complex modifications, such as global structural transformation [2]. As yet, two hypotheses could account for these allosteric phenomena, which are conformational selection model and induced fit model [3]. The former regards the molecular recognition as selection over population shift, where the ligand or receptor could adopt several conformations, including real bound-state conformation. The latter hypothesis assumes that the ligand would undergo a local structural change near the binding site triggered by its interaction with the receptor.

In antibodies, variations in the complementarity-determining region (CDR) contribute to the distinct antibody surface topographies necessary for epitope binding [4]. Conformational rearrangement in the CDR may occur during antigenic binding to fit with an optimal shape and facilitate appropriate physico-chemical properties. The best way to understand the underlying allosteric mechanism is to dynamically analyze the actual binding process rather than by comparing the antibody-bound and -unbound crystal structures [5]. Indeed, noncovalent and transient biological interactions dominate most of life, and insight into the dynamicity of such conformational changes as they take place could lead to improved peptide, drug, and antibody designs [6]. Such knowledge of the unique recognition and neutralization mechanisms triggered by conformational rearrangements is thus of great import [7].

Computer simulation techniques are increasingly used to predict the structural and thermodynamic properties of biological molecules and, sometimes, more complicated protein complexes, such as those with roles in molecular biological machinery [8]. Simulations using molecular dynamics (MD) are becoming increasingly popular as a tool to reproduce the behavior of a system and to investigate protein functions and properties [9]. The method employs Newton’s laws to solve the classical equations of motion for a group of atoms, and evaluates motion in a predefined periodic system; this includes water molecules, ions, small molecules, and biological molecules under specific conditions [10], [11].

In terms of single protein allostery, MD can provide atomic insight into the relation between allosteric process and functions, which thereafter facilitates interpretation of the biological mechanism [1]. For example, Yang and colleagues described how the acid-sensing ion channel 1 (ASIC1) undergoes a conformational change characterized by a “twist-to-open” motion [12]. However, it remains challenging to capture the dynamic conformational changes in recognition process of antibody–antigen and even protein–protein using general molecular simulation methods, due to the intrinsically global effect of complex dynamics and the use of potentially inappropriate sampling methods. To address these issues, we proposed a semi-biased stepwise docking molecular dynamics (SDMD) simulation in the light of induce-fit model to restore the dynamic structural changes of the ligand when approaching. The principle behind that was the differentiation of the docking process of two units. We carried out tests using our previously described flu antibody (13D4) in its binding to the head of hemagglutinin (HA), wherein we showed a significant conformational shift in the HCDR3 loop during binding by comparing the bound and unbound crystal structures. We reproduced the allosteric fit using different MD methods including Conventional MD (CMD), steered MD (SMD), and SDMD, hoping to mimic the actual process of antibody–antigen recognition. Our analysis of the trajectories suggests the SDMD best simulates the HCDR3 loop to its binding status. This work provides a dynamic view for the in-depth exploration of the conformational changes and insight into ways to use MD to perform flexible protein-protein docking.

Results

Structural analysis of the conformational changes of H5 13D4 antibody binding to the HA receptor binding site.

We previously identified a broad-spectrum neutralizing monoclonal antibody (nAb) 13D4 that could neutralize all representative H5 strains isolated from 1997 to 2009. We solved the crystal structures of the free 13D4 Fab and its co-crystal complex with the HA (strain VN1194) head region at resolutions of 2.3 Å and 2.33 Å, respectively [13]. Intriguingly, when the 13D4 Fab was superimposed with the 13D4:HA complex, we found distinctions between the HCD3 loops in the free Fab and in the Fab-bound complex, referred to F (free) and B (bound) states, respectively (Fig. 1a). In aligning the free Fab to the complex, HCDR3 in the F state showed partial clashing with HA (Fig. 1b), and the residues of antibody involved in the clash were A103H, V104H, E105H and R106H. Structurally, an approximately 8.2-Å swing at the tip of Fab in the F state as compared with the B state, indicates a conformational rearrangement where the L-bent loop adapt to fit into the grooved floor of the RBS in antibody–antigen recognition (Fig. 1c).

Fig. 1

The F (free) and B (bound) states of 13D4 HCDR3 in Fab-bound and unbound structures. a) Alignment of the unbound and bound 13D4 Fabs HCDR3, rendered in red and yellow ribbon representation, respectively, generates RMSD values of 1.123 Å and 4.063 Å for whole Fab and HCDR3, respectively. b) Residues in Fab-unbound state that collide with HA are shown as sticks (red) and independent surface (salmon) in the context of the HA surface. c) The swing motion of the Fab-unbound (red) and Fab-bound HCDR3 (yellow) indicates a maximum shifting distance of 8.2 Å at the tips. Abbreviations: HCDR3, heavy chain complementarity-determining region; VH, heavy chain variable region; VL, light chain variable region; HA, hemagglutinin; RMSD, root-mean-square deviation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Simulation of the conformational rearrangement of 13D4 Fab by CMD

We attempted to simulate this conformational change using several simulation approaches based on population shift model and induce fit model.

First, a 100 ns unbiased dynamics simulation of isolated 13D4 Fab was performed to inspect if the conformation of antibody would convert to the B state spontaneously. The stable curves of RMSD values for both antibody and HCDR3 indicated the convergence of structure (Fig. 2a). The simulated structure of HCDR3 still reveal distinction (RMSD value: 3.547) compared to co-crystal structure (Fig. 2b).

Fig. 2

Conventional molecular dynamics (CMD) of isolated 13D4 Fab and its complex with HA. a) RMSD values of Cα of the antibody (up) and HCDR3 loop (down) showing an equilibrium after 100 ns run time. b) Simulated structure of the HCDR3 loop of isolated 13D4 Fab (slate) is compared with that of crystal structure (yellow). c) Total energy and RMSD values of Cα of the complex showing an equilibrium after 600 ns run time. d) Conformations extracted at 200 ns (pink), 400 ns (magenta) and 600 ns (purple) are aligned by HA to the crystal structure (yellow) and are shown as a cartoon. e) Three simulated structures of the HCDR3 loop are compared with the F (free; red) and B (bound; yellow) states of 13D4 Fab and the minimized F-state Fab (deep teal), showing a negligible offset between the Fab-unbound HCDR3 loop with (red) and without (deep teal) energy minimization. f-h) CMD recovers some critical interactions including the pi-interaction between Val104H (purple) and Trp153 (cyan) (d), two hydrogen bonds between Glu105H (purple) and Ser137 (cyan) (e), and one electrostatic contact between Glu190 (cyan) and Arg106H (purple) (f). The corresponding amino acids in the co-crystal complex are rendered as yellow sticks. Abbreviations: HCDR3, heavy chain complementarity-determining region; VH, heavy chain variable region; VL, light chain variable region; HA, hemagglutinin; RMSD, root-mean-square deviation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Inspired by the “induced-fit” model, we initially wondered if the CDRs of antibody were reshaped and forced to rearrange by global interactions after binding to HA. We performed the simulation using an initial state where the free Fab was aligned to the complex structure and the bound Fab was evicted. The interatomic clashes were avoided by a rough energy minimization and the overall loop shape was unchanged (Fig. 2e). A regular MD simulation was performed for 600 ns, and the energy and RMSD changes over time suggested that the system was well equilibrated at the end of the procedure (Fig. 2c).

Among the 600-ns trajectories, the atomic models at 200 ns, 400 ns, and 600 ns were sampled for structural comparisons (Fig. 2d). As the simulation proceeded the HCDR3 loop underwent a conspicuous drift from its original position within a short timeframe, which may have been exerted by an integrated effect of the surrounding atoms, thereby allowing the loop to adapt to the groove region, as described in the co-crystal structure (Fig. 2e). This adaptation involved the establishment of some critical interactions involving Val104, Glu105, and Arg106 (Fig. 2f-h; Table 1).Yet, the equilibrium conformation of HCDR3 was distinctly different from that noted for the 13D4:HA complex, and, crucially, when the simulation time was extended to 600 ns, no other interactions found in the crystal structure could be established as compared with 200 ns and 400 ns models, and the structure of HCDR3 remained stable as the RMSD for the loop at 200 ns was 2.44 Å, which is almost the same as that at 400 ns (2.31 Å) and 600 ns (2.33 Å) (Fig. 2e). It can be inferred that the conformational rearrangement might not proceed if 13D4 Fab and HA were put together manually to form a complex at the beginning of the simulation, and this prolonged CMD approach was thus not suitable for the simulation.

Table 1

Interactions between HCDR3 of antibody and HA.

Key residues on HCDR3	Key residues on HA	Interaction mode	CMD	ASDMD
Val104	Trp153	Pi-Alkyl	✓	✓
Glu105	Ser137	Hydrogen Bond	✓	✓
		Hydrogen Bond	✓
Arg106	Asp187	Electrostatic		✓
		Hydrogen Bond		✓
	Glu190	Electrostatic	✓	✓
		Hydrogen Bond	✓	✓
Asp107	Lys222	Electrostatic		✓
Trp108	Ala189	Pi-Alkyl		✓
	Lys193	Pi-Cation
		Pi-Cation
		Pi-Alkyl		✓
		Pi-Alkyl

We next generated a detached complex model by pulling the antibody away from the binding interface at a distance of 5 Å along the direction passing the HA center of mass (COM). In the expectation, the Fab would gradually dock to HA while driven by the interaction between Fab and HA during MD simulation as observed in other studies [14], and meanwhile conformational changes might accompany this association.

The 5-Å separated Fab:HA system was simulated by CMD with similar parameter settings as mentioned above. Two atoms—HCDR3:Val104:Cα and HA:Trp153:Cα—were selected as the distance indicator to measure the potential approaching of the antibody; the original distance was 9.56 Å in the co-crystal structure. We observed no potential approaching of the antibody during the 100 ns of the simulation. The distance of two indicator atoms in the MD trajectory remained steady at around 12.0 Å when the simulation reported reaching an equilibrium; this is far from the proposed distance of 9.56 Å (Fig. 3a). Yet, an obvious deflection of the HCDR3 toward the co-crystal position was detected (Fig. 3b). The Fab failed to land in HA due to the noncomplementary surface contours of the binding site and the HCDR3 seemed to be trapped in a low-energy conformation. Theoretically, accurate free-energy calculations that recover Boltzmann distribution should cross the high-energy barriers that separate low-energy states, and, in final, cross them several times to obtain converged statistics. For specific biological processes, including conformational rearrangements, they are often random events and may occur on a time scale of microseconds or even milliseconds [15]. Thus, such a system may need additional time to cross the barrier.

Fig. 3

The 100-ns CMD of a pre-detached 13D4f:HA complex. a) The curve of the distance of the referenced atoms (HCDR3:Val100:Cα and HA:Trp153:Cα) in the simulation trajectory. b) Positions of the HCDR3 loops by HA superimposition: the pull-away loop in the initial state is colored blue whereas the loop extracted from a 100-ns simulation is orange; the unbound (red) and bound (yellow) states are also depicted as references. Abbreviations: HA, hemagglutinin; CMD, conventional molecular dynamics. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Simulation of the conformational adaption of 13D4 HCDR3-loop by SMD

To facilitate the antibody approaching, we next explored SMD to simulate the Fab:HA detached system. In SMD, an external force is applied to one or more atoms of the system, pulling the forced part in a predetermined direction. Again, we generated a separated 13D4f:HA model with a distance of 5 Å. To simplify and accelerate the simulation, only the HCDR3 loop (aa. 98–111) was used. Three rounds of individual SMD with forces of different values were performed to determine the most suitable force to drive the HCDR3 loop to anchor to the RBS.

A series of forces ranging from 10 pN to 10,000 pN were first applied to the structure for a 4-ns simulation. The loop moved away from its initial position when the force was set to 10 pN or 100 pN, and remained near the binding site at 1,000 pN or 10,000 pN (Fig. 4a). In the second round, a 2-ns SMD was performed with forces ranging from 1,000 pN to 10,000 pN. We found that a 4,000 pN or 6,000 pN force could push the two atoms close to the RBS (Fig. 4b). We then narrowed down the external force values to 4,000 pN to 6,000 pN in 500-pN increments and the third round of SMD showed that the HCDR3 loop stay engaged with HA only in 4,000 pN and 5,500 pN simulations, with the indicator atoms at a distance of 9.9 Å to 11.3 Å (4,000 pN) and 9.7 Å to 11.5 Å (5,500 pN). There was no linear correlation between the applied force and the resultant distance of the indicator atoms (Fig. 4c).

Fig. 4

Steered molecular dynamics (SMD) simulation with different forces applied to the detached HCDR3 loop. a) The final conformation of the HCDR3 loop from the stage I SMD simulations with various forces applied: 10 pN (black), 100 pN (magenta), 1,000 pN (green), and 10,000 pN (purple). The loops applied with a 10 pN or 100 pN force failed to approach the surface of HA as compared with the reference B (bound; yellow) and F (free; red) states. b) The distance (Å) of the referenced atoms over the SMD simulation time following the application of different external forces (1,000 pN to 10,000 pN). The deep gray columns indicate the average distance during the 2,000 ps simulation, which suggests that forces of 4,000 pN (grey line) and 6,000 pN (gold line) applied to the HCDR3 loop can guarantee the approach of the loop to HA. c) The final distance between HCDR3:Val100:CA and Trp153A:CA under different applied forces. Only when the force was set to 4,000 pN or 5,000 pN could the loop approach closely to HA. d) Conformational comparisons of the HCDR3 loops of the co-crystal structure (yellow), the aligned unbound structure (red), the SMD with an applied force of 4,000 pN (orange), and the SMD with an applied force of 5,500 pN (blue). Abbreviations: HCDR3, heavy chain complementarity-determining region; HA, hemagglutinin. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

When the SMD simulation finished, the force was released, and an additional 100-ns simulation was undertaken. The resultant RMSD and energy trajectory suggested that the simulation system reached a relative equilibrium by the end of the simulation (Supplemental Figure S2a and b). However, the steered loop under either 4,000 pN or 5,500 pN force could not bind deep within the RBS floor, unlike in the co-crystal structure. Additionally, the simulated loop shape exhibited no significant resemblance to that in the co-crystal structure (Fig. 4d), suggesting that SMD using a fixed force was not a suitable strategy for antibody recognition simulation.

Simulation of conformational adaptation of 13D4 HCDR3-loop by SDMD

SMD simulation demonstrated that the HCDR3 loop under external forces moved toward HA and rapidly reached the RBS within sub-picoseconds, ranging from 40 ps to 500 ps (equivalent to 20,000 steps to 250,000 steps), likely insufficient at time scale for antigen-induced structural rearrangement in silico. Based on this inference, we proposed a stepwise docking molecular dynamics (SDMD) simulation method to satisfy a gradual structural rearrangement during the 13D4 Fab stepwise docking to HA, given sufficient MD equilibrium is implemented for every step. The overall procedure for the simulation was summarized in Fig. 5. Similar to the SMD simulation, the antibody was first pulled away from HA, then subjected to a regular MD simulation with residues unrelated to the interface being fixed; this kept the antibody from drifting away but allowed the loop to adapt itself to the surrounding contours. When the MD trajectory reported reaching an equilibrium, the loop was moved 0.1 Å toward the RBS. The stepwise docking simulation was reiterated and the loop structure at each step was inspected until the distance dropped to zero. As compared with SMD, our SDMD approach incorporating differentiation thought into the binding process allowed us to acquire the most stable and rational conformation in each small step; slowing down the binding simulation allowed the interface to find its reasonable conformation. The model incorporating the whole 13D4 Fab and HA was submitted to an automatic SDMD (ASDMD) simulation.

Fig. 5

Flowchart of the stepwise docking simulation for 13D4 and HA. Abbreviations: HCDR3, heavy chain complementarity-determining region; HA, hemagglutinin; CMD, conventional molecular dynamics.

To view the holistic process of conformational change and underlying mechanism of 13D4 HCDR3 loop, major intermediate conformation scenarios were extracted from the whole SDMD trajectories and interaction variation were analyzed (Fig. 6a). As the simulation proceeded, electrostatic interactions were formed by HCDR3:Arg106-HA:Glu190, HCDR3:Arg106-HA:Asp187 and HCDR3:Arg105-HA:Lys222, of which the first two are critical for 13D4 binding in the crystal structure (Fig. 6b; Supplemental Figure S3). HCDR3:Ala103:O oriented toward HA:Gln:226:N to form hydrogen bonds, resulting in disruption of the intermolecular hydrogen bonds between Ala103:O and Arg106:N, following by a structural relaxation of the tip of HCDR3 loop (Fig. 6b). Distinct conformational rearrangement of HCDR3 loop firstly occurred when the distance reduced to 3.7 Å. The loop deflected towards the groove center where the Val104 moved against the floor concomitantly. HCDR3:Glu105 then established a stable hydrogen bond with HA:Ser137 (Observed in co-crystal structure), while lost the salt bridge interaction with HA:Lys222 (Fig. 6c; Supplemental Figure S3). The allosteric changes leaded to the conformational rearrangements of both side-chains and backbones of residues near Thr101, presenting a loop bulge (Fig. 6c). With further decreases in distances, the Val104 accompanied with the tip of HCDR3 stretched deeply into the hydrophobic groove. The loop bulge showed a shift away from the LCDR1 side (Fig. 6d). The next important turning point was that when the distance dropped to 1.6 Å, HCDR3 underwent a slight conformational rearrangement, with Asp:107 moving towards HA:Lys222 and gradually replacing Glu:105 to form hydrogen bond interactions with HA:Lys222 (Fig. 6e; Supplemental Figure S3). The movement of Asp:107 caused a further conformational adjustment of the adjacent Thr101, then at the simulation step of 1.1 Å, Thr101:O formed a hydrogen bond with Trp108:O, which was also observed in the co-crystal structure (Fig. 6f). Up to this point, the conformational rearrangement of HCDR3 loop was basically completed, showing a L-shape loop fit into the grooved floor. The ASDMD lasted 854 ns. In the final structure, the HCDR3 loop tended to superimpose the position of that in the co-crystal structure (Supplemental Movie S1). Before the SDMD, the angle formed by Arg106:Cα and Val104:Cα of 13D4:HA model and Val104:Cα of 13D4f:HA was 75.7°. After the simulation, the detached antibody returned to its original position and the angle decreased accordingly, falling to nearly 31.6° (Fig. 6g). Most importantly, the HCDR3 loop in the automated SDMD approach could best matched with that of the crystal complex, with an RMSD value of 0.926 Å. More interactions were recovered with the SDMD (Table 1), including electrostatic interactions and hydrogen bonds between HCDR3:Arg106 and HA:Asp187, and the Pi-Alkyl interaction formed by HCDR3:Trp108–HA:Ala189 (Fig. 6h and i). Overall, SDMD of the whole 13D4f:HA model in auto-supervised manner affords an MD approach to simulate the conformational changes of the HCDR3 loop during antibody binding, and the method may provide way to explore the allosteric mechanism behind the antibody-antigen or even protein–protein interactions that incorporated with complicated conformational rearrangement.

Fig. 6

The stepwise docking molecular dynamics (SDMD) simulation on the 13D4:HA complex, using a script with auto-equilibrium judgement (ASDMD). a) The conformational change of 13D4 HCDR3 loop during SDMD simulation. A serial of equilibrium trajectories at 3.9 Å, 3.7 Å, 1.6 Å, 1.7 Å and 1.1 Å intermediate distance of 13D4 from HA were rendered in cartoon for 13D4 and surface mode for HA while the complexes being superimposed against HA. The final docking model of SDMD and the co-crystal structure were shown in same mode with color pink and yellow, respectively. The complete picture showed that the 13D4 HCDR3 loop gradually resembles the shape as co-crystal structure while approaching step-wise to HA RBS during SDMD. b-g) The conformation variation occurred at HCDR3 loop for every SDMD steps showed in a). h-i) The structure comparison of HCDR3 loop for CMD (purple), SDMD (pink) and co-crystal structure (yellow). Abbreviations: HA, hemagglutinin; SDMD, stepwise docking molecular dynamics; CMD, conventional molecular dynamics. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Discussion

Structural changes play significant roles in receptor-ligand recognition. Despite this, the fundamental allosteric mechanisms involved in mediating these bindings are not fully understood. MD simulations have been widely used to trace the trajectory of atoms in predefined periodic systems, including target proteins, water molecules, ions and other macromolecules, the knowledge of which has afforded a better understanding of many specific biological processes.

In this study, we tested several simulation approaches based on the principle of conformational selection model and induce fit model to trace the allosteric trajectory of 13D4 HCDR3. Unbiased simulation of the free 13D4 Fab and aligned 13D4f:HA complex showed that no B-state conformation of HCDR3 could be observed. The antibody was then detached from the antigen to emulate its binding to HA, but we failed to capture any tendency that the antibody would approach the antigen within an acceptable calculation time. The SMD simulation based on a detached complex did not yield better simulation results; the HCDR3 loop could not fit into the RBS, presumably due to an inadequate equilibrium time and fast movement of the detached loop. To solve these issues, we proposed and tested a semi-biased stepwise docking MD where binding process of receptor and ligand was split into multiple steps. Unlike biased umbrella sampling, where conformational dynamics was sampled along several predefined reaction coordinates with simple harmonic potential added to the system’s Hamiltonian [16]. We performed the simulation confined with approximate binding direction and surmised that this unbiased MD performed near the binding interface may simulate a realistic mode of docking. The loop made contacts with the RBD floor and adopted a shape to optimally interact with the surrounding residues according to the surface contours of the site. This essentially provided an adequate sampling time for inducing the antibody to fit into the antigenic binding site. In our example, we submitted the whole 13D4f:HA model for SDMD simulation using an automated process, and obtained the best simulation result, with the lowest RMSD value (0.926) and the highest recovery of interactions (Table 1).

Previous work shows conformational rearrangement in bnAb H5.3 against HA, which undergoes a ∼90° twist in the HCDR3 when binding, as well as a ∼5-Å shift in the LCDR1 [17]. Our simulation method could be used to explain why a broad-spectrum antibody like bnAb H5.3 tended to form a “clever” HCDR3 capable of adapting to the shape of multiple types of antigens. Such knowledge of this and other intrinsic allosteric mechanisms may help to enhance peptide or biological drug design for RBS-oriented antiviral prophylaxis and/or therapeutics against influenza. To expand on this, we retrieved all the complex structure items from the PDB and screened for those that show discrepancies between the bound and unbound states. We found that antibody–antigen complexes comprise a large proportion of the screened items (data not shown). Thus, our method might pave the way for research into structurally dependent biological function of antibody-antigen interaction.

To check if traditional docking methods could provide a reasonable prediction of 13D4-HA docking, we used more docking methods to find the docked conformation, including rigid docking algorithm-based methods (ZDock [18], PatchDock [19], GRAMM [20], HDock [21]) and even flexible docking algorithm-based methods (SwarmDock [22], FiberDock [23], FireDock [24], HadDock [25]). Interestingly, all the methods could find roughly binding sites of 13D4 at HA RBM, only FireDock and HadDock provided a similar interaction orientation with that of crystal structure, however none of them can recapitulate the changed conformation of HCDR3 loop as our SDMD does (Supplemental Figure S4). As compared to these traditional docking methods, and MD-based docking methods such as Relaxed Complex Scheme (RCS) [26], Dynamic Molecular Docking (DMD) [27], Targeted molecular dynamics (tMD) [28] and Steered molecular dynamics (SMD), our proposed stepwise docking molecular dynamics (SDMD) allows the conformational extended HCDR3 loop to match the surface contour of HA receptor binding motif site by two-fold strategy, one is the stepwise approach preventing the potential clash, the other is the very sufficient MD per step driving the docking with a gradually but substantially conformational change to recapitulate the real interaction scenario as observed in the 13D4:HA co-crystal structure.

In terms of the success in recapitulating the conformational rearrangement in the study, the parameter setting for SDMD could be elaborated with four aspects of consideration when applied to other cases. Firstly, the binding orientation should be set as the exact opposite against the separating direction to avoid any unfavorable contact between long flexible loop and other unintended regions at the receptor (Supplemental Figure S5). The docking route could be determined with the suggested approach by Kang et al. [29]. Secondly, the distance separated away should be minimal while avoid from potential clash in terms of Lennard-Jones and short-range electrostatic interactions. Regarding to step size, more steps for SDMD may allow a more precise sampling of the conformational changes, but would tremendously increase the computational time. Thirdly, the regions distal from interface could indirectly affect structural changes by global influence during MD, and should be fixed with reasonable force strength. In this regard, pre-simulation using CMD method could be conducted for receptor and ligand individually to estimate the degree of structural flexibility of protein. Relatively rigid regions are more suitable to be fixed. Finally, the endpoint criteria for the system equilibrium could be optimized to avoid insufficient equilibrium or pointless longer calculation. We anticipate that the method could be further improved as a high-efficient MD-based docking method and be applied in simulating intricate allosteric mechanisms among a range of biological functions.

Conclusion

Our work demonstrates the utility of a stepwise docking simulation method to provide a realistic representation of the binding process and offers a final model that most closely resembles the crystal structure in terms of RMSD and key interactions. Binding between two biological macromolecules will inevitably involve some degree of conformational rearrangement at the interface, which poses difficulties for researchers. We propose that our method will help to improve our understanding of allosteric mechanism through more accurate simulations [2], [30]. In particular, we anticipate that this simulation approach can be used in antibody engineering and biological drug design. Furthermore, if possible, the algorithm and workflow could be optimized to identify highly efficient, flexible protein–protein docking methods.

Materials and methods

Structures acquisition and preprocessing

The crystal structures of individual 13D4 Fab (13D4f) and 13D4:HA complex were checked for missing atoms, and hydrogen atoms were added by Discovery Studio 2017 R2 (DS) software. Water molecules and 2-acetamido-2-deoxy-beta-D-glucopyranose (NAG) bound to HA were removed. The structure of the 1:1 complex of 13D4f:HA was prepared by superimposing the free-state 13D4 Fab with the 13D4:HA complex.

Simulation system setup

The Autopsf plugin in the Visual Molecular Dynamics software (VMD, version 1.9.2) [31] provides a streamlined process for generating a dynamics-ready protein atom coordinate file (pdb) and a protein structure file (psf). The prepared structure was submitted to Autopsf with the default settings and the CHARMM27 force field topology files. To keep the disulfide bridges stable, patches were implemented between pairs of cysteine residues. The structure was then embedded in the explicit solvent (water) box encompassing 10 Å from the protein boundary using the TIP3P (transferable intermolecular potential with three interaction sites) water potential model [32]. The system was neutralized by adding counter chloride ions to achieve zero charge followed by additional sodium and chloride ions to a final physiological concentration of 0.15 m. We applied CHARMM27 force field in our simulation, and periodic boundary conditions were employed to avoid an edge effect.

Molecular dynamics simulation

All MD simulations were performed using the NAMD version 2.12 MD package [33] compiled with cuda support under an InfiniBand-based cluster. In the default simulation conditions, the integration timestep of the simulation was set to 2 fs and the position coordinates (DCD file) were saved every 4 ps for analysis. The simulations were performed at 310 K and constant temperature was controlled by Langevin dynamics [34] under a pressure of 1 atm [35] kept constant using the Nose-Hoover thermostat. Long-range periodic electrostatic interactions were evaluated using the smooth Particle-Mesh Ewald (PME) [36] method, with a real cut-off radius of 10 Å. Lengths of all chemical bonds involving hydrogen bonds were constrained by the SHAKE algorithm [37]. The simulation was under NPT [33] ensemble conditions, generated by the Langevin equation. Before production, the system was energy-minimized by 2000 conjugate gradients steps to reduce steric conflicts between water molecules and the protein. The dynamic results were analyzed using the VMD program.

Conventional molecular dynamics (CMD) simulation

The 13D4 Fab in free state was simulated with default conditions. For the antibody-antigen complex, two independent systems were employed. One contained the aligned 13D4f:HA complex that underwent a 200-step energy minimization. Optimization was performed with no constraints in a GBSW (Generalize Born With a simple Switching) implicit solvent model by DS software and released almost all of the clashes between the aligned 13D4 Fab and HA. In the other system, the complex structure was separated by pulling the aligned 13D4 Fab away from the binding interface with a distance of 5 Å using PyMol [38], along the direction passing the HA COM. Both of the two systems were simulated by CMD with initial simulation conditions.

Steered molecular dynamics (SMD) simulation

SMD simulations were employed to accelerate the simulation procedure while ensuring the binding between the antibody and HA. To achieve this, the 13D4 HCDR3 loop was pulled away from HA at a distance of 5 Å, along the direction passing the HA COM. The distance between the antibody and HA—defined as the distance between the Val104:Cα on the forced loop and Trp153:Cα on HA—served as a reference for measuring the driving effect of the force. The system was submitted to SMD simulation with the constant temperature control switched off to decrease the freedom of the atoms as much as possible; the other parameters remained the same as those in the CMD. A constant force was then applied to atoms in the loop in a reverse direction and removed when the distance reached its minimum. To find a rational force for simulation, multiple independent SMD processes with various values of applied force were carried out, divided into three stages: stage I (10 pN, 100 pN, 1,000 pN and 10,000 pN), stage II (2,000 pN, 4,000 pN, 6,000 pN and 8,000 pN), and stage III (4,500 pN, 5,000 pN and 5,500 pN). The force was released and an additional 100 ns CMD simulation was employed after the SMD simulation reached an equilibrium while the indicator distance remaining relatively steady.

Stepwise-docking molecular dynamics (SDMD) simulation

The principle behind the SDMD simulation is differentiating the recognition of the receptor and ligand protein. The ligand protein was first pulled away from the receptor along the possible binding direction using VMD software [31]. (NB: Initially, the direction could be perpendicular to the binding interface). The region unrelated to the interface was constrained in subsequent calculations. The structural constraints confined the model from undesired swing or movement; the less the residues were restricted, the smaller the interference on the sampling of conformation during simulation. Approaching simulations was completed using a series of discrete but stepwise CMDs, with the ligand artificially moving closer to the receptor in a step-by-step fashion. In each step, when the current simulated system came to an equilibrium, a small shift was applied and the next-step simulation was initiated.

For this article, the system containing the whole 13D4 Fab and HA was submitted to SDMD simulation using an automatic procedure, which we referred to as “automatic stepwise docking molecular dynamics” (ASDMD) (Supplementary Figure S1). The separation of 13D4f and HA was done along the direction parallel to the 3-fold axis of HA, with a transition height of 4 Å. Regions far from the interface were fixed (Supplemental Table S1). The parameters setting in the configuration file was the same for all simulations as that in the aforementioned CMD. When the system reached an equilibrium, the structure with the lower energy was extracted and the distance was decreased by 0.1 Å for the next step run. After 40 steps, the simulated loop reached its binding site and the final conformation was recorded. A flowchart of the whole procedure is summarized (Fig. 5, S1). To our experience, the upper limit of the run time was set to 40 ns to prevent the simulation from becoming trapped in an unstable state.

Supplementary figure 1

During the runtime, equilibrium judgement was executed via an algorithm, as described in the next section. The decline of the ligand was calculated using a 4 × 4 matrix in VMD,where x , y, and z determine the coordinates of the vector used to set the direction of movement of the ligand protein.

Assessment of equilibrium during MD

Structural variation during dynamics was calculated in terms of the differences between the simulating structure and the initial model. The root mean square deviation (RMSD) of the alpha carbons (Cα) of the structure during the simulation were monitored. When the trajectory of the RMSD value exhibited a smooth and stable curve and the average value remained constant, this meant that the system had converted to a relative state of equilibrium.

Equilibrium in the simulation was estimated with a 2-step verification process. For the first step, the total energy of the system, Et, in the last 5 ns was collected to develop one-dimensional linear regression equations. When the slope of the line, defined as b1, was less than 2E-05, the system was regarded as being tentatively in a metastable state. The second step incorporated the RMSD values from the trajectory of the last 5 ns and the slope of the linear equation was defined as b2. When the slope dropped to 3E-05, the conformation was regarded as stable. The data were then evenly divided into 2 groups and 5 groups, respectively, for a deeper evaluation. When split in half, the data were used to judge for the presence of grossly asymmetric peaks and troughs in the two data sets, with each representing the independent lower slope; i.e., ba and bb. When split among 5 groups, the data were used to calculate the arithmetic mean and SD for each group (From Xavg1 to Xavg5) to ensure that the global RMSD value was maintained on average. The system was regarded as stable when a combination of the following conditions was satisfied: the b2 value was less than 3E-05; the ba and bb values were both less than 6E-05; the SD was less than 0.01. These empirical values were tested and confirmed using the RMSD data sets in the previous SDMD simulation.

Methods for result analysis

Dynamic allostery was evaluated using the non-bonded interaction and the corresponding binding energy between the 13D4 Fab and HA, the deflection angle, and the RMSD value of the HCDR3 with reference to the co-crystal structure.

Interaction analysis

Interactions between 13D4f and HA were identified using the DS Analyze Protein Interface module with default parameters. The maximum distance cutoffs for hydrogen bonds and salt bridges were set to 3.4 Å and 4.0 Å, respectively. The tool allows for calculations of the interaction between two group of atoms in a given structure or specified trajectory. Hydrogen bonding interactions throughout the trajectory were calculated using VMD HBonds Plugin version 1.2 with default parameters.

Measurement of deflection angle and RMSD

The deflection angle of HCDR3 was measured using PyMOL. For each simulated model extracted from the trajectory, the structure was aligned to the 13D4:HA co-crystal structure to calculate the RMSD value. Three atoms in HCDR3 were selected as points for defining the deflection angle: Arg106:Cα, Val104:Cα from the crystal complex, and Val104:Cα from the simulated complex.

Molecular simulation docking

Several docking methods was used to simulate the 13D4 Fab docking to HA (Supplemental Table S2). All methods were run on their respective servers with default parameters. Partial binding site residues were specified to narrow down the search to the actual binding pocket. For FiberDock and FireDock methods, they process a flexible refinement and scoring of the rigid-docking candidates from PatchDock. Finally, the binding model with the highest score for each docking result was selected for comparison.

Data accessibility

Supporting Information is available from the Wiley Online Library or from the author.

Authors’ contributions

Yang Huang and Zizhen Li contributed to formulating the idea of the research and performed simulations and analysis. All authors contributed to the writing of the paper, gave final approval for publication and agreed to be held accountable for the work performed herein.

Funding

This work was supported by grants from the National Natural Science Foundation (Grant No. U1705283), Xiamen Major Scientific and Technological Special Project (Grant No. 3502Z20193009, 3502Z20203023) and CAMS Innovation Fund for Medical Sciences (Grant No. 2019RU022).

CRediT authorship contribution statement

Yang Huang: Conceptualization, Methodology, Writing – original draft. Zizhen Li: Writing – original draft, Methodology, Data curation. Qiyang Hong: Data curation. Lizhi Zhou: Data curation. Yue Ma: Visualization. Yisha Hu: Visualization. Jiabao Xin: Formal analysis. Tingting Li: Formal analysis. Zhibo Kong: . Qingbing Zheng: Visualization, Funding acquisition. Yixin Chen: Visualization. Qinjian Zhao: Visualization. Ying Gu: Supervision. Jun Zhang: Resources, Visualization. Yingbin Wang: Visualization. Hai Yu: Conceptualization. Shaowei Li: Conceptualization, Writing – review & editing, Project administration, Resources, Funding acquisition. Ningshao Xia: Resources, Supervision, Funding acquisition.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.