Catalytic Mechanism of ATP Hydrolysis in the ATPase Domain of Human DNA Topoisomerase IIα.

Human DNA topoisomerase IIα is a biological nanomachine that regulates the topological changes of the DNA molecule and is considered a prime target for anticancer drugs. Despite intensive research, many atomic details about its mechanism of action remain unknown. We investigated the ATPase domain, a segment of the human DNA topoisomerase IIα, using all-atom molecular simulations, multiscale quantum mechanics/molecular mechanics (QM/MM) calculations, and a point mutation study. The results suggested that the binding of ATP affects the overall dynamics of the ATPase dimer. Reaction modeling revealed that ATP hydrolysis favors the dissociative substrate-assisted reaction mechanism with the catalytic Glu87 serving to properly position and polarize the lytic water molecule. The point mutation study complemented our computational results, demonstrating that Lys378, part of the important QTK loop, acts as a stabilizing residue. The work aims to pave the way to a deeper understanding of these important molecular motors and to advance the development of new therapeutics.

Introduction

DNA topoisomerases encompass a family of complex biological nanomachines that can catalyze the induction of topological changes in the DNA molecules.[1,2] By modulating DNA topology, they play an essential role in several crucial processes in the cell such as transcription, replication, and chromosome segregation.[3−5] They are divided into two groups: type I topoisomerases, which cleave a single strand of the double-stranded DNA, and type II topoisomerases, which cleave both strands of the DNA.[6,7] Both types of DNA topoisomerases are essential for DNA replication and transcription in all living cells as they maintain steady-state distributions of topological states.[8]

The human DNA topoisomerase IIα (htIIα) is a stable homodimer protein[9] and its structure can be divided into three separate domains.[10] The N-terminal or the ATPase domain is a part of the GHKL (Gyrase, Hsp90, histidinKinase, MutL) family of proteins and binds the T-segment of the DNA (Figure ).[11,12] The central domain binds and cleaves the DNA G-segment. Finally, the C-terminal section contains residues important for phosphorylation as well as interactions with other proteins.[6,7] The enzyme also undergoes acetylation, with several possible modification sites on the ATPase and central domains; in particular, the eukaryotic-specific Lys168 in the ATP binding pocket was found to be key for the dimerization of the ATPase domain.[13] The proposed catalytic cycle of a type II topoisomerase is characterized by a two-gate mechanism, where the homodimer shows remarkable flexibility as the molecular gates undergo substantial opening and closing motions.[14−17]

Figure 1

Human topoisomerase IIα ATPase dimer (PDB: 1ZXM), (Left): Green—GHKL domain of chain A, red—transducer domain of chain A, orange—chain B, violet—QTK loop, and blue—ATP. Right: Zoom in on the ATP binding site with the bound AMP–PNP ligand.

The role of the ATP molecule in type II topoisomerases has been a subject of intensive research.[14,18−21] Biochemical studies revealed that type II topoisomerases hydrolyze two ATP molecules per reaction cycle. The role of the first ATP molecule is assumed to support a unidirectional transfer of the T-segment of DNA from its binding at the N-terminal domain to its exit via the C-gate of the enzyme after the G-segment is cleaved at the DNA gate.[22] The role of the second ATP hydrolysis is even more ambiguous and some research suggests to be linked with the communication between the ATPase and the cleavage/ligation domains to support catalysis by resetting the enzyme conformation.[23]

The studies of the reaction mechanism of ATP hydrolysis in type II topoisomerases show that after binding of two ATP molecules, the enzyme stochastically hydrolyzes the first ATP before hydrolyzing the second ATP molecule. More precisely, following the hydrolysis of the first ATP, at least one slow step in the mechanism occurs before products of the second ATP hydrolysis can be detected. The sequential hydrolysis of ATP and product release steps indicate the complexity via which type II topoisomerase operate.[14] In human topoisomerase IIα, it was further shown that only one ATP molecule is required for the enzyme to fully complete its catalytic cycle, while nonhydrolyzable ATP analogues allow the enzyme to complete only a single turnover of DNA transport and leave it unable to reset for a new cycle.[14,23] ATP hydrolysis and ADP dissociation are important for rapid kinetics and opening of the N-gate, allowing the enzyme to enter into a succeeding catalytic cycle.[14,24] Biochemical investigations also studied the role of the QTK loop (376QTK378), which extends from its transducer domain to the ATP-binding pocket in the GHKL domain, and concluded it to be an important feature for interdomain communication (Figure ). It keeps the N-terminal gate open in the absence of nucleotides and facilitates the N-gate dimerization once a nucleotide binds. Deletion causes increased baseline cleavage levels but sharply reduces strand passage and virtually abolishes stimulation of cleavage by nucleotides. Furthermore, while the low baseline ATPase activity is retained, DNA stimulation of ATPase activity is abolished. Lastly, the operation of the N-gate becomes disordered; it can trap DNA not only in the presence of non-hydrolysable ATP analogues but also in the absence of nucleotides.[25]

Many structural studies investigated the mechanism of type II topoisomerases and predominantly focused on their DNA-binding/cleavage core of the enzyme.[26−28] The ATPase domain of htIIα was also studied separately, revealing a rigid-body movement of the structural modules within this domain.[29] Recently, an important breakthrough occurred with the determination of two conformations of htIIα bound to DNA trapped by the topo II poison etoposide. The study provided valuable information on how the ATPase domain is spatially connected to the DNA-binding/cleavage domain conformations and that the bound ATP analogue enables the rotation of the N-gate and opening of the DNA gate.[10]

The ATPase domain of htIIα is also a promising new target site for an emerging class of potential anticancer drugs, called catalytic topo II inhibitors.[30,31] These molecules include a class that targets the ATP binding site of htIIα and predominantly mimics the adenine moiety of ATP, with a variety of scaffolds already discovered and validated.[32−35] Such inhibition of topo II would, in principle, not result in an induction of excessive DNA damage, characterized by the DNA double strand breaks that are one of the main triggers of severe adverse effects observed when administering topo II poisons (e.g., etoposide and doxorubicin). These harmful effects comprise incidences of secondary malignancies after chemotherapy and cardiotoxicity.[36−38] Furthermore, the increased occurrence drug resistance to clinically used topo II poisons further fuels the need for revisiting this established anticancer target.[39]

With the aim to unravel some of the intricacies linked with the mechanism of action type II topoisomerases, we focused on the ATPase dimer of htIIα, and employed all-atom molecular simulations, multiscale quantum mechanics (QM)/molecular mechanics (MM) reaction modeling, and a point mutation study. The collected data offer new atomistic insights into this fascinating biological nanomachine and could support the development of new cancer therapeutics targeting this domain.

Results and Discussion

Dynamics of Human Topo IIα ATPase Dimer

To study the dynamics of the of the ATPase domain dimer and assess the influence of the ATP ligand, we performed molecular dynamics (MD) simulations of the ATPase topo II dimer with bound ATP (holo system), and its apo form where both ATPs were removed. Initially, we assessed the metrics connected with the overall system stability during the simulations. RMSD values for the apo dimer and the system with bound ATP were 2.9 and 3.4 Å, respectively. Visualizing the trajectories, we observed that, the elevated RMSD in the holo structure can be attributed to the opening of the transducer domain, which is not as prevalent in the apo structure. Interestingly, the resulting conformation of the holo system resembles the solved topoisomerase II-ADP ATPase dimer complex structure (PDB: 1ZXN); that is, the likely structure after hydrolysis of the first ATP molecule has occurred.[29] (Supporting Information, Sections S1 and S2).

The bound ATP molecules were stable in their binding sites with the RMSD of about 0.8 Å in both protomers. This was further reflected by the stability of other bonds than anchor the ATP to its binding pocket. The bond between the amide oxygen of Asn120 and the amino group of the ATP averaged 3 Å. Similarly, the average distance of the H-bond between the side chain nitrogen of Lys168 and the oxygen on the α-phosphate was stable at 2.7 Å. The Mg2+ ion also showed little spatial movement; RMSD was 0.5 Å. Additionally, we measured the interactions observed in the X-ray structure between the Mg2+ ion and its interacting partners—ATP, water molecules, and relevant amino acid residues determining that they all remained stable and within the observed distances completing the octahedral coordination of this ion (Table S1). With the crude system stability established, we proceeded to focus on catalysis-relevant parts of the htIIα ATPase domain in more detail.

First, we examined the Glu87 residue that was shown to be an essential residue for the catalysis.[40] In the crystal structure, Glu87 interacts with a water molecule that is already optimally positioned for a nucleophilic attack on the γ-phosphate of the ATP. During the simulation, the candidate for the catalytic water molecule was firmly nested in the ATP binding site and present throughout the simulation. The catalytic Glu87 also interacted with the ATP molecule. We monitored the bond lengths between the side chain oxygens of the Glu87 carboxylic group and the oxygen on the γ-phosphates. We noticed that in chain A, the carboxylic side chain group of Glu87 exhibits rotations, alternating which of its oxygen atoms was closer (∼3.5 Å) and which farther (∼5.5 Å) from the γ-phosphate, whereas in chain B, no such rotation was observed (Figure S2B). This marks the only substantial difference observed between chains A and B during our simulation. Next, we examined the QTK loop, which is completely conserved[41] and is critical for the interdomain communication and potentially also plays a role in the ATP catalytic reaction.[25] In both the holo and the apo system, the loop remains stable in both chains/monomers of the ATPase domain throughout the simulation (RMSD = 2 Å and 1.6 Å). Additional graphs for all analyzed distances and RMSD values can be found in Section S2 of the Supporting Information.

To further quantify the observed dimer dynamics, we performed a principal component analysis (PCA) and calculated per-residue root mean square fluctuation (RMSF) values. The RMSF protein heatmaps in Figure again reflect a higher flexibility of the transducer domain versus the GHKL domain. The different levels of movement between the GHKL and the transducer domain are also evident on the per-residue RMSF graphs. It is interesting to note that the apo system exhibited less mobility in the transducer domain and that after 100 ns, the holo system was closer to the structure of 1ZXN, that is, the supposed structure of topo IIα after hydrolysis has occurred (see Section S1).

Figure 2

RMSF heatmaps and per residue RMSF graphs of the simulated htIIα dimers: (A) the holo system and (B) the apo system without ATP. Blue on the graph represents the A chain and red the B chain of the ATPase dimer. Our model is missing the first 28 residues. The dotted line shows a divide between the GHKL domain (residues 29–264) and the transducer domain (residues 265–428).

In PCA for the holo system, the obtained eigenvectors of the first principal component (PC1) indicate, as expected, the opening/closing of the transducer domain. The vectors of the second principal component (PC2) are almost perpendicular to the first, signifying a twisting motion (Figure A). The system without ATP is less ordered in its movements, reflected both in the vector and per residue graphs (Figure B) and (Animation 1). This is in agreement with the observation that ATP binding promotes concerted movement of the ATPase domain.[18]

Figure 3

(Left) Eigenvectors and per-residue contribution for the first and second principal component for the (A) holo system and (B) apo system. On graphs the dotted line shows a divide between chain A and B. The GHKL domain (residues 29–264) is unshaded, and the transducer domain (residues 265–428) is shaded in green for the holo system. (Right) Residual correlative motions visualized by the dynamic cross-correlation maps for the following: (A) holo system and (B) apo system. The color scale in the matrices goes from blue (for values ranging between −1 and −0.25), through white (−0.25 to 0.25) to red (0.25 to 1). Negative values depict residue pairs having the opposite motion, namely, anti-correlated motions, while positive values indicate pairwise movement in the same direction. Above each matrix, the three-dimensional cartoon model shows the strongest anti-correlated motions in each structure (i.e., blue lines connecting the residue pairs).

We also determined the correlation of the displacements of all residue pairs and plotted the correlations in the dynamic cross-correlation matrices (DCCM) outlined in Figure . With DCCM analysis, we wanted to confirm and additionally characterize the domain dynamics previously studied in PCA. Already at first glance, we can detect differences in the degree of (anti)correlated motions that each system exhibits. Comparing the two matrices of the htIIα dimers, we can observe a presence of clearer patterns in the holo system, whereas the apo system exhibits no larger concentrated groups of either positively or negatively correlated motions. Visualizing the strongest anti-correlations in three-dimensional representation confirmed our observations, as we can detect a clearer cluster of anti-correlations in the holo transducer domain compared to the apo system. All results support the notion that binding of the ATP influences the dynamics of the ATPase dimer.

It should be noted that our truncated topo IIα systems lack the DNA molecule, the absence of which likely influences global protein shifts. Therefore, the conformational motions of the protein cannot be fully captured within the current simulations. More accurate evaluation of the dynamics would require a more comprehensive conformational space sampling of the system[42] and the inclusion of other components of the topo II molecular motor. Therefore, these simulations serve only as an initial assessment of the dimer’s dynamics.

Reaction Mechanism of ATP Hydrolysis

Next, we investigated the catalytic process of ATP hydrolysis on a single ATPase domain. Molecular modeling of such processes requires the application of quantum mechanics.[43,44] A standard multiscale QM/MM approach involves dividing the system into a reaction region described by QM, while the remaining part is treated with classical MM.[45,46] The determination of the course of events in the transformation from the initial (reactants) to the final (products) configuration of the system is a task solved by the reaction pathway methods.[47] One group of methods comprises the restrained coordinate driving (RCD)[47] where the reaction is forced to occur by restraining the system to a predefined chemically relevant reaction coordinate. Unfortunately, RCD methods suffer from hysteresis difficulties due to the constraints imposed on the reaction coordinate.[48] Another approach to determine the reaction pathways that can circumvent this concern is the replica path method (RPATh).[48−51] It is an extension of the self-penalty walk methods,[52,53] which involves a simultaneous optimization of a set of geometries of the system corresponding to a set of points along the reaction path, resulting in an approximate minimum energy pathway (MEP). RPATh removes the bias in choosing the reaction coordinate because the method uses the root mean square distance (RMSD) parameter to define the distance between two points on the approximate MEP. This allows the “global pathway” movement to define a reaction path rather than a combination of distances used by RCD.[54−56]

The literature lists two major mechanisms regarding the ATP hydrolysis.[57] One is a general base catalysis (GBC) where the catalytic base abstracts a proton from the lytic water molecule prior to or at the transition state (TS). An alternative to GBC is the substrate-assisted catalysis (SAC), in which the γ-phosphate group of the ATP molecule acts as the catalytic base for the cleavage of the lytic water.[58,59] The mechanism of ATP hydrolysis has been computationally investigated in several biological systems, including ATP-dependent molecular motors such as myosin, F1-ATPase, and kinesin, proposing diverse reaction paths, especially concerning the lytic water proton transfer.[60−63]

Again, the co-crystal structure of the htIIα ATPase domain with bound AMP–PNP offers an excellent modeling starting point as an approximation of the ATPase state before the hydrolysis takes place (Figure A). It further provides an optimal experimentally determined position of the lytic water molecule for the nucleophilic attack. In our study, we evaluated the GBC and SAC reaction mechanisms of ATP hydrolysis (Figure B). For modeling the QM region, we chose a well-established DFT hybrid functional B3LYP in conjunction with a standard 6-31G* basis set. Such a methodological setup has been shown to work well in previous computational studies of the catalytic mechanisms of ATP hydrolysis of various molecular motors.[62]

Figure 4

(A) ATP binding site on the htIIα and bound ATP molecule generated from the experimentally present AMP–PNP analogue (PDB: 1ZXM). Marked are also distances that we monitored during the modeling of the reaction pathways. (B) Schematics of the SAC and GBC reaction paths of ATP hydrolysis that were investigated in this study.

The QM region included three water molecules near the γ-phosphate, the Mg2+ ion, the methyl triphosphate fragment of ATP, and Glu87 and Lys378, that is, atoms/side chains possibly involved in H+ and OH– transfer (Figure S5). Lys378 was posited as a catalytic base and modeled as protonated, as Lys residues near phosphates have been shown to be protonated, with deprotonation only occurring in highly apolar (hydrophobic) environments.[64−66] We also examined potential roles of Gln376, which is also a part of the QTK loop. Considering the position of the side chain of this residue in relation to other partners, we could construct no feasible way for this residue to participate directly in proton transfer, for example, via a proton wire. Therefore, Gln376 was left in the region described by classic molecular mechanics, where its H-bonding effects on stabilizing ATP could still be fully included in the calculations.

In the obtained SAC reaction mechanism, the energy progression, changes of key distances, and key frames/replicas as determined by the RPATh method are depicted in Figure A (Animation 2). We defined four key distances to describe the reaction course. The R1 distance represents the breakage of the bond between the oxygen on β-phosphate (Pβ) and γ-phosphate (Pγ), R2 monitors the transfer of the OH– from the lytic water to the γ-phosphate. Finally, to describe the proton transfer of the Hw atom from the OW oxygen of the lytic water to the OE γ-phosphate oxygen, we measured RO and RO distances.

Figure 5

RPATh-determined reaction pathways with monitored key distances and reaction frames/replicas. Dashed lines denote coordination bonds with Mg2+ (green sphere) or formed hydrogen bonds. Distances are labeled in angstroms (Å). (A) SAC with key replicas 2, 9, and 16 of the reaction pathway and monitored R1, R2, RO and RO distances. (B) GBC with key replicas 4, 10, and 16 of the reaction pathway and monitored R1, R2, R3, R4, RN, and RO distances.

Between replicas 1 and 4, the lytic water molecule stabilized by two hydrogen bonds with Glu87 residue readjusted itself. Until step 4, the R1 distance remained consistent (1.75 vs 1.76 Å), while the R2 increased from 2.98 to 3.47 Å. Then, the system moved toward the TS and the Pγ–OS bond breaks in replica 8 (R1 = 1.94 Å). The γ-phosphate, OS, and water Ow oxygen atoms form a trigonal bipyramidal structure. The dissociation of the Pγ–OS bond thus occurs before the proton and OH– transfer. In reaction step 9, OH– and Hw+ comprising the lytic water have already moved to OE and Pγ atoms of ATP via the concerted reaction mechanism (RO = 0.99 Å, RO = 1.93 Å). In the succeeding replicas, a product state is reached (PS) in which R1 = 2.90 Å and R2 = 1.63 Å. The resulting energy profile has a structure approximating the TS with an activation energy of about 22.0 kcal/mol, which is reasonably close to previous theoretical and experimental studies of comparable molecular motors.[62]

Subsequently, we explored the GBC reaction mechanism of ATP hydrolysis. For this reaction mechanism, we measured six key distances to describe the reaction events and plotted them in Figure B (Animation 3). The R1 distance characterizes the breakage of the bond between the oxygen on β-phosphate and Pγ. R2 monitors transfer of the OH– nucleophile from the lytic water to the Pγ. R3 and Ow–Hw monitor the transfer of Hw from the lytic water to Lys378. With R4 and NE–HN distances, we follow the proton transfer from the QTK Lys378 to the γ-phosphate species.

In the first four reactions steps, the system slightly readjusts itself; R1 and R4 remained stable, while R2 increases from 3.66 to 3.80 Å, and R3 from 3.06 to 3.26 Å. Then, the system shifts toward the TS. In step 9, the proton transfers from the NE nitrogen of Lys378 to the OE of the γ-phosphate, signified by the RN increasing to 1.73 Å and R4 decreasing to 1.05 Å. ATP and the lytic water stabilized via the interaction with Glu87 continue to reorient themselves in replica 10. The lytic water approaches both Pγ and NE, causing a decrease in R2 and R3 to 2.31 Å and 1.79 Å, respectively. In this replica, R1 reaches 2.03 Å and the phosphodiester bond breaks. We again find Pγ forming a trigonal bipyramidal TS together with OS and Ow atoms. In the next reaction step, there is a simultaneous transfer of lytic water’s Hw+ proton and OH– nucleophile to the NE atom on Lys378 and γ-phosphorus of ATP. The products state is finally reached in the succeeding replicas with R1 distance of 3.33 Å and R2 of 1.65 Å (replica 16). With the activation energy of 52 kcal/mol, the GBC reaction mechanism is energetically much less favorable compared to the pathway obtained for SAC.

Comparing the energy of the products, the GBC configuration displays lower energy. The observed higher energy of the product structure in the SAC is probably related to the final position of the proton transferred from the lytic water molecule to the oxygen of the γ-phosphate. Unlike in GBC, it does not enable the formation of a H-bond between the γ-phosphate and ADP. This suggests that in SAC, some additional readjusting of the hydrolyzed γ-phosphate is likely to occur to lead to more favorable product energy. Although not presently modeled, such a readjustment would not change the main feature—the dissociative reaction mechanism.

In both investigated variations of the ATP hydrolysis, we observed a dissociative reaction mechanism with a simultaneous Hw+ proton, and OH– transfer of the lytic water in and could not identify a TS with only the OH– species being present. Such a dissociative mechanism has already been previously observed in multiscale simulations of other relevant biomolecular systems[60] including molecular motors,[62,67−69] which depend on the chemical energy of ATP to execute motility processes. Observed TS formed by the Pγ, OS, and Ow atoms resembled a trigonal bipyramidal structure. Our computational results heavily favored SAC and indeed, the need for a GB could be considered ambiguous, as the catalytic benefit of the water deprotonation is estimated to be rather meagre.[58]

Reaction pathways also showed that both Glu87 and Lys378 favorably stabilize the TS but have a different role in the catalytic reaction. Lys378 primarily stabilizes the reactant state and TS via H-bonding with the non-bridging oxygen of the γ-phosphate, rather than interacting with Glu87 or with the lytic water molecule. For the Glu87 residue, earlier studies, on htIIα and DNA gyrase pinned it as a catalytic residue.[40,70] During our study, we performed GBC calculations with Glu87 acting as a catalytic base, but our attempts were unsuccessful as no stable structure with a protonated Glu87 could be obtained. Instead our QM/MM simulations suggested Glu87 plays its catalytic role by properly positioning the crucial lytic water molecule. This would enable substantial catalysis of the reaction rendering the successful attack of Pγ possible.[58] In addition, simulations of ATP hydrolysis in other molecular motors already established the polarization of involved waters by a glutamate residue.[62]

Performed multiscale simulations of ATP hydrolysis catalyzed by molecular motors often suggested that the reaction takes place via a proton-wire reaction mechanism involving uncoordinated water molecules facilitating the lytic H+ transfer. In the active site of htIIα, there are two waters in the vicinity of the Pγ, atom that could theoretically support a similar proton transfer mechanism via further involvement of the Asp94 residue. However, both of these water molecules are tightly coordinated with the Mg2+ ion and remained so during the multiscale reaction modeling where they were included in the QM region. As mentioned, molecular simulations also demonstrated that the interaction between Mg2+ and these water molecules is stable. All this taken into consideration makes the wire transfer mechanism an unlikely catalytic strategy for the htIIα system.

The used QM DFT setup worked well in our case; nevertheless, it should be noted that the use of newly developed DFT functionals such as M06-2X and basis sets might provide even superior description of the studied reaction. Finally, the application of the RPATh method allowed the generation of an approximate MEP that is not biased toward a prechosen reaction coordinate. However, a more detailed representation of the reaction free energy surface would require more advanced methods with more comprehensive sampling. This could be achieved, for example, by metadynamics simulations[71] or QM/MM methods based on free energy perturbations to provide more precise insights into the reaction energetics and potentially reinvestigate and reassess the steps and roles in the reaction mechanism.

Experimental Study of the Lys378 Role

To further verify the role of Lys378, we sought to experimentally characterize ATP hydrolysis in a K378A htIIα mutant. In all expression experiments, the cells appeared to be infected by appearance. However, the protein appeared to be expressed at quite low levels or not at all on the sodium dodecyl sulfate-polyacrylamide gel electrophoresis gels and in the western blots. The latter showed that some intact protein was present 24 h after infection, with some degradation products also visible. After 48 h, no intact protein was seen. This suggests that the mutant is degraded upon expression, and we confirmed this by western blot analyzes (see Supporting Information, Section S3). Interestingly, when the analogous Lys337 residue in bacterial DNA gyrase was mutated to Glu, stability problems also occurred; activity began to slowly decline over time.[72] The mutant protein behaved as expected during purification in a larger scale infection, that is, it behaved as the wild-type (WT) protein However, the yield was still low because the amount of protein present in the starting material was low. It did not appear to degrade during the purification process. The final protein concentration of the K378A htIIα mutant was 25 μg/mL. The K378A mutant was initially tested for its activity in the topo II-mediated decatenation assay and was found to possess activity (data not shown). The WT used as a control was at 200 μg/mL.

Having verified that the K378A htIIα mutant retained its function via the decatenation study, we moved to examine the impact of this point mutation on the ATPase activity. The data from the runs after the addition of the ATP were plotted and lines were then fitted to this data by a linear regression (Figure ). The rates for each ATP concentration were then calculated as μM ATP hydrolyzed/min and rates for both WT and K378A htIIα mutant were then plotted against the ATP concentration in Figure and curves were fitted to the obtained data using a hyperbolic equation. Individual and average calculated values are further depicted in Table .

Figure 6

Results of the ATPase assay: (A) WT htIIα and (B) K378A mutant. The rates were plotted against the ATP concentration and curves fitted using the hyperbolic equation y = ax/(b + x), where y = rate, x = [ATP], a = Vmax, and b = Km. Detailed data in Section S4 of the Supporting Information.

Table 1

ATPase Assay; Calculated Km (mM ATP) and Vmax (μM ATP/min) Values

	WT run 1	WT run 2	av. WT	K378 run 1	K378 run 2	K378 run 3	av. K378
Vmax (μM/min)	1.22	0.96	1.09	0.23	0.29	0.25	0.26
Km (mM)	0.24	0.22	0.23	0.15	0.12	0.096	0.122
kcat (s–1)	1.27	1.00	1.14	0.55	0.69	0.6	0.65

In the K378A mutant, the Km and kcat were reduced to about half of the WT htIIα, while the Vmax was a mere 25% of the original. When the analogous Lys337 was mutated to Glu in bacterial DNA gyrase B, the loss of ATPase activity was much more profound. The mutant domain had a 103-fold decrease in kcat, and for the full-length enzyme, ATPase rates were more than 10 times lower.[72]

The results of the performed point-mutation study nicely complemented and substantiated our multiscale computational results of reaction modeling. Lys378 is not essential for the catalysis, and the most likely mechanism of ATP hydrolysis in htIIα is SAC. However, Lys378 does seem to play an important role in enzymatic efficiency. We deduce it helps in the stabilization of the reactants and the TS so that the ATP hydrolysis can proceed smoothly; in the reaction profile, we observed stable H-bonding of the side chain of Lys378 with the non-bridging oxygen of the γ-phosphate of the ATP. Additionally, mutations of the analogous Lys359 in Drosophila showed dramatically reduced but not completely abolished ATPase activity, severe reduction of the stimulating effect of DNA on ATP hydrolysis, and also disruption of the N-gate operating mechanism; the enzyme could trap DNA even without ATP or ATP analogues present.[73] Same was true for a mutant of htIIα with a deleted QTK loop.[25] This positions the QTK loop as essential for interdomain communication and Lys378 as having the additional role of a sensor amino acid required for the correct operation of the clamp. Broadly, the proposed role of Lys378 is consistent with observations in several molecular motors, such as myosin, kinesin, and F1-ATPase, in which a stabilizing contribution has been attributed to the Lys residue located near the transferred γ-phosphate of ATP.[62] Moreover, such a stabilizing function appears to be observed in other systems in which phosphate group transfer occurs,[74−76] further confirming the proposed catalytic mechanism in this biomolecular system.

Conclusions

Human DNA topoisomerase IIα is a biological nanomachine that regulates DNA topological changes. It is also an established anticancer target that is currently being revisited in the development of new catalytic inhibitors, some of which target the ATP binding site. To further unravel some of the mechanistic aspects of the htIIα mode of action, we have combined all-atom simulations, multiscale QM/MM calculations, and a point mutation study to investigate the htIIα ATPase domain dimer. All-atom simulations of the holo and apo systems of the htIIα ATPase domain dimer showed that the binding of ATP influences its dynamics. Furthermore, using multiscale QM/MM methods, we investigated two reaction mechanisms of ATP hydrolysis, substrate-assisted and GBC. The computational results indicated that ATP hydrolysis proceeds via a dissociative mechanism, consistent with the proposed pattern observed when studying ATP-driven molecular motors, and heavily favored the SAC reaction pathway. To supplement our computational results and further asses the probability of the two reaction mechanisms, we constructed a K378A htIIα point mutant, which revealed that Lys378 plays a role in stabilizing the reacting complex but does not act as a catalytic base. The obtained results thus positioned the determined SAC mechanism of ATP hydrolysis as more probable. The catalytic Glu87 residue was not observed to act as a catalytic base, but instead served to properly position and polarize the lytic water molecule, which is also consistent with previous studies. We hope that this work will contribute to a deeper understanding of these intriguing molecular motors and provide new information to advance the development of new therapeutics.

Experimental Section

All-Atom Simulations of the Human Topo ΙΙα Dimer

The X-ray structure of ATPase domain dimer of htIIα with bound non-hydrolysable ATP analogues AMP–PNP is publicly available (PDB: 1ZXM). MD calculations were performed using the CHARMM molecular modeling suite.[77] The missing side chains and protein loop residues not provided in the 1ZXM PDB file were generated using the PDB Hydro web server (http://lorentz.dynstr.pasteur.fr/pdb/index.php) in a two-stage procedure.[78] First, the 1ZXM crystal structure was submitted to the PDB Hydro web server to construct the missing protein loop. The server automatically detected the missing loop residues and added them into the structure.[79] The upgraded structure was retrieved and visually inspected to assess the generated conformation of the loop. Another module of the Hydro server was then used to detect the missing side chain and screen possible rotamers of these side chains selecting those with the lowest van der Waals energy to be added to the structure. The process was done consecutively for each residue resulting in the final structure with all protein atoms included.[78] Also, the AMP–PNP ligand present in the 1ZXM structure was changed to ATP. CHARMM-GUI[80] was used to generate the solvated complexes and two systems of the human topo IIα dimer were subsequently prepared for MD simulations, (1) htIIα devoid of ATP ligands (apo structure) and (2) a topo II dimer with two ATP and two Mg2+ ions bound to their corresponding binding sites (holo structure).

CHARMM parameter and topology files (version 26) specified the force field parameters of the amino acid residues of the topo II ATPase dimers.[81,82] CHARMM general force field (CGenFF) modeled the partial charges and atom types of the ATP molecule.[83] The system was immersed into TIP3 water molecules[84] with truncated octahedral shape and edge distance of 10 Å. Chlorine and potassium ions were added to make the system electroneutral with further ions added to the final KCl concentration of 0.15 M. Ion placement was performed using the Monte Carlo method. The periodic boundary conditions were applied and grid information for the particle-mesh Ewald fast Fourier transform was generated automatically. The apo and holo topo II structures comprised 78,948 and 79,078 atoms, respectively. Short energy minimizations were then carried out to remove bad contacts. Both topo II systems were subsequently minimized by the steepest descent method, followed by the modified adopted basis Newton–Raphson (ABNR) method (both for 2000 steps) and an equilibration MD of 4 ns using 1 fs simulation step. The production MD trajectories were generated by leap-frog integration scheme and 2 fs simulation step coupled with SHAKE algorithm. A 0.10 μs long MD simulation production runs were performed for each structure. Sampling of conformations occurred every 100th step extracting 50,000 conformations for analysis. Visualization and analysis of the geometry parameters of the production MD trajectories were performed using Visual MD (VMD) program.[85]

Analysis of the MD Trajectories

MD Trajectories were examined by employing the Cpptaj module of AmberTools 18[86] to calculate root mean square deviations (RMSDs), RMSFs, and PCA. Bio3D library (version 2.3-0)[87] within the R environment enabled the calculations of the dynamic cross-correlation maps (DCCMs). Results were visualized by VMD,[85] R software environment[88] and PyMOL.[89]

Principal Component Analysis

Covariance matrices of positional fluctuations were built and diagonalized to obtain the eigenvalues, and their corresponding eigenvectors. The eigenvectors possessing the largest eigenvalues correspond to the most relevant motions sampled in the simulation (known as principal components, PCs). We calculated the cumulative variances accounted by the PCs for both systems. First two PCs were considered in further explanation of differences in motions between simulations. In the holo system and apo system, first two PCs captured approximately 40% of overall position fluctuations. We used the normal mode Wizard[90] plugin of VMD program to visualize the large-scale collective motions observed in our simulations.

Dynamical Cross-Correlation Map (DCCM) Calculations

The cross-correlation maps were calculated using the DCCM function in the Bio3D package, which determines the covariance matrices and computes the Pearson correlation coefficient (C) on the Cα atom pairs, i and j via eq

In the equation, c is defined as , and Δr represents a displacement vector of atom i and Δr of atom j with the brackets signifying an ensemble average. Cross-correlation coefficient C was obtained by normalizing the covariances. C values range from −1 to 1, corresponding to type of correlated motion between pairs. The positive vales indicate positively correlated motions (movement in the same direction), while the negative values indicate the negatively correlated motions, anticorrelations (movement in the opposite direction). Visualization of the anticorrelated motions corresponding to the residues in the matrices was completed in PyMOL.

Multiscale QM/MM Study of the ATP Hydrolysis

Preparation of the Reactant and Product Structures

Again, we used coordinates of one protomer of the experimentally determined ATPase dimer of the htIIα complexed with the bound AMP–PNP ligand (PDB: 1ZXM). We left crystal waters inside the system as they were considered a part of the reaction mechanism. Multiscale QM/MM calculations[46] were also performed within the CHARMM environment.[91] HBUILD command was used to add hydrogen atoms. Then, the system was solvated with a TIP3P waters cubic box, and 20 potassium and 13 chloride ions were added as counterions, making the system neutral. The topo II system was divided to the QM region which was treated with quantum mechanics, while the rest of the system was described with molecular mechanics. The QM region encompassed the methyl-triphosphate fragment of the ATP molecule, lytic water, Mg2+ ion, two Mg2+-coordinated water molecules, and side chains of residues Lys378 and Glu87, that is, atoms that could feasibly take part in the H+ and OH– transfer (Figure S5). For the MM part of the calculations, we used CHARMM (version 36). QM DFT calculations were performed using GAMESS (general atomic and molecular electronic structure system),[92] interfaced with CHARMM. The quantum regions were treated with the B3LYP DFT method and 6-31G* basis set. ABNR minimization algorithm was used for geometry optimization. MM calculations were performed with a dielectric constant of ε = 1 with a classical force shift method and a cutoff distance of 12 Å. The empty valances of the QM Glu87, Lys378, and ATP moieties were filled with a link atoms (hydrogen).

We used CHARMM restrain distance (RESD) methodology (eq ) to form required covalent bonds of the product structures where the ATP hydrolysis was completed. This was achieved by adding restrains to the CHARMM energy function, thereby directing the movement along the chemically relevant reaction coordinate and forcing the breakage/formation of the anticipated chemical bond (water hydrolysis, proton transfer).[47]k denotes an applied force constant, d0 is the current distance between a pair of selected atoms, and deq is the targeted/equilibrium distance.

For the SAC mechanism, the product structure was obtained by first restraining a distance between the γ-phosphorus atom Pγ and the Ow of the lytic water molecule to 1.6 Å in 100 steps of RESD QM/MM minimization and then minimizing the obtained structure with the previous constraint removed until the root-mean-squared gradient was smaller than 0.0001 kcal/mol/Å. For the product structure resulting from the GBC reaction mechanism, a similar procedure was applied only more restrains were introduced. Here, we restrained (1) the γ-phosphorus atom Pγ and the Ow of the lytic water molecule to 1.6 Å, (2) the Hw proton of the lytic molecule and NE of Lys378 to 1 Å, and (3) the proton from the NE of Lys378 and OE oxygen of Pγ to 1 Å. The obtained structure was then minimized without the present restrains.

Generation of Reaction Pathways

The replica path (RPATh) method[49] is an extension of the self-penalty walk methods that can model reaction pathways.[52,53] The 16 initial replicas of both investigated reaction pathways were created by a linear interpolation of coordinates between atoms of the QM/MM minimized starting and the corresponding product structures. The method further uses a combined minimization involving the sum of the configurational energies and two RPATh penalty energy terms. First penalty term restrains (eq ) the distances between adjacent replica points. This safeguards that the pathway is regularly spaced and smooth.

NREP is the number of generated replicas and rms(i,i + 1) is the best fit RMS value (eq ) between successive replicas⟨rms⟩ denotes the average RMS value (eq ), while w is the atom weighting factor (eq )

To restrain the angle between an adjacent and the next adjacent pathway points, we used a second RPATh force constant, which uses the law of cosines, and with that, we avoid paths that double back on themselves (eqs –9). In this manner, the RPATh method ensures that the optimized replicas represent the reaction pathway.[48,49]

QM/MM RPATh optimizations of the initially generated two pathways were performed by applying an ABNR method.[50] Parameters of the RPATh penalty terms were as follows: Krms = 2000.0 kcal/mol/Å2, Kangle = 100.0 kcal/mol/Å, and COSMAX = 0.95 radian. To evade the overestimation of contributions of those atoms not involved in the reaction, only the QM region was weighted in the RPATh RMS calculation.[50] Pathways were optimized for 10,000 steps or until for at least 30 consecutive steps, the total pathway root mean square gradient was less than 0.01 kcal/mol/Å and the change in total pathway energy was less than 1.0 kcal/mol. After initial reaction pathways were obtained, we further optimized them by selecting the new initial and end coordinates closer to the observed TS. We then recalculated the pathways by generating again 16 replicas using the same settings for RPATh penalty terms. This process could be subject to multiple iterations to obtain appropriate initial and end reaction coordinates. Finally, the energy profiles were calculated for both reaction mechanisms under investigation.

Construction and Expression of the Human Topo IIα K378A Mutant

The K378A mutant was made in the WT human htIIα gene in the pFastbac vector. This was used to transform DH10bac cells and positive clones were identified by sequencing of the mutant gene. The bacmid was used to transfect Spodoptera frugiperda (SF21) insect cells. Virus from the transfection was used to infect more insect cells to amplify the virus stock. A larger (300 mL) culture of insect cells was infected with the P4 stock of virus, grown 24 h and the human topo IIα purified by nuclear preparation, followed by column chromatography. Expression of the protein was confirmed by western blot using polyclonal anti-human topo II alpha antibodies.

ATPase Assay of the WT and htIIα K378A Mutant

The ATPase activity of both systems was examined using a linked assay. The hydrolysis of ATP produces ADP which coupled with pyruvate kinase/lactate dehydrogenase mix leads to a conversion of NADH into NAD. The reduction of NADH is then monitored at 340 nm. A mix of assay buffer (5 μL of 10× buffer per assay: final conc. 20 mM Tris-HCl, 5 mM magnesium acetate, 125 mM potassium acetate, 2 mM dithiothreitol, pH = 7.9), linear pBR322 (1.5 μL of 1 mg/mL per assay), phosphoenol pyruvate (0.5 μL of 80 mM per assay), pyruvate kinase/lactate dehydrogenase mix (0.75 μL per assay), NADH (1 μL of 20 mM per assay), dimethyl sulfoxide (1 μL per assay), and water (32.85 μL per assay) was first prepared. 41.1 μL was aliquoted into the wells of a 384-well microtiter plate. 5 μL of the dilution buffer or human topo II (WT; 170 nM stock giving 17 nM final concentration: K378; 70 nM stock giving 7 nM final) was then added and mixed. The OD340 was monitored for 10 min and then the reaction started by adding 3.4 μL of the fitting concentration of ATP and the OD340 was monitored for up to 40 min. Assays were performed at 37 °C and two negative controls (dilution buffer but no enzyme) were also run.

Data and Software Availability

All molecular simulations, analysis, and visualization were performed with widely used programs available freely for academic institutions: CHARMM (version 40), GAMESS (May 2013 R1 release), AmberTools 18, PyMol 2.0, VMD 1.9.2, and R software environment with the Bio3D library 2.3-0. The starting structure was obtained from the Protein Data Bank public structure database. All procedures and workflows are described in the Methods section. Structures used for MD and QM/MM replica path simulations in pdb format are provided in the Supporting Information.