First Principles Calculation of Protein–Protein Dimer Affinities of ALS-Associated SOD1 Mutants

Cu,Zn superoxide dismutase (SOD1) is a 32 kDa homodimer that converts toxic oxygen radicals in neurons to less harmful species. The dimerization of SOD1 is essential to the stability of the protein. Monomerization increases the likelihood of SOD1 misfolding into conformations associated with aggregation, cellular toxicity, and neuronal death in familial amyotrophic lateral sclerosis (fALS). The ubiquity of disease-associated mutations throughout the primary sequence of SOD1 suggests an important role of physicochemical processes, including monomerization of SOD1, in the pathology of the disease. Herein, we use a first-principles statistical mechanics method to systematically calculate the free energy of dimer binding for SOD1 using molecular dynamics, which involves sequentially computing conformational, orientational, and separation distance contributions to the binding free energy. We consider the effects of two ALS-associated mutations in SOD1 protein on dimer stability, A4V and D101N, as well as the role of metal binding and disulfide bond formation. We find that the penalty for dimer formation arising from the conformational entropy of disordered loops in SOD1 is significantly larger than that for other protein–protein interactions previously considered. In the case of the disulfide-reduced protein, this leads to a bound complex whose formation is energetically disfavored. Somewhat surprisingly, the loop free energy penalty upon dimerization is still significant for the holoprotein, despite the increased structural order induced by the bound metal cations. This resulted in a surprisingly modest increase in dimer binding free energy of only about 1.5 kcal/mol upon metalation of the protein, suggesting that the most significant stabilizing effects of metalation are on folding stability rather than dimer binding stability. The mutant A4V has an unstable dimer due to weakened monomer-monomer interactions, which are manifested in the calculation by a separation free energy surface with a lower barrier. The mutant D101N has a stable dimer partially due to an unusually rigid β-barrel in the free monomer. D101N also exhibits anticooperativity in loop folding upon dimerization. These computational calculations are, to our knowledge, the most quantitatively accurate calculations of dimer binding stability in SOD1 to date.

S1 Error analysis for each energy contribution S1.1 S * and ∆G free o Since S * and ∆G free o are calculated analytically, there are no error associate with the two terms. S1.2 ∆G restr dist+a The integral of the PMF in Eq. (7) (main text) is dominated by the immediate vicinity of the distance in the bound state. Since the PMF error at r * is typically much larger than the PMF error at the bound state ( Fig. 3 (main text) right most column), the error of I * could be safely approximated by the error propagation of W (r * ) alone (Eq. (S1)).
where the error σ W (r * ) is taken from the larger value of 1) statistical error of MBAR estimator calculated from the last 50% REMD-US trajectories, or 2) the difference of W (r * ) calculated from the last 75% and the last 50% REMD-US trajectories. The error of ∆G restr dist+a is then given by the simple expression:

S1.3 Other contributions to Eq. (11) (main text)
The error of the other free energy contributions in Eq. (11) (main text) are approximated by the larger of two values, as obtained from either of the following two methods: 1) The standard deviation of free energy perturbation calculation from 10,000 re-samplings of the PMF. The re-sampled PMF at certain reaction coordinate value, ξ i , is sampled from a Gaussian probability function where F 0 (ξ i ) and σ F 0 (ξ i ) are the MBAR estimated PMF value and error at ξ i .
2) The difference of the free energy contribution (e.g. ∆G bound LA,c ) calculated by the last 75% and the last 50% REMD-US trajectories.

S2 Data sets
All reaction coordinate trajectories in biasing umbrella potentials (files named Umbrella.txt) and umbrella bias potential files (plumed.dat), are available in a compressed .tar.gz file at https://phas.ubc.ca/~steve/Dimer/ These files are sufficient to construct the PMFs in the present work. The five folders named EESS, EESH, EESS A4V, EESS D101N, and CuZnSS contain the data for five respective SOD1 variants (in the same order as in Table 5 (main text)). For each variant, there are 16 subfolders: • 6 subfolders are for conformational contributions in the dimer (C LA, C LB, C BA, C BB, C IA, and C IB), • 5 subfolders are for orientational/angular contribution in the dimer (O THETA, O PHI, O PSI, A theta, and A phi), • 3 subfolders are for conformational contributions in the monomer (C L, C B, and C I), • 1 subfolder is for the separation PMF (I), • and 1 subfolder is for the equilibration (equil). In each folder that represents free energy contribution (folders other than equil), there is a folder serial and a folder RE, which store the data for initial configuration preparation and REMD-US. The reaction coordinate trajectory of each umbrella in REMD-US simulation is stored in a folder named starting with "Z" followed by the bias center (units are in either nm or rad) of the reaction coordinate (e.g. Z0.52 in C L corresponds to an umbrella with bias center at RMSD= 0.52 nm). Figure S1: The PMFs for different fractions of the MD trajectories for WT E,E(SS), to illustrate convergence. PMFs are constructed with either 100%, the last 75%, or the last 50% of the REMD-US trajectories, to show convergence. The PMFs shown in each panel are ordered from the top-left to the lower-right in the same order that the restraints are sequentially applied or released, following the arrow in Fig. 2 (main text).       The initial configurations are steered from the RMSD∼ 4Å equilibrated structure. Right panel: The initial configurations are steered from the RMSD∼ 8Å equilibrated structure. The schematic arrows below the x-axis in each panel indicate the process of RMSD steering from equilibrated structures (blue dots) where the initial configurations are sampled. The two PMFs are significantly different and do not converge on the time scales of our simulations. For this reason, we did not use either of these seeding methods, but rather mixed initial conditions from both small and large values of RMSD.