To explore how the gate-opening mechanism occurs when an intermediate bond forms in the residence of dsDNA, the two experimental coordinates of dTopo-IA and dTopo-I, obtained via X-ray crystallography, were taken from the Protein Data Bank (PDB ID: 1MW8 and 3PX7) (Perry and Mondragón, 2003; Zhang et al., 2011). Missing residues in the protein part were modeled using the Chimera MODELLER program, and the best-scoring models of 3PX7.pdb (zDOPE = −0.98) and 1MW8.pdb (zDOPE = −1.23) were selected for further simulations (Šali and Blundell, 1993) (Supplementary Tables S1A,B). The template of both enzymes is given in Supplementary Table S2A. The DNA sequence was converted into double-stranded DNA using UCSF Chimera Tools (Couch et al., 2006). To convert active-site tyrosine into phosphorylated tyrosine (PTR) in dTopo-IA, the necessary parameter for this residue (Supplementary Table S2B) was generated with the Gaussian 09 package (Frisch, 1995), using the density functional theory (DFT) method with the B3LYP/6-31++G (2d, 2p) level of theory. For the artificial cleavage in the double-stranded DNA, one strand of DNA was broken and the corresponding nucleotide was modified as thymine (MDT5) (Pavlin et al., 2023; Tiwari, 2018) (MDT5.lib file given in Supplementary Table S2C). A bond between PTR residue and the DNA nucleotide was made using the xleap bond command. The dTopo-IA complex with covalent intermediate and dTopo-I without covalent bond was then taken to compare the action in the enzyme with and without the intermediate covalent bond (Bates et al., 2011) (Figures 2A,B).The entire modified complex molecule was input into the xleap program of the AMBER18 package. The protein was treated with the ff14SB force field (Raguette et al., 2020), and the DNA with DNA.OL15 (Raguette et al., 2024). For covalent bond phosphorylated tyrosine, leaprc.phosaa10 of the Generalized Amber Force Field (GAFF2) was used (Wang et al., 2004). Partial atomic charges and missing parameters for PTR were obtained with the RESP charge fitting method (Wang et al., 2004) at the B3LYP/6-31++G (2d, 2p) level of theory (Bayly et al., 1993; Cornell et al., 2002) using an optimized geometry methodology level of theory with the Gaussian 09 program (Becke, 1988).

Since the covalent PTR linkage is not parameterized in standard AMBER force fields, custom charge derivation was required. Although ff14SB and DNA.OL15 are historically calibrated using RESP charges at the HF/6-31G level, this approach can over-polarize highly charged phosphate groups. Therefore, RESP charges derived from B3LYP/6-31++G (2d, 2p) electrostatic potentials were employed to obtain a balanced description of phosphate charge distribution and P–O bond polarity relevant to the covalent intermediate while maintaining compatibility with the non-polarizable AMBER framework through RESP charge restraints.

To neutralize the protein–DNA complex, 14 Na⁺ ions in dTopo-I and 19 Na⁺ ions in dTopo-IA were added according to their charges. The whole system was solvated in an octahedral box of the TIP3P model, extended to a minimum cutoff of 12.0 Å from the protein boundary (Dubey and Shaik, 2019).

After generating the required coordinate and topology files, minimization was performed up to 10,000 steps in two parts using the steepest descent (5000 steps) and coupled gradient (5,000 steps) methods (Muralidhar and Tiwari, 2024; Sahoo et al., 2021). In the first case, the complex molecules were restrained at a weight of 500 kcal/mol-Ų, while the water and neutralizing ions were free (Xiong et al., 2008). In the second step, the entire system was minimized without any restraint.

The system was gradually heated to reach 300 K under the NVT ensemble over a period of 50 ps (Dubey et al., 2013). This step was followed by a 1 ns density equilibration under the NPT ensemble, maintaining a temperature of 300 K and a pressure of 1 atm, regulated using a Langevin thermostat (Ka et al., 2005) and a Berendsen barostat (Izaguirre et al., 2001), with a collision frequency of 2 ps and pressure relaxation time of 1 ps (Berendsen et al., 1984). Subsequently, a 3 ns equilibration phase without restraints was performed prior to the 200 ns molecular dynamics (MD) production run. The production run utilized a Monte Carlo barostat (Case et al., 2005) for all complex systems. Long-range electrostatic interactions were treated using the particle mesh Ewald (PME) method (Åqvist et al., 2004) with a cutoff distance of 10 Å, and the SHAKE (Darden et al., 1993) algorithm was applied to constrain bonds involving hydrogen atoms. All MD simulations were executed using the AMBER18 package (Ryckaert et al., 1977).

All MD simulations were performed using the CPU version of AMBER18 version. For each system, three independent replicas of 200 ns were generated for each system, and convergence of key structural and dynamical metrics (RMSD, RMSF, hydrogen-bonding, and PCA) was observed within the simulated timescale.

Trajectory analyses were performed using the CPPTRAJ module (Salomon-Ferrer et al., 2013), which is available in AMBER Tools. The analysis included calculating various essential parameters, such as root mean square deviation (RMSD), root mean square fluctuation (RMSF), number of hydrogen bonds (Hbond), hydrogen bond distance variation, radius of gyration (RoG), and solvent-accessible surface area (SASA). To visualize and render the figure, VMD and PyMOL were employed (Roe et al., 2013). We utilized VMD’s visualization capabilities to study the system’s hydrogen bond dynamics in addition to its structural properties (Humphrey et al., 1996). This required computing polar bond formation and occupancies using stringent criteria, such as a donor–acceptor distance of less than 3.0 Å and an angle cutoff of 20°. The results demonstrated the critical intermolecular interactions that underlie the system’s behavior, thus exposing the intricate network of the hydrogen bonds that support its stability and functionality. Furthermore, we explored the system’s molecular architecture by examining decomposition-free energy per residue, providing deeper insights into the conformational transitions and structural transformations that were occurring.

After completing a 200 ns molecular dynamics (MD) simulation, an in-depth evaluation of the structural stability and dynamic flexibility of the protein-DNA complex was performed using several key parameters. The RMSD of residues 1 to 621, including the DNA strands, was meticulously monitored, revealing valuable insights into the conformational stability of the system. Additionally, the RMSF, RoG, and SASA were calculated using the CPPTRAJ module of the AMBER package, providing a thorough analysis of the system’s stability and flexibility throughout the simulation.

From the RMSD calculation for enzymes residues as well as dsDNA residues, we found that all systems exhibited an initial equilibration phase during the first 20–30 ns, followed by a clear stabilization period for the remainder of the simulation.

For the dTopo-IA complex, the enzyme backbone RMSD (black line) gradually increased during equilibration and stabilized at approximately 3.5–4.0 Å, indicating that the overall fold of the enzyme remained stable throughout the 200 ns trajectory. The corresponding DNA RMSD (red line) remained lower, fluctuating between 2.0 and 2.8 Å, suggesting that the DNA duplex experienced only moderate conformational adjustments while maintaining structural integrity. After 30 ns, the RMSD of the complex molecule stabilized, indicating that after the formation of the PTR intermediate bond, the DNA strands and topoisomerase enzyme molecules adjusted to a stable orientation, with a minimum deviation of approximately 3.2 Å, while for the system without cleavage, this deviation was approximately 4.5 Å. A detailed inspection of the DNA RMSD indicated that its variation was mainly due to the contribution of the double-stranded DNA bases (Figure 3). The RMSD of the system with covalent complex oscillated at approximately 3.2 Å. In contrast, the DNA bases showed a relatively high RMSD (nearly 3.0 Å), suggesting that the deviations observed in the crystal study were consistent with our findings. The corresponding snapshot is shown in Supplementary Figure S1A-C.

Both enzyme systems attained stable conformations overall, according to the RMSD analysis; however, the dTopo-IA complex showed significantly lower RMSD values for both protein and DNA, which is consistent with improved stabilization after PTR formation. The non-cleavage (dTopo-I) system’s higher RMSD indicates that the DNA–enzyme interface is less rigid and more flexible.

To ensure the stability and reproducibility of the trajectories, three independent MD simulations were performed for each system and their corresponding RMSD plots (Supplementary Figure S2A–D). Figure 3 depicts a representative MD replica chosen for early convergence to illustrate the conformational behavior.

The RMSF versus residue number of both complexes is depicted in Figure 4A. The highest fluctuated residue ranges of different domains were from domains D-I (27–126), D-II (205–305), and D-III (306–500), while D-IV (501-621) remained stable.

The highest fluctuations were found in the DNA residues (1–26) of dTopo-IA. These indicate that during gate opening, when the active domain III (306–430) allows the passage of DNA T-strands, these residues experience repulsion, resulting in such fluctuations. This phenomenon originates from the dynamic motion of domain I (residues 27–179) within the complex dTopo-IA (Supplementary Figure S3). Upon domains II and III aligning, domain I exhibits a downward displacement of approximately 4.5 Å (black dotted line), followed by a backward shift of approximately 5.5 Å (black dotted line) (Figures 4B,C). This movement closely correlates with structural observations from studies (Bugreev and Nevinsky, 2009). However, the study protocols elucidate the specific conformations of key residues such as Pro-255, Ser-256, Ala-259, and Glu-296, which may contribute to this domain mobility (Schoeffler and Berger, 2008; Gunn et al., 2017). The conformational interactions of the residues are illustrated in Figures 4D,E. For more clarity, the domain-wise deviations and corresponding snapshots are shown in Supplementary Figures S4A–D.

Figure 5A provides information about the formation of a number of hydrogen bonds during the simulation time. Figure 5B depicts the variation in donor–acceptor distances vs. simulation time. Initially, the distance between the O3P of PTR-345 and the NH₂ group of Arg-347 extends up to 8.5 Å, and the distance involving the OH of PTR-345 and NH₂ of Arg-347 reaches nearly 8.0 Å. After 10 ns, when the enzyme–dsDNA complex became stable and the gate closes, PTR-345 and Arg-347 moved closer, and their interaction distance stabilized at 3.7 Å, forming weak donor–acceptor interactions. Concurrently, the weak hydrogen bond between the oxygen of PTR-345 and the nitrogen of Thr-348 stabilized at approximately 2.6 Å. A notable shift occurred as the hydrogen bond between the OH of PTR-345 and the NH2 of Arg-347 contracted to 2.5 Å. Throughout the simulation, the weak hydrogen bond between the oxygen in PTR-345 and the nitrogen in Thr-348 remained steady at 2.8 Å. However, after 170 ns, the donor–acceptor distance between O3P in PTR-345 and NH2 of Arg-347 re-stabilized at 4.5 Å, exceeding the typical hydrogen bond cutoff; this should be interpreted as a donor–acceptor contact rather than a hydrogen bond. The molecular interaction of PTR-345 to other molecules is shown in Figure 5C. The highly favorable hydrogen bond occurrence for the two complexes is given in Supplementary Tables S3A,B.

In our MD simulation of the protein–DNA complex involving dTopo-IA, the RoG served as a crucial metric for understanding conformational dynamics. The plot of RoG against simulation time showcased intriguing dynamics (Figure 6A). The RoG exhibited significant fluctuations, reaching up to 34.0 Å for the complex without PTR residue, while for the system with PTR residues with covalent bonds, the RoG value was 31.5 Å. After 50 ns, at the average value (32.3±5.1) Å of dTopo-I while at average values (61.9±11.3) of dTopo-IA, both complexes became stable. This variability is attributed to the collective interaction of non-hydrogen atoms aligning along rotational axes. After this period, RoG stabilizes, indicating the system relaxing under the applied conditions. These observations underscore the complex interplay of molecular forces that govern conformational dynamics in the simulated environment, providing valuable insights into the stability and behavior of the protein–DNA complex (Dasgupta et al., 2020; Liu et al., 2017).

The SASA measures the surface area of an enzyme structure accessible to solvent molecules. This metric provides insights into conformational changes and protein–DNA interactions. Figure 6B illustrates SASA versus simulation time for both complexes. Here the average values for the complexes first were 365.8±23.5) nm² while, for the other, it was 304.7±17.3) nm². The stabilization of the accessible area is indicative of conformational adjustments in the complex, which may create a more favorable environment for DNA relegation, as shown in Figure 6C taken as a snapshot during trajectory analysis. However, SASA itself does not drive this process.

The MM/GBSA calculation provides the binding free energy of the receptor ligand (DNA). The highly negative values of binding free energy between the receptor DNA provide the notion of stronger interactions. Here we calculated a binding free energy molecule for two different simulating systems for last 200 frames (after 175 ns) for the generalized Born model, with igb = 5. In the first simulating system without intermediate bond, the value of binding free energy was −442.1 kcal/mol. Furthermore, in the case of dsDNA with the dTopo-IA enzyme, the binding free energy was −108.6 kcal/mol. This loss of binding energy (B.E) of the complex with the covalent intermediate bond indicates that after the intermediate bond formation, the T-strand of dsDNA became unstable, and thus because of this high energy, the strand passage of T-strands occurs.

The molecular mechanics energy function is a pair-wise additive function that can be broken into per-residue components. Per-residue binding energy decomposition analysis revealed the contribution of various amino acids toward total binding energy (Tiwari et al., 2024). From the decomposition analysis (Figure 7), the contributions of the consistently interacting amino acids were extracted, and their correlation with experimental activity was deduced. Arg-519, Gln-223, Arg-533, Ala-216, Lys-212, Arg-90, Thr-65, Arg-96, Val-201, Gly-202, PTR-345, Lys-41, and Lys-117 showed contributions of energy (in kcal/mol) of −15.9, −17.2, −5.7, −5.9, −5.6, −7.8, −13.5, −11.3, −11.1, −6.8, −24.0, −6.5, and −5.1, respectively (Supplementary Table S6). The energies of Arg-519, Gln-223, Arg-533, Lys-212, Val-201, and Gly-202 provide favorable interactions that show that their energy contributions correlate with stabilizing enzyme−DNA interactions, facilitating gate opening and the passage of T-strands of DNA helix. Among the amino acids, Gln-223, Arg-533, Val-201, and Lys-212 interact with the DNA strands for gate opening. The DNA interaction energy map with amino acids of dTopo-IA enzymes is shown in Figure 7. The remnant MM/GBSA is used to compute the interaction map, which is quite helpful because it identifies non-polar and hydrophobic interactions. It is impossible to examine these interactions in the same way as for salt bridges and hydrogen bonds. Among the amino acids, Gln-223, Arg-533, Val-201, and Lys-212 interact with the DNA strands for gate opening.

The DNA strand interacts favorably with Val-201, Glu-223, and Arg-533 through hydrophobic interactions in the presence of an aqueous environment but with PTR-345 covalently bonded. These non-polar and hydrophobic residues make a big difference in the entire positive shift in enthalpy that occurs during the interaction between enzyme residues and DNA interacting residues (Genheden and Ryde, 2015) (Figure 8).

To find the global motions of the systems, principal component analysis (PCA) was calculated on the MD trajectories of dTopo-I and dTopo-IA. The resulting projections along the first two principal components (PC1 and PC2) (Figures 9A,B) reveal that the covalent intermediate complex (dTopo-IA) exhibits broader conformational sampling, indicating enhanced structural flexibility compared to the non-covalent dTopo-I system.