The crystal structure 2V0A (Strange et al., 2007) was recruited. This structure did not contain any mutations and captured in the dimeric form which represent the active form the wild-type SOD1 (Fridovich, 1975). We have utilized the dimeric form in our studies, which holds two metal ions. The metal ions, zinc and copper, were kept in the holo structures as Z n 2 + and C u 2 + , representing the second ground dimer Prior to simulations, multiple occupancies were eliminated from the structure. In silico mutagenesis was applied to obtain five different mutations, namely H63A, H63R, K136A, G37R and H46R/H48D. For every wild-type and mutant dimers, four distinct structures were generated: 1) apo dimer (non-metallated) with an oxidized disulfide bond (C57-C146), 2) apo dimer (non-metallated) with a reduced disulfide bond, 3) holo dimer (metallated) with an oxidized disulfide bond and 4) holo dimer (metallated) with a reduced disulfide bond. Overall, 24 different SOD1 dimers were generated. All of the structures were protonated according to pH 7.0 by PROPKA (Li et al., 2005; Søndergaard et al., 2011). Protonated structures were solvated in a cubic box of water modeled by TIP3P (Jorgensen et al., 1983). Solvated systems were neutralized by sodium and chloride counter-ions.

Resulting systems containing wild-type and SOD1 mutant enzymes were initially energy minimized by the steepest descent approach without any constraints. Molecular dynamics simulations were performed for all of the 24 dimeric SOD1 by using Gromacs 5.1.4 (Berendsen et al., 1995; Van Der Spoel et al., 2005; Pronk et al., 2013; Abraham et al., 2015) and adopting the CHARMM36 force field Phillips et al. (2005); MacKerell and Bashford, (1998); Brooks et al. (2009); Huang and MacKerell (2013). CHARMM was previously used for the SOD1 dimer producing results that were well aligned with experiments (Strange et al., 2007). Particle Mesh Ewald (PME) method was used for the electrostatic energy calculation (Darden et al., 1993). The time step was kept at 2 fs for all simulations. A cut-off distance of 10 Å was implemented for prediction the short-range interactions. Initially, the systems were slowly heated to the temperature of 298 K over 100 ps After stabilization of the system temperature at 298 K, an additional 100 ps of simulation was performed by using NPT ensembles at the constant pressure of 1 bar. The productions simulations were also run using NPT ensembles and at two different temperatures 298 and 400 K for 100 ns Systems were directly heated to 400 K. Overall 48 different MD simulations were performed. Trajectory analysis were done by monitoring backbone displacements and fluctuations, radius of gyration (R _G ) of the dimer and metal coordination states of both monomers. VMD was used for visual inspections of the structures and trajectories (Humphrey et al., 1996).

To date, 217 mutations in the coding region of SOD1 gene have been identified to be associated with ALS (https://alsod.ac.uk/) (Lill et al., 2011). A large portion of these mutations (184/217) are missense mutations which have been scattered throughout the SOD1 sequence (Figure 1A) and structure (Figure 1B). Many studies on SOD1 variants have converged that demetallation of SOD1 structure can significantly alter stability of the enzyme inducing formation of misfolded neurotoxic aggregates (Tiwari et al., 2009; Banci et al., 2007; Ding and Dokholyan, 2008; Wright et al., 2019; Wang et al., 2002, 2003, 2007). From this respect, metal-binding SOD1 mutations that directly affect the stability of the enzyme have been particularly well-studied in the aim of understanding the SOD1 toxicity in ALS (Wright et al., 2019).

Except for the H63 which acts as the bridging ion between metal ions, all of the amino acids that are directly involved in metal coordination; H46, H48 and H120 for Cu²⁺ and H71, H80 and D120 for Zn²⁺, were mutated at least in one ALS case (Figure 2A). H63 in the imidazolate form coordinates both of the metals, Cu²⁺ or Zn²⁺, through its δN and ϵN respectively. Despite its potential contribution to the integrity of both of the metallic centers, H63 was not mutated in any of the ALS cases documented so far (Abel et al., 2012). Furthermore, K136 which is not directly involved in metal coordination but found in the first shell (<5Å) was not linked to ALS. Similar to H63, K136 also is located at the center of the bimetallic center (Figure 2A). Along with H63, another bridge in the second shell, D124 exists forming hydrogen bonds with the H46 of Cu-shell and H71 of Zn-shell. D124 was previously observed to be mutated in ALS (Figure 1A) (Mera-Adasme et al., 2013). While all amino acids that were directly involved in metal coordination of the Cu/Zn SOD1 have been reported to be somehow associated with ALS, two amino acids, H63 and K136, that are located in the fist shell of the bimetallic center have not been linked to this disease. Hence despite the absence of epidemiological evidence bridging these positions with ALS, we appraise that any substitutions at these positions have the potential to affect the metallation of the enzyme leading to misfolding, in the same manner other metal binding mutants did. A similar prediction has been made for another metal binding mutant, H80R, for which compelling in vitro evidence suggesting a dramatic effect on the metalation and stability of SOD1 came first and clinical findings bridging this mutation to ALS was established, then (Wright et al., 2019).

Here, we have investigated these two non-ALS linked positions in three variants by MD simulations. We have also included a double mutant of metal-binding, H46R/H48D and a wild-type-like mutant, G37R, in our simulations. Particularly, the mutant, H46R/H48D that was mutated in two of the Cu-binding positions, has been expected to lose its metals, particularly Cu, parallel to previous findings on the Cu-binding variants (Winkler et al., 2009; Lelie et al., 2011). G37R, on the other hand, was not classified as a metal-binding mutant as G37 has been located far away from the bimetallic center (Lynch et al., 2004). Nevertheless, investigations of G37R SOD1 suggested altered metallation in this variant (Milardi et al., 2010). Hence, alongside with the wild-type SOD1, overall six SOD1 variants were recruited to our study to test whether the putative variants show similar structural dynamics to any of these ALS-linked mutants.

Cu/Zn SOD1 catalyzes the conversion of two superoxide radical ions, O 2 . − to hydrogen peroxide, H ₂ O ₂ and molecular oxygen, O ₂ (Figure 2B). This reaction is a disproportionation (dismutation) reaction in which simultaneous oxidation and reduction of the O 2 . − radical take place by the active site of the Cu-center (Pelmenschikov and Siegbahn, 2005). According to the proposed mechanism (Figure 2B), the first enzyme-substrate (ES) complex is transformed into the first intermediate by the binding of one of O 2 . − to the Cu(II) center. This intermediate is then broken down to the second ES complex through a homolytic cleavage of the Cu-O bond, Cu(II) center is reduced to Cu(I), O 2 . − is oxidized to O ₂ and H63 loses its contact with the Cu-center gaining a proton at its δN. The second ES complex is formed through binding of the reduced Cu (I) center to a second O 2 . − . In the final step, the Cu (I) center catalyzes the reduction of O 2 . − to H ₂ O ₂ transforming back into the first ES complex.

Given this reaction mechanism (Figure 2B), we stress the importance of the accurate force field parameters for the unique imidazolate-bridged bimetallic center. For this purpose, a number of studies have reported ab initio parameters of the Cu/Zn center by quantum mechanical (QM) calculations (Shen et al., 1990; Strange et al., 2007; Peng et al., 2018). Otherwise standard force fields do not provide the parameters for this bimetallic center particularly for the deprotonated H63 (Vanommeslaeghe et al., 2010). Here we defined the SOD1 variants through the CHARMM force field (Strange et al., 2007) but did not include the ab initio parameters for the metallic center (Shen et al., 1990; Strange et al., 2007; Peng et al., 2018). Because the QM parameters were developed for the wild-type metallic center (Shen et al., 1990; Strange et al., 2007; Peng et al., 2018), whilst half of the variants to be analyzed here were H63 mutants in which the imidazolate ion was deleted. Thus, we considered that usage of the first ground state of the SOD1 for all variants would lead to biased results due to incorporation of different parameters for their metallic centers, i.e. QM for the wild-type-like centers while standard parameters for the H63 variants. Given this caveat, we opted for the second ground state with a reduced Cu-center and a protonated H63 which can be easily defined by conventional force fields. Previous MD studies of SOD1 either did not disclose any details for the metal center parameters (Srinivasan and Rajasekaran, 2017; Pereira et al., 2021) or similar to this study modeled the SOD1 from the second ES complex (Mera-Adasme et al., 2013).

For every SOD1 variant including the wild-type, we have generated four different structures which differed by means of the metallation and presence of the intramolecular disulfide bridge. RMSD plots show that none of the systems underwent a large structural movement regardless of the variants tested (Figure 3A). At 400 K the backbone mobility increased regardless of the metallation or disulfide bond or mutation. A similar pattern was spotted in the R _G of the dimer (Figure 3B). In line with RMSD measurements, the compactness of the dimer was slightly increased at 400 K for every system analyzed. On the hand, the variants, including the metal-binding double mutant, did not significantly affect the R _G of dimer. Overall, we have further analyzed these simulations to compare the dynamics of the hypothesised mutants with the ALS-linked variants and to elaborate the impact of metallation and disulfide bond on the dynamics of SOD1 variants.

Fluctuations of the Cα were monitored to compare the flexibility of the putative SOD1 variants with metal-binding and/or wild-type variants (Figure 4). One notable observation is the fact that the metal binding loops of the SOD1 structure which ligate the metals were the most flexible portion of the SOD1 structure regardless of mutation, temperature and/or disulfide bond. As such, temperature increase led to apparent increase in the flexibility of these loops which were highlighted in Figure 4. Otherwise, fluctuations of the β-barrel structure and termini did not significantly higher. For the wild-type enzyme we noted that the fluctuations were reduced as the structure became metalated and/or possessed an intact disulfide bond. This outcome is in line with the previous experimental observation that the stability of SOD1 structure is contributed by the disulfide bridge (Goodsell and Olson, 2000; Furukawa and O’Halloran, 2005), and the imidazolate-bridged metallic center (Arnesano et al., 2004). Among the SOD1 variants tested we noted an increase in the flexibility of the relatively shorter electrostatic loop (Figure 4). At 298 K, for the loops in the chain A (Figure 4A), the highest flexibility was observed for the H63R variant for the holo-oxidized form while for chain B (Figure 4B) the metal-binding double variant and K136A. These observations were also spotted in the R _G analyses of the loops (Supplementary Figure S1, 2). For the K136A variant, both of the loops showed reduced compactness for both the chains in the holo-oxidized form (Supplementary Figure S1), the observation which was parallel to its enhanced fluctuations (Figure 4A). Similarly, the variant H63R led to a loosened electrostatic loop particularly for the chain B (Supplementary Figure S2). On the other hand, the other variants H63A and G37R including the wild-type enzyme did not show particularly high flexibility and altered compactness compared to these three; H63R, K136A and the double mutant. When the temperature was increased to 373 K, we observed a consistent increase in the fluctuations of the metal-binding loops for every variant including the wild-type enzyme (Figures 4A,B bottom panels). Consistently the dimer compactness and flexibility have been affected implying reduced stability at high temperature. Notably, all variants including the wild-type enzyme showed enhanced fluctuations for the shorter loop reaching as high as 15 Å. Together with the compactness and fluctuation analyses, we surmise that the putative SOD1 variants could lead to altered conformations implying their potential toxicity in ALS.

To investigate whether the SOD1 variants analyzed here restored an intact bimetallic center or not, we visually monitored the metals along with the metal-binding loops (Figure 5) and measured the average coordination distance for all systems (Supplementary Figure S3). Our results promptly suggested that the Cu-center was less tolerant to increased kinetic energy than the Zn-center (Supplementary Figure S3). This observation was aligned with the fact that Zn-center or explicitly the Zn ion has been attributed to structural stability whilst Cu ion was linked to enzymatic activity (Wright et al., 2019). At high temperature, almost all of the variants including the wild-type SOD1 showed an increase in the Cu-coordination distance with the exception of G37R (Supplementary Figure S3). This wild-type-like variant maintained intact Zn and Cu centers particularly when the disulfide bond was intact and less intact metal centers when the intramolecular disulfide bond was reduced. Taken together with these analyses, we pointed out that even the wild-type SOD1 and the wild-type-like variant G37R failed to maintain a fully intact bimetallic center, confirming that demetallation could still be the case for the wild-type or wild-type-like variant (Crow et al., 1997). From this point, our results were in line with the paradigm explaining altered metallation as one of the unifying toxicity mechanisms behind all ALS cases (Hilton et al., 2015).

We have visualized all of the variants in the holo-oxidized state which was shown to be have the least flexible for the wild-type (Figure 4). At 298 K we observed that Cu was released from the H63R and K136A variants while Zn was kept within the metallic center (Figure 5). Despite the increase in the average coordination distance (Supplementary Figure S3), the metal-binding double mutant did not particularly leave the structure. Given the increase in the average coordination distance along with the particularly increase fluctuations of the metal binding loops, we stressed that this double mutant also had an disrupted Cu-center even if Cu was not released completely. On the other hand, none of the other variants showed a particular disruption of the metal binding sites including the variant of H63A.

Overall the putative variants of H63R and K136A showed similar to or even more extreme dynamics than the metal-binding double mutant by means of fluctuations and metallation status. Therefore, we consider that these SOD1 variants, although they have not yet linked to ALS, could exhibit similar toxicity through demetallation. One possible explanation of why neither of these positions was not reported in any ALS cases could be due to their highly toxic nature. A similar case was noted for the metal-binding mutant H80R which was first generated in an animal model and its structure was shown to have a distorted Zn center due to clashes from H71 and K136 Seetharaman et al. (2010), albeit the absence of clinical evidence linking it to ALS (Wright et al., 2019). Later, a single patient who had displayed ALS symptoms at an extremely early age with a rapidly progressive disease was reported to have this mutation (Alexander et al., 2002). Hence in line with the case of H80R which has now been linked to ALS (Figure 1A), the bridging H63 and its neighbouring K136 could lead to toxic SOD1 variants leading to ALS.

Fully demetallated SOD1 variants could be targeted by the proteosomal system due to their extremely low structural stability while partially metallated variants could escape from degradation due to the fact that partial demetallation, i.e. losing one of the metals, leads to meta-stable structures (Hilton et al., 2015). Having promising therapeutic effects from both of metal-chelating and metal-delivering agents can in fact be used to support this hypothesis (Hilton et al., 2015). Accordingly, the destructive impacts of H63R and K136A on both of the metal centers resulting in the release of Cu from the structure (Figure 5) and loosing of the Zn coordination (Supplementary Figure S3) could lead to a dramatic destabilization of the structure. Thus, another plausible explanation behind lacking epidemiological data between H63R and K136A variants would lie at the other extreme. As such, these variants could be degraded by the proteosaomal system and did not include SOD1-related toxicity.