The amino acid sequence of hSOD1 was achieved from the UniProt database (ID: P00441) (Consortium, 2015). The structure for the WT-SOD1 protein was acquired through the Protein Data Bank [PDB ID: 2C9V] (Berman et al., 2000). The PDB ID: 2C9V structure was employed because it has the highest atomic resolution (1.07Å) compared to other SOD1 structures (Pereira et al., 2021).

The functional effects of E49K and R115G mutants on hSOD1 protein were predicted using the following algorithms: Panther (Mi et al., 2005), PolyPhen-1 (Ramensky et al., 2002), PolyPhen-2 (Adzhubei et al., 2013), MAPP (Stone and Sidow, 2005), Phd Snp (Capriotti et al., 2006), SIFT (Ng and Henikoff, 2003), SNAP (Bromberg and Rost, 2007), and PredictSNP for accurate predictions along with the prediction score. The majority of these tools are designed to predict whether a certain replacement is neutral or deleterious by using a variety of factors obtained from the evolutionary, physicochemical, or structural characteristics (Bendl et al., 2014). In addition, the effects of mutants on hSOD1 stability were analyzed using the different prediction web servers, viz., DUET (Pires et al., 2014a), I-Mutant2.0 (Capriotti et al., 2005), mCSM (Pires et al., 2014b), DDGun (Montanucci et al., 2019), SAAFEQ-SEQ (Li et al., 2021), DynaMut2 (Rodrigues et al., 2021), and i-Stable. These algorithms evaluated the effect of mutation on hSOD1 stability by using free energy change (ΔΔG) calculations, which could report positive and negative values for stabilizing and destabilizing mutations, respectively (Chen et al., 2013).

Using Chimera software (Pettersen et al., 2004) version 1.16 (PDB ID: 2C9V), we have observed changes in the interactions of neighboring amino acids and structural alterations in E49K and R115G mutants compared to WT-hSOD1. The software was used to scan the three-dimensional structure of SOD1 proteins and then replace the native amino acid with the candidate to demonstrate the impact.

The molecular net charge of WT-hSOD1, E49K, and R115G mutants was computed at pH 7.4 using the Prot pi Protein Tool (https://www.protpi.ch/Calculator/ProteinTool). In addition, using Chimera software (PDB ID: 2C9V), we have visualized the electrostatic potential map onto the molecular surface of WT-hSOD1 and E49K and R115G mutants.

Molecular Dynamics (MD) simulations of the WT-hSOD1 protein and its mutants (E49K and R115G) were performed using the GROMACS 2022.6 software package (Van Der Spoel et al., 2005). The program Visual Molecular Dynamics 1.9.3 was utilized to induce the mutations (E49K and R115G) on the crystallographic structure of the WT-hSOD1 dimer (PDB ID: 2C9V). The force field CHARMM36-JUL2022 was selected for the simulations. The molecules were solvated in a triclinic box with TIP3P water molecules, and to maintain overall charge neutrality, sodium (Na⁺) and chloride (Cl⁻) ions were introduced to substitute some solvent molecules. The simulations were performed in four distinct stages. In the initial stage, the energy in the system was minimized using the steepest descent algorithm to eliminate additional forces acting on the atoms, thereby relaxing the systems. The subsequent second and third stages maintained a temperature of 300 K and a pressure of 1 bar, achieved through two 1-nanosecond (ns) simulations conducted in the canonical (NVT) and isothermal-isobaric (NPT) ensembles, with position restraints applied to the heavy atoms. The final phase removed position restraints on all atoms and included a 300 ns simulation with a time step of 2 femtoseconds (fs), allowing for trajectory recording to extract various physical properties. Copper (Cu) and zinc (Zn) ions were restrained to their designated binding sites throughout the simulations. Subsequently, the MD trajectories were analyzed using the following GROMACS package distribution programs: gmx rms, gmx hbond, gmx gyrate, gmx sasa, gmx rmsf, and gmx dssp to compute root-mean-square deviation (RMSD), hydrogen bonds (H-bond), radius of gyration (Rg), solvent-accessible surface area (SASA), root-mean-square fluctuation (RMSF), and secondary structure (SS), respectively. PCA is referred to as essential dynamics (ED) when used for the analysis of molecular dynamics (MD) simulations (Pereira et al., 2020). Theoretical aspects of PCA have been thoroughly discussed in prior studies (Zaji et al., 2023). The GROMACS tools performed PCA of the WT-SOD1, E49K, and R115G mutants in each system. Additionally, the popular conformations were examined utilizing the clustering scheme of GROMACS with a cutoff of 0.2 nm.

The free energy landscape of a protein was obtained using the conformational sampling approach, which produces the near-native structural conformation. The free energy landscape (FEL) for conformational sampling of the WT-SOD1, E49K, and R115G mutants was constructed using simulated trajectory data. Consequently, FEL was derived by mixing RMSD and Rg trajectories with free energy according to the following Equation 1: Δ G   p 1 , p 2 = − k B T ⁡ ln ⁡ ρ p 1 , p   2 (1) Where in, ΔG depicts the Gibbs free energy of state, k _B represents the Boltzmann constant and T corresponds to the temperature of the simulation. The FEL’s reaction coordinates p1 and p2, are shown by joint probability distribution ρ (p1, p2) (Papaleo et al., 2009; Noorbakhsh Varnosfaderani et al., 2023).

The g_mmpbsa tool was employed to calculate the binding free energy of the WT-hSOD1 and its mutants (E49K and R115G). g_mmpbsa computed the binding free energy of the complex structure using the Molecular Mechanics Poisson–Boltzmann Surface Area (MM/PBSA) methodology (Basith et al., 2022). The following Equation 2 is used to determine the binding free energy ΔG_binding of a complex: Δ G binding = Δ G complex   –   Δ G protein 1 + Δ G protein 2 (2)

Here, Δ_Gcomplex demonstrates the total MMPBSA energy of the protein-protein complex; ΔG_protein1 and ΔG_protein2 are solution-free energies of the individual protein₁ and protein₂, respectively. The free energy of the individual existence can be expressed as follows Equation 3: Δ G = E MM + G solvation − TS (3)

E_MM shows the average molecular mechanical potential energy in the vacuum; G_solvation defines the free energy of solvation. T and S explain the temperature and entropy, respectively, and together TS exhibits the entropic contribution to the free energy in vacuum. In addition, the E_MM contains both bonded and nonbonded interactions of the molecules, including the bond angle, torsion, and electrostatic (E_elec) and van der Waals (E_vdw) interactions. Last, the free energy of solvation and G_solvation include both electrostatic and nonelectrostatic (G_polar and G_nonpolar) components (Kumari et al., 2014). To gain more insights into the residue contributions, the residues of SOD1 responsible for the interface surface were computed using the MM-PBSA per residue free energy decomposition method (mm_pbsa module) of GROMACS-CHARMM36. These calculations were based on the same snapshots as those used in the binding free energy calculations.

A plasmid extraction kit was purchased from GeneAll (South Korea). Pfu polymerase and Ni-NTA agarose resin were obtained from Agarose Bead Technologies (ABT) and Wizbiosolutions, respectively. Dpn I was purchased from Thermo Scientific. Phenylmethylsulfonyl fluoride (PMSF), Kanamycin, Dialysis tubing cellulose membrane (12,400 Da cutoff), 1-Anilinonaphthalene-8-sulfonate (ANS), and Isopropyl-D-thiogalactopyranoside (IPTG) were acquired from Sigma. Dithiothreitol (DTT), Yeast Extract, Tryptone, and Agarose were obtained from Bio Basic. Copper sulfate (CuSO₄) and zinc sulfate (ZnSO₄) were purchased from Merck. All other reagents were of analytical reagent grade and were used without further purification.

The recombinant pET28a-hSOD1 served as the template to create E49K and R115G mutations using Pfu DNA polymerase and the QuickChange site-directed mutagenesis method. Two mutagenesis primers for E49K were designed: forward primer (5′TGGATTCCATGTTCATAAGTTTGGAGATA ATACAGCAG-3′) with mutation site (highlighted) and reverse primer (5′-AGCCTGCTGTATTAT CTCCAAACTTATGAACATGG-3′) with mutation site (highlighted). Additionally, two primers were designed for the R115G mutation: the forward primer (5′-ACCATTGCATCATTGGCGGCA CACTGGTGG-3′) with mutation site (highlighted) and the reverse primer (5′-TCATGGACCACCA GTGTGCCGCCAATGATGC-3′) with mutation site (highlighted) were designed to amplify pET28a-hSOD1 as a template. The product PCR is digested with Dpn I to destroy methylated parental plasmid DNA. The recombinant expression plasmid was transformed into E. coli BL21 (DE3) competent cells for protein production by heat shock. The protein was expressed, purified using a Ni-NTA affinity column, and characterized. The purified SOD1 proteins were reconstituted with CuSO₄ and ZnSO₄ by dialysis in 3 main steps: metal removal, metal charging, and removal of unbound metals at 4°C, as previously described (Mohammadi et al., 2024). The concentration of purified WT-hSOD1, E49K, and R115G mutant proteins was determined by the Bradford method. Using a 12.5% SDS-PAGE electrophoresis, a single band of about 16 kDa was observed with high purity (data not shown).

Enzyme activity assay was performed in a volume of 300 μL composed of Tris-HCl buffer (50 mM) containing EDTA (1 mM; pH 8.2), 0.1 mM pyrogallol (dissolved in 10 mM HCl), and 20 μL of enzyme purified by the Marklund method (Marklund and Marklund, 1974) in a microplate spectrophotometer reader (Epoch, BioTek, USA). The rate of pyrogallol autoxidation depends on the availability of superoxide and is monitored spectrophotometrically (at 420 nm, 5 min, 30 s intervals) by a color change. One unit of SOD1 activity is defined as the enzyme amount inhibiting pyrogallol autoxidation by 50% per minute, with specific activity measured in U/mg protein (Tiwari et al., 2019).

Intrinsic (Trp) fluorescence in proteins is a widely used tool for monitoring protein changes, providing detailed information on dynamics and conformational changes (Vivian and Callis, 2001; Ghisaidoobe and Chung, 2014). Intrinsic fluorescence measurements were carried out using a spectrofluorometer JASCO fluorescence spectrophotometer (FP-8300) with excitation fixed at 295 nm, and emission spectra were measured between 300 and 500 nm (Exc/Em slit widths = 5/10 nm, respectively). Extrinsic fluorescence studies using ANS are applied to monitor exposed hydrophobic pockets (Hawe et al., 2008). ANS fluorescence emission spectra were measured from 400 to 700 nm, with excited at 370 nm. The purified enzyme concentration was 20 μg/mL for each assay in phosphate buffer (20 mM, pH 7.4). Additionally, intrinsic fluorescence and extrinsic fluorescence emission spectra were investigated under amyloidogenic conditions (DTT 50 mM, Tris-HCl 50 mM, and KSCN 0.16 M). The concentration of ANS was 30 μM with a protein-to-ANS ratio of 1:30.

Using potassium bromide (KBr) tablets, FTIR spectroscopy (BRUKER TENSOR 27, Germany) was used to examine the protein’s secondary structure of the WT-SOD1, E49K, and R115G mutants. The infrared spectra were recorded at a resolution of 4 cm⁻¹, spanning the range of 400–4,000 cm⁻¹. The protein samples were investigated under amyloidogenic conditions (50 mM Tris-HCl, 50 mM DTT, 0.16 M KSCN) at a concentration of 25 μM. Subsequently, 10 μL of each sample was used for FTIR evaluations. The buffer baseline was then subtracted before taking each spectrum (Chowdhury et al., 2019). The peak assignments were accordingly made by the previously described spectral components associated with various secondary structure elements (Gregoire et al., 2012). The raw spectra in the amide I region (1,600–1,700 cm⁻¹) were further processed, using the OriginPro 2021 software package. In the FTIR curve fitting method for analyzing the amide I band components, Fourier self-dissociation (FSD) and second derivative spectra were used to elucidate the overlapping peaks in the amide I region.

The fibrillation kinetics was investigated by thioflavin-T (ThT) fluorescence assay based on the rise in the fluorescence quantum efficiency of ThT after binding to the amyloid fibers. The aggregated SOD1 protein samples, under the conditions inducing amyloid aggregate formation (viz., 0.16 M KSCN, 50 mM DTT, 50 mM Tris-HCl) at 37°C (pH 7.4) and 190 rpm were then incubated. The concentration of the SOD1 variants fibrils (30 μM dimer, 20 mM phosphate buffer, pH 7.4) was also obtained. ThT (20 μM) was then added to the incubated samples at different time intervals. Afterward, fluorescence was recorded with the excitation and emission at 444 nm and 485 nm, respectively. The excitation and emission slit widths were subsequently set at 5 and 10 nm, in that order. The reported data were thus the mean measurements in three replicates. As previously reported (Esmaeili et al., 2020), the fibrillation kinetic parameters were acquired using the following Equation 4: F = F m a x   1 + exp − t − t m / τ (4) where F and F _max represent the fluorescence intensity at time t and the end time of ThT signal recording, respectively. As well, t shows the characteristic time constant and t _m stands for the time through which the half of the amyloid aggregates are formed. F _max and τ are the floating parameters, according to the Best-Fit (BF) method. Furthermore, k _app denotes the apparent rate constant, given by 1/τ for the increase in the length of fibrils, and the formula of t _D = t _m - 2τ is then practiced to calculate the delay time (t _D ).

The analysis of single nucleotide variations and their prioritization for experimental characterization heavily relies on computational methods to predict the effects of mutations on protein function (Bendl et al., 2014). The effect of mutations on the hSOD1 protein function was analyzed by eight algorithms: PANTHER, SIFT, SNAP, PhD-SNP, MAPP, PolyPhen-1, PolyPhen-2, and Predict-SNP were combined into a consensus classifier, PredictSNP. The results from all those tools predicted that the mutations are deleterious except Polyphen-2 for the E49K mutant, indicating the neutral effect of the mutation on hSOD1 function with a confidence score of 0.68. Using Predict-SNP, the confidence score predicted for E49K and R115G mutants was found to be 0.72 and 0.87, respectively, which suggested that the mutations have a deleterious effect on hSOD1 protein function. The confidence score that was calculated from each tool was tabulated in Table 1. These results demonstrate the importance of using different algorithms to ensure optimal reliability. Protein stability is crucial in thermodynamics, determined by the change in free energy between the folded and unfolded states (Srinivasan and Rajasekaran, 2018). The stability hSOD1 concerning mutations was examined through ΔΔG calculations using several predictive tools available on the web, including i-Stable, I-mutant2.0, DUET, mCSM, DDGun, DynaMut2, and SAAFEC-SEQ. The results of integrated servers computed in this line for E49K and R115G mutations as presented in Table 2 showed destabilizing effects of mutation on the hSOD1 structure. Previous reports have established the role of some mutations on the function and stability of SOD1 using various computational tools (Srinivasan and Rajasekaran, 2017; 2018; Da Silva et al., 2019; Srinivasan and Rajasekaran, 2019).

To elucidate the effects of amino acid substitutions, we utilized the UCSF Chimera software for evaluating the impacts of single-point mutations on conformational alterations of the E49K and R115G mutants compared to WT-hSOD1 (Figure 2A), thereby enhancing our understanding of their contribution to disease development. As shown in Figure 2B1, the glutamic acid residue at position 49 (E49) is vital for enabling significant interactions inside protein structures by forming hydrogen bonds and van der Waals interactions. Before the E49K mutation, E49 establishes significant hydrogen bonds with the residues Arg115 and Pro62. These hydrogen bonds promote the three-dimensional structure of the protein to be stabilized and enable important chemical interactions required for appropriate biological function. Apart from its hydrogen bonding features, E49 interacts by van der Waals with Pro62. This connection increases the stability of the protein structure by supporting the integrity of the protein’s core and enhancing the hydrophobic environment. Substituting Lys with Glu at position 49 destabilizes the surrounding area. The absence of van der Waals interaction with Pro62 may modify the spatial arrangement of this residue inside the protein, despite the hydrogen bonds between Arg115 and Pro62 being intact, the alteration in the spatial arrangement of this residue may affect the protein’s overall stability and functionality (Figure 2B2). As illustrated in Figure 2C1, the original residue R115 is crucial for the structural stability and functional interaction of proteins since its various bonds can endow the protein’s integrity. The ability of R115 to form hydrogen bonds with key residues (Glu49, Cys111, Ilu112, Ilu149) of the protein is conspicuous. Besides hydrogen bonding, R115 establishes van der Waals interactions with Phe64 and His48, along with hydrophobic interaction with Val47. Moreover, Arg115 forms a van der Waals interaction with Ilu151 in chain F. These types of interactions contribute significantly to the overall stability of the protein structure. The R115G mutation significantly alters the hydrogen bonding networks and interaction profiles by substituting R115 with G. After this mutation, the hydrogen bonds that were engaged in Glu49 and Cys111 are eliminated, resulting in the local structure of the protein being disturbed. Specifically, Glu49 located in the E49-T54 loop region is crucial for the robust binding of monomers at the dimer interface; any abnormal fluctuations in this loop region might adversely affect the stability of the SOD1 dimer (Ghosh et al., 2020). Notably, Cys111, located on the surface adjacent to the dimer interface of the SOD1 subunit, is oxidized by redox substrates in human wild-type SOD1. These oxidative changes of Cys111 may alter the structure of the SOD1 dimer, leading to monomerization and subsequent aggregation formation of the SOD1 protein (Nagano et al., 2015). Furthermore, the R115G mutation reduces the van der Waals interactions with Phe64, His48, and Ilu151 (chain F) and the hydrophobic contact with Val47, which were formerly sustained by R115. In addition to disrupting these crucial connections that enhance the protein’s dimer stability, the absence of R115 may also reveal hydrophobic pockets typically buried within the protein’s core (Figure 2C2). This exposure may result in increased aggregation tendency.

The charged and polar side chains of SOD1, similar to other soluble proteins, are located on the protein surface, protrude freely into the solvent, or participate in solvent-accessible salt bridges and hydrogen bonding (Byström et al., 2010). The active site is positively charged and accounts for around 11% of the exposed surface, whereas the remainder is negatively charged (Rakhit and Chakrabartty, 2006). A protein’s net charge is a key physical characteristic that affects its solubility and aggregation (Gitlin et al., 2006). Biomolecules often possess a net negative charge to ensure solubility inside the cellular interior. Each SOD1 monomer normally has a net charge of −6 (Sandelin et al., 2007). Point mutations in SOD1 that cause ALS often reduce the net negative of SOD1 (Prudencio et al., 2009). At pH 7.4 using the Prot pi Protein Tool, the WT-hSOD1, E49K, and R115G mutants have the molecular net charge (−6.52), (−4.53), and (−7.52), respectively. The results obtained from this tool demonstrated a significant alteration in the molecular net charge of both mutants. Notably, for mutant E49K, there was a reduction in net negative charge compared to the wild-type protein. As illustrated in Figure 3, the electrostatic potential was mapped onto the surfaces of WT-hSOD1. Glu49 is located on the outer surface of metal-binding loop IV, surrounded by positively charged His48 and Arg115, which likely serves to neutralize these positive charges and maintain the surface electrostatic potential of SOD1. Mutation of Glu49 to a positively charged residue (Lys) was found to alter the net negative charge (reduce net negative charge) on the protein’s surface (Figure 3), which could reduce the stability of the SOD1 dimer, promote aberrant interactions and misfolding, and/or nucleate aggregation. Arg115 is located on the outer surface of Greek key loop VI, surrounded by negative charge residue (Glu49). On the other hand, an ALS-associated mutation of Arg115 to a neutral residue (Gly) was found to increase the net repulsive charge of the SOD1 (Figure 3). A previous study suggested that two ALS-associated point mutations at the surface residue Glu100 (E100G and E100K) in SOD1, which decrease the net negative charge of the protein by eliminating negatively charged or introducing positively charged amino acids, respectively, could reduce the stability of SOD1 and promote aggregation (Tompa and Kadhirvel, 2020).

MD simulations capture a diverse range of crucial biomolecular processes that include structural changes and protein folding, alongside forecasting the impacts of mutations on the structure of the protein (Idrees et al., 2017; Hollingsworth and Dror, 2018; Idrees et al., 2022). Therefore, the simulations were conducted on WT-hSOD1, E49K, and R115G mutants (holo dimer forms) to assess the structural dynamic changes. The root-mean-square deviation (RMSD) was applied to quantify the difference between the initial structural conformation and the final positions of the Cα atoms in the protein backbone. The RMSD calculation was carried out to evaluate the changes in the structure and dynamics of the modeled protein (Aier et al., 2016). The backbone RMSD of WT-hSOD1, E49K, and R115G mutants of the proteins were calculated over 300 ns simulations (Figure 4A) and pointed to an unstable behavior at the beginning of simulations. After approximately 150 ns, E49K, and R115G mutants, the structures were floating in a stable state. In contrast, the RMSD value for the WT-SOD1 demonstrated fluctuation for the remainder of the MD trajectory. The average root-mean-square deviation (RMSD) values for WT-SOD1, E49K, and R115G mutants were (0.27 nm), (0.24 nm), and (0.22 nm), respectively. These changes in mutants demonstrated the loss of protein conformational deviation as compared to the WT-SOD1. Our results concerning the RMSD of SOD1 mutants align with previous research on different mutants, which exhibited a relative conformational deviation compared to the WT-SOD1 (Srinivasan and Rajasekaran, 2017; 2019). Besides, hydrogen bonds, one of the most common noncovalent interactions in proteins, play an extremely important role in protein structure and function (Adhav and Saikrishnan, 2023). They were also the first bonds to react to the structural alterations caused by various protein mutations (Alemasov et al., 2017). To determine the potential impact of mutations on the conformational changes in SOD1, the number of hydrogen bonds was computed. The mean number of H-bonds for the holo dimeric WT-SOD1, E49K, and R115G mutations were 204, 202, and 199, respectively (Figure 4B). The number of hydrogen bonds in the mutants is significantly decreased compared to the WT-SOD1, which may influence their structural stability and function, potentially leading to misfolding and aggregation. Previous studies indicate that the removal of hydrogen bonds considerably decreases the stability of SOD1 and promotes aggregation (Tompa et al., 2020; Noorbakhsh Varnosfaderani et al., 2023). Correspondingly, protein stability and the number of intramolecular interactions are both affected by the relative compactness of the protein structure, which is shown by the protein’s radius of gyration (Rg) (Wise-Scira et al., 2013). According to the simulation data presented in Figure 4C, the Rg value computed for the E49K mutant indicated a phase of structural instability, which was observed in the first half of the simulation. After approximately 150 ns, the Rg value exhibited stable behavior for the remainder of the MD trajectory. In comparison, the Rg values for both the WT-SOD1 and mutant R115G displayed initial structural instability, which persisted throughout the entirety of their respective simulations. This analysis also indicated that the average Rg values for the WT-SOD1 (2.07 nm) are relatively similar to those of mutants E49K (2.02 nm) and R115G (2.04 nm). This result suggests that there are slightly increased compactness alterations in the analyzed mutants compared to WT-SOD1. Our findings are in agreement with previous studies (Kalia et al., 2023). In this context, the solvent-accessible surface area (SASA) refers to the exposed regions of protein structures that are accessible to solvent molecules, a key factor in studying protein folding and stability. Amino acid residues in a protein may be categorized as exposed or buried based on their solvent-accessible surface area (SASA) values (Ausaf Ali et al., 2014). The SASA values acquired from all simulations (Figure 4D) showed an unstable behavior at the beginning of all simulations. Following approximately 150 ns, the SASA values exhibited stable behavior until the end of the MD trajectories. This analysis indicated that the average SASA values for the WT-SOD1 (152.78 nm²) were similar to those of the mutants E49K (152.27 nm²) and R115G (153.21 nm²). This finding implies that no significant changes in SASA occurred among the examined mutants compared to WT-SOD1. The findings of our research are in agreement with the previous study, which indicated that no SASA alterations in the analyzed mutants when compared to the WT-SOD1 (Pereira et al., 2021). Next, root-mean-square fluctuation (RMSF) was used to assess the residual flexibility of the WT-SOD1 E49K and R115G mutants (holo-dimer forms). Of note, RMSF, a critical parameter, provides information about the structural flexibility of the Cα atoms in each residue of the corresponding system through MD simulation (Haq et al., 2017). As seen in Figures 5A, B, both monomers displayed asymmetric fluctuation patterns. Specifically, the WT-SOD1 has notable fluctuations in chain A at residues 62–77, which are located in the zinc-binding region of loop IV. Notably, in the E49K and R115G mutants, all atoms contribute almost equally to their relatively stable fluctuations. The average residual RMSF for the WT-SOD1 and mutants in the A and F chains are presented in Table 3. These results suggested that WT-SOD1 is more flexible than the mutants in the holo-dimer form. We find our data to correspond with the results of the previous study, which indicated several mutants have less fluctuation relative to WT-SOD1 (Srinivasan and Rajasekaran, 2019). Secondary structure is a vital component in understanding the conformational stability and functionality of proteins (Kumar et al., 2021). The composition of secondary structures is pivotal in the aggregation process of proteins linked to neurodegenerative diseases (Zaji et al., 2023). Moreover, the transition of the alpha helix into beta-sheets is regarded as a primary issue in the formation of amyloid fibrils, which directly correlates with an increase in aggregation propensity (Dong et al., 2016; Srinivasan and Rajasekaran, 2019). Given this context, the analysis of the secondary structure content for the WT-SOD1, E49K, and R115G mutants was performed using the Dictionary of Secondary Structure in Proteins (DSSP) algorithm to monitor the secondary structural changes throughout the simulation (Table 4). The analysis of the results revealed that the alpha helix arrangement in E49K and R115G mutants was lower than that of WT-SOD1. Conversely, we found that the propensity of beta-sheets was slightly higher in E49K and R115G mutants compared to WT-SOD1. Consequently, we could suggest that E49K and R115G mutants compared to the WT-SOD1 have unequal proportions of different secondary structures, which might be detrimental not only in altering protein stability and function but also in inducing misfolding in SOD1, ultimately leading to protein aggregation. We found that our findings are in harmony with those identified in prior studies (Zaji et al., 2023; Ashkaran et al., 2024). A common multivariate statistical method for minimizing the number of dimensions describing the dominant protein domain motion is Principal Component Analysis (PCA) (Kalia et al., 2023). The projections for the MD trajectories of WT-SOD1, E49K, and R115G mutants into the subspace spanned by PC1 and PC2. As shown in Figures 6A, B, the E49K and R115G mutants occupied a smaller area in conformational space. We also observed changes in cluster shape in the conformational space of both mutants compared to the WT-SOD1. The analysis thus pointed to alternations in the overall essential dynamics of both mutants. Our findings are in agreement with previous studies (Pereira et al., 2021). Cluster analysis of molecular simulation trajectories for molecular systems where the conformation and orientation of the sampled states are important parameters requires a different approach from the standard clustering methods. These systems require cluster analysis methods that discriminate based on molecular orientation and conformation. New clustering methods that account for these parameters are presented and demonstrated by the analysis of trajectories produced in protein simulations for validation (Abramyan et al., 2016). To better visualize the effect of mutation on SOD1 structure, we analyzed the popular conformations obtained using the clustering method with a 0.2 nm cutoff. To further compare the WT-SOD1 with the variant systems, we performed a hierarchical clustering analysis, which allowed us to compare each system with all the others. Different clustering was observed in the wild-type and mutant proteins, with cluster classifications ranging from 5 to 10. The number of clusters for the WT-SOD1, E49K, and R115G mutants was 10, 5, and 6, respectively. The fewer number of clusters may imply a faster aggregation process. On the other hand, the highest conformations for the WT-SOD1, E49K, and R115G mutants was 1916, 2,889, and 2,688, respectively. The most populated conformations are shown as cluster 1 for the WT-SOD1 and mutants in Figure 7.

Protein aggregation can be well evaluated through the multiple conformational states occupied in the free energy landscape (FEL) (Zheng et al., 2013). Misfolding of protopathic proteins is a key biomolecular trigger for amyloid formation, ultimately leading to amyloid fibrils (Adamcik and Mezzenga, 2018). Accordingly, we separately evaluated the FEL for the WT-SOD1 and mutants. Similarly, a FEL was constructed to enumerate the occurrence of global energy minima constructed to investigate the pathology of the E49K and R115G mutants during the simulation period (Yang et al., 2013). RMSD and Rg were the two coordinates used for FEL analysis to show the compaction and structural stability of the protein during the molecular dynamics period. The Gibbs free energy of the generated FEL of WT-SOD1, E49K, and R115G mutants ranged from 0 to 10 kcal/mol, as shown in Figure 8. The free energy ranges from purple to yellow, challenging the purple region for the minimum global energy composition and the yellow for metastable states. We concluded that the substitution mutation in SOD1 significantly altered the folding pattern, possibly leading to increased RMSD values toward the multiple energy minimum achieved by the mutant conformers. The minimum global energy number for a native protein is usually limited to one. The free energy landscape for the WT-SOD1 structure shows a unique favorable region (Figure 8A) that lies between Rg and RMSD values of ∼2.08 nm and ∼0.38 nm, respectively. The results for the E49K mutant showed a global energy minimum that could occur in more compact conformational states, each of which could cause misfolding and lead to aggregation. The formation of favorable free energy domains (Figure 8B) was 2.02 and ∼0.23 nm for Rg and RMSD, respectively. However, results for the R115G mutant protein indicated multiple global energy minima that could occur in various thermodynamically conceivable conformational states, any of which can lead to aggregation. The effect of substitution mutation on this process led to the formation of numerous favorable free energy basins (Figure 8C), which were ∼2.03 and ∼0.25 nm for Rg and RMSD, respectively. Hence, as the free energy profile showed, the presence of R115 in the wild-type plays a prominent role in the presence of the limited free energy basin along RMSD and Rg, which by replacing Gly led to the change of SOD1 proteins and reaching several global free energy minima. Furthermore, the free energy domain obtained by WT-SOD1 was smaller than the conformational structures of the R115G mutant, indicating a more favorable conformation adopted by the wild-type conformers than the mutant. Altogether, the results support the idea that an increased percentage of multiple conformers with lower free energy in the mutant than the wild-type are directed toward forming predominantly unfolded states. Furthermore, in confirmation of our results indirectly compared to previous studies, it suggests that the aggregated proteins acquire the minimum multiple energy for the conformational structures that indicate the formation of toxic aggregates in mutant SOD1 (Maitham Qassim et al., 2023; Noorbakhsh Varnosfaderani et al., 2023; Abdullah Waheed et al., 2024).

The MM-PBSA strategy for binding free energy estimation has become one of the most adopted techniques to measure interaction energies in biomolecular studies (Wang et al., 2019). To quantify SOD1 dimerization in terms of energy, an analysis of MM-PBSA was performed. The results identified the WT-SOD1 system as the most stable form with the highest binding energy, −105.17 kcal/mol. The results further suggested that E49K and R115G mutants compromised the dimerization potential between monomers, leading to reduced binding energies of −32.1 and −79.47 kcal/mol, respectively. The binding energies in the mutants are significantly decreased compared to the WT-SOD1, which may influence their dimer structural stability, potentially leading to aggregation. Our results concerning the MM-PBSA of SOD1 mutants align with previous research on several mutants, which exhibited mutations that compromised the dimerization potential between monomers, leading to reduced binding energies compared to WT-SOD1 (Basith et al., 2022). The technique of per-residue binding free energy decomposition can reveal the contributions of the key residues responsible for the protein-protein interactions at the interface. A total of 1,000 frames extracted from the last 50 ns of the MD trajectories were decomposed by the MM-PBSA method. The analysis of free energy decomposition was performed on the WT-SOD1 and mutants. Free energy decomposition analysis helps to find the contribution of a single residue by summing its interactions over the entire residues. The results of this study showed that in the wild-type state (Figure 9A), the binding free energy for residues E49 and R115 in chain-A was 0.64 and −4.27 kJ/mol, respectively. For chain-F, the binding energy values for residues E49 and R115 were 0.06 and 0.26 kJ/mol, respectively. The E49K mutation (Figure 9B) changed the binding free energy value, which was 9.31 and 7.25 kJ/mol in the A and F chains, respectively. On the other hand, the binding free energy values for the R115G mutation (Figure 9C) were −3.26 and −3.43 kJ/mol in chains A and F, respectively. The analysis thus pointed to alternations in the per-residue free energy decomposition of both mutants compared to the WT-SOD1, which may influence their dimer structural stability, potentially leading to aggregation.

Almost all mutations in SOD1-amyotrophic lateral sclerosis result in a significant decrease in SOD1 enzyme activity, according to published data from individuals with this disease (Saccon et al., 2013). The active site of the SOD1 protein contains a Cu²⁺ ion coordinated with four histidine residues: His 46, His 48, His 120, and His 63 (Keerthana and Kolandaivel, 2015). The enzymatic activity of WT-hSOD1, E49K, and R115G mutants was 7,031, 5,536, and 4305 U/mg, respectively. The results of our study revealed reduced enzymatic activity for E49K and R115G compared to holo WT-hSOD1. At position 49, glutamic acid has a negative charge and is vital for setting up the electrostatic environment that His48 needs to coordinate the copper ion. This contact is essential for the enzyme’s performance. The E49K mutation substitutes essential glutamic acid with positively charged lysine. This change may substantially disturb the delicate balance necessary for the correct coordination of the his48 with the copper ion. On the other hand, the R115G mutation diminishes the van der Waals interaction with His48. The lack of this interaction may further affect the structural integrity and activity of the enzyme, especially considering the significance of copper coordination at the His48 location. Initially, researchers believed that ALS-causing SOD1 mutations altered the protein’s structure, resulting in an unstable or misfolded enzyme that lowered dismutase activity and heightened oxidative stress. However, different FALS-linked SOD1 mutations had different amounts of dismutase activity. This suggests that mutant SOD1’s toxicity could not be explained by a simple loss of function. Additionally, there are no FALS patients known to have mutations resulting in complete loss of SOD1 (Barber et al., 2006). Our results correspond with previous investigations indicating that the specific activity of several SOD1 mutants is decreased in comparison to the WT-SOD1 (Mavadat et al., 2023b; Ashkaran et al., 2024).

Among the aromatic amino acids (Trp, Tyr, and Phe), intrinsic Trp fluorescence is the most often used method to track local structural and dynamic changes in proteins (Raghuraman et al., 2019). WT-hSOD1 only possesses one Trp at residue 32 (W32), which is highly solvent-exposed (Grad et al., 2011). Figure 10A, A shows the intrinsic fluorescence of the E49K mutant is lower than that of WT-hSOD1, indicating that the environment of the single solvent-exposed Trp residue 32 is not buried. In other words, the polarity of the Trp residue microenvironment increased and was more exposed to the solvent. The R115G mutant fluorescence spectrum shows an increase in fluorescence intensity with a change of λmax in the emission spectra and a blue shift compared to the wild-type, showing that the alterations in the circumambient environment, the tryptophan (W32), toward a more hydrophobic environment. The ability to measure the hydrophobicity of a protein may be advantageous in comprehending and forecasting the implications of alterations to the sequence of its structure. The increased ANS fluorescence signified that ANS bound to the exposed contiguous hydrophobic regions (Cardamone and Puri, 1992; Idrees et al., 2016). The hydrophobic inner core of the monomer SOD1 contains roughly half of the protein backbone (Wright et al., 2019), and protein dimerization provides approximately 640 Å of hydrophobic surface area at the dimer interface through a complex network of hydrogen bonds and hydrophobic interactions, resulting in a compact mature homodimer with minimal solvent accessibility. Accordingly, holo-SOD1 is one of the most stable proteins known (Trist et al., 2021). To evaluate the impact of mutations on changes in SOD1’s hydrophobicity, ANS fluorescence was utilized. As shown in Figure 10B, the ANS fluorescence of E49K and R115G indicates a significant increase compared to the WT-SOD1, indicating the role of mutations in conformational changes and increased amounts of exposed hydrophobic surfaces. Our findings are in agreement with previous studies (Namadyan et al., 2023). FTIR, which is known as one of the accurate techniques for reflecting the secondary structure of proteins, was used to study misfolding and aggregate formation (Miller et al., 2013). Accordingly, the formation of β-sheet structures in WT-SOD1 and the mutants (under destabilizing conditions) was confirmed by FTIR spectroscopy (Figures 10C–E). The FTIR spectra were also employed to distinguish between types of β-sheet structures and the amyloidogenic conformers (amide I bands) were then determined. According to Table 5, the peak in 1,610–1,630 cm⁻¹, indicating the formation of fibrils and β-oligomers with the intermolecular β-sheet structure due to protein aggregation, was reported for WT-SOD1 and mutants. The native β-sheets could produce an absorption peak in the 1,623–1,641 cm⁻¹ range. Examining the amide (I) bands, a parallel arrangement could be thus marked in the anti-parallel arrangement of the β-strands in the protein aggregates. In this line, the parallel β-sheets could display only an elevated component at 1,623–1,641 cm⁻¹ (Hoffner et al., 2014). The disordered conformation (random coil) is usually associated with the IR band between 1,640 and 1,648 cm⁻¹. A distinctive feature of a random coil is that it is non-repetitive. Such an assignment is supported by the position of a prominent band in the spectrum of apparently disordered proteins (Byler and Susi, 1986). A component centered at approximately 1,650–1,660 cm⁻¹ is assigned to the α-helix, which is consistent with theoretical calculations and the observation of bands in the spectrum of α-helix proteins (Sadat and Joye, 2020). α-Helix structure was observed for WT-SOD1 at peak position 1,655, cm⁻¹. On the other hand, α-helical structure (3₁₀-helix) was observed for the E49K and R115G mutants at peak positions of 1,666 and 1,664 cm⁻¹, respectively. Bands near 1,663 ± 3 cm⁻¹ are assigned to 3₁₀ helices, although this structure is rarely found in proteins (Dong et al., 1990). The 3₁₀-helix structure is classified as a subset of the α-helix and appears in proteins when the α-helix is distorted due to the presence of (mutation) an undesirable residue (Peters, 2005). The assignment of bands around 1,670–1,699 cm⁻¹ to β-turns has been suggested (Krimm and Bandekar, 1986). Hence, the β-turn structure was observed for the E49K and R115G mutants in the peak region of 1,682 and 1,679 cm⁻¹, respectively. Therefore, the FTIR spectroscopy results showed that the mutants have different changes in the natural β-sheet structure of the protein compared to the wild-type. Therefore, the changes in the secondary structure of the SOD1 variants indicate a structural rearrangement of SOD1, potentially leading to aggregation. The results are consistent with what has been previously reported (Ashkaran et al., 2024; Farrokhzad et al., 2024).

To obtain a deeper understanding about the relationship between aggregate formation and changes in the Trp environment, we looked at the intrinsic fluorescence of pathological aggregates formed from wild-type and mutants under destabilizing circumstances. Obtained data indicated that fluorescence intensity in both mutants during aggregation compared to wild-type shows that Trp-32 is not buried. In other words, the polarity of the Trp residue microenvironment increased and was more exposed to the solvent, as seen in Figure 11A. To evaluate the impact of mutations on changes in SOD1’s hydrophobicity, as well as to observe the exposure of hydrophobic pockets in β-sheet structures during aggregation formation, ANS fluorescence was utilized. Several factors, such as formulation conditions, protein sequences, and structure, influence protein aggregation. The classical hydrophobic effect drives the interaction of exposed hydrophobic pockets on the protein surface, leading to the formation of protein aggregates (Tosstorff et al., 2019). According to Figure 11B, ANS fluorescence emission of mutants under amyloid induction conditions showed lower emission intensity than WT-SOD1. The decrease in ANS fluorescence intensity is probably due to the coverage of large parts of the exposed hydrophobic surfaces, which results from the interaction of hydrophobic segments between protein molecules during the amyloid process. These results are aligned with observations noted in earlier studies (Abbasabadi et al., 2013; Salehi et al., 2015; Baziyar et al., 2022; Mavadat et al., 2023b).

The end-point level of ThT fluorescence intensity has been on mature fibrils, as generally reported. A growing number of structural models for amyloid fibrils thus strongly suggest that ThT binds to the long axis of fibril by aligning parallel, intercalating to the repeating side-chain interactions running across β-strands within the β-sheet layer (Wu et al., 2009). In this study, WT-SOD1 and mutants were incubated under inducing conditions for time intervals of 0–72 h with and without agitation. Once NaCl was substituted for KSCN, no amyloid was formed. This implied that the effect of KSCN was related to its chaotropic properties but not its impact on ionic strength. The representative fluorescence time courses for the WT-SOD1 and mutants are illustrated in Figure 12 and a complete set of the kinetic parameters are then given in Table 6. The results showed that the WT-SOD1 and mutants formed amyloid aggregates when there were high concentrations of DTT. Furthermore, experiments showed that SOD1 fibrillation preceded the growth in ThT fluorescence, which was different in mutant and WT-SOD1. This delay was different in mutants and WT-SOD1. Comparison of the lag phases showed that the mutants tended to have shorter lag phases and faster fibrillation than holoWT-SOD1. This step increased thioflavin T fluorescence to different levels for WT-SOD1 and mutants. However, this change in lag was shown to be compatible with other reported results that have been studied under different settings (Abbasabadi et al., 2013; Salehi et al., 2015). In particular, FALS-causing mutations wouldn't always reduce the latency and/or increase the aggregation propensity of SOD1. In addition, the kinetics of amyloid fibril aggregation is delayed, with an initial increase in ThT fluorescence intensity, which was more pronounced for mutants than for WT-SOD1, followed by a decrease in ThT intensity. This may be due to the formation of mat-like aggregates of fibrils, resulting in a number of modified or unavailable binding sites for ThT. This difference in the kinetics of aggregate formation in the mutant forms may be due to the amino acid substitution. Finally, when the ThT fluorescence intensity reaches its maximum, the plateau phase begins. However, this stage does not last long because amyloid fibers tend to form web-like network structures. The findings of our study may support the idea that SOD1 aggregation occurs through a pathway including weakening of the dimer interface, dimer dissociation, reducing the binding/affinity of metal ions from monomers, and aggregation of monomers (Tiwari and Hayward, 2003; Hough et al., 2004). Furthermore, the results show that SOD1 forms amyloid aggregates under different conditions including agitation, ionic strength, temperature and physiological pH. We were able to demonstrate the potential of WT-SOD1 and mutants to form amyloid structures under conditions inducing amyloid aggregation. These results are consistent with previous studies showing the importance of fibrillation or oligomerization of SOD1 in the pathogenesis of ALS (Oztug Durer et al., 2009; Abbasabadi et al., 2013; Salehi et al., 2015; Baziyar et al., 2022).