Limitations of field-theory simulation for exploring phase separation: The role of repulsion in a lattice protein model.

Field-theory simulation by the complex Langevin method offers an alternative to conventional sampling techniques for exploring the forces driving biomolecular liquid-liquid phase separation. Such simulations have recently been used to study several polyampholyte systems. Here, we formulate a field theory corresponding to the hydrophobic/polar (HP) lattice protein model, with finite same-site repulsion and nearest-neighbor attraction between HH bead pairs. By direct comparison with particle-based Monte Carlo simulations, we show that complex Langevin sampling of the field theory reproduces the thermodynamic properties of the HP model only if the same-site repulsion is not too strong. Unfortunately, the repulsion has to be taken weaker than what is needed to prevent condensed droplets from assuming an artificially compact shape. Analysis of a minimal and analytically solvable toy model hints that the sampling problems caused by repulsive interaction may stem from loss of ergodicity.

Finite-size shifts in simulated protein droplet phase diagrams.

Computer simulation can provide valuable insight into the forces driving biomolecular liquid-liquid phase separation. However, the simulated systems have a limited size, which makes it important to minimize and control finite-size effects. Here, using a phenomenological free-energy ansatz, we investigate how the single-phase densities observed in a canonical system under coexistence conditions depend on the system size and the total density. We compare the theoretical expectations with results from Monte Carlo simulations based on a simple hydrophobic/polar protein model. We consider both cubic systems with spherical droplets and elongated systems with slab-like droplets. The results presented suggest that the slab simulation method greatly facilitates the estimation of the coexistence densities in the large-system limit.

When a foreign gene meets its native counterpart: computational biophysics analysis of two PgiC loci in the grass Festuca ovina.

Duplicative horizontal gene transfer may bring two previously separated homologous genes together, which may raise questions about the interplay between the gene products. One such gene pair is the "native" PgiC1 and "foreign" PgiC2 in the perennial grass Festuca ovina. Both PgiC1 and PgiC2 encode cytosolic phosphoglucose isomerase, a dimeric enzyme whose proper binding is functionally essential. Here, we use biophysical simulations to explore the inter-monomer binding of the two homodimers and the heterodimer that can be produced by PgiC1 and PgiC2 in F. ovina. Using simulated native-state ensembles, we examine the structural properties and binding tightness of the dimers. In addition, we investigate their ability to withstand dissociation when pulled by a force. Our results suggest that the inter-monomer binding is tighter in the PgiC2 than the PgiC1 homodimer, which could explain the more frequent occurrence of the foreign PgiC2 homodimer in dry habitats. We further find that the PgiC1 and PgiC2 monomers are compatible with heterodimer formation; the computed binding tightness is comparable to that of the PgiC1 homodimer. Enhanced homodimer stability and capability of heterodimer formation with PgiC1 are properties of PgiC2 that may contribute to the retaining of the otherwise redundant PgiC2 gene.

Finite-size scaling analysis of protein droplet formation.

The formation of biomolecular condensates inside cells often involve intrinsically disordered proteins (IDPs), and several of these IDPs are also capable of forming dropletlike dense assemblies on their own, through liquid-liquid phase separation. When modeling thermodynamic phase changes, it is well known that finite-size scaling analysis can be a valuable tool. However, to our knowledge, this approach has not been applied before to the computationally challenging problem of modeling sequence-dependent biomolecular phase separation. Here we implement finite-size scaling methods to investigate the phase behavior of two 10-bead sequences in a continuous hydrophobic-polar protein model. Combined with reversible explicit-chain Monte Carlo simulations of these sequences, finite-size scaling analysis turns out to be both feasible and rewarding, despite relying on theoretical results for asymptotically large systems. While both sequences form dense clusters at low temperature, this analysis shows that only one of them undergoes liquid-liquid phase separation. Furthermore, the transition temperature at which droplet formation sets in is observed to converge slowly with system size, so that even for our largest systems the transition is shifted by about 8%. Using finite-size scaling analysis, this shift can be estimated and corrected for.

Stability and Local Unfolding of SOD1 in the Presence of Protein Crowders.

Using NMR and Monte Carlo (MC) methods, we investigate the stability and dynamics of superoxide dismutase 1 (SOD1) in homogeneous crowding environments, where either bovine pancreatic trypsin inhibitor (BPTI) or the B1 domain of streptococcal protein G (PGB1) serves as a crowding agent. By NMR, we show that both crowders, and especially BPTI, cause a drastic loss in the overall stability of SOD1 in its apo monomeric form. Additionally, we determine chemical shift perturbations indicating that SOD1 interacts with the crowder proteins in a residue-specific manner that further depends on the identity of the crowding protein. Furthermore, the specificity of SOD1-crowder interactions is reciprocal: chemical shift perturbations on BPTI and PGB1 identify regions that interact preferentially with SOD1. By MC simulations, we investigate the local unfolding of SOD1 in the absence and presence of the crowders. We find that the crowders primarily interact with the long flexible loops of the folded SOD1 monomer. The basic mechanisms by which the SOD1 ß-barrel core unfolds remain unchanged when adding the crowders. In particular, both with and without the crowders, the second ß-sheet of the barrel is more dynamic and unfolding-prone than the first. Notably, the MC simulations (exploring the early stages of SOD1 unfolding) and the NMR experiments (under equilibrium conditions) identify largely the same set of PGB1 and BPTI residues as prone to form SOD1 contacts. Thus, contacts stabilizing the unfolded state of SOD1 in many cases appear to form early in the unfolding reaction.

Fitting a function to time-dependent ensemble averaged data.

Time-dependent ensemble averages, i.e., trajectory-based averages of some observable, are of importance in many fields of science. A crucial objective when interpreting such data is to fit these averages (for instance, squared displacements) with a function and extract parameters (such as diffusion constants). A commonly overlooked challenge in such function fitting procedures is that fluctuations around mean values, by construction, exhibit temporal correlations. We show that the only available general purpose function fitting methods, correlated chi-square method and the weighted least squares method (which neglects correlation), fail at either robust parameter estimation or accurate error estimation. We remedy this by deriving a new closed-form error estimation formula for weighted least square fitting. The new formula uses the full covariance matrix, i.e., rigorously includes temporal correlations, but is free of the robustness issues, inherent to the correlated chi-square method. We demonstrate its accuracy in four examples of importance in many fields: Brownian motion, damped harmonic oscillation, fractional Brownian motion and continuous time random walks. We also successfully apply our method, weighted least squares including correlation in error estimation (WLS-ICE), to particle tracking data. The WLS-ICE method is applicable to arbitrary fit functions, and we provide a publically available WLS-ICE software.

Markov modeling of peptide folding in the presence of protein crowders.

We use Markov state models (MSMs) to analyze the dynamics of a ß-hairpin-forming peptide in Monte Carlo (MC) simulations with interacting protein crowders, for two different types of crowder proteins [bovine pancreatic trypsin inhibitor (BPTI) and GB1]. In these systems, at the temperature used, the peptide can be folded or unfolded and bound or unbound to crowder molecules. Four or five major free-energy minima can be identified. To estimate the dominant MC relaxation times of the peptide, we build MSMs using a range of different time resolutions or lag times. We show that stable relaxation-time estimates can be obtained from the MSM eigenfunctions through fits to autocorrelation data. The eigenfunctions remain sufficiently accurate to permit stable relaxation-time estimation down to small lag times, at which point simple estimates based on the corresponding eigenvalues have large systematic uncertainties. The presence of the crowders has a stabilizing effect on the peptide, especially with BPTI crowders, which can be attributed to a reduced unfolding rate ku, while the folding rate kf is left largely unchanged.

Peptide folding in the presence of interacting protein crowders.

Using Monte Carlo methods, we explore and compare the effects of two protein crowders, BPTI and GB1, on the folding thermodynamics of two peptides, the compact helical trp-cage and the ß-hairpin-forming GB1m3. The thermally highly stable crowder proteins are modeled using a fixed backbone and rotatable side-chains, whereas the peptides are free to fold and unfold. In the simulations, the crowder proteins tend to distort the trp-cage fold, while having a stabilizing effect on GB1m3. The extent of the effects on a given peptide depends on the crowder type. Due to a sticky patch on its surface, BPTI causes larger changes than GB1 in the melting properties of the peptides. The observed effects on the peptides stem largely from attractive and specific interactions with the crowder surfaces, and differ from those seen in reference simulations with purely steric crowder particles.

Equilibrium simulation of trp-cage in the presence of protein crowders.

While steric crowders tend to stabilize globular proteins, it has been found that protein crowders can have an either stabilizing or destabilizing effect, where a destabilization may arise from nonspecific attractive interactions between the test protein and the crowders. Here, we use Monte Carlo replica-exchange methods to explore the equilibrium behavior of the miniprotein trp-cage in the presence of protein crowders. Our results suggest that the surrounding crowders prevent trp-cage from adopting its global native fold, while giving rise to a stabilization of its main secondary-structure element, an α-helix. With the crowding agent used (bovine pancreatic trypsin inhibitor), the trp-cage-crowder interactions are found to be specific, involving a few key residues, most of which are prolines. The effects of these crowders are contrasted with those of hard-sphere crowders.

Thermodynamics of amyloid formation and the role of intersheet interactions.

The self-assembly of proteins into ß-sheet-rich amyloid fibrils has been observed to occur with sigmoidal kinetics, indicating that the system initially is trapped in a metastable state. Here, we use a minimal lattice-based model to explore the thermodynamic forces driving amyloid formation in a finite canonical (NVT) system. By means of generalized-ensemble Monte Carlo techniques and a semi-analytical method, the thermodynamic properties of this model are investigated for different sets of intersheet interaction parameters. When the interactions support lateral growth into multi-layered fibrillar structures, an evaporation/condensation transition is observed, between a supersaturated solution state and a thermodynamically distinct state where small and large fibril-like species exist in equilibrium. Intermediate-size aggregates are statistically suppressed. These properties do not hold if aggregate growth is one-dimensional.

Conformational and aggregation properties of the 1-93 fragment of apolipoprotein A-I.

Several disease-linked mutations of apolipoprotein A-I, the major protein in high-density lipoprotein (HDL), are known to be amyloidogenic, and the fibrils often contain N-terminal fragments of the protein. Here, we present a combined computational and experimental study of the fibril-associated disordered 1-93 fragment of this protein, in wild-type and mutated (G26R, S36A, K40L, W50R) forms. In atomic-level Monte Carlo simulations of the free monomer, validated by circular dichroism spectroscopy, we observe changes in the position-dependent ß-strand probability induced by mutations. We find that these conformational shifts match well with the effects of these mutations in thioflavin T fluorescence and transmission electron microscopy experiments. Together, our results point to molecular mechanisms that may have a key role in disease-linked aggregation of apolipoprotein A-I.

Hybrid Monte Carlo with non-uniform step size.

The Hybrid Monte Carlo method offers a rigorous and potentially efficient approach to the simulation of dense systems, by combining numerical integration of Newton's equations of motion with a Metropolis accept-or-reject step. The Metropolis step corrects for sampling errors caused by the discretization of the equations of motion. The integration is usually performed using a uniform step size. Here, we present simulations of the Lennard-Jones system showing that the use of smaller time steps in the tails of each integration trajectory can reduce errors in energy. The acceptance rate is 10-15 percentage points higher in these runs, compared to simulations with the same trajectory length and the same number of integration steps but a uniform step size. We observe similar effects for the harmonic oscillator and a coarse-grained peptide model, indicating generality of the approach.

Robust Estimation of Diffusion-Optimized Ensembles for Enhanced Sampling.

The multicanonical, or flat-histogram, method is a common technique to improve the sampling efficiency of molecular simulations. The idea is that free-energy barriers in a simulation can be removed by simulating from a distribution where all values of a reaction coordinate are equally likely, and subsequently reweight the obtained statistics to recover the Boltzmann distribution at the temperature of interest. While this method has been successful in practice, the choice of a flat distribution is not necessarily optimal. Recently, it was proposed that additional performance gains could be obtained by taking the position-dependent diffusion coefficient into account, thus placing greater emphasis on regions diffusing slowly. Although some promising examples of applications of this approach exist, the practical usefulness of the method has been hindered by the difficulty in obtaining sufficiently accurate estimates of the diffusion coefficient. Here, we present a simple, yet robust solution to this problem. Compared to current state-of-the-art procedures, the new estimation method requires an order of magnitude fewer data to obtain reliable estimates, thus broadening the potential scope in which this technique can be applied in practice.

Local unfolding and aggregation mechanisms of SOD1: a Monte Carlo exploration.

Copper, zinc superoxide dismutase 1 (SOD1) is a ubiquitous homodimeric enzyme, whose misfolding and aggregation play a potentially key role in the neurodegenerative disease amyotrophic lateral sclerosis (ALS). SOD1 aggregation is thought to be preceded by dimer dissociation and metal loss, but the mechanisms by which the metal-free monomer aggregates remain incompletely understood. Here we use implicit solvent all-atom Monte Carlo (MC) methods to investigate the local unfolding dynamics of the ß-barrel-forming SOD1 monomer. Although event-to-event variations are large, on average, we find clear differences in dynamics among the eight strands forming the ß-barrel. Most dynamic is the eighth strand, ß8, which is located in the dimer interface of native SOD1. For the four strands in or near the dimer interface (ß1, ß2, ß7, and ß8), we perform aggregation simulations to assess the propensity of these chain segments to self-associate. We find that ß1 and ß2 readily self-associate to form intermolecular parallel ß-sheets, whereas ß8 shows a very low aggregation propensity.

Mechanical resistance in unstructured proteins.

Single-molecule pulling experiments on unstructured proteins linked to neurodegenerative diseases have measured rupture forces comparable to those for stable folded proteins. To investigate the structural mechanisms of this unexpected force resistance, we perform pulling simulations of the amyloid ß-peptide (Aß) and α-synuclein (αS), starting from simulated conformational ensembles for the free monomers. For both proteins, the simulations yield a set of rupture events that agree well with the experimental data. By analyzing the conformations occurring shortly before rupture in each event, we find that the mechanically resistant structures share a common architecture, with similarities to the folds adopted by Aß and αS in amyloid fibrils. The disease-linked Arctic mutation of Aß is found to increase the occurrence of highly force-resistant structures. Our study suggests that the high rupture forces observed in Aß and αS pulling experiments are caused by structures that might have a key role in amyloid formation.

Aggregate geometry in amyloid fibril nucleation.

We present and study a minimal structure-based model for the self-assembly of peptides into ordered ß-sheet-rich fibrils. The peptides are represented by unit-length sticks on a cubic lattice and interact by hydrogen bonding and hydrophobicity forces. Using Monte Carlo simulations with >10(5) peptides, we show that fibril formation occurs with sigmoidal kinetics in the model. To determine the mechanism of fibril nucleation, we compute the joint distribution in length and width of the aggregates at equilibrium, using an efficient cluster move and flat-histogram techniques. This analysis, based on simulations with 256 peptides in which aggregates form and dissolve reversibly, shows that the main free-energy barriers that a nascent fibril has to overcome are associated with changes in width.

Distinct phases of free α-synuclein--a Monte Carlo study.

The α-synuclein protein (αS), implicated in Parkinson's disease, shows conformational versatility. It aggregates into ß-sheet-rich fibrils, occurs in helical membrane-bound forms, is disordered as a free monomer, and has recently been suggested to have a folded helical tetramer as its main physiological form. Here, we use implicit solvent all-atom Monte Carlo methods to explore the conformational ensemble sampled by the free αS monomer. We analyze secondary structure propensities, size, and topological properties and compare with existing experimental data. Our study suggests that free αS has two distinct phases. One phase has the expected disordered character. The other phase also shows large conformational variability. However, in this phase, the ß-strand content is substantial, and the backbone fold shows statistical similarities with that in αS fibrils. Presence of this phase is consistent with data from low-temperature experiments. Conversion of disordered αS to this fibril-like form requires the crossing of a rather large apparent free-energy barrier.

Mutation-induced fold switching among lattice proteins.

Recent experiments uncovered a mutational pathway between two proteins, along which a single mutation causes a switch in fold. Searching for such paths between real proteins remains, despite this achievement, a true challenge. Here, we analyze fold switching in the minimalistic hydrophobic/polar model on a square lattice. For this analysis, we generate a comprehensive sequence-structure database for chains of length ≤ 30, which exceeds previous work by five units. Single-mutation-induced fold switching turns out to be quite common in the model. The switches define a fold network, whose topology is roughly similar to what one would expect for a set of randomly connected nodes. In the combinatorially challenging search for fold switches between two proteins, a tempting strategy is to only consider paths containing the minimum number of mutations. Such a restricted search fails to correctly identify 40% of the single-mutation-linked fold pairs that we observe. The thermodynamic stability is correlated with mutational stability and is, on average, markedly reduced at the observed fold switches.

Accelerating atomic-level protein simulations by flat-histogram techniques.

Flat-histogram techniques provide a powerful approach to the simulation of first-order-like phase transitions and are potentially very useful for protein studies. Here, we test this approach by implicit solvent all-atom Monte Carlo (MC) simulations of peptide aggregation, for a 7-residue fragment (GIINFEQ) [corrected] of the Cu/Zn superoxide dismutase 1 protein (SOD1). In simulations with 8 chains, we observe two distinct aggregated/non-aggregated phases. At the midpoint temperature, these phases coexist, separated by a free-energy barrier of height 2.7 k(B)T. We show that this system can be successfully studied by carefully implemented flat-histogram techniques. The frequency of barrier crossing, which is low in conventional canonical simulations, can be increased by turning to a two-step procedure based on the Wang-Landau and multicanonical algorithms.