An Allosteric Pathway in Copper, Zinc Superoxide Dismutase Unravels the Molecular Mechanism of the G93A Amyotrophic Lateral Sclerosis-Linked Mutation.

Several different mutations of the protein copper, zinc superoxide dismutase (SOD1) produce the neurodegenerative disorder amyotrophic lateral sclerosis (ALS). The molecular mechanism by which the diverse mutations converge to a similar pathology is currently unknown. The electrostatic loop (EL) of SOD1 is known to be affected in all of the studied ALS-linked mutations of SOD1. In this work, we employ a multiscale simulation approach to show that this perturbation corresponds to an increased probability of the EL detaching from its native position, exposing the metal site of the protein to water. From extensive atomistic and coarse-grained molecular dynamics (MD) simulations, we identify an allosteric pathway that explains the action of the distant G93A mutation on the EL. Finally, we employ quantum mechanics/molecular mechanics MD simulations to show that the opening of the EL decreases the Zn(II) affinity of the protein. As the loss of Zn(II) is at the center of several proposed pathogenic mechanisms in SOD1-linked ALS, the structural effect identified here not only is in agreement with the experimental data but also places the opening of the electrostatic loop as the possible main pathogenic effect for a significant number of ALS-linked SOD1 mutations.

The protein copper, zinc superoxide dismutase (SOD1) catalyzes the dismutation of the superoxide radical into molecular oxygen and peroxide.[1] It protects the cell from the radical, which is naturally formed during cell respiration. Human SOD1 is a dimeric protein, which binds a Zn(II) and a Cu(I/II) ion per monomer (see Figure A,B).

Figure 1

Atomistic structure and dynamics of SOD1. (A) Ribbon view of the crystallographic structure of dimeric SOD1.[27] The EL is colored blue; the G93 residue is shown as red sticks, and Cu(I) and Zn(II) ions are shown as orange and gray spheres, respectively. (B) Details of the metal site. Oxygen atoms are colored red, and nitrogen atoms blue. (C and D) Superimposed frames colored from blue to white to red for the first eigenvector of the covariance matrix for WT and G93A mutant MD trajectories, respectively, where the EL opening can be seen.

Well more than 100 mutations, scattered along the whole sequence of SOD1, are linked to familial forms of the neurodegenerative disease amyotrophic lateral sclerosis (ALS).[2] The study of SOD1-linked cases of ALS is expected to shed light on the pathogenic mechanism of the most prevalent sporadic form of the disease. The pathogenic mechanisms of the ALS-linked mutations of SOD1 have not been elucidated in detail. The current data indicate that a gain of a pathological function is responsible for the toxicity of the mutant, rather than a loss of enzymatic activity.[3] Several hypotheses for the pathological function exist: SOD1 misfolding could interfere with cellular metabolism.[4−8] A distortion of the active site of the protein could cause aberrant redox catalysis and cellular oxidative stress.[9] Finally, a loss of metal ions, particularly Zn(II), from the SOD1 active site, could cause a harmful metal dyshomeostasis.[10−13]

A loss of Zn(II) from the active site of SOD1 appears to play a critical role in the hypotheses for SOD1 pathogenicity:[14] It destabilizes the protein’s structure[15−17] and alters its active site,[18,19] and mutant SOD1 has been shown to cause Zn(II) dyshomeostasis in mice.[20] Close to the Zn(II) site of SOD1 is the electrostatic loop (EL), which is known to be destabilized in several ALS-linked mutants.[21] From a physicochemical point of view, the diversity of the pathogenic mutations, regarding both the amino acid substitution and its location in the protein structure, is intriguing. Because demetalation has been shown to be important for the misfolding process,[14] computational investigation of pathogenic mutants of SOD1 has been preferentially performed on the zinc-deficient apoenzyme,[7,8,18] confirming that ALS-linked mutations also increase the level of oligomerization. Mera-Adasme and collaborators have investigated the effect of SOD1 mutations on the metal site by means of computational chemistry.[22,23] Their focus was exclusively on mutations in the proximity of the metal site, where altered Zn(II) binding or interactions with the EL determined the impact of the mutations. The effect of non-active site mutations of SOD1 on the metal affinity of the protein and the underlying molecular mechanisms are still unclear. To the best of our knowledge, neither aspect has been elucidated.

This work focuses on the G93A mutation of SOD1, a widely studied mutation, as the expression of its human form in mice constitutes a common animal model of ALS.[24] Because the mutation occurs away from the metal site (Figure , top), it must affect the metal site allosterically. While the G93A mutant is experimentally known to affect the EL and the metal site of the protein,[20,21,25] the effect is subtle, as only a methyl group is added, and the purified variant has also been shown to be metalated and catalytically active.[26] We use atomistic and coarse-grained (CG) molecular dynamics (MD) simulations of the SOD1 wild type (WT) and its G93A mutant to compare their EL loop opening. Moreover, to evaluate the impact of the EL loop opening on the metal site, we also perform hybrid quantum mechanical/molecular mechanical (QM/MM) MD simulations. The thorough statistics from the CG simulations allow us to identify for the first time the allosteric pathway from the distant G93A mutation to the metal site of SOD1.

We performed atomistic MD simulations on models for the monomer of the wild type (WT) and G93A variants of human SOD1, both derived from the high-resolution crystallographic structure for the WT enzyme.[27] Three 1 μs simulations were performed for each variant and show a stable protein fold. In one of the mutant trajectories and two of the WT ones, full opening of the EL could be observed (Figures S4 and S5). During the opening, weak interactions between one half of the loop and the rest of the protein are broken, and the detached loop adopts a distorted configuration (Figure C,D). The “open” half of the loop matches the amino acid sequence previously shown by H/D exchange to have enhanced hydration in ALS-linked mutants of SOD1.[21] However, the obtained atomistic simulation times are insufficient for judging whether there is an altered tendency for loop opening in G93A SOD1.

To obtain reliable statistics, we performed CG MD simulations with the Martini force field. To the best of our knowledge, this is not only the first exploration of the combined Martini 3 (open-beta version[28−30]) and Go̅-like models[31] but also among the first uses of Martini for the study of a conformational transition in a protein.[28,32] A total of 480 μs of simulation time (12 × 40 μs) was obtained for each variant (WT and G93A) of SOD1. The CG simulations revealed similar EL opening events compared to the atomistic simulations, indicating that the CG model reliably reproduces the protein flexibility. Figure summarizes the key findings of the CG simulations. The distribution of the average protein backbone root-mean-square deviation (RMSD) depicted in Figure A has a similar maximum position for the WT (blue) and the G93A mutant (red). However, above an RMSD of 2.3 Å, a clear shoulder is present in the mutant that is missing in the WT. Thus, the mutation results in an increased population at a higher RMSD. Figure B shows that the EL exhibits the highest RMSD of ∼4.5 Å in the case of the WT, indicating that the shoulder at the high RMSD might impact the EL to a large extent. Interestingly, the vicinity of the mutation site in the G93A mutant becomes rigid compared to that of the WT (see Figure C). In contrast, the EL exhibits an increased flexibility in the mutant, which is in agreement with the experimentally reported increased mobility of the SOD1 EL.[21] Moreover, the increased flexibility of the Zn loop (ZL) as well as the unaffected dimerization surface are in agreement with observations from H/D exchange experiments.[25]

Figure 2

Comparison of the SOD1 WT and G93A mutant based on the CG Martini simulations. (A) Distribution of the average protein backbone RMSD for the WT (blue) and the G93A mutant (red) obtained from 480 μs of simulation each. The inset shows the region between 1.75 and 2.55 Å including error bars. (B) Flexibility of the protein backbone of the WT. Snapshots were taken every 1 μs; the color scale represents the backbone RMSD. (C) Change in RMSD between the WT and G93A. Blue indicates rigidification in G93A; red indicates increased flexibility. (D and E) Free energy surfaces of WT and G93A, respectively, spanned by the EL RMSD and the water–active site contacts (top left). The underlying probability distribution of the EL RMSD (top right) and the number of water–active site contacts (bottom) are also depicted. The differences in the probability distributions are small, but meaningful as shown by the error bars in Figure S6. For the sake of clarity, the error bars have been omitted.

The G93A mutation has been shown to impair the Zn metabolism in the cell and to destabilize the metal site.[6,20,25] We had previously proposed an increased level of hydration of the metal site via EL opening, resulting in a decreased affinity of SOD1 for Zn(II), as a pathogenic mechanism for ALS-linked SOD1 mutants.[22,33] We analyzed the number of contacts between the active site of SOD1 and CG water beads. The G93A mutant shows more water–active site contacts than the WT (bottom panels in Figure D,E). Note that one CG water bead represents four water molecules. Thus, 1.5 CG water beads are sufficient to form a hexaaqua complex with the Cu(I/II) or Zn(II) ions. To resolve the connection between the water–active site contacts and the RMSD of the EL, we calculated the corresponding two-dimensional free energy surface for the WT (Figure D) and the G93A mutant (Figure E). In both cases, the global minimum is located at a low RMSD (∼2.5 Å) and ∼70 water–active site contacts. A closer look at the tail of the minimum at a higher RMSD reveals that the G93A mutant exhibits a wider and slightly deeper minimum at an RMSD of ∼6 Å and ∼80 water–active site contacts. Overall, our CG simulations show a stronger tendency for EL opening in the G93A mutant, in agreement with experimental observations, together with a higher level of hydration of the active site.

To determine whether the increased level of hydration of the metal site has an effect on the binding energy of Zn(II), implicated in proposed mechanisms for SOD1 pathogenic function, we performed QM/MM MD simulations starting from structures taken from the atomistic trajectories (see Figure A,B). The QM subsystem, treated at the composite HF-3c level of theory, contained 590 atoms, in line with what we previously showed to be adequate[34] (see the Supporting Information). Four simulations of 4 ps each were performed: with the EL in open and closed configurations and, for each of them, with and without Zn(II) bound to the protein. Subtraction of the potential energy averages for the trajectories can be assimilated to the difference in ΔΔH of Zn(II) binding between the protein in open and closed EL conformations. The closed EL conformation has a more favorable ΔH of Zn(II) binding by −99 kJ/mol. A frequency-based estimation for the ΔΔS for the process, employing the full protein plus the water molecules within ∼10 Å of the metal site (see the Supporting Information), gave a −TΔΔS value of 15 kJ/mol, close to an order of magnitude lower than the calculated ΔΔH.

Figure 3

(A and B) QM/MM-optimized structures for the G93A mutant in the closed and open EL conformations, respectively. The EL is colored blue; the G93 residue is shown as red sticks, and Cu(I) and Zn(II) ions are shown as orange and gray spheres, respectively. (C) Matrix representation of the integrated absolute difference in the distance distributions between all backbone beads. The bottom left triangle represents the full data set. In the top right triangle, only values of >0.3 are depicted; all other values are colored blue. (D) Residues exhibiting an integrated absolute difference between the distributions of the distance of the respective residue and the mutation site of WT and G93A of >0.3 are colored red. The red bonds between the backbone chain indicate transfer pathways based on the Go̅-like model.

To unravel the pathway transmitting the subtle mutation G93A toward the EL, we analyzed the distance distributions between all pairwise residue combinations in the CG trajectories (as described in Figures S3 and S8A,B). The integrated absolute differences between the distance distributions of the WT and the G93A mutant are depicted in Figure C. The largest changes are observed between the region of residues 90–95 (around the mutation site) and all other residues. In Figure D, the residues of the SOD1 monomer for which the integrated absolute difference between the distance distributions of mutation site 93 and all other residues is larger than 0.3 are colored red. Besides the vicinity of the mutation site, several residues connecting the mutation site with the joint for the EL opening (residues 122 and 141) exhibit significant differences in the distance distributions (see Figure S8C,D).

A closer look at the atomistic structure reveals a network comprising hydrogen and peptide bonds that are responsible for the mechanical communication of G93 with the EL. The backbone of L38 is closely packed with G93 in the WT structure and directly linked to the first residue in the network, T39. The packing is altered in the mutant by a steric clash between Cβ of A93 and the backbone oxygen of L38. This clash can also be observed in our atomistic simulations resulting in a change in the G/A93(Cβ3)–L38(O) distance distribution (see Figures S9 and S10). The identified network promotes this altered distance via three hydrogen and two peptide bonds toward the joint of the EL. Namely, residues T39, H43, K122, A123, A140, and G141 are involved in the network and all exhibit significant changes in their distributions (Figure S8). The allosteric pathway thus identified is comprised mainly of backbone-mediated interactions, which explains its lack of susceptibility to ALS-causing mutations, as previously discussed.[22] The exception is the H43 residue, which participates in the interactions with its side chain and is the target of at least one known ALS-causing mutation.[35]

In summary, our multiscale simulation study provides a molecular mechanism for the pathogenic effect of the widely studied G93A mutation of SOD1. Our results directly link the G93A mutation to a Zn(II) site destabilization via an increased level of electrostatic loop opening. It arises from a steric clash between G93 and L38 in the mutant. We identified a network connecting G93 with the main point of EL opening, which is responsible for transmitting the steric clash to the EL. In a more general view, our results back the general involvement of the EL in SOD1-linked ALS pathogenesis by providing molecular level insights into EL destabilization. Moreover, the results also mark a milestone in the combined use of the Martini 3 force field and Go̅-like models to study long-range structural communication and allosteric changes caused by single-point mutations.

Computational Methods

All MD simulations were performed with the MD package GROMACS (version 2018).[36,37] Atomistic simulations were realized with the AMBER99SB-ILDN force field for the protein,[38,39] while metal parameters were obtained in a previous work.[40] The open-beta version of the Martini 3 force field[28,30] in combination with a Go̅-like model[31] was used for the CG MD simulations. The beta release of Martini 3 is available online at the webpage of Marrink’s group.[29] QM/MM calculations were performed with the program packages pDynamo,[41] ORCA 4.0,[42] Turbomole version 7.3,[43,44] and xtb.[45] The HF-3c composite method and GFN0-xTB semiempirical method from Grimme[46,47] were employed for the QM part of the system, while the AMBER-ILDN force field[38,39] was employed for the MM subsystem. Further details are given in the Supporting Information.