Interfacial Dynamics and Growth Modes of β2-Microglobulin Dimers

Protein aggregation is a complex process, strongly dependent on environmental conditions and highly structurally heterogeneous, both at the final level of fibril structure and intermediate level of oligomerization. Since the first step in aggregation is the formation of a dimer, it is important to clarify how certain properties of the latter (e.g., stability or interface geometry) may play a role in self-association. Here, we report a simple model that represents the dimer’s interfacial region by two angles and combine it with a simple computational method to investigate how modulations of the interfacial region occurring on the ns−μs time scale change the dimer’s growth mode. To illustrate the proposed methodology, we consider 15 different dimer configurations of the β2m D76N mutant protein equilibrated with long Molecular Dynamics simulations and identify which interfaces lead to limited and unlimited growth modes, having, therefore, different aggregation profiles. We found that despite the highly dynamic nature of the starting configurations, most polymeric growth modes tend to be conserved within the studied time scale. The proposed methodology performs remarkably well taking into consideration the nonspherical morphology of the β2m dimers, which exhibit unstructured termini detached from the protein’s core, and the relatively weak binding affinities of their interfaces, which are stabilized by nonspecific apolar interactions. The proposed methodology is general and can be applied to any protein for which a dimer structure has been experimentally determined or computationally predicted.


■ INTRODUCTION
Protein aggregation is the process in which soluble protein conformations (monomers) with exposed hydrophobic patches self-associate into dimers and higher-order oligomers. 1 Aggregation may lead to amorphous aggregates with a granular appearance, protofibrils (including annular oligomeric aggregates), or other oligomeric aggregated states. Often, the end product of protein aggregation is amyloids, insoluble aggregates comprising long unbranched fibers, characterized by the crossbeta structure. 2 While the formation of amyloids is the hallmark of a plethora of conformational disorders, including Parkinson's and Alzheimer's diseases, 3 it is thought that the most cytotoxic species within the aggregation cascade are the early forming oligomers. Indeed, the latter not only cause an overload of the proteostasis machinery but also can seed the aggregation of identical proteins, cross-seed the aggregation of homologous proteins, and even disrupt biological membranes (reviewed in ref 4).
Since protein aggregation typically starts with the selfassociation of two monomers into a dimer, by blocking dimerization one could -at least in principle -hinder aggregation. Therefore, it is critical to have a molecular-level understanding of the mechanism of dimer formation and to establish how certain traits of the dimer (including structure, conformational dynamics, energetic and thermodynamic stability, etc.) may affect aggregation and its final outcome. It has been suggested that a properly formed dimer with increased stability may be protective against aggregation. 5 Additionally, it has been observed that a larger conformational variability at the dimer level appears in association with a faster aggregation process and higher toxicity. 6 The identification and structural characterization of early formed oligomers is particularly challenging because they exist in a complex dynamic equilibrium with each other and with insoluble higher order aggregates. 7,8 An additional difficulty is that they form transiently, in a time scale that is not compatible with the temporal resolution of commonly used biophysics apparatus. Furthermore, their formation is under kinetic control, being strongly dependent on environmental conditions. 9,10 Computational modeling covering different resolutions and time scales has been playing a pivotal role in helping understand protein aggregation 11,12 and in providing predictions for the structure of the oligomeric species formed along the aggregation pathway, 7,13,14 and for the understanding of fibril growth. 15 In recent years, some of us developed an integrated computational methodology (reviewed in ref 16) combining data from discrete molecular dynamics (DMD) of a structurebased model, 17−19 constant-pH MD, 20−23 and a rigid-body Monte Carlo ensemble docking method developed in house (MC-ED) 17 that uses a low-resolution cost function based on step potentials to model electrostatic, hydrophobic, and hydrogen bond interactions. 24 The MC-ED provides a probabilistic view of the dimerization phase, in particular, which regions of the monomer are more likely to associate, and which residues play a key role in aggregation acting as aggregation hot spots. We applied it to study the dimerization phase of protein beta-2-microglobulin (β 2 m) under different environmental conditions, 16 but the method is general and, in principle, can be used for any model system. β 2 m is a small (99 residues long) protein with biomedical interest. 16 It has a classical beta-sandwich fold in which the native structure is stabilized by a disulfide bridge between the sulfur atoms of the cysteine residues 25 (at the B strand) and 80 (at the F strand) 25 ( Figure 1A). The disulfide bridge has been considered critical for β 2 m fibrillogenesis. 26,27 The wild-type form (wt-β 2 m) is the canonical agent of dialysis-related amyloidosis, a disease that affects individuals with kidney impairment undergoing long-term hemodialysis, 28 while the single point mutant D76N is the causing agent of a familial amyloidosis affecting the visceral organs. 29 By exploring the folding transition of D76N with DMD simulations of a structure-based Go ̅ model, 30 we were able to predict an intermediate state with exposed hydrophobic patches (i.e., with the potential to trigger self-association), which is structurally characterized for having a well-preserved core and two unstructured and detached terminal regions. 30 For this reason, we termed this intermediate I 2 (Figure 1B). Afterward, a study that combined in crystallo with in silico analyses 31 reported a highly dynamic conformation for D76N where the loss of the beta structure at both the N-and C-terminal strands causes the exposition of very aggregation-prone regions, in consonance with our results. More recent DMD folding simulations indicate that the wt-β 2 m can also populate a similar conformer provided the disulfide bond is established. 32 While the role of the experimentally detected intermediate state in the aggregation pathway of D76N still lacks consensus, 34 we decided to use the topologically similar I 2 as a model system to explore the dimerization phase of β 2 m. 16 By conducting extensive protein−protein docking simulations, we found a prominent role for the DE-and EF-loops as adhesion regions, with the N-terminal region becoming also relevant at acidic pH. 30 Our analyses also predict that Trp60, located at the apex of the DE-loop, is no longer the leading hot spot at acidic pH, losing its role for Arg3 located at the N-terminus. 24 Given the coarse-grained nature of the deployed cost function, the MC-ED method should not be used to predict accurate structural dimer models. Motivated by this limitation, we recently developed a computational protocol that combines relatively short (100 ns) Molecular Dynamics (MD) simulations with MM-PBSA calculations. 14 The MD simulations correct structural errors (e.g., steric clashes), while the MM-PBSA method provides accurate interfacial energies for the relaxed dimers. Additionally, by using PyMol 35 and a simple geometric protocol to replicate the dimers interface, we could determine by visual inspection the binding modes (i.e., dimer interfaces) leading to self-limited growth (i.e., polymeric chains that close upon themselves) and unlimited growth (i.e., polymeric chains that can grow indefinitely) in the absence of structural rearrangement. By deploying this methodology to an ensemble of 212 dimers of the D76N mutant of β 2 m obtained with the MC-ED, we found that they are mainly stabilized by apolar interactions and that both the N-and C-termini are involved in the interface of the 10 most stable dimer interfaces. 14 We also found that some of the considered dimer interfaces, including the most energetically stable, lead to self-limited growth, while others can grow indefinitely in the absence of structural rearrangement.
Here, we go a step further by developing a method that allows the determination of a dimer's growth mode from the values of two angles that provide a minimal description of the dimer's interface. By using a simple computational protocol, one is able to generate the three-dimensional representation of the corresponding polymerized chain. To illustrate our method, we consider the ten most stable binding modes of the D76N mutant of β 2 m identified in ref 14, as well as five additional binding modes of intermediate and lower stability, and investigate how the growth mode is modulated by the geometric changes occurring at the dimer's interface in the ns−μs time scale. We observe that for most binding modes (including the less stable ones), the growth mode is conserved within the  ■ METHODS Molecular Dynamics Simulations. All molecular dynamics (MD) simulations were performed using GROMACS 2018.6, 36 the GROMOS 54A7 force field, 37,38 and the SPC water model. 39 The starting conformation of our simulations is a relaxed homodimer of the I 2 intermediate, ∼20k water molecules, and 2 Na + ions. We chose a set of 15 dimers previously investigated 14 that consisted of the following: the top 10 dimers with the lowest binding energies (best affinity) (∼−83 to ∼−67 kcal mol −1 ); 2 dimers with very high energies (weak affinities; ∼−10 kcal mol −1 ); and 3 dimers with average binding energies (∼−43 kcal mol −1 ). To improve our sampling, we performed 3 replicates of 500 ns per dimer. Starting from the dimer simulations previously published, 14 the first replicate went straight into production. The starting structures from the second and third replicates were obtained by performing a new velocity initialization (50 ps with a random seed) before proceeding to the long production runs. The production MD was conducted with a 2 fs time step, and the electrostatics were treated with the Particle-Mesh Ewald (PME) method, 40,41 with a verlet scheme cutoff of 1.4 nm, a Fourier grid spacing of 0.12 nm, and an interpolation order of 4 (cubic). van der Waals interactions were simply truncated above 1.4 nm. 42 The neighbor list was updated every 10 integrator steps. The LINCS algorithm 43 was used to constrain all the protein bonds, and SETTLE was used on the water molecules. 44 Structural Analysis. All structural analyses of the dimers were carried out using the GROMACS package and other inhouse tools. The periodic boundary conditions effects were corrected using the fixbox tool. 45 When dealing with equilibrium properties, the first 250 ns of the simulation were discarded. For the structural alignments and Root Mean Square Deviation (RMSD) calculations, only the Cα atoms from the protein core were considered. These correspond to residues 23−27, 36−39, 51−55, 62−66, and 78−82 of each monomer (strands B−F of Figure 1B). This approach discards spurious structural variability due to loops and unstructured N-and C-terminal regions. As the RMSD reference structure, we used the initial structure, which corresponded to the relaxed structure obtained in our previous work. 14 The interfacial area of the dimers was calculated using the difference between the SASA values of both monomers in the presence and absence of the partner. 46,47 This procedure was also applied to each residue individually, with the final area being normalized (converted into a percentage) by the maximum SAS area of each residue type observed in all simulations. Residues with an interfacial area higher than 10% were considered to be located at the interface.
The asphericity of each monomer was calculated using the radius of gyration (R g ) and the eigenvalues of the gyration tensor (λ i ), given by the square of the three principal radii of gyration (R x , R y , R z ) obtained using the GROMACS tool polystat. The asphericity of each monomer, A s , is given by eq 1 48, 49 A trT . Values of A s higher than 0 indicate a deviation from the ideal sphere shape of our protein. Since 0 ≤ A s ≤ 1, we can obtain the sphericity, ψ, in percentage, by taking To obtain the average structures of each equilibrated dimer replicate, the GROMACS clustering tool was employed with an RMSD cutoff of 2 nm, following the GROMOS clustering algorithm. 50 All error bars indicate the standard error of the mean. Images were rendered using PyMOL 35 and Blender. 51 Calculation of Binding Energies with PyBindE. All binding free energy (E bind ) calculations were performed using PyBindE, a new in-house implementation of the Molecular Mechanics Poisson−Boltzmann Surface Area (MM/PBSA) method, 52,53 written in the programming language Python (https://github.com/mms-fcul/pybinde). We used a dielectric constant of 4, a grid scaling factor (scale) of 2, and a convergence criterion (convergence) of 0.01 kT/e. 14 Additionally, the Poisson−Boltzmann (PB) calculations used 500 and 50 linear (nlit) and nonlinear (nonit) iterations, respectively. The grid size was calculated as double the maximum value of the 3 atomic spatial coordinates, plus 1 if the final number is even, while the center was located in the geometric center of all the atomic coordinates of the system. For each dimer, the MM-PBSA calculations were performed every 100 ps, which resulted in 3× 5000 frames per dimer.
A Simple Protocol for Predicting and Visualizing Dimer Replication. In our previous study, 14 we deployed a simple replication protocol that allows predicting the threedimensional structure of a polymerized chain of dimers forming in the absence of structural rearrangement. In particular, we used PyMOL to visualize the structure of a polymerized chain of β 2 m dimers that was created by using a simple geometric rule ( Figure  2). In doing so, we found that some dimer interfaces lead to chains that close upon themselves (self-limited growth), while other interfaces can sustain indefinite growth. 14 A Simple Model for Dimer Polymerization. In the present study, we develop a simple model for chain polymerization and combine it with the protocol for dimer replication. In doing so, we are able to predict in a systematic manner which dimer interfaces lead to self-limited or unlimited growth modes. To generate a polymerized chain of dimers, we consider a starting dimer conformation formed by monomers A and B. We subsequently create a duplicate dimer conformation where monomers A and B are respectively renamed A′ and B′. Finally, we use PyMOL to structurally align monomer A′ from the duplicate dimer with monomer B of the preceding dimer and repeat this operation a desired number of times.
In real conditions, the formation of the fibrils may differ from consecutive additions of conformationally similar monomers, highlighting how complex and vast the polymerization process is. In our model, we reduced the polymerization process to a simple propagation of growth following the monomer interfaces to reduce the dimension of the process and allow the study of this problem. The model represents a dimer by two spheres of unit radius (one sphere per monomer) and reduces the description of dimer polymerization to two angles: a polymerization angle (θ pol ), which determines the structure of the polymerized chain, and a polymerization dihedral angle (ϕ pol ), which establishes the relative orientation of the two monomers in the dimer. The polymerization angle is defined by where V ⃗ IA A is the vector between the geometric center of monomer A (GC A ) and the geometric center of the set of interfacial residues of A which are in contact with monomer B (GC IA ), and V ⃗ IA B is the vector between the geometric center of monomer B (GC B ) and GC IA mapped onto the surface of monomer B ( Figure 3A,B). On the other hand, the polymerization dihedral is given by Consider two initial spheres representing two monomers A and B. One of the spheres, say monomer A, has its center placed at the origin of a Cartesian reference frame, while the other has its center placed at (−2,0,0). A pair of angles (θ pol , ϕ pol ) is chosen for this dimer. By duplicating the dimer thus created, one can use the previously described geometric protocol that aligns sphere A′ from the duplicate dimer with monomer B of the preceding one ( Figure 2). By iterating this procedure a certain number of times, one obtains a sphere model for dimer growth ( Figure 4). By using this simplified model and the protocol for dimer replication, it becomes possible to map out self-limited and indefinite growth modes for a pair of parameters (θ pol , ϕ pol ) ( Figure 5). Considering several (θ pol , ϕ pol ) pairs separated by 5°w ithin the range 0°−180°, we were able to generate several models of a polymerized chain and identify its growth mode. A growth mode is classified as limited for the models exhibiting overlapping spheres. Within the 0°−60°range, the polymerization mode is, in general, of limited type, while the combination of θ pol > 60°with some dihedral values leads to unlimited growth. Based on this analysis, we created what we term by growth landscape, a two-dimensional surface generated by θ pol and ϕ pol that outlines the regions where growth is limited, uncertain, or unlimited ( Figure 5). The transition points between limited and unlimited growth modes are shown as dots, and we used an exponential fit to this data set to draw the delimitation region.
This simple approach allows one to predict the type of growth mode associated with any dimer for which a conformation has been experimentally determined or computationally predicted because it suffices to calculate the angles θ pol and ϕ pol for that particular conformation. When a dimer falls into the uncertain  region, it is necessary to visually confirm the type of growth mode by using the previously described protocol for dimer replication. We note that the observed uncertainty is mainly due to deviations of the monomer from the spherical shape considered in the simple model. Those deviations can be quantified by the monomer's sphericity described in the previous section.

■ RESULTS AND DISCUSSION
Binding Energy and Structural Analysis. In our previous work, 14 we performed 100 ns MD simulations to relax an ensemble of 212 β 2 m dimer configurations obtained with the MC-ED method (the Relax configurations). We evaluated the dimer's binding energies, investigated which regions of the monomer are most likely present in the dimer's interface, identified the so-called hot-spot residues (the residues that establish a higher number of intermolecular interactions), and determined which (of the most stable) interfaces could be propagated indefinitely allowing the formation of long, polymerized chains. Here, we will study long-term dimer stability and determine to which extent the latter is a determinant of the polymerization growth mode. This is an important question, whose answer will contribute toward a better understanding of structural heterogeneity in protein aggregation. In doing so, we will also address the methodological problem of gauging the quality of dimer conformation obtained with the MC-ED, which uses a low-resolution cost function to generate a dimer interface. Indeed, if a short MD simulation proves to be enough to relax the dimer conformation obtained with MC-ED and establish its polymerization growth mode, the dimer conformation can be considered as a good predictive model. This is of practical importance since the MC-ED provides a dimer interface at a low computational cost.
To explore these questions, we performed long (3 × 500 ns) MD simulations of a subset of 15 dimers (the Equil conformations) that includes the 10 dimers with the lowest binding energy, 3 dimers with average binding energy, and the 2 Figure 5. Growth landscape as a function of the θ pol angle and ϕ pol dihedral angle. The black dots represent the upper limit of angle/ dihedral combinations that yield a limited growth polymer as predicted by the simple model. The orange region represents the region where polymer growth is limited, the yellow color represents the uncertain region (one that requires visual inspection), and the region in white represents the unlimited growth. The Relax data was obtained from our previous work, 14 while the Equil data refers to the equilibrated segment (last 250 ns) of the long MD simulations carried out in the present study. The % shown in parentheses in the Equil E bind corresponds to the variation from the Relax value. dimers with the highest binding energy from the ensemble studied in ref 14 (Table 1). We generically refer to these dimers as binding modes (BMs). For each dimer, we evaluated the binding energies together with several structural properties, such as RMSD, interfacial area, and secondary structure content. Perhaps not surprisingly, the most stable BMs equilibrate relatively fast (100 ns), with the MD simulations starting from the most unstable BMs requiring more time to equilibrate (Figures S1−S4 of the Supporting Information). The RMSD ( Figure S1 of the Supporting Information) is highly sensitive to small structural and conformational rearrangements, making it challenging to achieve equilibration and/or replicate convergence. Nevertheless, even the most dynamic BMs/replicates roughly equilibrated within ∼250 ns as hinted even by their RMSD time evolution. Furthermore, to quantify the structural stability of the interfacial region in the long MD simulations, we calculated the percentage of residues, present in the relaxed interface, 14 that are conserved after the 500 ns equilibration simulation (Table 1).
With the exception of BM-3 and BM-8, whose binding energies increased by 10% and 18%, respectively, after the long MD equilibration, all the other dimer interfaces considered in this study became more stable, with the highest increase in stability (55%) being observed for H-2, followed by A-3 (37%) and BM-2 (23%). The increase in stability of the less stable BMs was expected as these were probably local minima that should rearrange into more stable configurations in our time scale. The small losses of binding affinities observed in BM-3 and BM-8 were attributed to changes in at least one C-terminus that moved away from the interface during the long simulations and were determinant for the overall stability gain. The energetic stability of the interfaces is also related to their size, with the most stable BMs typically exhibiting the highest interfacial area. This is consistent with our previous finding that the van der Waals interactions are a major source of stability for these dimer interfaces. 14 With the exception of BM-9, all the other most stable BMs were able to conserve a moderate (70%) to high (85%) fraction of their starting interface, which can thus be classified as being moderately stable from a structural point of view. Interestingly, although the two less stable dimers (H-1 and H-2) were able to decrease their binding energies after the full equilibration, the 500 ns time scale was not long enough to observe the structural rearrangement required to significantly decrease their binding energy to values comparable with BMs 1−10.
To gauge the quality of the data obtained from the previous MD relaxation step, 14 we calculated the correlation between the binding energies and the average binding energies during the long equilibration step (calculated in intervals of 10 ns) ( Figure  6). The high correlation values observed even after 500 ns of MD simulations indicate that, although the binding interfaces may have changed considerably, the ranking based on their binding energies is not significantly affected. Altogether, these results suggest that the relaxation protocol adopted in our previous study 14 is sufficiently accurate in ranking the different binding modes by their binding energies.
Interfacial Dynamics and Polymerization Modes. The results presented in the previous section show that the interfacial region of the considered dimers changed both in terms of energetic stability and residue composition during the 500 ns MD simulations. It is thus natural to ask if these modifications are able to drive different growth modes. In other words, does the growth mode change with time? In order to answer this question, one needs to be able to determine the growth mode of any dimer conformation extracted from a MD simulation. By evaluating the polymerization angle (θ pol ) and dihedral (ϕ pol ) for each dimer conformation, we can then map them onto the growth landscape in the region (limited, unlimited, or uncertain) that determines the growth mode ( Figure 5). We recall that the latter was constructed by reducing monomers to spheres. For many conformations, this approximation will work fine, but it breaks down when the dimer conformation deviates strongly from this geometric ideal, either due to rearrangements within the protein core or due to protrusions of the protein termini. In this case, visual inspection is needed to confidently classify the growth mode.
To characterize the interfacial dynamics and its relation with the growth mode, we have averaged the polymerization angle θ pol and polymerization dihedral ϕ pol of the equilibrated segment of each replicate. We have then identified the dimer conformation within the ensemble of sampled dimers whose θ pol and ϕ pol more closely match the average value. We term it the average conformation (or average dimer). With this approach, we can evaluate the dispersion of both angles during the simulation and compare the results of the present simulations with those obtained in the shorter MD relaxation step reported in ref 14. Despite a few exceptions, there is an overall consistency  (Table 2).
Our results show that, for the most stable binding modes within the analyzed set (BMs 1−10), the polymerization angle (θ pol ) tends to change within an ∼30°range ( Figure S5 of the Supporting Information). This is expected since significant changes in the residue composition of the dimer's interfacial region are usually required to trigger larger variations to the polymerization angle. An exception to this observation is illustrated by BM-3, which conserves only 61.2% of the original interface and shows a relatively high dispersion in the polymerization angles (large error values) for all replicates ( Figure 8A). Nevertheless, the large θ pol and ϕ pol polymerization angles ensure a low number of clashes between monomers, leading to an unambiguous assignment of unlimited growth.
A more detailed analysis of replicates 1 and 3 shows a high θ pol (∼120°) combined with a high ϕ pol (∼150°), which lead to a linearized growth mode. Interestingly, the smaller θ pol (∼80°) in Table 2. Average Polymerization Angle, θ pol , Dihedral Angle, ϕ pol , Sphericity of Monomers A (ψ A ) and B (ψ B ), and Type of Growth Mode (GM) -Limited (L) or Unlimited (U), for the Average Structure (Equil) a Within the L and U growth modes, we can further classify these structures as h2h (head-to-head), don (doughnut), hel (helix), and lin (linear) (see Figure 7 and Figure S6 of the Supporting Information). Also shown is the GM (Relax) reported in our previous work. 14

Journal of Chemical Information and Modeling
pubs.acs.org/jcim Article replicate 2 caused a more pronounced rotation/distortion in the linear structure. These findings indicate that larger θ pol angle values will favor less structured growth modes, eventually converging to fully linear and unlimited structures. The BM-2 is quite interesting owing to the structural sphericity of monomer B being significantly decreased (up to ∼30%) in replicates 1 and 3 ( Table 2). Since this BM has a θ pol ∼ 60°, it falls into the uncertain region of the polymerization landscape ( Figure 8B). However, since monomer B, which is the one that is being consecutively added, has its detached termini in the dimers' interfacial region, it deviates from the ideal spherical interface of the simple model. As a result, both replicates 1 and 3 have only a small structural overlap and are able to polymerize linearly, albeit being in a low θ pol region. However, this sphericity deviation is smaller in replicate 2, and its limited growth prediction is quite clear.
BM-9 also shows heterogeneity between replicates in their growth modes ( Figure 8C). Similar to BM-2, we observe some loss of sphericity in monomer B, specially in replicate 3 (Table  2), due to the detachment of one of its termini that became part of the interfacial region. This leads to a decrease in the θ pol polymerization angle, bringing some uncertainty to the growth mode viability. However, a visual inspection shows that, despite a few structural clashes that can be mitigated with small conformational rearrangements, this BM can easily grow into an unlimited helical pattern.

■ CONCLUSION
As we have observed in our previous work on β 2 m, 14 dimer structures can be quite dynamic. Hence, it is reasonable to propose that structural modulations of the dimers' interface (or binding mode) change the way it grows or polymerizes, with potential impact on aggregation (e.g., by contributing to its heterogeneity or actually blocking it). In order to explore this scenario, we developed a simple model that represents the geometry of the dimers' interface by two angles, a polymerization angle (θ pol ) and a polymerization dihedral (ϕ pol ), and calculated a (theoretical) growth landscape for idealized dimers (i.e., dimers formed by spherical monomers). The growth landscape is the space spanned by the two angles and contains Figure 8. Growth landscape of BM-3 (A), 2 (B), and 9 (C). The black, blue, and red dots represent the ensemble of 250 dimer conformations extracted from the equilibrated part of replicates 1, 2, and 3, respectively. The structure closer to the average θ pol and ϕ pol values of each replicate is shown in the corresponding color tone (gray, cyan, and dark red). These are the structures used to generate the polymerization modes represented on the right-hand side of each growth landscape. The green triangle represents the starting structure obtained from the initial MD relaxation step. three regions. A region corresponding to self-limited growth, a region corresponding to unlimited growth, and a section in between which we termed the uncertain region. Given any protein dimer conformation/configuration, it suffices to determine its polymerization angles and map them on the growth landscape. When they fall into the uncertain region, an additional visual inspection step based on the replication protocol becomes necessary to establish the growth mode.
In the present study, we employed the proposed methodology to dimers of the D76N mutant of protein β 2 m. In particular, we conducted long (3 × 500 ns) MD simulations of 15 β 2 m dimer binding modes exhibiting different interfacial stabilities, starting from the relaxed conformations, which were previously obtained. 14 With only two exceptions, all the considered dimers adopted more energetically stable interfaces, with stability increasing up to ∼23% for the top-10 dimers or ∼55% for the less stable ones. Despite these changes, the original ranking of BMs based on interfacial stability was not significantly altered. Since nonspecific apolar interactions are the main drivers for these BMs, 14 we continue to find a high correlation between the dimer stability and the interfacial area (r = 0.93). However, we could not establish a clear relation between the percentage of the conserved interfacial area (61−85%) and the corresponding change in interfacial stability (r = 0.25), indicating that all these BMs are highly dynamic in nature, even the most stable ones.
To investigate if the growth mode changes within the time span of the MD simulations, an ensemble of 250 dimer conformations was extracted from the equilibrated part of each MD replicate (1 conformation/ns). For each conformation, we measured the polymerization angle and dihedral and mapped them onto the growth landscape. Since it is not feasible to replicate the dimer interface for the total 11250 dimer conformations, we took a conservative approach and considered the conformation whose polymerization angle and dihedral more closely match the averaged ensemble angles per replicate. This methodology also allows establishing some level of comparison with the results obtained in our previous study. 14 We noticed that for some binding modes (e.g., BM-9) the dispersion around the average is considerably large with different regions of the growth landscape being populated for the same replicate. This is a clear indication that the growth mode of the considered dimers can change within the ns−μs time scale, contributing to the heterogeneity of the aggregation process. When considering the dimer conformation representing the ensemble average, we observe consistency with our previous results, 14 which can be taken as an indication that a methodology that combines MC-ED with a standard relaxation MD protocol was able to deliver proper dimer models. Altogether, our approach led to a robust classification of which dimer configurations have limited and unlimited growth modes, with different consequences to their aggregation potential and biological impact.

■ ASSOCIATED CONTENT Data Availability Statement
The GROMACS package is freely available software used to perform MD simulations and can be downloaded at https:// manual.gromacs.org/documentation/2018.6/download.html. PyMOL v2.5 is also free software for molecular visualization and generating high quality images. It can be downloaded from https://pymol.org/2. Blender is a free and open-source 3D creation suite, licensed as GNU GPL, and it can be downloaded from https://www.blender.org/.