Determination of the absolute binding free energy of dimer is carried out using a generalization of the method of Roux and colleagues (Woo and Roux, 2005; Gumbart et al., 2013a) in which a series of restraints are successively applied and released to divide the binding free energy calculation into separate calculations, each with manageable convergence.

The restraints in this study include conformational restraints, orientational restraints, and angular restraints. The restraining potentials are listed in Table 1, along with their values for the protein mutants and variants considered in this study. The conformational restraints restrain the backbone atoms (N, C _α , and C for each residue) of the entire protein and the sidechain atoms of the dimer interface residues (described below). SOD1 contains two long loops, from residues 48–83 (loop 4) and residues 121–143 (loop 7), which contain minimal secondary structure and undergo large conformational rearrangements in the metal-depleted dimer (Elam et al., 2003a; Elam et al., 2003b; Strange et al., 2003; Roberts et al., 2007; Cao et al., 2008; Ahmad et al., 2009; Kevin et al., 2015), apo monomer (Banci et al., 2003; Tom et al., 2009; Das and Plotkin, 2013a), and disulfide-reduced monomer (Das and Plotkin, 2013a; Kumar et al., 2018a) (turquoise in Figures 1A,B).

The potentials inducing conformational restraints on backbone atoms of the central barrel in chains A and B of the homodimer (u _BA,c , u _BB,c ) and the corresponding loops (u _LA,c , u _LB,c ) are applied separately. The free energies corresponding to applying these restraints to either the bound dimer state (e.g., Δ G L A , c bound ) or free state (e.g., G L A , c free ) are given in Table 1, in the order they are applied (bound state) or released (free state) (see Figure 2). In all cases in Table 1, except for the separation PMF, the free energy is given in terms of applying the restraint (a positive contribution) so that the terms may be simply added to find the binding free energy.

The central barrel consists of residues 1–48, 84–120, and 143–153 (dark blue in Figures 1A,B). The dimer interface sidechain restraint for chains A and B (u _IA,c , u _IB,c ) is imposed on side chains of residues 5, 7, 50–54, 114, 148, 150–153, which are at least partially buried in the dimer interface in our simulations and are consistent with interface residues reported in previous experimental and theoretical studies (Hough et al., 2004; Das and Plotkin, 2013b) (yellow licorice in Figures 1A,B). The orientational restraint potential (u _o ) is imposed on Θ, Φ, Ψ angles defined in Figure 1B. The restraint potential u _o (Θ, Φ, Ψ) ensures that the same faces of chain A and chain B point towards each other when calculating other contributions to the dimer binding free energy, such as the potential of mean force (PMF) as a function of separation distance. Similarly, the angular restraint (u _a ) is imposed on the polar and azimuthal angles θ and ϕ defined in Figure 1B. The potential u _a (θ, ϕ) obviates the need to sample the full 4π solid angle, whose phase space sampling contribution can be accounted for analytically.

The orientational angles shown in Figure 1B are defined as follows: three groups of atoms are chosen in chain A. P ₁ is at the centre of mass of the central beta-barrel structure, represented by residues 1–48, 84–120, and 143–153. P ₂ is at the centre of mass of residues 5–7, 17–19, and 32–34, creating a reference point on the surface of the beta sheet of the central barrel. P ₃ is at the centre of mass of residues 11–13, 40–42, and 120–122, describing the “lid” of the central beta barrel. Likewise, the same groups of atoms are chosen in chain B to define P 1 ′ , P 2 ′ , and P 3 ′ . The spherical coordinate system establishing the position of chain B relative to chain A is by the distance r ( P 1 P 1 ′ ̄ ) , angle θ ( ∠ P 1 ′ P 1 P 2 ) , and the dihedral angle ϕ ( P 1 ′ -P ₁-P ₂-P ₃). The Euler angles needed to define the orientation of chain B relative to chain A as the angle Θ ( ∠ P 1 P 1 ′ P 2 ′ ) , the dihedral angle Φ (P ₁- P 1 ′ - P 2 ′ - P 3 ′ ), and the dihedral angle Ψ (P ₂-P ₁- P 1 ′ - P 2 ′ ). The same definitions for the reference frame apply to all variants.

The absolute binding free energy can be defined in terms of equilibrium binding constant as Δ G bind ≡ − k B T ⁡ ln K eq c ° by assuming a standard state concentration c° of 1 mol/L (1 molecule/1661 ų).

The equilibrium constant may be written as a ratio of two integrals, one in the bound state and one in the free state: K e q = ∫ bound d B ∫ d A e − β U ∫ free d B δ ( r A B − r A B ∗ ) ∫ d A e − β U , (1) where A and B correspond to the degrees of freedom of each protein, along with the solvent degrees of freedom that equilibrate about each protein configuration. It is convenient in practice and theoretically justified (Boresch et al., 2003) to use relative coordinates and hold one protein (A) fixed while separating the other protein (B) from it when calculating the potential of mean force (PMF) and corresponding restraints, as described below.

The essence of the calculation is that the ratio in Eq. 1 may be split by several intermediate integrals involving restraining potentials that effectively multiply the expression by unity but make the thermodynamic averaging tractable (Hermans and Shankar, 1986; Boresch et al., 2003; Woo and Roux, 2005; Gumbart et al., 2013a). The restraining potentials bias the relative orientation, relative position in spherical coordinates, and conformation of each protein to be similar to that in the bound state, as described above. For example, K _eq in Eq. 1 may be written as K e q = ∫ bound d B ∫ d A e − β U ∫ bound d B ∫ d A e − β ( U + u L A , c ) × ∫ bound d B ∫ d A e − β ( U + u L A , c ) ∫ free d B δ ( r A B − r A B ∗ ) ∫ d A e − β U , (2) where the first term in curly brackets is a configurational integral equal to e − β u L A , c ( b o u n d , U ) − 1 , which is equal to a free energy difference e β Δ G L A , c bound that can be calculated using free energy perturbation techniques. With this approach, the equilibrium binding constant K _eq in Eq. 1 can be written as the product of the following free energetic terms: e β Δ G L A , c bound = ∫ bound d B ∫ d A e − β U ∫ bound d B ∫ d A e − β ( U + u L A , c ) (2a) e β Δ G L B , c bound = ∫ bound d B ∫ d A e − β ( U + u L A , c ) ∫ bound d B ∫ d A e − β ( U + u L A , c + u L B , c ) (2b) e β Δ G B A , c bound = ∫ bound d B ∫ d A e − β ( U + u L A , c + u L B , c ) ∫ bound d B ∫ d A e − β ( U + u L A , c + u L B , c + u B A , c ) (2c) e β Δ G B B , c bound = ∫ bound d B ∫ d A e − β ( U + u L A , c + u L B , c + u B A , c ) ∫ bound d B ∫ d A e − β ( U + u L A , c + u L B , c + u B A , c + u B B , c ) (2d) e β Δ G I A , c bound = ∫ bound d B ∫ d A e − β ( U + u L A , c + u L B , c + u B A , c + u B B , c ) ∫ bound d B ∫ d A e − β ( U + u L A , c + u L B , c + u B A , c + u B B , c + u I A , c ) (2e) e β Δ G I B , c bound = ∫ bound d B ∫ d A e − β ( U + u L A , c + u L B , c + u B A , c + u B B , c + u I A , c ) ∫ bound d B ∫ d A e − β ( U + u L A , c + u L B , c + u B A , c + u B B , c + u I A , c + u I B , c ) (2f) e β Δ G o bound = ∫ bound d B ∫ d A e − β ( U + u L A , c + u L B , c + u B A , c + u B B , c + u I A , c + u I B , c ) ∫ bound d B ∫ d A e − β ( U + u L A , c + u L B , c + u B A , c + u B B , c + u I A , c + u I B , c + u o ) (2g) e β Δ G a bound = ∫ bound d B ∫ d A e − β ( U + u L A , c + u L B , c + u B A , c + u B B , c + u I A , c + u I B , c + u o ) ∫ bound d B ∫ d A e − β ( U + u L A , c + u L B , c + u B A , c + u B B , c + u I A , c + u I B , c + u o + u a ) (2h) e − β Δ G dist + a restr = c ° ∫ bound d B ∫ d A e − β ( U + u L A , c + u L B , c + u B A , c + u B B , c + u I A , c + u I B , c + u o + u a ) ∫ free d B δ ( r − r ∗ ) ∫ d A e − β ( U + u L A , c + u L B , c + u B A , c + u B B , c + u I A , c + u I B , c + u o ) (2i) e − β Δ G o free = ∫ free d B δ ( r − r ∗ ) ∫ d A e − β U + u L A , c + u L B , c + u B A , c + u B B , c + u I A , c + u I B , c + u o ) ∫ free d B δ ( r − r ∗ ) ∫ d A e − β ( U + u L A , c + u L B , c + u B A , c + u B B , c + u I A , c + u I B , c ) (2j) e − β × 2 Δ G I , c free = ∫ free d B δ ( r − r ∗ ) ∫ d A e − β U + u L A , c + u L B , c + u B A , c + u B B , c + u I A , c + u I B , c ) ∫ free d B δ ( r − r ∗ ) ∫ d A e − β ( U + u L A , c + u L B , c + u B A , c + u B B , c ) (2k) e − β × 2 Δ G B , c free = ∫ free d B δ ( r − r ∗ ) ∫ d A e − β U + u L A , c + u L B , c + u B A , c + u B B , c ) ∫ free d B δ ( r − r ∗ ) ∫ d A e − β ( U + u L A , c + u L B , c ) (2l) e − β × 2 Δ G L , c free = ∫ free d B δ ( r − r ∗ ) ∫ d A e − β U + u L A , c + u L B , c ) ∫ free d B δ ( r − r ∗ ) ∫ d A e − β U (2m)

Note that the numerator of a given equation in Eq. 2 is generally the denominator of the previous equation, and Eqs 2k–2m contain contributions for both monomers, thus including a factor of 2. Eq. 2 may be written as a product of averages: K e q = e − β u L A , c ( b o u n d , U ) − 1 e − β u L B , c ( b o u n d , U + u L A , c ) − 1 e − β u B A , c ( b o u n d , U + u L , c ) − 1 e − β u B B , c ( b o u n d , U + u L , c + u B A , c ) − 1 × e − β u I A , c ( b o u n d , U + u L , c + u B , c ) − 1 e − β u I B , c ( b o u n d , U + u L , c + u B , c + u I A , c ) − 1 e − β u o ( b o u n d , U + u c ) − 1 e − β u a ( b o u n d , U + u c + u o ) − 1 × e − β u o ( f r e e , U + u c ) e − β u I , c ( f r e e , U + u L , c + u B , c ) e − β u B , c ( f r e e , U + u L , c ) e − β u L , c ( f r e e , U ) × c ° ∫ bound d B ∫ d A e − β U + u c + u o + u a ∫ free d B δ ( r A B − r A B * ) ∫ d A e − β U + u c + u o . (3)

In the above equation, we use the shorthand notation u _L,c = u _LA,c + u _LB,c (conformational restraint on the backbone atoms for loops 4 and 7 in both monomers A and B), u _B,c = u _BA,c + u _BB,c (conformational restraint on barrel backbone for both monomers), u _I,c = u _IA,c + u _IB,c (conformational restraint on interface residue sidechains for both monomers), and u _c = u _L,c + u _B,c + u _I,c (total conformational restraint). The last term written as a ratio of two integrals will be treated separately below.

Each average corresponds to a free energy change, for example, e − β Δ G L A , c bound = ⟨ e − β u L A , c ⟩ ( b o u n d , U ) , where G L A , c bound is the free energy change due to the addition of the restraining potential u _LA,c on the loops of chain A, in the bound state with no other restraints applied. Because u _c is intended to apply conformational restrictions that may be already partially restricted due to binding itself, the corresponding free energy contributions are expected to be smaller in the bound state than in the free state: G c bound is thus expected to have smaller magnitude than G c free . A similar expectation holds for the other applied potentials. The numbers in Table 1 do not always follow this expectation; however, we discuss this further below.

In terms of free energies, Eq. 3 can be written in the following form by pairing terms containing the same restraining potentials in the bound and free states: K e q = exp − β 2 Δ G L , c free − Δ G L A , c bound − Δ G L B , c bound + 2 Δ G B , c free − Δ G B A , c bound − Δ G B B , c bound + 2 Δ G I , c free − Δ G I A , c bound − Δ G I B , c bound + Δ G o free − Δ G o bound + Δ G dist + a restr − Δ G a bound (4)

The 2nd to last terms in Eq. 4, defined as e − β Δ G dist + a restr in Eq. 2i, can be written as e − β Δ G dist + a restr ≡ c ° S ∗ I ∗ (5) where the term S ^∗ addresses the removal of the relative angular restraints: S ∗ = r ∗ 2 ∫ 0 π d θ sin ( θ ) ∫ 0 2 π d ϕ e − β u a (6) and the term I ^∗ can be recast as the difference in the potential of mean force (PMF) W(r) between the bound and free states, in the presence of the configurational, orientational, and axial restraints: I ∗ = ∫ bound d r A B ′ e − β W r A B ′ − W r A B ∗ . (7)

In the above equations, r _AB is the scalar distance between the centers of masses of proteins A and B (r in Figure 1B), r A B ∗ is an arbitrary fixed location (r ^∗, θ ^∗, ϕ ^∗) in the unbound region far from the other protein, U is the total potential energy of the system in the absence of restraints, potentials with lower case u are restraint potentials, and W(r) is the separation PMF in the presence of all restraints.

Several terms in Eq. 4 can be calculated analytically. For example, when protein B is in the unbound state sufficiently far away from its binding partner, the potential U + u _c is isotropic with respect to rotations about the angles Θ, Φ, and Ψ. This allows the term e − β u o ( f r e e , U + u c ) in Eq. 4 to be calculated as e − β G o free = 1 8 π 2 ∫ 0 π ⁡ sin ( Θ ) d Θ ∫ 0 2 π d Φ ∫ 0 2 π d Ψ e − β u o ( Θ , Φ , Ψ ) , (8) where u _o is a parabolic potential described by u o = 1 2 k Θ ( Θ − Θ 0 ) 2 + 1 2 k Φ ( Φ − Φ 0 ) 2 + 1 2 k Ψ ( Ψ − Ψ 0 ) 2 . (9)

Likewise, S ^∗ in Eq. 6 can be calculated analytically, where u _a (θ, ϕ) in Eq. 6 is a parabolic potential given by u a = 1 2 k θ ( θ − θ 0 ) 2 + 1 2 k ϕ ( ϕ − ϕ 0 ) 2 (10)

In Eqs 2a–2h, the numerators and denominators are partition functions before and after imposing restraints. Thus, the ratio of the partition function gives the free energy cost to impose the restraint. As shown schematically in Figure 2 (terms prior to separation), the restraints are imposed serially in the following order: u _LA,c → u _LB,c → u _BA,c → u _BB,c → u _IA,c → u _IB,c → u _o → u _a .

The term in Eq. 2i includes both the free energy cost due to monomer separation in the presence of all the restraints (related to I ^∗ in Eq. 7) and the free energy gain to release the axial angle restraints (related to S ^∗ in Eq. 6). These terms combine to yield Δ G dist + a restr in Eq. 5 and Figure 2.

In Eqs 2j–2m, the numerators and denominators are partition functions before and after releasing restraints. Thus, the partition function ratio gives the free energy gain (lowering) upon release of the restraint. As shown schematically in Figure 2 (terms after monomer separation), the restraints are released serially in the reverse order of which they were applied: u _o → u _IB,c → u _IA,c → u _BB,c → u _BA,c → u _LB,c → u _LA,c . The free energy change to release the restraints is gained for each monomer in the dimer. Thus, there is a coefficient of 2 in the exponents of Eqs 2k–2m. The exponent of Eq. 4 gives the binding free energy, written now in the order in which the terms are calculated: Δ G bind = − Δ G L A , c bound − Δ G L B , c bound − Δ G B A , c bound − Δ G B B , c bound − Δ G I A , c bound − Δ G I B , c bound − Δ G o bound − Δ G a bound + Δ G dist + a restr + G o free + 2 Δ G I , c free + 2 Δ G B , c free + 2 Δ G L , c free (11)

Each term in Eq. 11 is calculated by the free energy perturbation method. In other words, the potential of mean force (PMF) is calculated as a function of an order parameter, and the effects of imposing or releasing a given restraint are calculated by averaging the Boltzmann factor corresponding to that restraint over the unperturbed potential, as described further below.

Reference structures are prepared as follows:

1) WT E,E (SS): We started from the E,Zn (SS) dimer of chains A and B in PDB structure 1HL4 (Strange et al., 2003), and we removed the Zn ion. The N-terminal acetyl-modification was also removed. Although eukaryotic SOD1 was N-terminally acetylated, the bacterial-expressed SOD1 used in most in vitro biophysical and structural studies was not acetylated (Stathopulos et al., 2003; Ahl et al., 2004; Arnesano et al., 2004; Lindberg et al., 2004; Vassall et al., 2006; Hörnberg et al., 2007; Svensson et al., 2010; Broom et al., 2015b).

For all four apo variants above, the histidine protonation states are 43HSP, 46HSD, 48HSD, 63HSE, 71HSE, 80HSE, 110HSD, and 120HSD, respectively, where HSD is protonated on the delta nitrogen, HSE is protonated on the epsilon nitrogen, and HSP is doubly protonated. These are the states observed in NMR structures [PDB 2AF2 (Banci et al., 2006) and 1L3N (Banci et al., 2002)]. For holo SOD1, the histidine protonation states are 43HSP, 46HSEM, 48HSDM, 63HSN, 71HSEM, 80HSEM, 110HSDM, and 120HSDM in which HSN, HSEM, and HSDM are histidine side chains that must be reparametrized from the putative CHARMM36m force field to facilitate metal binding (Section 2.3). The histidine protonation state is determined by the metal coordination in structure 1HL5 after building hydrogens using the GROMACS module pdb2gmx (Abraham et al., 2015). For all five variants, N- and C-termini are charged (NH3+ and COO−). The monomer reference structure for each SOD1 variant is taken as chain A of the respective dimer reference structure.

To model WT Cu,Zn (SS) SOD1, we reparametrized the coulomb partial charges of all histidines in coordination with metal ions, including histidines 46, 48, 63, 71, 80, 110, and 120. Histidine 63 bridges the Cu and Zn ions in the native structure and is doubly deprotonated (Banci et al., 2002). Such a residue is not present in the putative CHARMM force field. Reparametrized histidines have been used in previous studies for CHARMM27 (Peng et al., 2018) and AMBER and OPLSAA force fields (Wells et al., 2021) but not for the CHARMM36m force field used in this study. Histidines 46, 48, 120 interact with Cu, and histidines 71 and 80 interact with Zn. Reparametrizing these histidines corrects the partial charges due to the charge polarization in the electric fields of the metal ions and allows for proper metal-coordinating geometry consistent with experimental structures. The force field reparametrization is performed using a hybrid approach in which energy gradients (Maple et al., 1988; Waldher et al., 2010) and interaction energy with metal ions (Peng et al., 2018) are constrained to a target value. The partial charges of all atoms in the aromatic rings of the sidechains of metal-coordinating histidines are allowed to relax by minimizing the following loss function: l o s s = w U ( U − U ° ) 2 + w g ∑ α , i ( g α , i − g α , i ° ) 2 , (12) where w _U and w _g are weighting factors set to 1 and 10, respectively. U° is the quantum mechanical metal interaction energy of the subsystem consisting of the metal-coordinating histidine rings H46, H48, H63, H71, H80, and H120 and all atoms in the side chain of D83 in chain A of the holo SOD1 structure 1HL5. The quantum mechanical interaction energy was calculated in GAUSSIAN09 (Frisch et al., 2009) in a previous study (Peng et al., 2018). g _α,i ° is the potential energy gradient at the position of (or equivalently the force on) the Cu and Zn ions and should be zero in all directions i = 1, 2, 3 in the experimentally resolved structure. The subscript α = 1-4 includes the Cu and Zn ions in both chains A and H of structure 1HL5. As mentioned above, the calculation of interaction energy and gradient involves a subsystem in the vicinity of the metals. This includes Cu²⁺, Zn²⁺, all atoms in the histidine rings of the metal-coordinating amino acids (H46, H48, H63, H71, H80, and H120), and all atoms in the side chain of D83.

The partial charges in each reparametrized histidine are allowed to relax within a constrained range relative to the initial charge. The initial charges of H48, H110, and H120 are set to those in HSD in the CHARMM36m force field, and the initial charges of H46 and H71 are set to those of HSE in the CHARMM36m force field. The assignment of HSE or HSD is determined by metal coordination in 1HL5. The initial charge of the doubly deprotonated histidine H63 is designed from HSD in the CHARMM36m force field as follows: the hydrogen on the epsilon nitrogen (HE2) is removed, and then the surplus negative charge and the two previous nitrogen charges (ND1 and NE2) are redistributed evenly on the two nitrogens. Partial charges belonging to histidine 63 are restrained from being within ±0.5e of their initial charges, and the partial charges of the other atoms in the histidine side chains are constrained to be within ±0.1e of their respective initial charges.

Two additional constraints are also applied: 1) the charges of the nitrogens in H63 cannot be lower than −0.7, which is the partial charge of the deprotonated nitrogen in a neutral histidine, and 2) the partial charge of the CG atom in H63 must remain within ±0.1e of its initial charge. This latter constraint prevents the Zn ion from being shielded by the aromatic ring of H63, which prevents the SOD1 electrostatic loop from detaching from the beta barrel. This procedure is implemented in scipy constr-trust minimizer (Conn et al., 2000), and it gives the reproducible parameters in Table 2. With the reparametrized atomic partial charges in Table 2, the dimer and monomer of holo SOD1 remain in correct metal coordination for the full duration of our 200 ns MD simulations. From the last 160 ns of the monomer equilibrium trajectories, we calculated the all-atom root mean squared fluctuations (RMSF) for each amino acid coordinating either metal (H46, H48, H63, H71, H80, D83, and H120). Metalation structurally stabilized the coordinating amino acids, reducing the RMSF from 1.29 ± 1.28 A ̊ to 0.42 ± 0.17 A ̊ , a 67% decrease.

Before carrying out the potential of mean force (PMF) calculations, we obtained properly equilibrated initial structures as follows: three distinct systems required equilibration: bound dimer, interacting monomers during dimer separation, and isolated monomers. Special treatment was also applied to equilibrate the long disordered loops 4 and 7 of SOD1 for both dimer and monomer. All simulations were carried out using the CHARMM36m potential (Huang et al., 2017) in an explicit TIP3P solvent (Jorgensen et al., 1983), using GROMACS 2019.2 (Abraham et al., 2015) patched with PLUMED 2.5.2 (Bonomi, 2019) unless otherwise stated (e.g., Section 4). All simulations in this study were performed on the Sockeye computing cluster (UBC ARC Sockeye, 2019) using NVIDIA Tesla V100 GPU and Intel Xeon Silver 4216 CPU.

1) Dimer: A dodecahedron unit cell with box boundary 1.2 nm distance away from the closest atom on the protein was used for dimer simulations, wherein each variant was solvated with explicit TIP3P water, and K⁺ and CL⁻ ions were added to neutralize the system charge and maintain an aqueous salt concentration of 150 mM. System energy was then minimized through steepest descent until a maximum force < 100 kJ/mol/nm, followed by 300 ps NVT thermostat through the V-rescale method, with 1,000 kJ/mol/nm positional restraints on the heavy atoms. Protein and solvent thermostats had a coupling time of 0.1 ps. A time step of 2 fs was used in all simulations followed by a 300 ps NPT thermostat using the Parrinello−Rahman and V-rescale method with 1,000 kJ/mol/nm positional restraints on heavy atoms. The pressure coupling was isotropic with a coupling time of 2 ps and compressibility of 4.5 × 10⁻⁵ bar⁻¹. Electrostatics was calculated by the PME method with order 4 and Fourier spacing of 0.16. The electrostatics cutoff and van der Waals cutoff were both 1.2 nm. LINCS constraints method of order 4 was applied on heavy atom-H bonds, with iteration set to 1. The temperature and pressure were maintained at 300 K and 1.0 bar, respectively.

Following the NPT thermostat, each dimer variant is relaxed for 200 ns using conventional MD, using the same method as the NPT thermostat, except that positional restraints are no longer applied. Because dimers are restrained progressively throughout the binding free energy calculation using the additional restraint potentials in Eq. 2, the PMF at different stages must be calculated with all the previous restraints present. As a result, equilibrated structures must be prepared with restraints successively applied. A 50 ns MD equilibration with the conformational restraint on loop backbone of chain A (u _LA,c ) was first implemented, followed by another 50 ns equilibration with both u _LA,c and the loop backbone restraint of chain B (u _LB,c ). In addition to the two loop backbone restraints, 10 ns MD equilibrations with the other restraints are then applied successively in the following order: u _BA,c → u _BB,c → u _IA,c → u _IB,c → u _Θ,o → u _Φ,o → u _Ψ,o → u _θ,a → u _ϕ,a .

2) Interacting monomers during dimer separation: The initial protein structures used for constructing the separation PMF were taken from the final structure of the dimer equilibration simulation, with all restraints applied. The simulation unit cell was a dodecahedron with a box boundary of 3 nm from the closest atom on the protein to prevent the protein complex from interacting with its image during umbrella sampling on the separation distance. The procedure of solvation, ionization, energy minimization, and NVT/NPT equilibration followed the same procedure of the dimer. Following this, conventional MD with all restraints applied was run for 40 ns. The protein structures for chains A and B were then translated to various separation distances (Table 3) to be used as initial conditions for replica-exchange molecular dynamics umbrella sampling (REMD-US), as described in Section 2.5.

The calculations of all of the PMFs resulting from the applied restraints used replica-exchange MD combined with umbrella sampling (REMD-US). Distance, RMSD, or angle information obtained using PLUMED is analyzed, and a potential of mean force (PMF) for each reaction coordinate is obtained using the multistate Bennett acceptance ratio (MBAR) algorithm implemented through pymbar (Shirts and Chodera, 2008). After obtaining the PMF for a reaction coordinate, the cost of restraining that reaction coordinate can be computed using one of Eqs 2a–2m.

To implement REMD-US, a series of configurations along the reaction coordinate are generated by biased simulations. These configurations will serve as the starting configurations for the REMD-US windows. Starting from an equilibrated structure, a simulation with an increasing bias center and a simulation with a decreasing bias center are conducted in parallel. In other words, simulations are performed two at a time, with bias centers moving outwards from the original equilibrated structure. This ensures that the conformational changes across reaction coordinates are smooth and continuous.

Special treatment is applied to preparing the REMD-US initial configurations for the loop PMFs, to calculate Δ G L A , c bound , Δ G L B , c bound , and G L , c bound . For loop terms, the two equilibrated structures (one with loops 4 and 7 structured and one with loops 4 and 7 disordered) generate four sets of biased simulations, using the above method of moving bias centers. The two biased simulations that start from the conformation with structured loops generate initial configurations with loop RMSD 0.2–7.0 Å, and the other two biased simulations that start from the conformation with unstructured loops generate initial configurations with loop RMSD 5.2–15.2 Å. The overlapped region, 5.2–7.0 Å, ensures sufficient exchange between umbrellas starting from different initial conformations.

The parameters in each REMD-US are listed in Table 3. Each window in REMD-US was initially run for 20 ns and was extended until the corresponding free energy contribution had converged. In other words, it did not change significantly as simulation time was increased (Supplementary Figure S7). Exchanges between neighboring windows are attempted every 1 ps and are accepted or rejected according to a Metropolis energy criterion.

Because the separation PMF W(r) has a steeper slope at the initiation of dissociation, higher umbrella spring constants are used for this region of rapidly changing PMF.

For all PMFs, the initial 50% of the REMD-US trajectory is discarded in constructing the PMF for calculating the free energy contribution.

All the restraint potentials are harmonic. The same conformational restraints are applied to all four apo SOD1 variants, with parameters given in Table 4. The loop restraints of the holo SOD1 variant bias more closely to the native structure, with restraint center at 1.2 Å instead of 3.5 Å. Otherwise, the conformational restraint parameters are the same as apo ones. Although the conformational restraints for different variants take the same formula, they differ because the RMSD is calculated against different reference structures. The changes in PMFs after applying restraints are shown in Supplementary Figure S6.

The orientation restraint potentials confine the angles Θ, Φ, and Ψ to a variant-specific value, by three separate harmonic potentials u _Θ,o , u _Φ,o , and u _Ψ,o (Table 5). Likewise, the angular restraint harmonic potentials u _θ,a and u _ϕ,a restrain the angles θ and ϕ to values near those given in Table 5. Again, because the reference structure for each variant is different, each variant has a slightly different restraint center.

Given the coordinate system used in Figure 1, another factor that weakly affects the final binding free energy is the choice of r ^∗ in Eqs 6, 7, where r ^∗ is determined as the last point in the separation PMF of each variant. The values used for r ^∗ for the five variants (in the same order as in Table 5) are 38.8, 37.8, 38.8, 38.8, and 38.8 Å.

The accumulated simulation time in this work spent on each process includes equilibration (3.33 µs), serial umbrella construction (2.12 µs), REMD-US (175.46 µs), and force field reparametrization (3.4 µs). The total cumulative simulation time is thus 184.31 µs divided into 111.01 μs of simulation time on the dimer system and 73.3 μs on the monomer system.

The binding free energy ΔG _bind error comes from the error propagation from each component term in Eq. 11. The error could come from two sources: the statistical error from the MBAR estimator (Shirts and Chodera, 2008) or the systematical error caused by insufficient convergence of REMD-US. In this study, the larger error of the two is used, so sufficient convergence of the REMD-US must be achieved to acquire a low enough ΔG _bind error. The detailed error calculation of each term is described in Supplementary Section S1.

The largest free energy contribution is the cost of monomer separation ( Δ G dist + a restr ) . This term also contains a contribution from the relaxation of the axial angle restraints. Δ G dist + a restr is much larger than the experimental value of the binding free energy because of the numerous restraints that minimize unfavorable entropic factors in the binding process, which make sampling appreciably more efficient.

A4V E,E (SS) shows less binding free energy due to Δ G dist + a restr than WT E,E (SS). This is sensible because A4 is adjacent to the dimer interface, and its mutation, depending on the sidechain, could either remove dimer stabilizing interactions or stereochemically disrupt the dimer interface. WT E,E (SH) also has a smaller value of Δ G dist + a restr than that of WT E,E (SS). This is also sensible because residues 50–54 in loop 4 form part of the dimer interface. These are disordered in the disulfide-reduced state, as disulfide bond reduction removes the stable anchor between the loop to the β barrel.

Interestingly, although WT Cu,Zn(SS) SOD1 has more stable structure due to metal coordination and disulfide-bonded loop constraints, it has similar separation contribution Δ G dist + a restr to WT E,E (SS) and D101N E,E (SS) (Table 1). The separation contribution between monomers was similar in our calculations for all variants without mutations in the dimer interface due to the various restraints present in the calculation that minimize differences arising from entropic factors (Zhang et al., 2016).

All the REMD-US simulations converge within 20 ns per window (last row of Supplementary Figure S7). This fast convergence may be attributed to the convex binding interface for SOD1 dimers, in which no entangled loops are present, and the interface is mainly composed of β-sheets. A concave binding pocket, or a binding interface with entangled loops, often leads to artificially high unbinding free energies due to the long-time relaxations required for the structures on the dissociation pathway (Joshi and Lin, 2019; Walther Perthold and Oostenbrink, 2019).

We applied conformational restraints to the flexible loop region (loops 4 and 7) first because their relaxation time is otherwise very slow. We observed that several loop PMFs have a double-well structure (the first column of Figure 3; Supplementary Figure S6), in which the left well consists mainly of well-structured conformations and the right well consists mainly of entropically driven disordered structures. The double-well free energy surface suggests weak two-state-like transitions between ordered and disordered states of the long loops.

For WT E,E (SS) and A4V E,E (SS), the loop contributions from chain A are larger than chain B (i.e., Δ G L A , c bound > Δ G L B , c bound , see Table 1), indicating that the folding of chain A facilitates the folding of chain B in the dimer. This cooperative folding effect is also reflected in the shape of PMF, where chain B has a deeper dip in the left well (RMSD ∼ 4 Å ) than chain A. Moreover, the loop PMF for the free monomer has the widest and deepest free energy for the right well (RMSD ∼ 8 Å ), suggesting that loop regions are further disordered in the free state.

For WT E,E (SH), restraining loops in the free monomer ( Δ G L , c free ) takes the largest energy among all the variants in this study, consistent with this variant’s lack of loop stabilization by the disulfide bond. Surprisingly however, its Δ G L A , c bound and Δ G L B , c bound are the lowest among the apo variants, suggesting that, in the apoprotein, the lack of disulfide bond may lower the free energy of stable structures or the disulfide bond may strain the apoprotein (but not the holoprotein). A similar effect has been observed by us previously (Das and Plotkin, 2013c) using simulated mechanical force spectroscopy probes. In this previous study, the formation of the disulfide bond in the apoprotein weakened the mechanical coupling between the disulfide bonding residues 57/146 and the rest of the protein. In contrast, for holo-SOD1, the presence of the disulfide bond mechanically stabilized those residues. One caveat in interpreting the results here is that the reference structure for WT E,E (SH) is already less compact than that of the other disulfide-bonded variants (Hörnberg et al., 2007), making a direct comparison of the free energy cost to restrain loops less straightforward.

For D101N E,E (SS) mutant SOD1, the free energy cost to restrain the loops on chain B is much larger than the free energy cost to restraint loop on chain A (Table 1), which indicates that this mutation induces anticooperativity in the folding of loops, as opposed to apo WT and apo A4V. The loss of cooperative folding due to this mutation is discussed further in Section 4.1.

For WT Cu,Zn(SS), the loops are greatly stabilized by the metal cations, so a loop conformational restraint with a smaller bias center, 1.2 Å, is imposed (see bottom-left subpanel in Supplementary Figure S6). Despite this tighter restraint, both dimer and monomer loop free energy contributions were the smallest among the variants studied. The shift in free energy surface upon metalation for SOD1 towards a more structured free energy minimum (Figure 3, lower left panel) is consistent with previous experimental results showing that the presence of Zn facilitated the folding of disordered loops (Kayatekin et al., 2008).

The net free energy contribution from loops to the binding free energy is given by 2 Δ G L , c free − Δ G L A , c bound − Δ G L B , c bound (Table 6). Each Δ in the equation is the cost to constrain the loops, so a smaller constraining cost in the dimer means a relatively larger free energy decrease due to conformational relaxation in the monomer versus the dimer. We thus found that loop free energy has a destabilizing effect upon dimerization for all variants studied except for WT E,E (SS). Consistent with previous experimental results (Hörnberg et al., 2007), the loop stability penalty is the largest for the disulfide-reduced variant WT E,E (SH).

Perhaps surprisingly, for WT E,E (SS), the loop free energy does not disfavor dimerization in our calculation, and dimerization has almost no effect on loop stability. Moreover, the loops do not appear to have reduced conformational freedom in the dimer. The conformational fluctuations of the loops in the dimer have an average RMSF of 3.59 Å, while the monomers have a slightly smaller average RMSF of 3.26 Å. Rather than supporting a mechanism of entropy-enthalpy compensation upon dimerization acting on loop conformations, this supports significant conformational freedom of the loops in the WT E,E (SS) dimer.

The β-barrel backbone is the next structural region restrained after the loops, before monomer separation. By comparing the free energy cost to restrain the barrel backbone in monomers versus dimer ( 2 Δ G B , c free − Δ G B A , c bound − Δ G B B , c bound , Table 6), we sensibly found that the barrel backbone generally had larger flexibility in unbound the monomer than in the bound dimer and thus opposed dimerization. Interestingly, the magnitude of this effect was often as large as that of the loops.

D101N E,E (SS) is again exceptional in having a rigid barrel backbone in the unbound monomer, which approaches the stability of the WT Cu,Zn(SS) β-barrel (Table 1). The barrel in the monomer is actually more stable than in the dimer, indicating that the β-barrel conformational free energy favors rather than opposes dimer binding.

A4V E,E (SS) has the least stable β-barrel of the variants in this study (Table 1), and the backbone conformational free energy of A4V E,E (SS) destabilizes the dimer and opposes its formation most strongly of all the variants.

To ensure that the above effects on D101N E,E (SS) and A4V E,E (SS) were not artifacts of insufficient sampling, we used a larger spring constant, k = 50 kcal/mol/Ų, and longer simulation time for the REMD-US method for constructing the PMF for Δ G B , c free (Table 3; Supplementary Figure S7).

The final conformational restraint is applied to all (sidechain and backbone) heavy atoms of the dimer binding interface residues. We found the sensible result that the bound states always had increased interface stability relative to the free state (the third column in Figure 3 and 2 Δ G I , c free − Δ G I A , c bound − Δ G I B , c bound in Table 6), meaning that ordering of the interface sidechains strongly opposes dimer binding. This effect is largely entropic, as the energetic terms mediated by these sidechains (as well as other atoms) that favor dimer binding are accounted for in the PMF calculation for monomer separation (Section 3.1). The enhanced structural order of loops upon metalation also reduces the conformational disorder of interfacial residues as roughly five interface residues reside in loop 4.

The remaining PMFs for orientational and angular restraints in bound states all converge rapidly, with almost perfect overlap between PMFs constructed using either the last 75% or last 50% of the sampling trajectories (Supplementary Figures S1–S5). As expected, all the orientational free energy contributions in the bound state are negligibly small (Table 1). On the contrary, their contributions in the free monomer state, determined analytically and similar for all, are significant. In the free monomer state, the orientational restraints cost 7.62–7.63 kcal/mol. The opposition to dimer binding by restriction of rotational freedom of independent monomers is simple universal free energy cost that is variant and independent.

Because the equilibrium constant in Eq. 1 has dimensions of volume, the calculation involves a volume that is contained in terms S ^∗ and I ^∗ (Eqs 6, 7). The angular restraints are manifested in S ^∗, representing the area available to chain B on the sphere of r ∗ ≈ ( 37.8 − 38.8 ) A ̊ surrounding chain A, and are ≈ ( 5.2 − 5.5 ) A ̊ 2 for all the SOD1 variants. In other words, the solid angle is restrained to be ∼ 1 270 of the entire 4π during the separation.

The binding free energies calculated in this study (Table 1) are systematically weaker than the experimentally determined values, which is an issue reported before (Siebenmorgen and Zacharias, 2019) for the method we have used here. It may be rooted in inaccuracies of the non-polarizable force field we have used here for molecular dynamics simulations (CHARMM36m), particularly when used to evaluate protein binding free energies (Hazel et al., 2018). The mechanically induced unfolding of SOD1 has been observed to have a mechanism that is robust to force field and coarse-grained model for early events but sensitive to the force field and model for late stage unfolding events when the protein is more significantly disordered (Habibi et al., 2016). Differences in binding free energy between SOD1 variants (i.e., ΔΔG _bind) may be more robust to the force field, and we compare these here as well.

The dimer binding free energy ΔG _bind of WT E,E (SS) has been experimentally measured by several research groups. Published values include −12 kcal/mol⁵², −11.0 ± 0.4 kcal/mol at 23°C (Broom et al., 2015b), and −10.3 ± 0.5 kcal/mol at 37°C (Broom et al., 2015b). As mentioned above, this is much larger in magnitude than our calculated value of ≈ − 3.5 kcal/mol.

The conformational entropy increase associated with disulfide reduction for WT E,E (SH) was reported to lead to the dissociation of the apo SOD1 dimer in physiological concentration (Lindberg et al., 2004) and has also been reported to have at least 4 orders of decrease in the association constant (Arnesano et al., 2004), corresponding to about a 5.5 kcal/mol shift towards weaker ΔG _bind. In our calculation, the ΔΔG _bind between WT E,E (SH) and WT E,E (SS) is 4.5 ± 3, which is quite close to this experimental value. The dissociation constant for WT E,E (SH) has been reported to be approximately 85 ± 50 mM based on measured transient populations (Sekhar et al., 2015), which correspond to the binding free energy between −1 kcal/mol and −2 kcal/mol. While these experiments have measured marginal stability for the E,E (SH) dimer, our calculations have yielded a marginal instability for the E,E (SH) dimer of +1.04 ± 0.94 kcal/mol.

The experimental binding free energy ΔG _bind for A4V E,E (SS) SOD1 has been reported as − 7.2 ± 0.2 at 23°C (Broom et al., 2015b), −7.9 ± 0.7 at 25°C (Broom et al., 2015b), and −6.4 ± 0.3 kcal/mol at 37°C (Broom et al., 2015b). It has also been reported that the dissociation constant (K _d ) is in the mM range (corresponding to G _bind ≈ − 4 kcal/mol) (Hörnberg et al., 2007). The ΔΔG _bind of A4V E,E (SS) relative to WT E,E (SS) based on these experiments is 3.9 ± 0.6 kcal/mol at 37°C or 3.8 ± 0.5 kcal/mol at 23°C. In our calculations, the ΔΔG _bind of A4V E,E (SS) relative to WT E,E (SS) is about 5.7 ± 3.3, which is higher than the value from these experiments but still in the approximate experimental range.

The collective increase in the loop, barrel, and interface conformational free energy upon monomerization for all variants studied (Table 6) is consistent with experimental observations of extensive disruption of native structure upon apo SOD1 dimer dissociation (Broom et al., 2015b). These experimental measurements are based on the overall heat capacity and enthalpy changes upon dissociation, so they are not structurally resolved.

We may compare the total conformational free energy change between a variant and WT in both the free monomer and bound dimer states (Table 6). This calculation shows that the D101N mutation on the apoprotein [comparing to WT E,E (SS)] has a stabilizing effect on the free monomer ( Δ ∑ X = L , B , I 2 Δ G X , c free = − 3.55 ± 2.97 ) , and almost no effect in the bound dimer Δ ∑ X = L , B , I ( Δ G X A , c bound + Δ G X B , c bound ) = 0.43 ± 2.23 ).

The thermodynamic effects of D101N have been experimentally resolved for the unfolding of the apo monomer (ΔΔG _D-M = −0.80 kcal/mol) and unfolding of the apo dimer Δ Δ G D - M 2 = − 0.75 kcal/mol) (Byström et al., 2010). However, these numbers couple in the unfolding free energy and rely on a linear extrapolation from 5.8 M urea to 0 M urea (Byström et al., 2010), which does not permit a direct comparison to our values for the dimer binding free energy. Under the additional assumption that the dimer association rate is the same for WT and D101N, the experimental value of the difference in dimer binding free energy mutant to WT is − 2.3 R T ⁡ log 10 ( k d W T / k d m u t ) giving a value of ΔΔG ≈ 0.01 kcal/mol for D101N, or essentially equal to the WT dimer binding free energy.

The difference in the binding free energy of D101N-WT heterodimer to the WT and mutant homodimers, ΔG _het = 2ΔG _WT-mut − ΔG _WT-WT − ΔG _mut-mut, has also been experimentally resolved by Shi et al. (2016) to be −0.71 kcal/mol. Based on our dimer binding free energies, this gives a predicted value for the heterodimer binding free energy of −5.4 kcal/mol for the D101N-WT apo heterodimer.

Our calculations also show that the A4V mutation conformationally destabilizes both bound dimer ( Δ ∑ X = L , B , I ( Δ G X A , c bound + Δ G X B , c bound ) = 0.84 ± 2.16 kcal/mol) and free monomer ( Δ ∑ X = L , B , I 2 Δ G X , c free = 1.48 ± 3.14 kcal/mol). Because the free monomer has larger required constraining free energy than the bound dimer, this conformational disruption further opposes dimer binding. To our knowledge, there is no direct experimental measurement of these free energies. However, it has been shown that A4V is one of the most destabilizing mutants for both monomer and dimer unfolding (ΔΔG _D-M = 1.62 kcal/mol and Δ Δ G 2 D - M 2 = 4.31 kcal/mol) (Lindberg et al., 2005).