The Impairment of TorsinA's Binding to and Interactions With Its Activator: An Atomistic Molecular Dynamics Study of Primary Dystonia

Primary dystonia's prolonged muscle contractions and the associated abnormal postures and twisting movements remain incurable. Genetic mutation/deletion of GAG from TorsonA's gene resulting in ΔE303 (which weakens the binding between TorsinA and its activator, such as LULL1) primarily cause this neurodegenerative disorder. We studied TorsinA-LULL1 (or TorsinAΔE303-LULL1) bindings and interactions. For the first time, we show the atomic details of TorsinA-LULL1 dynamic interactions and TorsinAΔE303-LULL1 dynamic interactions and their binding affinities. Our results show extensive effects of ΔE303 on TorsinAΔE303-LULL1 interactions, and suggest that the differences between TorsinA-LULL1 interactions and TorsinAΔE303-LULL1 interactions are non-subtle. ΔE303 significantly weakens TorsinAΔE303-LULL1's binding affinity. We present pieces of evidence proving that the effects of ΔE303 (on the differences between TorsinA-LULL1 interactions and TorsinAΔE303-LULL1 interactions) are more pronounced than previously suggested, and that the nanobody used for achieving the X-ray crystallization in the previous study attenuated the differences between TorsinA-LULL1 and TorsinAΔE303-LULL1 interactions. Our accounts of the dynamic interactions between “TorsinA and LULL1” and between “TorsinAΔE303 and LULL1” and the detailed effects of ΔE303 on TorsinA-/TorsinAΔE303-LULL1 build on previous findings and offer new insights for a better understanding of the molecular basis of Primary Dystonia. Our results have long-term potentials of guiding the development of medications for the disease.

Introduction

Primary dystonia is a neurodegenerative disease characterized by prolonged involuntary muscle contractions, twisting movements, and abnormal postures (Fahn, 1988; Breakefield et al., 2008; Demircioglu et al., 2016). The disease affects 2-7320 individuals per million people (Defazio et al., 2004) with prevalence varying very widely across populations (Dystonia in Europe (ESDE) Collaborative Group and others, 2000; Müller et al., 2002; Defazio et al., 2004; Steeves et al., 2012). Dystonia as a whole (i.e. primary dystonia, secondary dystonia, and dystonia plus) are essentially motor symptoms of one or more underlying pathophysiological states and could be triggered by a number of insults (Kojovic et al., 2013). Since the symptoms of dystonia resembles those from basal ganglia lesions, researchers initially thought that primary dystonia was a basal ganglia disease (Marsden et al., 1985). Although cerebellar deficit may contribute to primary dystonia's symptoms (Sadnicka et al., 2012), studies have shown that genetic mutation in DYT1 gene (which is often hereditary and seldom sporadic) resulting in the deletion of glutamic acid 303 (hereafter referred to as ΔE303, Figure 1) from TorsinA protein (primarily found in nerve cells) is the most common cause of primary dystonia (Ozelius et al., 1997; Leung et al., 2001; Breakefield et al., 2008. Since primary dystonia results from a dominant allele, a copy of the mutated DYT1 gene is sufficient for the disease to manifest such that a child could have the disease if at least one of his/her parents carry the gene for the disease. However, the genetic defect exhibits incomplete penetrance such that only about 30% of the carriers show clinically detectable symptoms of dystonia (Sharma et al., 2005).

FIGURE 1

Figure 1. The four molecular systems studied are shown in (A) TorsinA+LULL1, (B) TorsinAΔE303+LULL1, (C) TorsinA+LULL1+VHH-BS2, and (D) TorsinAΔE303+LULL1+VHH-BS2. The inset in (A) shows the position of E303. VHH-BS2 (shown in C,D) is a nanobody used in the crystallization. The atomic coordinates obtained from the nanobody-aided X-ray crystallography have been described in a separate study (Demircioglu et al., 2016) and are in the PDB with accession identifications 5J1S and 5J1T for the wild type and the mutant respectively.

Being a member of the AAA+ ATPase enzymes family, a normal and functional TorsinA (like other family members) essentially hydrolyzes ATP and uses the energy released to drive cellular activities (Hanson and Whiteheart, 2005), such as modification of the structures of other molecules. However, unlike the other self-activating AAA+ ATPase enzymes, TorsinA must bind to and interact with either “Lamina-Associated Protein 1 (LAP1)” or “LUminal Domain Like LAP1 (LULL1, a paralog of LAP1, Figure 1)” to be activated (Zhao et al., 2013; Brown et al., 2014; Rose et al., 2015). The TorsinA's ΔE303 reduces its binding to (and its interactions with) LULL1 or LAP1. Unable to be properly activated (due to its weak binding to or its dissociation from its activator), the defective TorsinA ultimately causes disruptions in neuronal communications and muscular controls, which results in primary dystonia (Breakefield et al., 2008; Naismith et al., 2009).

The atomic structure of TorsinA and TorsinAΔE303 were recently solved (Demircioglu et al., 2016) by nanobody-aided X-ray crystallography. The observed differences in the structures of TorsinA and TorsinAΔE303 made the researchers believe that there is a subtle difference in TorsinA-LULL1 interactions and TorsinAΔE303-LULL1 interactions. The researchers claimed that “a comparison of these structures shows, in atomic detail, the subtle differences in activator interactions that separate the healthy from the diseased state” (Demircioglu et al., 2016). While it is correct that a subtle difference was observed in the structure of TorsinA and TorsinAΔE303 from the nanobody-aided X-ray crystallography structures, the differences between TorsinA-LULL1 interactions and TorsinAΔE303-LULL1 interactions may not necessarily be subtle in the real biological conditions involving dynamic interactions, which are different from the nanobody-aided X-ray crystallography conditions.

Here, we show that the previous account of the effects of ΔE303 on the differences between TorsinA-LULL1 interactions and TorsinAΔE303-LULL1 interactions may have underestimated the effects of ΔE303, and lack atomic details of the dynamic interactions between TorsinA (or TorsinAΔE303) and LULL1. Using 4.8 μs trajectories from our unbiased explicitly-solvated all-atom molecular dynamics (MD) simulations (of the four systems shown in Figures 1A–D), 3.0 μs of Gaussian accelerated MD (GaMD) simulations, and the results from our binding free energy change calculations using Adaptively Biased MD (ABMD) simulations, we show for the first time (1) the atomic details of the dynamic interactions between TorsinA (or TorsinAΔE303) and LULL1; (2) their binding free energy change; (3) the extensive effects of ΔE303 on TorsinAΔE303-LULL1 interactions; (4) that the differences between TorsinA-LULL1 interactions and TorsinAΔE303-LULL1 interactions are non-subtle but extensive/prominent; and (5) that the nanobody-aided X-ray crystallization conditions attenuated the differences between TorsinA and TorsinAΔE303; such that the effects of the ΔE303 are more pronounced than previously suggested (by nanobody-aided X-ray crystallography results). Nonetheless, we acknowledge that the previous work (Demircioglu et al., 2016), which makes the atomic structure of TorsinA-LULL1 complex and TorsinAΔE303-LULL1 complex available, remains a remarkable and very important contribution to this field of science, and the current work would be impossible without theirs.

Results and Discussions

Based on our four (A:TorsinA+LULL1, B:TorsinAΔE303+LULL1, C:TorsinA+LULL1+VHH-BS2, and D: TorsinAΔE303+LULL1+VHH-BS2, each with three replicates, Figure 1A) unbiased explicit-solvent molecular dynamics (MD) simulations resulting in 4.8 μs production MD trajectories, our 3.0 μs Gaussian accelerated MD (GaMD) simulations, and our binding free energy calculations (all with AMBER16 Pearlman et al., 1995; Case et al., 2005; Salomon-Ferrer et al., 2013 and ff14SB (Maier et al., 2015 force-fields), it is evident that the deletion of glutamic acid 303 from TorsinA, thus TorsinAΔE303, weakens its binding to LULL1, which is in agreement with the findings of previous studies (Németh, 2002; Sharma et al., 2005) which have shown that the ΔE303 mutation disrupts TorsinA's interaction with its activator (Demircioglu et al., 2016), TorsinA's ATPase activities (Konakova and Pulst, 2005), and TorsinA's ability to form multimeric complex (Pham et al., 2006). Furthermore, our results, as discussed in the following subsections, show that important intermolecular interactions are lost or weakened in TorsinAΔE303-LULL1 complex compared to TorsinA-LULL1 complex; that ΔE303 has more severe effects on TorsinA-/TorsinAΔE303-LULL1 interactions and binding than previously suggested in a nanobody-aided X-ray crystallography study (Demircioglu et al., 2016); and that nanobody-aided X-ray crystallization conditions attenuated the differences between TorsinA and TorsinAΔE303, and the differences between their interactions with LULL1. We provide the details and the respective explanations and discussions in the following subsections.

The Deletion of Glutamic Acid 303 From TorsinA, Thus TorsinAΔE303, Weakens Its Binding to LULL1

We calculated the difference between the binding free energy change for TorsinAΔE303-LULL1 and the wild type, WT, TorsinA-LULL1 (i.e. ΔΔG_ΔE303−WT, shown in Table 1) using both Adaptively Biased Molecular Dynamics (ABMD) simulations (Huber et al., 1994; Darve and Pohorille, 2001; Wang and Landau, 2001; Babin et al., 2008) and Molecular Mechanics Poisson-Boltzmann Surface Area (Kollman et al., 2000; Wang et al., 2001, 2006) (MM-PBSA). Our results show that TorsinAΔE303 has a weaker binding to and fewer interactions with LULL1 compared to the WT TorsinA. This is in agreement with the results of previous studies (Naismith et al., 2009; Zhu et al., 2010; Zhao et al., 2013) which also show that ΔE303 weakens TorsinAΔE303-LULL1's or TorsinAΔE303-LAP1's binding in comparison to the WT TorsinA. Quantitatively, the results of the binding free energy change calculations (via both ABMD and MM-PBSA) show that ΔE303 weakens the binding of TorsinAΔE303 to LULL1 by a ΔΔG_ΔE303−WT ≥ 4.7 ± 1.3 kcal/mol (Table 1), which corresponds to TorsinAΔE303-LULL1 having a binding affinity that is significantly lower than the binding affinity of TorsinA-LULL1 (Table 1). In addition to the agreement between our results and the findings of the previous studies, which showed that the ΔE303 weakens the interactions between TorsinAΔE303 and its activator (Naismith et al., 2009; Zhu et al., 2010; Zhao et al., 2013), our results show a quantitative estimate of the differences between the TorsinAΔE303-LULL1's and TorsinA-LULL1's binding free energy change (Table 1) thereby providing new quantitative insights on the effects of ΔE303 in TorsinAΔE303-LULL1's binding. At this point, we must acknowledge that ABMD is more rigorous than MM-PBSA and that the results from the ABMD (Table 1) are better estimates of the free energy changes than the results from MM-PBSA (Table 1) which has been shown to often overestimates free energy changes (Genheden et al., 2011; Homeyer and Gohlke, 2012).

TABLE 1

Table 1. Binding Gibbs free energy change, ΔG, for TorsinA-LULL1, and TorsinAΔE303-LULL1.

Furthermore, the ΔG_ΔE303 of −9.45 ± 0.57 kcal/mol (Table 1) suggests that ΔE303 may not guarantee complete and spontaneous dissociation of LULL1 from TorsinAΔE303 at all times, even though the ΔE303 significantly weakens the binding of LULL1 to TorsinAΔE303 (ΔΔG_ΔE303−WT = 4.67 ± 1.26 kcal/mol, Table 1). We believe that this is the reason why it was possible to obtain the X-ray crystallography structure of LULL1 bound to TorsinAΔE303 (Demircioglu et al., 2016). However, it is also very important to note that binding ΔG involving protein-protein interactions (such as in the current case) are prone to overestimation (Noskov and Lim, 2001; Gohlke et al., 2003; Kiel et al., 2004; Andberg, 2011). Therefore, even though the obtained binding ΔΔG_ΔE303−WT from the ABMD (Table 1) could be seen as a reliable estimate of the effects of the mutation in weakening TorsinAΔE303-LULL1 binding as compared to TorsinA-LULL1 binding, the individual binding ΔG (namely ΔG_WT and ΔG_ΔE303, Table 1) may have been overestimated and should be interpreted with caution. The effects of potential overestimation of the binding ΔG involving protein-protein interactions should have canceled out in the relative binding energy change, ΔΔG_ΔE303−WT, making the interpretations of the ΔΔG_ΔE303−WT more reliable than the interpretations of the individual binding ΔG (namely ΔG_WT and ΔG_ΔE303).

Important Intermolecular Interactions Are Lost or Weakened in TorsinAΔE303-LULL1 Complex

Using a cut-off distance of 6 Å between residues' centers of mass, we carried out residue-residue contact analysis for both the TorsinA-LULL1 and the TorsinAΔE303-LULL1 molecular systems over the unbiased production MD trajectory. Furthermore, we decomposed the binding enthalpy obtained from the MM-BPSA over the unbiased production MD trajectory into the residue-level contributions so as to be able to assess the contributions of each of the residues toward TorsinA-/ TorsinAΔE303-LULL1 binding. The results from the residue-residue contact analysis and from the MM-PBSA show that Y328, K317, T135, L136, H140, K320, T321, T104, T305, N208, D173, A209, A211, D327, K174, G210, E302, are the amino acids of TorsinA playing the most vital roles in the dynamic interactions between (as well as in the favorable binding of) TorsinA and LULL1 (Figure 2, Supplementary Figure S1). The favorable TorsinA-LULL1 interactions/binding that are mediated by these residues are highly compromised in the mutant, TorsinAΔE303-LULL1 (Figures 2, 3, Supplementary Figure S1) as shown by the red bins in Figure 2C, and longer residue-reside distances in the TorsinAΔE303-LULL1 systems shown in Figures 2D, 3.

FIGURE 2

Figure 2. ΔE303 significantly and persistently reduces TorsinAΔE303-LULL1 interactions compared to TorsinA-LULL1 interactions. We show the most active residues (A) in TorsinA-LULL1 interactions and (B) in TorsinAΔE303-LULL1 interactions, and (C) the interactions changes due to ΔE303 (i.e., B − A). The red cells in C indicate the residue-residue interactions that are lower in TorsinAΔE303-LULL1 compared to TorsinA-LULL1. The details of the weakened TorsinAΔE303-LULL1 interactions as compared to TorsinA-LULL1 interactions from a representative frame from the MD trajectory are shown in (D). The first image in (D) show the overall differences in the TorsinAΔE303-LULL1 interactions and the TorsinA-LULL1 interactions. The position of the E303 in the wild type TorsinA (which is deleted from the TorsinAΔE303) is shown at the top of the image. The comparisons of each pair of residues are shown in the rest of panel D with TorsinA-LULL1 interactions shown in blue and green, and TorsinAΔE303-LULL1 interactions shown in orange and red. Only the top 11 most weakened interactions are shown in details.

FIGURE 3

Figure 3. Detailed account of the effects of ΔE303 on the distributions of contacts/distances between TorsinA's residues and LULL1's residues. We show a more detailed account of the effects of ΔE303 on residue-residue (“amino acid of TorsinA/TorsinAΔE303”–“amino acid of LULL1” e.g., K317-W447) distances. A red bounding-box shows that (for the given TorsinA-/TorsinAΔE303-LULL1 residue pair) TorsinA-LULL1 interactions are better than TorsinAΔE303-LULL1 based on a cut-off distance of 6 Å, otherwise a blue bounding-box. In each group of four boxes, the top row shows the distributions of the positions of the amino acid of TorsinA-/TorsinAΔE303 and that of LULL1 in form of occupancy probability ranging from 0.0 (white) to 1.0 (dark blue), while the bottom row contains histograms showing the distributions of TorsinA-/TorsinAΔE303-LULL1 distances. The residue pairs average distances were sorted in descending order based on the magnitude of the differences between TorsinA-LULL1 and TorsinAΔE303-LULL1 such that those with the largest absolute differences are listed first. Only 16 pairs are shown in this figure. Additional examples are presented in Supplementary Figure S1.

Furthermore, we found that the differences between TorsinA-LULL1 interactions and TorsinAΔE303-LULL1 interactions do not involve only a few amino acids pairs but rather involves many amino acids pairs. For example, we observed significant and persistent weakening of the interactions between the following 21 pairs of amino acids of “TorsinAΔE303 and LULL1” compared to the corresponding interactions between the amino acids of “TorsinA and LULL1” sorted in the descending order of the effect of ΔE303 on their interactions: K317&W447, A141&K344, A141&L345, A209&S444, K320&S448, T321&H452, G105&S444, E212&D441, K174&K442, H137&K380, L136&Y379, A211&D441, T305&E416, G210&D441, A209&K442, L136&K380, K317&S448, G210&K442, T321&S451, T305&E415, and K320&H384 (Figures 2, 3, Supplementary Figure S1). However, it appears that the adverse effects of ΔE303 were mildly ameliorated by increases in Y328&K279, G105&S448, L136&L376, D327&K279, and E303&R412 interactions for TorsinAΔE303-LULL1 (Figures 2C, 3, Supplementary Figure S1).

ΔE303 Has More Severe Effects on TorsinA/TorsinAΔE303-LULL1 Interactions and Binding Than Previously Suggested

While investigating the TorsinA-/TorsinAΔE303-LULL1 intermolecular interactions, we assessed the differences between TorsinA-LULL1 and TorsinAΔE303-LULL1 hydrogen-bonding frequencies. Our detailed investigations of these basic intermolecular interactions, show significant differences between TorsinA-LULL1 hydrogen bonds frequencies and TorsinAΔE303-LULL1 hydrogen bonds frequencies (Figure 4, Supplementary Figure S2).

FIGURE 4

Figure 4. Hydrogen bonding between TorsinA and LULL1 is highly compromised by ΔE303. (A) Compared to TorsinA, TorsinAΔE303 forms fewer hydrogen bonds with LULL1. The difference is statistically significant (p < 0.001). Panels (B–L) show examples of hydrogen bonds that are frequently seen in TorsinA-LULL1 interactions but are never (or seldom) seen in TorsinAΔE303-LULL1 interactions. The hydrogen bonds were assessed in 4,500 frames out of 45,000 frames from the MD trajectories by choosing every tenth frame. Please, note that one frame may contain two or more hydrogen bonds between a given residue pairs, hence the possibility of having the overall hydrogen bonds frequency for a given residue pairs (such as 8,999 for e302-R412) greater than the number of frames. The hydrogen bonds are represented by yellow broken lines. A black triangle is placed next to each of the hydrogen bonds to guide the readers' eyes. The amino acids of TorsinA are shown in blue, while those of LULL1 are shown in green. The ratio of the frequency of a hydrogen bonding pattern in TorsinA-LULL1 interactions to its frequency in TorsinAΔE303-LULL1 interactions are presented as the numbers below each panel (e.g. TorsinA-LULL1:TorsinAΔE303-LULL1 = 8999:0 in panel B). Additional examples are provided in Supplementary Figure S2.

We observed that hydrogen bonding between TorsinAΔE303 and LULL1 is compromised by the deletion of glutamic acid 303, ΔE303, in TorsinAΔE303-LULL1. Compared to TorsinA, TorsinAΔE303 forms statistically significant (p < 0.001) fewer hydrogen bonds with LULL1 (Figure 4A). Some important TorsinA-LULL1 hydrogen bonds are completely absent (e.g. E302-R412, R288-E292, F306-E416, etc.) or seldom present (e.g. K325-H452, F306-R419, D331-R278, etc.) in TorsinAΔE303-LULL1 interactions (Figures 4B–L, Supplementary Figure S2). In other words, for example, the hydrogen bonds between TorsinAΔE303's E302 and LULL1's R412 are abolished in TorsinAΔE303-LULL1 (Figure 4B), while the hydrogen bonds between TorsinAΔE303's K325 and LULL1's H452 (Figure 4J) are reduced by about 100 folds, etc. These further show that the differences between TorsinA-LULL1 interactions and TorsinAΔE303-LULL1 interactions are not subtle. We think that the previous (Demircioglu et al., 2016) observation of subtle difference may have resulted from the effects of the nanobody used in the nanobody-aided X-ray crystallization of the complex and/or because X-ray crystallography cannot offer information on the dynamic interactions between the proteins. We further explain this in the next section.

It is surprising, but interesting, that the deletion of glutamic acid 303, ΔE303, from TorsinAΔE303-LULL1 would result in the complete loss of the hydrogen bonds between TorsinAΔE303's E302 and LULL1's R412 in TorsinAΔE303-LULL1 (TorsinA:TorsinAΔE303 = 8999:0, Figure 4B), the complete loss of the hydrogen bonds between TorsinAΔE303's R288 and LULL1's E292 in TorsinAΔE303-LULL1 (TorsinA:TorsinAΔE303 = 3720:0, Figure 4C), etc. These results reinforce the findings of previous studies where it was shown that both E302 (Ozelius et al., 1997) and R288 (Zirn et al., 2008) are important for TorsinA's activation and that their deletion or mutation results in weaker bindings of TorsinA to its activators (Ozelius et al., 1997; Zirn et al., 2008). Our findings (Figures 4B–L, Supplementary Figure S2) suggest that ΔE303 has severe effects on TorsinAΔE303-LULL1 interactions, to the extent of harming other residues of TorsinAΔE303 that would normally have favorable interactions with TorsinA's activator in the wild type with no deletion of E303.

Generally, hydrogen bonds often offer strong intermolecular binding forces, and their loss/absence could, in some contexts, quickly degrade intermolecular binding and interactions. Therefore, the compromised hydrogen bonds in TorsinAΔE303-LULL1 complex (Figure 4, Supplementary Figure S2) can considerably explain the weaker binding between TorsinAΔE303 and LULL1 (compared to the binding between TorsinA and LULL1) (Table 1, Figures 2, 3).

Nanobody-Aided X-Ray Crystallization Attenuated the Differences Between TorsinA and TorsinAΔE303 and the Differences Between Their Interactions With LULL1

The well-pronounced differences we have found between TorsinA-LULL1 interactions and TorsinAΔE303-LULL1 interactions (Figures 2–4, Supplementary Figures S1, S2) made us to carry out further investigations on the dynamic behaviors of TorsinA-LULL1 complex and TorsinAΔE303-LULL1 complex but now in a condition that is a closer to that of the previous nanobody-aided X-ray crystallography study (Demircioglu et al., 2016) by including VHH-BS2 in the studied molecular systems. At this point, we must emphasize that the new molecular system (wherein the VHH-BS2 nanobody is in complex with TorsinA-LULL1 or TorsinAΔE303-LULL1) is essentially a less physiological condition than our original molecular systems' setup (without the nanobody) because TorsinA-LULL1 or TorsinAΔE303-LULL1 is not in complex with VHH-BS2 in the real biological environment. However, this setup allowed us to investigate the potential influence of the nanobody (i.e., VHH-BS2 used in the X-ray crystallography) on the TorsinA-/TorsinAΔE303-LULL1 interactions.

We observed that the addition of VHH-BS2 to the molecular systems reduced the differences between TorsinA-LULL1 interactions (Figure 5A) and TorsinAΔE303-LULL1 interactions (Figure 5B) as compared to their differences in the absence of VHH-BS2 (i.e., the light red bins in Figure 5C compared to the dark red bins in Figure 5D). These findings (regarding the reduction in the differences between TorsinA-LULL1 and TorsinAΔE303-LULL1 systems in the presence of VHH-BS2) are consistent when we assessed the Root Mean Square (RMS) deviations from the X-ray structure (not shown, because they are comparable to those from the energy-minimized structure) or from the energy-minimized structure {TorsinA: 2.168 ± 0.004; TorsinAΔE303: 2.634 ± 0.009; TorsinA(VHH-BS2): 1.756 ± 0.004Å; TorsinAΔE303(VHH-BS2): 1.772 ± 0.006}, the RMS fluctuations (Figure 5E), and the secondary structure elements (by RaFoSA Salawu, 2016, Supplementary Table S1) of TorsinA and TorsinAΔE303 in the presence and in the absence of VHH-BS2. Each of these sets of results show that the X-ray crystallization conditions (e.g., the present of the VHH-BS2 nanobody used in the nanobody-aided X-ray crystallography) has the potentials of attenuating the differences between TorsinA and TorsinAΔE303 and the differences between their interactions with LULL1. The RMS deviations and RMS fluctuations essentially reflect the flexibility of a molecular system with larger values signifying more change/flexibility. On the other hand, for a given molecular system, higher proportions of sheets and helixes and lower proportion of coils or loops may indicate less flexibility (Salawu, 2016).

FIGURE 5

Figure 5. Reduced differences between TorsinA-LULL1 interactions and TorsinAΔE303-LULL1 interactions in the presence of VHH-BS2 nanobody. We show the most active residues (A) in TorsinA-LULL1 interactions and (B) in TorsinAΔE303-LULL1 interactions, and (C) the reduced interactions (i.e., B-A) due to ΔE303. (D) Shows similar information as C, but for the TorsinA-/TorsinAΔE303-LULL1 system without VHH-BS2. Additional pieces of evidence that the crystallization conditions (e.g., the presence of VHH-BS2) reduce the difference between TorsinA and TorsinAΔE303 is shown in (E) wherein the differences in the RMSF between the amino acids of TorsinA and the amino acids of TorsinAΔE303 is attenuated by the presence of VHH-BS2. The E303 deletion from TorsinA (thus TorsinAΔE303) significantly and persistently reduces TorsinAΔE303-LULL1 interactions compared to TorsinA-LULL1, but the effect of the mutation is attenuated by the crystallization environment/conditions (such as the presence of VHH-BS2). A comparison of (C,D) shows (by the fainter shades of blue and red in the current (C), but the darker shades in D) that the environment/conditions wherein the X-ray structures for TorsinA-LULL1 and TorsinAΔE303-LULL1 were solved (such as the presence of VHH-BS2) attenuated the differences between TorsinA-LULL1 interactions and TorsinAΔE303-LULL1 interactions.

Furthermore, we investigated the effects of the nanobody, VHH-BS2, on the differences between the hydrogen bonding patterns in TorsinA-LULL1 complex and TorsinAΔE303-LULL1 complex. Our comparison of the two complexes in the presence and in the absence of the VHH-BS2 shows that, VHH-BS2 makes TorsinAΔE303-LULL1 unnecessarily stable and makes both Torsin-LULL1 complex and TorsinAΔE303-LULL1 complex to have a few new TorsinA-/TorsinAΔE303-LULL1 hydrogen bonds (such as between D277 and S279, H426 and K266, etc.) that were not present in the absence of VHH-BS2 (Table 2). The effects of VHH-BS2 on TorsinA-/TorsinAΔE303-LULL1 hydrogen bonds is easier to assess by comparing the ratio of the hydrogen bonds in the presence and in the absence of VHH-BS2 as shown in Table 2 as HbR1 compared to HbR2. The hydrogen bonds ratio 1, HbR1, is obtained by dividing the number of hydrogen bonds in TorsinAΔE303-LULL1 by the number of hydrogen bonds in TorsinA-LULL1 for the molecular systems with VHH-BS2. HbR2 is obtained in a similar way but from the molecular systems without VHH-BS2. For a given pair of residues, a value of HbR1 (or HbR2) greater than 1.0 signifies that there are more hydrogen bonds in TorsinAΔE303-LULL1 complex than in the TorsinA-LULL1 complex for that pair of residues (Table 2). The comparison of HbR1 and HbR2 helps in objectively assessing the effects of VHH-BS2 on TorsinA-/TorsinAΔE303-LULL1 hydrogen bonds such that HbR1 > HbR2 suggests that TorsinAΔE303-LULL1 complex has unreasonably high hydrogen bonds relative to TorsinA-LULL1 complex in the presence of VHH-BS2 than in the absence of VHH-BS2, which is the case for a number of the residue pairs shown in Tables 2. This may explain why the weakening of TorsinAΔE303-LULL1 interactions appears to be attenuated in the presence of VHH-BS2.

TABLE 2

Table 2. In the presence of VHH-BS2 nanobody, TorsinAΔE303-LULL1 complex has atypically high number of hydrogen bonds and show atypical stability^a.

In addition, we carried out enhanced sampling for TorsinA-/TorsinAΔE303-LULL1 complexes with or without VHH-BS2 using Gaussian accelerated Molecular Dynamics (GaMD) simulations (Miao et al., 2015) and constructed the energy landscape for each of the molecular systems studied using TorsinA-/TorsinAΔE303-LULL1's native contacts and TorsinA-/TorsinAΔE303-LULL1's hydrogen bonds as the reaction coordinates (Figure 6). The results show that VHH-BS2 attenuated the difference between TorsinA-LULL1 interactions and TorsinAΔE303-LULL1 interactions. This is evident from the comparable energy profiles in Figures 6A,B. On the other hand, the true difference between TorsinA-LULL1 interactions and TorsinAΔE303-LULL1 interactions becomes more observable in the absence of VHH-BS2 as shown in Figures 6C,D.

FIGURE 6

Figure 6. Comparisons of the energy landscapes of the interactions between TorsinA-LULL1 (left column) and TorsinAΔE303-LULL1 (right column) in the presence (top row) of and in the absence (bottom row) of VHH-BS2 nanobody. The energy landscapes were obtained from accelerated sampling using Gaussian accelerated Molecular Dynamics (GaMD) simulations. The energy profiles in (A,B) are very comparable because VHH-BS2 attenuated the difference between TorsinA-LULL1 interactions and TorsinAΔE303-LULL1 interactions. The true difference between TorsinA-LULL1 interactions and TorsinAΔE303-LULL1 interactions are more evident in (C,D) in the absence of VHH-BS2.

From the foregoing, it is apparent that the attenuation of TorsinA-TorsinAΔE303 differences by the nanobody-aided X-ray crystallization conditions (such as the co-crystallization of VHH-BS2 with TorsinA-/TorsinAΔE303- LULL1 complex) might have led to the conclusion of the previous study on the subtle nature of the differences between the healthy TorsinA from the diseased TorsinAΔE303 (Demircioglu et al., 2016). Our results suggest that the nanobody-aided crystallization condition for the X-ray crystallography could not allow the capturing of all the essential interactions between TorsinA (or TorsinAΔE303) and LULL1, and more importantly, the biologically-relevant dynamic interactions between the molecules could not possibly be accounted for, correctly and/or completely, by the X-ray crystallography's results alone as the conditions necessary for the crystallization and/or the presence of the nanobody stripped off some of the proteins' true biological properties and virtually all of their dynamic behaviors. In general, it may be very important that one exercises cautions when making biological/functional interpretations of structural data whenever the crystallization conditions (such as in some nanobody-aided X-ray crystallography) differ sharply from biological conditions because such sharp contractions in the conditions may influence the true biophysical and biochemical nature of the molecule of interest and the true intermolecular interactions between molecular components of the molecular complex of interest.

We acknowledge that a small structural change may signify a considerable functional change or a function loss and that this possibly has motivated structural biology studies. Nonetheless, it is important to recognize the potential limitations of the structural methods such as the possibility that the crystallization environment (for instance, the use of VHH-BS2 nanobody in the current example) may not accurately reflect the biological environment of the molecule(s) of interest and that such difference may affect the observed properties of the molecule(s) of interest, in addition to the inability to observe the dynamic behaviors of the molecule(s). We suggest that the observations from this study regarding the effects of the nanobody on the proteins of interest and their interactions be further investigated, for example on other protein-protein systems, protein-ligand systems, etc. whose structures are only solvable with nanobody-aided crystallography and would not crystalize in the absence of the nanobodies.

Conclusions

Genetic mutations resulting in the deletion of glutamic acid 303, ΔE303, from TorsinA (thereby weakening the binding between TorsinA and its activator, such as LULL1) primarily cause the (currently) incurable primary dystonia which is a neurodegenerative disorder characterized by prolonged muscle contractions, abnormal postures, and twisting movements. Building on previous works, we have shown the atomic details of TorsinA-LULL1 dynamic interactions and TorsinAΔE303-LULL1 dynamic interactions and their binding free energy changes/binding affinities. Our account of the quantitative value of how much the ΔE303 weakens the binding of TorsinAΔE303 to LULL1 as compared to the binding of TorsinA to LULL1 (ΔΔG_ΔE303−WT) has never been previously documented and constitutes one of the important contributions of this study to this field. Furthermore, we have presented a detailed account of TorsinAΔE303-LULL1 residue-residue interactions that are compromised by the ΔE303. The presented extensive effects of ΔE303 on TorsinAΔE303-LULL1 interactions show that the differences between TorsinA-LULL1 interactions and TorsinAΔE303-LULL1 interactions are non-subtle but well-pronounced. The complete loss or weakening of vital hydrogen bonds in TorsinAΔE303-LULL1 interactions considerably explains the weakening of the TorsinAΔE303-LULL1 binding (compared to TorsinAΔE303-LULL1 binding). We have also shown that nanobody-aided X-ray crystallization environment in the previous study attenuated the differences between TorsinA-LULL1 and TorsinAΔE303-LULL1 interactions. Pieces of evidence leading to the conclusion that the effects of ΔE303 (on the differences between TorsinAΔE303-LULL1 interactions and TorsinA-LULL1 interactions) are more pronounced than previously suggested have also been presented. Overall, our accounts of the dynamic interactions between “TorsinA and LULL1” and between “TorsinAΔE303 and LULL1” (such as, but not limited to, the patterns and the frequencies of hydrogen bonds in TorsinA-/TorsinAΔE303-LULL1 complexes) offer new insights for a better understanding of the molecular basis of Primary Dystonia and have long-term potentials of guiding the development of medications for the disease.

Materials and Methods

Structure of TorsinA/TorsinAΔE303-LULL1 Complex

We obtained X-ray crystallographic structures of TorsinA/TorsinAΔE303 and its activator (LULL1) in complex with the nanobody (VHH-BS2) from the Protein Data Bank (Berman et al., 2000), PDB IDs: 5j1s, 5j1t (Demircioglu et al., 2016).

Creation of the Initial Molecular Systems

Using tLeap/AmberTools16 (Pearlman et al., 1995; Case et al., 2005; Salomon-Ferrer et al., 2013) and ff14SB (Maier et al., 2015) force-fields, we created four explicitly solvated initial molecular systems for Molecular Dynamics (MD) simulations with AMBER16 (Pearlman et al., 1995; Case et al., 2005; Salomon-Ferrer et al., 2013): (A) TorsinA+LULL1, (B) TorsinAΔE303+LULL1, (C) TorsinA+LULL1+VHH-BS2, and (D) TorsinAΔE303+LULL1+VHH-BS2 (Figure 1). VHH-BS2 is a nanobody co-crystallized with TorsinA-LULL1 complex and TorsinAΔE303-LULL1 complex. The crystallization was not possible without VHH-BS2 (Demircioglu et al., 2016).

Energy Minimization

Each of the systems was energy-minimized using AMBER16 (Pearlman et al., 1995; Case et al., 2005; Salomon-Ferrer et al., 2013). The energy minimizations were done in three stages – weakly (0.5 kcal/mol/Ų) restraining all non-water atoms in the first stage, all alpha carbon atoms in the second stage, and without any restraints in the third stage. With weak restraints (0.5 kcal/mol/Ų) on alpha carbon atoms, each of the systems was gradually heated to 310K in canonical ensemble.

Molecular Dynamics (MD) Simulations

The systems' temperatures (Langevin thermostat Pastor et al., 1988 with a collision frequency of 2 ps⁻¹) and pressures (Berendsen barostat Berendsen et al., 1984) were controlled during the equilibration and production runs. Full electrostatic interactions energies were calculated using Particle Mesh Ewald method (Darden et al., 1993). A cutoff distance of 10 Å and a cubic spline switch function were used when calculating nonbonded interactions.

Overall, more than 8 μs of all-atom explicit solvent MD simulations were performed. Out of the entire MD simulations performed, 4.8 μs are unbiased production MD–TorsinA+LULL1: 500 ns * 3 replicates; TorsinAΔE303+LULL1: 500 ns * 3; TorsinA+LULL1+VHH-BS2: 300 ns * 3; and TorsinAΔE303+LULL1+VHH-BS2: 300 ns * 3. Thus, 300 ns *3 + 300 ns*3 + 500 ns*3 + 500 ns*3 = 4,800 ns = 4.8 μs. Gaussian accelerated MD (GaMD) simulations (Miao et al., 2015) for enhanced sampling made up 3.0 μs (i.e., 750 ns for each of the four systems), and Adaptively Biased MD (ABMD) simulations(Babin et al., 2008) made up 240 ns (120 ns for each of “TorsinA + LULL1” and “TorsinAΔE303 + LULL1” systems).

Binding Free Energy Change Calculations by ABMD

Binding free energy change calculations were done using both Adaptively Biased Molecular Dynamics (ABMD) simulations (Huber et al., 1994; Darve and Pohorille, 2001; Wang and Landau, 2001; Babin et al., 2008) and Molecular Mechanics Poisson-Boltzmann Surface Area (Kollman et al., 2000; Wang et al., 2001, 2006) (MM-PBSA).

For the adaptively biased molecular dynamics (ABMD) simulations (Huber et al., 1994; Darve and Pohorille, 2001; Wang and Landau, 2001; Babin et al., 2008), three sets of reaction coordinates/collective variables (CVs, Figure 7A) were used. They are (1) CV1 and CV2, (2) CV1 and CV3, and (3) CV2, and CV3. Where CV1 is a distance-based reaction coordinate (Figures 7A,B), CV2 is an angle-based reaction coordinate (Figures 7A,C), and CV3 is a torsion-angle-based reaction coordinate (Figures 7A,D). Each of the ABMD simulations was run until the observed Gibbs free energy change converged.

FIGURE 7

Figure 7. Reaction coordinates for calculating binding Gibbs free energy change through ABMD. In each of the panels/sub-figures, the smaller stripped structure represents the ligand, while the bigger unstripped structure represents the receptor. (A) The three reaction coordinates/collective variables (CV) used to define three pairs of CVs are shown as CV1, CV2, and CV3. CV1 is the distance between V (which is the center of mass, CoM, of T and U) and Y (which is the CoM of W and X). CV2 is the angle formed by T, U, and Y, while CV3 is the torsion angle formed by T, U, W, and X. Examples of changes in the position of the ligand relative to the receptor when CV1 is explored/sampled is shown in (B); when CV2 is explored is shown in (C); when CV3 is explored is shown in (D). D is made up by views from the top so as to easily demonstrate the effects of changes in the torsional angle (CV3) on the relative orientations of the ligand and the receptor. In each of the (A–D), the setup numbered “#1” depicts the initial setup.

Binding Free Energy Change Calculations by MM-PBSA

In addition to using ABMD, we calculated the binding free energy for each frame from the production MD simulations trajectory using Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) (Kollman et al., 2000; Wang et al., 2001, 2006). Although less rigorous than our ABMD approach, the MM-PBSA method we used is advantageous in allowing us to be able to isolate the contributions of each residue to the overall binding energy change. Analytical linearized Poisson-Boltzmann model (Luo et al., 2002; Sigalov et al., 2005, 2006; Wang and Luo, 2010) with the default parameters (except otherwise stated) was used through the MMPBSA.py (Miller et al., 2012) implemented for AmberTools16. The atomic radii specified in the ff14SB force-fields (Maier et al., 2015) and in the ion parameters for Ewald and TIP3P water (Joung and Cheatham III, 2008) were automatically read from the parameter and topology file. The solvent probe radius was also based on the force-fields parameters (Joung and Cheatham III, 2008; Maier et al., 2015). A surface tension of 0.0072 kcal/mol/Ų was used for calculating the nonpolar contribution to solvation free energy. An ionic strength of 0.150 nM was used for the PB solvent.

The binding free energy (ΔG_binding) accompanying the complexation of a receptor and a ligand to form a complex can be estimated from MM-PBSA as

$\begin{array}{ll}\Delta G_{binding}=\Delta H-T\Delta S\approx \Delta E_{MM}+\Delta G_{sol}-T\Delta & (1)\end{array}$

where ΔE_MM, ΔG_sol and –TΔS are respectively gas-phase Molecular Mechanics (MM) energy change, solvation free energy change, and conformational entropy change upon binding. While –TΔS can be computed by normal-mode analysis, ΔE_MM and ΔG_sol can be expressed as follows.

$\begin{array}{ll}\Delta E_{\text{MM}}=\Delta E_{\text{internal}}+\Delta E_{\text{electrostatics}}+\Delta E_{\text{vdw}}& (2)\end{array}$

$\begin{array}{ll}\Delta G_{\text{sol}}=\Delta G_{\text{PB}}+\Delta G_{\text{SA}}& (3)\end{array}$

where ΔG_PB is (the polar contribution to) electrostatic solvation energy calculated using Poisson Boltzmann (PB) model. ΔG_SA is (the non-electrostatic contributions to) the solvation component estimated by solvent accessible surface area.

Enhanced Sampling by Gaussian Accelerated MD (GaMD) Simulations

By adding harmonic boost potentials that smoothen the potential energy surface of a molecular system of interest, GaMD enhances the sampling of the conformational space of the molecular system (Miao et al., 2015). Given a molecular system with N atoms, such that the atomic coordinates is $r^{3^{*}N}=r_1^3,r_2^3,\dots ,r_N^3$, the GaMD algorithm modifies the molecular system's potential energy, $V_{original}(r^{3^{*}N})$, to $V_{updated}(r^{3^{*}N})$ by adding a boost potential, $V_{boost}(r^{3*N})$, whenever the original system's potential, $V_{original}(r^{3^{*}N})$, is lower than a given threshold, E_threshold. Mathematically,

$V_{updated}(r^{3\ast N})=\{\begin{array}{ll}{V}_{original}(r^{3\ast N})+V_{boost}(r^{3\ast N})\hfill & if{V}_{original}(r^{3\ast N})\hfill \\ \hfill & <E_{threshold}\hfill \\ V_{original}(r^{3\ast N})\hfill & else\hfill & (4)\end{array}$

where

$\begin{array}{ll}{V}_{boost}(r^{3*N})=\frac{1}{2}k{(E_{threshold}-V_{original}(r^{3*N}))}^2& (5)\end{array}$

where k is a harmonic force constant.

Both E_threshold and k can be adjusted while their optimal values can be automatically estimated from the regular unbiased MD simulation based on three criteria.

Criterion 1. $V_{boost}(r^{3*N})$ must be a monotonic function that guarantees that the relative rank/order of the biased potentials is the same as the relative rank/order of the original unbiased potentials, such that, given any two arbitrary potentials $V_{original\text{\_}1}(r^{3*N})$ and $V_{original\text{\_}2}(r^{3*N})$ on the original energy surface, if $V_{original\text{\_}1}(r^{3*N})<V_{original\text{\_}2}(r^{3*N})$ the function $V_{boost}(r^{3*N})$ must guarantee that $V_{updated\text{\_}1}(r^{3*N})<V_{updated\text{\_}2}(r^{3*N})$.

Criterion 2. The difference between the two updated potential energies must be less than the difference between the two original potential energies, such that $|V_{update{d}_2}(r^{3*N})-V_{update{d}_1}(r^{3*N})|<|V_{origina{l}_2}(r^{3*N})-V_{origina{l}_1}(r^{3*N})|$. The combination of criteria 1 and 2, and equations (4) and (5) suggest that equation (6) must be valid, which in turn requires that equation (7) must be valid and that 0 ≤ k₀ ≤ 1.

$\begin{array}{ll}{V}_{max}\le E_{threshold}\le V_{min}+\frac{1}{k}& (6)\end{array}$

$\begin{array}{ll}k=k_0*\frac{1}{V_{max}-V_{min}}\le \frac{1}{V_{max}-V_{min}}& (7)\end{array}$

Criterion 3: $V_{boost}(r^{3*N})$ must have a narrow, thus a small standard deviation (equation 8) to make the reweighing using cummulant expansion to the second order possible.

$\begin{array}{ll}{\sigma }_{V_{boost}}=k(E_{threshold}-V_{origina{l}_{mean}}(r^{3*N})){\sigma }_{V_{original}}\le {\sigma }_0& (8)\end{array}$

Where V_originalmean and σ_Voriginal are the mean and the standard deviation of the unbiased potential energies, σ_Vboost is the standard deviation of the boosted potential energy, and σ₀ is an upper limit specified by the user.

The k₀ in equation 7 can be computed as $k_0=min\left(1.0,\frac{{\sigma }_0}{{\sigma }_{V_{original}}}*\frac{V_{max}-V_{min}}{V_{max}-V_{origina{l}_{mean}}}\right)$.

GaMD Reweighing Using Cummulant Expansion

The probability along a reaction coordinate, RC(r), given as p′(RC), where r is the atomic coordinate $r_1^3,r_2^3,\dots ,r_N^3$, from the boost potential $V_{boost}(r^{3^{*}N})$ of each frame, p′(RC) can be reweighted to obtain the probability distribution, p(RC) of the canonical ensemble as shown in equation 9.

$\begin{array}{ll}p(R{C}_j)=p^{\prime }(R{C}_j)*\frac{{\langle e^{k_{B}T{V}_{boost}(r^{3*N})}\rangle }_j}{{\displaystyle {\sum }_{j=1}^M}{\langle e^{k_{B}T{V}_{boost}(r^{3*N})}\rangle }_j}forj=1,\dots ,M& (9)\end{array}$

The reweighted free energy can then be obtained as $F(R{C}_j)=-\frac{1}{k_{B}T}lnp(R{C}_j)$.

Residue-Residue Contact

Detail interactions between TorsinA (or TorsinAΔE303) and LULL1 were examined by residue-residue contact analysis wherein an amino acid of TorsinA/TorsinAΔE303 interacts with an amino acid of LULL1 if the distance between their centers of mass is within 6 Å.

Hydrogen Bonds Analysis

A hydrogen bond is recorded between two units whenever a hydrogen-bond acceptor (HBA) of one of the units is within 3.0 Å of a hydrogen-bond donor (HBD) of the other unit and the angle formed by HBA-Hydrogen-HBD is greater than or equal to 135°.

Data Availability

The initial structure of TorsinA/TorsinAΔE303 and its activator (LULL1) as well as the nanobody (VHH-BS2) are available in the Protein Data Bank (PDB IDs: 5j1s, 5j1t). Other dataset (such as Trajectories from Molecular Dynamics Simulations, etc.) upon which the results presented in the paper are based can be requested by sending an email to the corresponding author or to Bioinformatics Center at tools@Bioinformatics.Center.

Author Contributions

ES designed the study, carried out the study, analyzed and interpreted the data, and wrote the manuscript.

Conflict of Interest Statement

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The author acknowledges financial supports from Academia Sinica Taipei, Taiwan, and financial supports from National Tsing Hua University, Hsinhcu, Taiwan. The author thanks colleagues for their feedback on the work and its manuscript, and Professor David A. Case (of the Department of Chemistry & Chemical Biology, Rutgers University, United States) for providing an academic license of AMBER16.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2018.00064/full#supplementary-material

References

Andberg, T. A. H. (2011). Protein-protein Binding Affinities Calculated using the LIE Method. Master's thesis in chemistry, Faculty of science and technology, University of Tromsø.

Babin, V., Roland, C., and Sagui, C. (2008). Adaptively biased molecular dynamics for free energy calculations. J. Chem. Phys. 128:134101. doi: 10.1063/1.2844595

PubMed Abstract | CrossRef Full Text | Google Scholar

Berendsen, H. J. C., Postma, J. P. M., van Gunsteren, W. F., DiNola, A., and Haak, J. R. (1984). Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81, 3684–3690. doi: 10.1063/1.448118

PubMed Abstract | CrossRef Full Text | Google Scholar

Berman, H. M., Westbrook, J., Feng, Z., Gilliland, G., Bhat, T. N., Weissig, H., et al. (2000). The protein data bank. Nucleic Acids Res. 28, 235–242. doi: 10.1093/nar/28.1.235

CrossRef Full Text | Google Scholar

Breakefield, X. O., Blood, A. J., Li, Y., Hallett, M., Hanson, P. I., and Standaert, D. G. (2008). The pathophysiological basis of dystonias. Nat. Rev. Neurosci. 9, 222–234. doi: 10.1038/nrn2337

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown, R. S. H., Zhao, C., Chase, A. R., Wang, J., and Schlieker, C. (2014). The mechanism of torsin ATPase activation. Proc. Natl. Acad. Sci. U.S.A. 111, E4822–E4831. doi: 10.1073/pnas.1415271111

PubMed Abstract | CrossRef Full Text | Google Scholar

Case, D. A., Cheatham, T. E., Darden, T., Gohlke, H., Luo, R., Merz, K. M., et al. (2005). The amber biomolecular simulation programs. J. Comput. Chem. 26, 1668–1688. doi: 10.1002/jcc.20290

PubMed Abstract | CrossRef Full Text | Google Scholar

Darden, T., York, D., and Pedersen, L. (1993). Particle mesh Ewald: an N· log (N) method for Ewald sums in large systems. J. Chem. Phys. 98, 10089–10092. doi: 10.1063/1.464397

PubMed Abstract | CrossRef Full Text | Google Scholar

Darve, E., and Pohorille, A. (2001). Calculating free energies using average force. J. Chem. Phys. 115, 9169–9183. doi: 10.1063/1.1410978

PubMed Abstract | CrossRef Full Text | Google Scholar

Defazio, G., Abbruzzese, G., Livrea, P., and Berardelli, A. (2004). Epidemiology of primary dystonia. Lancet Neurol. 3, 673–678. doi: 10.1016/S1474-4422(04)00907-X

CrossRef Full Text | Google Scholar

Demircioglu, F. E., Sosa, B. A., Ingram, J., Ploegh, H. L., and Schwartz, T. U. (2016). Structures of TorsinA and its disease-mutant complexed with an activator reveal the molecular basis for primary dystonia. Elife 5:e17983. doi: 10.7554/eLife.17983

PubMed Abstract | CrossRef Full Text | Google Scholar

Dystonia in Europe (ESDE) Collaborative Group and others, (2000). A prevalence study of primary dystonia in eight European countries. J. Neurol. 247, 787–792. doi: 10.1007/s004150070094

PubMed Abstract | CrossRef Full Text

Fahn, S. (1988). Concept and classification of dystonia. Adv. Neurol. 50, 1–8. doi: 10.1212/WNL.50.5_Suppl_5.S1

PubMed Abstract | CrossRef Full Text | Google Scholar

Genheden, S., Nilsson, I., and Ryde, U. (2011). Binding affinities of factor Xa inhibitors estimated by thermodynamic integration and MM/GBSA. J. Chem. Inf. Model. 51, 947–958. doi: 10.1021/ci100458f

PubMed Abstract | CrossRef Full Text | Google Scholar

Gohlke, H., Kiel, C., and Case, D. A. (2003). Insights into protein–protein binding by binding free energy calculation and free energy decomposition for the Ras–Raf and Ras–RalGDS complexes. J. Mol. Biol. 330, 891–913. doi: 10.1016/S0022-2836(03)00610-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanson, P. I., and Whiteheart, S. W. (2005). AAA+ proteins: have engine, will work. Nat. Rev. Mol. Cell Biol. 6, 519–529. doi: 10.1038/nrm1684

PubMed Abstract | CrossRef Full Text | Google Scholar

Homeyer, N., and Gohlke, H. (2012). Free energy calculations by the molecular mechanics poisson- boltzmann surface area method. Mol. Inform. 31, 114–122. doi: 10.1002/minf.201100135

PubMed Abstract | CrossRef Full Text | Google Scholar

Huber, T., Torda, A. E., and Gunsteren, W. F. (1994). Local elevation: a method for improving the searching properties of molecular dynamics simulation. J. Comput. Aided Mol. Des. 8, 695–708. doi: 10.1007/BF00124016

PubMed Abstract | CrossRef Full Text | Google Scholar

Joung, I. S., and Cheatham III, T. E. (2008). Determination of alkali and halide monovalent ion parameters for use in explicitly solvated biomolecular simulations. J. Phys. Chem. B 112, 9020–9041. doi: 10.1021/jp8001614

PubMed Abstract | CrossRef Full Text | Google Scholar

Kiel, C., Selzer, T., Shaul, Y., Schreiber, G., and Herrmann, C. (2004). Electrostatically optimized ras-binding ral guanine dissociation stimulator mutants increase the rate of association by stabilizing the encounter complex. Proc. Natl. Acad. Sci. U.S.A. 101, 9223–9228. doi: 10.1073/pnas.0401160101

PubMed Abstract | CrossRef Full Text | Google Scholar

Kojovic, M., Pareés, I., Kassavetis, P., Palomar, F. J., Mir, P., Teo, J. T., et al. (2013). Secondary and primary dystonia: pathophysiological differences. Brain 136, 2038–2049. doi: 10.1093/brain/awt150

PubMed Abstract | CrossRef Full Text | Google Scholar

Kollman, P. A., Massova, I., Reyes, C., Kuhn, B., Huo, S., Chong, L., et al. (2000). Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models. Acc. Chem. Res. 33, 889–897. doi: 10.1021/ar000033j

PubMed Abstract | CrossRef Full Text | Google Scholar

Konakova, M., and Pulst, S. M. (2005). Dystonia-associated forms of torsinA are deficient in ATPase activity. J. Mol. Neurosci. 25, 105–117. doi: 10.1385/JMN:25:1:105

PubMed Abstract | CrossRef Full Text | Google Scholar

Leung, J., Klein, C., Friedman, J., Vieregge, P., Jacobs, H., Doheny, D., et al. (2001). Novel mutation in the TOR1A (DYT1) gene in atypical, early onset dystonia and polymorphisms in dystonia and early onset parkinsonism. Neurogenetics 3, 133–143. doi: 10.1007/s100480100111

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, R., David, L., and Gilson, M. K. (2002). Accelerated Poisson–Boltzmann calculations for static and dynamic systems. J. Comput. Chem. 23, 1244–1253. doi: 10.1002/jcc.10120

PubMed Abstract | CrossRef Full Text | Google Scholar

Maier, J. A., Martinez, C., Kasavajhala, K., Wickstrom, L., Hauser, K. E., and Simmerling, C. (2015). ff14SB: improving the accuracy of protein side chain and backbone parameters from ff99SB. J. Chem. Theory Comput. 11, 3696–3713. doi: 10.1021/acs.jctc.5b00255

PubMed Abstract | CrossRef Full Text | Google Scholar

Marsden, C. D., Obeso, J. A., Zarranz, J. J., and Lang, A. E. (1985). The anatomical basis of symptomatic hemidystonia. Brain 108, 463–483. doi: 10.1093/brain/108.2.463

PubMed Abstract | CrossRef Full Text | Google Scholar

Miao, Y., Feher, V. A., and McCammon, J. A. (2015). Gaussian accelerated molecular dynamics: unconstrained enhanced sampling and free energy calculation. J. Chem. Theory Comput. 11, 3584–3595. doi: 10.1021/acs.jctc.5b00436

PubMed Abstract | CrossRef Full Text | Google Scholar

Miller, B. R. III., McGee Jr, T. D., Swails, J. M., Homeyer, N., Gohlke, H., and Roitberg, A. E. (2012). MMPBSA. py: an efficient program for end-state free energy calculations. J. Chem. Theory Comput. 8, 3314–3321. doi: 10.1021/ct300418h

PubMed Abstract | CrossRef Full Text | Google Scholar

Müller, J., Kiechl, S., Wenning, G. K., Seppi, K., Willeit, J., Gasperi, A., et al. (2002). The prevalence of primary dystonia in the general community. Neurology 59, 941–943. doi: 10.1212/01.WNL.0000026474.12594.0D

PubMed Abstract | CrossRef Full Text | Google Scholar

Naismith, T. V., Dalal, S., and Hanson, P. I. (2009). Interaction of torsinA with its major binding partners is impaired by the dystonia-associated ΔGAG deletion. J. Biol. Chem. 284, 27866–27874. doi: 10.1074/jbc.M109.020164

PubMed Abstract | CrossRef Full Text | Google Scholar

Németh, A. H. (2002). The genetics of primary dystonias and related disorders. Brain 125, 695–721. doi: 10.1093/brain/awf090

PubMed Abstract | CrossRef Full Text | Google Scholar

Noskov, S. Y., and Lim, C. (2001). Free energy decomposition of protein-protein interactions. Biophys. J. 81, 737–750. doi: 10.1016/S0006-3495(01)75738-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Ozelius, L. J., Hewett, J. W., Page, C. E., Bressman, S. B., Kramer, P. L., Shalish, C., et al. (1997). The early-onset torsion dystonia gene (DYT1) encodes an ATP-binding protein. Nat. Genet. 17, 40–48. doi: 10.1038/ng0997-40

PubMed Abstract | CrossRef Full Text | Google Scholar

Pastor, R. W., Brooks, B. R., and Szabo, A. (1988). An analysis of the accuracy of Langevin and molecular dynamics algorithms. Mol. Phys. 65, 1409–1419. doi: 10.1080/00268978800101881

CrossRef Full Text | Google Scholar

Pearlman, D. A., Case, D. A., Caldwell, J. W., Ross, W. S., Cheatham, T. E., DeBolt, S., et al. (1995). AMBER, a package of computer programs for applying molecular mechanics, normal mode analysis, molecular dynamics and free energy calculations to simulate the structural and energetic properties of molecules. Comput. Phys. Commun. 91, 1–41. doi: 10.1016/0010-4655(95)00041-D

CrossRef Full Text | Google Scholar

Pham, P., Frei, K. P., Woo, W., and Truong, D. D. (2006). Molecular defects of the dystonia-causing torsinA mutation. Neuroreport 17, 1725–1728. doi: 10.1097/WNR.0b013e3280101220

PubMed Abstract | CrossRef Full Text | Google Scholar

Rose, A. E., Brown, R. S. H., and Schlieker, C. (2015). Torsins: not your typical AAA+ ATPases. Crit. Rev. Biochem. Mol. Biol. 50, 532–549. doi: 10.3109/10409238.2015.1091804

PubMed Abstract | CrossRef Full Text | Google Scholar

Sadnicka, A., Hoffland, B. S., Bhatia, K. P., de Warrenburg, B. P., and Edwards, M. J. (2012). The cerebellum in dystonia–help or hindrance? Clin. Neurophysiol. 123, 65–70. doi: 10.1016/j.clinph.2011.04.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Salawu, E. O. (2016). RaFoSA: random forests secondary structure assignment for coarse-grained and all-atom protein systems. Cogent Biol. 2:1214061. doi: 10.1080/23312025.2016.1214061

CrossRef Full Text | Google Scholar

Salomon-Ferrer, R., Case, D. A., and Walker, R. C. (2013). An overview of the Amber biomolecular simulation package. Wiley Interdiscip. Rev. Comput. Mol. Sci. 3, 198–210. doi: 10.1002/wcms.1121

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharma, N., Baxter, M. G., Petravicz, J., Bragg, D. C., Schienda, A., Standaert, D. G., et al. (2005). Impaired motor learning in mice expressing torsinA with the DYT1 dystonia mutation. J. Neurosci. 25, 5351–5355. doi: 10.1523/JNEUROSCI.0855-05.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Sigalov, G., Fenley, A., and Onufriev, A. (2006). Analytical electrostatics for biomolecules: beyond the generalized born approximation. J. Chem. Phys. 124:124902. doi: 10.1063/1.2177251

PubMed Abstract | CrossRef Full Text | Google Scholar

Sigalov, G., Scheffel, P., and Onufriev, A. (2005). Incorporating variable dielectric environments into the generalized Born model. J. Chem. Phys. 122:94511. doi: 10.1063/1.1857811

PubMed Abstract | CrossRef Full Text | Google Scholar

Steeves, T. D., Day, L., Dykeman, J., Jette, N., and Pringsheim, T. (2012). The prevalence of primary dystonia: a systematic review and meta-analysis. Mov. Disord. 27, 1789–1796. doi: 10.1002/mds.25244

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, F., and Landau, D. P. (2001). Efficient, multiple-range random walk algorithm to calculate the density of states. Phys. Rev. Lett. 86, 2050–2053. doi: 10.1103/PhysRevLett.86.2050

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Hou, T., and Xu, X. (2006). Recent advances in free energy calculations with a combination of molecular mechanics and continuum models. Curr. Comput. Aided Drug Des. 2, 287–306. doi: 10.2174/157340906778226454

CrossRef Full Text | Google Scholar

Wang, J., and Luo, R. (2010). Assessment of linear finite-difference Poisson-Boltzmann solvers. J. Comput. Chem. 31, 1689–1698. doi: 10.1002/jcc.21456

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W., Donini, O., Reyes, C. M., and Kollman, P. A. (2001). Biomolecular simulations: recent developments in force fields, simulations of enzyme catalysis, protein-ligand, protein-protein, and protein-nucleic acid noncovalent interactions. Annu. Rev. Biophys. Biomol. Struct. 30, 211–243. doi: 10.1146/annurev.biophys.30.1.211

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, C., Brown, R. S. H., Chase, A. R., Eisele, M. R., and Schlieker, C. (2013). Regulation of Torsin ATPases by LAP1 and LULL1. Proc. Natl. Acad. Sci. U.S.A. 110, E1545–E1554. doi: 10.1073/pnas.1300676110

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, L., Millen, L., Mendoza, J. L., and Thomas, P. J. (2010). A unique redox-sensing sensor II motif in TorsinA plays a critical role in nucleotide and partner binding. J. Biol. Chem. 285, 37271–37280. doi: 10.1074/jbc.M110.123471

PubMed Abstract | CrossRef Full Text | Google Scholar

Zirn, B., Grundmann, K., Huppke, P., Puthenparampil, J., Wolburg, H., Riess, O., et al. (2008). Novel TOR1A mutation p. Arg288Gln in early-onset dystonia (DYT1). J. Neurol. Neurosurg. Psychiatry 79, 1327–1330. doi: 10.1136/jnnp.2008.148270

PubMed Abstract | CrossRef Full Text | Google Scholar