1 Introduction
Biomolecular condensates may form via phase separation into coexisting solvent-rich and macromolecule-rich phases. Phase separation (PS) is driven by multiple, often transient, interactions which are in many cases engendered by intrinsically disordered proteins (IDPs) and low-complexity domains (LCDs) of multi-domain proteins1–5. The amino acid sequence dictates both the propensity of IDPs to phase separate and the viscoelastic properties of the condensates. Moreover, condensates of some IDPs reconstituted in vitro tend to undergo a transition to a dynamically arrested state, in which oligomeric species can nucleate and ultimately aggregate into fibrils2,3,6–14. As accumulating evidence suggests that these processes may be involved in neurodegeneration and cancer15–18, understanding how the amino acid sequence governs PS and rheological properties of condensates is a current research focus. Due to the transient nature of the protein-protein interactions that underpin PS, quantitative characterization of biomolecular condensates via biophysical experimental methods is challenging, and hence molecular simulations have played an important role in aiding the interpretation of experimental data on condensates reconstituted in vitro19. However, molecular simulations of the PS of IDPs require a minimal system size of ∼100 chains and long simulation times to sample the equilibrium properties of the two-phase system. To enhance the computational efficiency of these simulations, the atomistic representation of the phase-separating protein is typically coarse-grained to fewer interaction sites while the solvent is modelled as a continuum.
A widely used class of coarse-grained models of IDPs describes each residue as a single site centered at the Cα atom. Charged residues interact via salt-screened electrostatic interactions whereas the remaining nonionic nonbonded interactions are incorporated in a single short-range potential characterized by a set of “stickiness” parameters. The stickiness parameters are specific to either the single amino acid or pairs of residues and were originally derived by Dignon et al. from a hydrophobicity scale20. Other models, based on the lattice simulation engines LaSSI21 and PIMMS22, classify the amino acids into a reduced number of residue types with distinct stickiness, ranging from binary categorizations into stickers and spacers4,23,24 to more detailed descriptions and parameterizations25,26. Recently, the accuracy of the stickiness parameters has been considerably improved. This has been achieved by leveraging (i) experimental data on single-chain properties, (ii) statistical analyses of protein structures, and (iii) residue-residue free energy profiles calculated from all-atom simulations26–31. Notably, automated optimization procedures to improve of the stickiness parameters have been proposed by us and others26,28,29,32. In our previous work32, the procedure maximized the accuracy of the model with respect to experimental data reporting on conformational properties of IDPs, namely, small-angle X-ray scattering (SAXS) and paramagnetic relaxation enhancement (PRE) NMR data. To ensure the transferability of the model across sequence space, we trained the model on a large experimental data set and employed a Bayesian regularization approach32,33. As the regularization term, we defined the prior knowledge on the stickiness parameters in terms of 87 hydrophobicity scales reported in the literature. The resulting optimal parameters (previously referred to as M132) were shown to capture the relative propensities to phase separate of a wide range of IDP sequences. However, we also observed that a systematic increase in simulation temperature of about 30 °C is needed to quantitatively reproduce the experimental concentration of the dilute phase coexisting with the condensate on an absolute scale. Herein, we refer to this model as the first version of the CALVADOS (Coarse-graining Approach to Liquid-liquid phase separation Via an Automated Data-driven Optimisation Scheme) model (CALVADOS 1).
In this class of coarse-grained models of IDPs, nonionic interactions are modelled via a Lennard-Jones-like potential, which decays to zero only at infinite residue-residue distances. For computational efficiency, the potential is typically calculated up to a cutoff distance, rc, and interactions between particles that are farther apart are ignored. Although this truncation may introduce severe artefacts, in the different implementations of the models, the value of rc has varied considerably between 1 and 4 nm20,29–32,34,35. Here, we systematically investigate the effect of the cutoff of nonionic interactions on single-chain compaction and PS propensity. We find that decreasing the cutoff from 4 to 2 nm results in a small increase in the radius of gyration whereas the PS propensity significantly decreases. We exploit this effect to improve the temperature-dependence of the CALVADOS 1 model by tuning the cutoff of the nonionic potential. Further, we perform a Bayesian optimization of the stickiness parameters using a cutoff of 2.4 nm and an augmented training set. We show that the updated model (CALVADOS 2) has improved predictive accuracy.
2 Methods
2.1 Molecular simulations
Molecular dynamics simulations are conducted in the NVT ensemble using the Langevin integrator with a time step of 10 fs and friction coefficient of 0.01 ps−1. Non-bonded interactions between residues separated by one bond are excluded from the energy calculations. Functional forms and parameters for bonded and nonbonded interactions are reported in the “Bonded and nonbonded interactions” subsection. Single chains of N residues are simulated using HOOMD-blue v2.9.336 in a cubic box of side length 0.38 × (N − 1) + 4 nm under periodic boundary conditions, starting from the fully extended linear conformation. Conformations are saved every ∆t ≈ 3 × N2 fs if N > 100 and ∆t = 30 ps otherwise. Each chain is simulated in ten replicas for a simulation time of 600 × ∆t. The initial 100 frames of each replica are discarded, so as to obtain 5,000 weakly correlated conformations for each protein (Figure 1). The functional form of ∆t was inferred from calculations of the autocorrelation function of the radius of gyration, Rg, for proteins of various N. Direct-coexistence simulations are performed using openMM v7.537 in a cuboidal box of side lengths [Lx,Ly,Lz] = [25, 25, 300], [17, 17, 300] and [15, 15, 150] nm for Tau 2N4R, Ddx4 LCD, and for the remainder of the proteins, respectively. In the starting configuration, 100 chains are aligned along the z-axis and with their middle beads placed in the xy-plane at random (x, y) positions which are more than 0.7 nm apart. Multi-chain simulations are carried out for at least 2 µs, saving frames every 0.5 ns (Figure S1, S2, and S3). After discarding the initial 0.6 µs, the slab is centered in the box at each frame as previously described32 and the equilibrium density profile, ρ(z), is calculated by averaging over the trajectory of the system at equilibrium. The densities of the dilute and protein-rich phases are estimated as the average densities in the regions |z| < zDS − t/2 and |z| > zDS + 6t nm, where zDS and t are the position of the dividing surface and the thickness of the interface, respectively. zDS and t are obtained by fitting the semi-profiles in z > 0 and z < 0 to ρ(z) = (ρa + ρb)/2 + (ρb − ρa)/2 × tanh [(|z| − zDS )/t], where ρa and ρb are the densities of the protein-rich and dilute phases, respectively. The uncertainty of the density values is estimated as the standard error obtained from the blocking approach38 implemented in the BLOCKING software (github.com/fpesceKU/BLOCKING).
Figure 1. Values of the autocorrelation function (ACF) of the Rg for a lag time of one, two and three frames as a function of sequence length, N.
The autocorrelation is calculated for the sequences of Table 1 and Table 2 simulated using (A) CALVADOS 1 and (B) CALVADOS 2 for ∼ 6 × 0.3 × N2 ps if N > 100 and for 18 ns otherwise. 600 simulation frames are saved every ∼ 0.003 × N2 ps if N > 100 and every 30 ps otherwise. The initial 100 frames are discarded from each simulation.
2.2 Bonded and nonbonded interactions
In this study, we used the following truncated and shifted Ashbaugh-Hatch potential39,
where ϵ = 0.8368 kJ mol−1, rc = 2 or 4 nm, and uLJ is the Lennard-Jones potential:
σ and λ are arithmetic averages of amino acid specific parameters quantifying size and hydropathy, respectively. For σ, we use the values calculated from van der Waals volumes by Kim and Hummer40 whereas, for λ, we use the recently proposed M1 parameters32 and the values optimized in this work.
Salt-screened electrostatic interactions are modeled via the Debye-Hückel potential,
where q is the average amino acid charge number, e is the elementary charge, is the Debye length of an electrolyte solution of ionic strength cs and B(ϵr) is the Bjerrum length. Electrostatic interactions are truncated and shifted at the cutoff distance rc = 4 nm, irrespective of the value of rc used for Eq. 1. We use the following empirical relationship41
to model the temperature-dependent dielectric constant of the implicit aqueous solution. As previously observed31, accounting for the temperature-dependence of ϵr has a small effect on the predictions of the model. Indeed, the relative change in D upon an increase in temperature from 4 to 50 °C is only −3%. Similarly, at cs = 150 mM, the Debye-Hückel energy between like-charged residues at the cutoff distance, uDH(r = 4 nm), is 2.6 J mol−1 at 4 °C and 2.8 J mol−1 at 50 °C. The Henderson–Hasselbalch equation is used to estimate the average charge of the histidine residues, assuming a pKa value of 642.
The amino acid beads are connected by harmonic potentials,
of force constant k = 8033 kJ mol–1 nm–2 and equilibrium distance r0 = 0.38 nm.
2.3 Optimization of the stickiness scale
The optimization of the stickiness parameters, λ, is carried out to minimize the cost function using an algorithm which is analogous to the one we previously described32. and quantify the discrepancy between model predictions and experimental data, and are defined as and respectively, where σexp is the error on the experimental values, Y is either PRE rates or intensity ratios and Nlabels is the number of spin-labeled mutants used for the NMR PRE data. In the expression for the cost function, the coefficients are set to η = 0.1 and θ = 0.05. The prior is the distribution of λ, P(λ), derived from a subset of the hydrophobicity scales reported in Table 3 and 4 of Simm et al.43. Specifically, only the 70 scales that are unique after min-max normalization (Figure S4) are used for the calculation of P(λ), namely Wimley, BULDG reverse, MANP780101, VHEG790101, JANIN, JANJ790102, WOLR790101, PONP800101–6, Wilson, FAUCH, ENGEL, ROSEM, JACWH, CowanWhittacker, ROSM880101 reverse, ROSM880102 reverse, COWR900101, BLAS910101, CASSI, CIDH920101, CIDH920105, CIDBB, CIDA+, CIDAB, PONG1–3, WILM950101–2, WILM950104, Bishop reverse, NADH010101–7, ZIMJ680101, NOZY710101, JONES, LEVIT, KYTJ820101, SWER830101, SWEET, EISEN, ROSEF, GUYFE, COHEN, NNEIG, MDK0, MDK1, JURD980101, SET1–3, CHOTA, CHOTH, Sweet & Eisenberg, KIDER, ROSEB, Welling reverse, Rao & Argos, GIBRA, and WOLR810101 reverse. P(λ) is obtained via multivariate kernel density estimation, as implemented in scikit-learn44, using a Gaussian kernel with bandwidth of 0.05. This prior is 20-dimensional and contains information on the λ-distribution of the single amino acid as well as on the covariance matrix inferred from our selection of 70 hydrophobicity scales (Figure 2).
Figure 2.
(A–T) Probability distributions of the stickiness parameters, P(λ), obtained from 70 min-max normalized hydrophobicity scales selected from the set by Simm et al.43. Blue bars are histograms with bin width of 0.1. Black lines are obtained as 1D projections of the multivariate kernel density estimation implemented in scikit-learn44, using a Gaussian kernel with bandwidth of 0.05. (U) Covariance matrix of the 70 min-max normalized hydrophobicity scales selected from the set by Simm et al.43. The upper triangle of the matrix shows the covariance calculated directly from the 70 min-max normalized hydrophobicity scales whereas the lower triangle of the matrix shows the covariance calculated from the multivariate kernel density estimation averaging over the 70 min-max normalized hydrophobicity scales.
In the first step of the optimization procedure, the λ values for all the amino acids are set to 0.5, λ0 = 0.5, and these parameters are used to simulate the proteins of the training set (Table 1 and Table 3). We proceed with the first optimization cycle, wherein, at each k-th iteration, the λ values of a random selection of five amino acids are nudged by random numbers picked from a normal distribution of standard deviation 0.05 to generate a trial λk set. For each ith frame, we calculate the Boltzmann weight as wi = exp {−[U(ri,λk) − U(ri,λ0)]/kBT}, where U is the nonionic potential. The trial λk is discarded if the effective fraction of frames, , is lower than 60%. Otherwise, the acceptance probability follows the Metropolis criterion, where ξk is a unitless control parameter. Each optimization cycle is divided into ten micro-cycles, wherein the control parameter, ξ, is initially set to ξ0 = 0.1 and scaled down by 1% at each iteration until ξ < 10−8. From the complete optimization cycle, we select the λ set yielding the lowest estimate of ℒ. Consecutive optimization cycles are performed from simulations runs carried out with the intermediate optimal λ set. To show that the procedure is reproducible and that the final λ set is relatively independent of the initial conditions, we performed an additional optimization procedure starting from the M1 model, λ0 =M132. The optimization performed in this work differs from our previous implementation32 also for the following details: (i) nine additional sequences have been included in the training set (Table 1 and Table 3); (ii) single chains are simulated as detailed in the “Molecular simulations” Subsection; (v) the average radius of gyration is calculated as 〈Rg〉 instead of .
Table 1. Solution conditions and experimental radii of gyration of proteins included in the training set for the Bayesian parameter-learning procedure.
| Protein | N | Rg (nm) | T (K) | cs (M) | pH | Ref. |
|---|
| Hst5 | 24 | 1.38 ± 0.05 | 293 | 0.15 | 7.5 | 45 |
| (Hst5)2 | 48 | 1.87 ± 0.05 | 298 | 0.15 | 7.0 | 46 |
| p53 (20–70) | 62 | 2.39 ± 0.05 | 277 | 0.1 | 7.0 | 47 |
| ACTR | 71 | 2.6 ± 0.1 | 278 | 0.2 | 7.4 | 48 |
| Ash1 | 81 | 2.9 ± 0.05 | 293 | 0.15 | 7.5 | 49,50 |
| CTD2 | 83 | 2.61 ± 0.05 | 293 | 0.12 | 7.5 | 50,51 |
| Sic1 | 92 | 3.0 ± 0.4 | 293 | 0.2 | 7.5 | 52 |
| SH4UD | 95 | 2.7 ± 0.1 | 293 | 0.2 | 8.0 | 53 |
| ColNT | 98 | 2.8 ± 0.1 | 277 | 0.4 | 7.6 | 54 |
| p15PAF | 111 | 2.8 ± 0.1 | 298 | 0.15 | 7.0 | 55 |
| hNL3cyt | 119 | 3.2 ± 0.2 | 293 | 0.3 | 8.5 | 56 |
| RNaseA | 124 | 3.4 ± 0.1 | 298 | 0.15 | 7.5 | 57 |
| A1 | 137 | 2.76 ± 0.02 | 298 | 0.15 | 7.0 | 58 |
| -10R | 137 | 2.67 ± 0.01 | 298 | 0.15 | 7.0 | 58 |
| -6R | 137 | 2.57 ± 0.01 | 298 | 0.15 | 7.0 | 58 |
| +2R | 137 | 2.62 ± 0.02 | 298 | 0.15 | 7.0 | 58 |
| +7R | 137 | 2.71 ± 0.01 | 298 | 0.15 | 7.0 | 58 |
| -3R+3K | 137 | 2.63 ± 0.02 | 298 | 0.15 | 7.0 | 58 |
| -6R+6K | 137 | 2.79 ± 0.01 | 298 | 0.15 | 7.0 | 58 |
| -10R+10K | 137 | 2.85 ± 0.01 | 298 | 0.15 | 7.0 | 58 |
| +12D | 137 | 2.80 ± 0.01 | 298 | 0.15 | 7.0 | 58 |
| +4D | 137 | 2.72 ± 0.03 | 298 | 0.15 | 7.0 | 58 |
| +8D | 137 | 2.69 ± 0.01 | 298 | 0.15 | 7.0 | 58 |
| -9F+3Y | 137 | 2.68 ± 0.01 | 298 | 0.15 | 7.0 | 58 |
| +12E | 137 | 2.85 ± 0.01 | 298 | 0.15 | 7.0 | 58 |
| +7K+12D | 137 | 2.92 ± 0.01 | 298 | 0.15 | 7.0 | 58 |
| +7K+12D blocky | 137 | 2.56 ± 0.01 | 298 | 0.15 | 7.0 | 58 |
| -4D | 137 | 2.64 ± 0.01 | 298 | 0.15 | 7.0 | 58 |
| -8F+4Y | 137 | 2.71 ± 0.01 | 298 | 0.15 | 7.0 | 58 |
| -10F+7R+12D | 137 | 2.86 ± 0.01 | 298 | 0.15 | 7.0 | 58 |
| +7F-7Y | 137 | 2.72 ± 0.01 | 298 | 0.15 | 7.0 | 58 |
| -12F+12Y | 137 | 2.60 ± 0.02 | 298 | 0.15 | 7.0 | 58 |
| -12F+12Y-10R | 137 | 2.61 ± 0.02 | 298 | 0.15 | 7.0 | 58 |
| -9F+6Y | 137 | 2.65 ± 0.01 | 298 | 0.15 | 7.0 | 58 |
| αSyn | 140 | 3.55 ± 0.1 | 293 | 0.2 | 7.4 | 59 |
| FhuA | 144 | 3.34 ± 0.1 | 298 | 0.15 | 7.5 | 57 |
| K27 | 167 | 3.70 ± 0.2 | 288 | 0.15 | 7.4 | 60 |
| K10 | 168 | 4.00 ± 0.1 | 288 | 0.15 | 7.4 | 60 |
| K25 | 185 | 4.10 ± 0.2 | 288 | 0.15 | 7.4 | 60 |
| K32 | 198 | 4.20 ± 0.3 | 288 | 0.15 | 7.4 | 60 |
| CAHSD | 227 | 4.8 ± 0.2 | 293 | 0.07 | 7.0 | 61 |
| K23 | 254 | 4.9 ± 0.2 | 288 | 0.15 | 7.4 | 60 |
| Tau35 | 255 | 4.7 ± 0.1 | 298 | 0.15 | 7.4 | 62 |
| CoRNID | 271 | 4.7 ± 0.2 | 293 | 0.2 | 7.5 | 63 |
| K44 | 283 | 5.2 ± 0.2 | 288 | 0.15 | 7.4 | 60 |
| PNt | 334 | 5.1 ± 0.1 | 298 | 0.15 | 7.5 | 57,64 |
| PNt Swap1 | 334 | 4.9 ± 0.1 | 298 | 0.15 | 7.5 | 64 |
| PNt Swap4 | 334 | 5.3 ± 0.1 | 298 | 0.15 | 7.5 | 64 |
| PNt Swap5 | 334 | 4.9 ± 0.1 | 298 | 0.15 | 7.5 | 64 |
| PNt Swap6 | 334 | 5.3 ± 0.1 | 298 | 0.15 | 7.5 | 64 |
| GHRICD | 351 | 6.0 ± 0.5 | 298 | 0.35 | 7.3 | 65,66 |
3 Results and discussion
When applying a cutoff scheme, we neglect the interactions of residues separated by a distance, r, larger than the cutoff, rc. For the most strongly interacting residue pair (between two tryptophans), the nonionic potential of the CALVADOS 1 model at rc = 2 nm takes the value of -5 J mol–1, that is, only a small fraction of the thermal energy (Figure 3A). However, the Lennard-Jones potential falls off slowly whereas the number of interacting partners increases quadratically with increasing r. Therefore, in a simulation of a protein-rich phase, decreasing the cutoff from 4 to 2 nm (Figure 3A) can imply ignoring a total interaction energy per protein of several times the thermal energy. We first look into the effect of the choice of cutoff on the conformational ensembles of isolated proteins. We simulated single IDPs of different sequence length, N = 71–441, and average hydropathy, 〈λ〉 = 0.33–0.63. The average radii of gyration, 〈Rg〉, calculated from simulation trajectories are systematically larger when we use rc = 2 nm, compared to the values obtained using rc = 4 nm. CALVADOS 1 was optimized using the longer rc and estimating the ensemble average Rg values as the root-mean-square Rg, . Since is systematically larger than 〈Rg〉, decreasing rc to 2 nm results in a slight improvement of the agreement between the calculated 〈Rg〉 values and the experimental data (Figure 3B).
Figure 3. Effect of cutoff size on predictions of radii of gyration, Rg, and saturation concentration, csat, from simulations performed using the CALVADOS 1 parameters.
(A) Nonionic Ashbaugh-Hatch potentials between two W residues with cutoff, rc, of 4 (blue solid line) and 2 nm (orange dashed line). The inset highlights differences between the potentials for rs ≤ r ≤ rc. (B) Relative difference between experimental and predicted radii of gyration, 〈Rg〉, for rc = 4 (blue) and 2 nm (orange). values reported in the legend are calculated for all the sequences in Table 1. Error bars represent the experimental error relative to . (C) 〈Rg〉 of α-Synuclein, hnRNPA1 LCD, PNt and human full-length tau (Table 1 and Table 2) from simulations performed with increasing cutoff size, rc, and normalized by the value at rc = 2 nm. (D) Saturation concentration, csat, for hnRNPA1 LCD, the randomly shuffled sequence of LAF-1 RGG domain, LAF-1 RGG domain and human full-length tau for increasing values of rc and normalized by the csat at rc = 2 nm. (E–G) Correlation between csat from simulations and experiments for (E) A1 LCD variants, (F) A1 LCD∗ WT at [NaCl] = 0.15, 0.2, 0.3 and 0.5 M and (G) variants of LAF-1 RGG domain (Table 4).
To gain further insight into the effect of the cutoff, we performed simulations of single chains of α-Synuclein, hnRNPA1 LCD, PNt and Tau 2N4R (Table 1 and Table 2) using rc = 2, 2.5, 3 and 4 nm. Irrespective of the sequence, 〈Rg〉 decreases monotonically with increasing rc. However, the effect on compaction appears to be larger for long sequences and high content of hydrophobic residues, both of which result in an increased number of shorter intramolecular distances. For example, upon increasing the rc from 2 to 4 nm, the 〈Rg〉 of α-Synuclein (N = 140, 〈λ〉 = 0.33) decreases by 2.3% whereas the effect is more pronounced for hnRNPA1 LCD (N = 137, 〈λ〉 = 0.61) and Tau 2N4R (N = 441, 〈λ〉 = 0.38), with a decrease in 〈Rg〉 of 4.0% and 7.7%, respectively.
Table 2. Solution conditions and experimental radii of gyration of proteins simulated in this study but not included in the training set for the Bayesian parameter-learning procedure.
| Protein | N | Rg (nm) | T (K) | cs (M) | pH | Ref. |
|---|
| DSS1 | 71 | 2.5 ± 0.1 | 288 | 0.17 | 7.4 | 66 |
| p27Cv14 | 107 | 2.936 ± 0.13 | 293 | 0.095 | 7.2 | 67 |
| p27Cv15 | 107 | 2.915 ± 0.10 | 293 | 0.095 | 7.2 | 67 |
| p27Cv31 | 107 | 2.81 ± 0.18 | 293 | 0.095 | 7.2 | 67 |
| p27Cv44 | 107 | 2.492 ± 0.13 | 293 | 0.095 | 7.2 | 67 |
| p27Cv56 | 107 | 2.328 ± 0.10 | 293 | 0.095 | 7.2 | 67 |
| p27Cv78 | 107 | 2.211 ± 0.03 | 293 | 0.095 | 7.2 | 67 |
| PTMA | 111 | 3.7 ± 0.2 | 288 | 0.16 | 7.4 | 66 |
| NHE6cmdd | 116 | 3.2 ± 0.2 | 288 | 0.17 | 7.4 | 66 |
| A1 LCD∗ | 131 | 2.645 ± 0.02 | 293 | 0.05 | 7.5 | 68 |
| A1 LCD∗ | 131 | 2.65 ± 0.02 | 293 | 0.15 | 7.5 | 68 |
| A1 LCD∗ | 131 | 2.62 ± 0.02 | 293 | 0.3 | 7.5 | 68 |
| A1 LCD∗ | 131 | 2.528 ± 0.02 | 293 | 0.5 | 7.5 | 68 |
| ANAC046 | 167 | 3.6 ± 0.3 | 298 | 0.14 | 7.0 | 66 |
| Tau 2N3R | 410 | 6.3 ± 0.3 | 298 | 0.15 | 7.4 | 62 |
| Tau 2N4R | 441 | 6.7 ± 0.3 | 298 | 0.15 | 7.4 | 62 |
Table 3. Protein and conditions related to the intramolecular PRE data included in the training set.
| Protein | N | Nlabels | ωI /2π (MHz) | T (K) | cs (M) | pH | Ref. |
|---|
| FUS | 163 | 3 | 850 | 298 | 0.15 | 5.5 | 2 |
| FUS12E | 164 | 3 | 850 | 298 | 0.15 | 5.5 | 2 |
| OPN | 220 | 10 | 800 | 298 | 0.15 | 6.5 | 69 |
| αSyn | 140 | 5 | 700 | 283 | 0.2 | 7.4 | 70 |
| A2 | 155 | 2 | 850 | 298 | 0.005 | 5.5 | 3 |
To investigate the effect of the cutoff distance on PS propensity, we performed direct-coexistence simulations of 100 chains of hnRNPA1 LCD, LAF-1 RGG domain (WT and shuffled sequence with higher charge segregation), and Tau 2N4R (Table 4). From the simulation trajectories of the two-phase system at equilibrium, we calculate csat, i.e. the protein concentration in the dilute phase coexisting with the condensate. The higher the csat value, the lower the propensity of the IDP to undergo PS. As expected from the increased contact density in the condensed phase, the choice of cutoff has a considerably larger impact on csat than on chain compaction: decreasing rc from 4 to 2 nm results in an increase in csat of over one order of magnitude. In contrast to what we observed for the 〈Rg〉, the decrease in csat does not show a clear dependence on sequence length and average hydropathy.
Table 4. Proteins and conditions used for the direct-coexistence simulations performed in this study and references to the experimental data.
Shaded rows highlight systems which are not included in the correlation plot of Figure 7C.
| Protein | N | cs (mM) | pH | Ref. | T (K) |
|---|
| 4 nm | 2 nm | Figure 1C |
|---|
| 6His-TEV-Lge11−80-StrepII WT | 114 | 100 | 7.5 | 71 | - | 293 | - |
| 6His-TEV-Lge11−80-StrepII -11R+11K | 114 | 100 | 7.5 | 71 | - | 293 | - |
| 6His-TEV-Lge11−80-StrepII -14Y+14A | 114 | 100 | 7.5 | 71 | - | 293 | - |
| A1 LCD WT | 137 | 150 | 7.0 | 4,
58 | 310 & 323 | 277 & 293 | 310 |
| A1 LCD +7F-7Y | 137 | 150 | 7.0 | 58 | 310 & 323 | 277 & 293 | - |
| A1 LCD -12F+12Y | 137 | 150 | 7.0 | 58 | 310 & 323 | 277 & 293 | - |
| A1 LCD -23S+23T | 137 | 150 | 7.0 | 58 | 310 & 323 | 277 & 293 | - |
| A1 LCD -14N+14Q | 137 | 150 | 7.0 | 58 | 310 & 323 | 277 & 293 | - |
| A1 LCD -10G+10S | 137 | 150 | 7.0 | 58 | 310 & 323 | 277 & 293 | - |
| A1 LCD -20G+20S | 137 | 150 | 7.0 | 58 | 310 & 323 | 277 & 293 | - |
| A1 LCD -30G+30S | 137 | 150 | 7.0 | 58 | 323 | 293 | - |
| A1 LCD +23G-23S | 137 | 150 | 7.0 | 58 | 323 | 293 | - |
| A1 LCD +23G-23S+7F-7Y | 137 | 150 | 7.0 | 58 | 323 | 293 | - |
| A1 LCD +23G-23S-12F+12Y | 137 | 150 | 7.0 | 58 | 323 | 293 | - |
| A1 LCD -9F+3Y | 137 | 150 | 7.0 | 58 | 310 | 277 | - |
| A1 LCD -8F+4Y | 137 | 150 | 7.0 | 58 | 310 | 277 | - |
| A1 LCD -3R+3K | 137 | 150 | 7.0 | 58 | 310 | 277 | - |
| A1 LCD -6R | 137 | 150 | 7.0 | 58 | 310 | 277 | - |
| A1 LCD -4D | 137 | 150 | 7.0 | 58 | 310 | 277 | - |
| A1 LCD +4D | 137 | 150 | 7.0 | 58 | 310 | 277 | - |
| A1 LCD +8D | 137 | 150 | 7.0 | 58 | 310 | 277 | - |
| A1 LCD +2R | 137 | 150 | 7.0 | 58 | 310 | 277 | - |
| A1 LCD∗ WT | 131 | 150 | 7.0 | 72 | 323 | 293 | - |
| A1 LCD∗ WT | 131 | 200 | 7.0 | 72 | 323 | 293 | - |
| A1 LCD∗ WT | 131 | 300 | 7.0 | 72 | 323 | 293 | - |
| A1 LCD∗ WT | 131 | 500 | 7.0 | 72 | 323 | 293 | - |
| LAF-1 RGG Domain | 176 | 150 | 7.5 | 73 | 323 | 293 | 293 |
| LAF-1 RGG Domain Shuffled | 176 | 150 | 7.5 | 73 | 323 | 293 | 323 |
| LAF-1 RGG Domain ∆21–30 | 166 | 150 | 7.5 | 73 | 323 | 293 | - |
| A2 LCD | 155 | 10 | 5.5 | 74 | - | 297 | - |
| FUS LCD | 163 | 150 | 7.4 | 75 | - | 297 | - |
| Ddx4 LCD | 236 | 130 | 6.5 | 76 | - | 297 | - |
| Human Full-Length Tau (2N4R) | 441 | 70 | 7.4 | - | - | - | 277 |
From the multi-chain trajectories of hnRNPA1 LCD, LAF-1 RGG domain (WT and shuffled sequence) and Tau 2N4R obtained using rc = 4 nm, we estimate that the increase in nonionic energy per protein upon decreasing the cutoff from 4 to 2 nm is 13±1 kJ mol–1 (mean±standard deviation), respectively (Figure 4A). Assuming that the number of interactions neglected by the shorter cutoff is proportional to the sequence length and the amino acid concentration in the condensate, the small variance in the energy increase across the different IDPs finds explanation in the fact that the simulated systems display similar values of N2 × ccon (Figure 4A), where ccon is the protein concentration in the condensate. The ratio U (rc = 2 nm)/U (rc = 4 nm) of the nonionic energies for rc = 2 and 4 nm is also largely system independent (Figure 4B). Moreover, decreasing the temperature by ∼30 K in the range between 310 and 323 K has a rather small impact on the relative strength of the electrostatic interactions with respect to the thermal energy, due to the decrease in the dielectric constant of water with increasing temperature (Figure 4C). Therefore, we speculate that the effect of decreasing rc can be compensated by simulating the system at a lower temperature (Figure 4B).
Figure 4.
(A) Comparison between nonionic energy difference per protein (∆U = U (rc = 2 nm) − U (rc = 4 nm), hatched) and N2 × ccon (orange), where N is the sequence length and ccon is the molar protein concentration in the condensate. (B) Ratio between nonionic energies calculated with rc = 2 and 4 nm (open circles) compared to the ratio of the thermal energy, , at T′ = T − 20 K and at T (orange), where R is the gas constant. (C) Increase in electrostatic energy relative to the thermal energy upon decreasing the temperature by 30 (black) and 20 K (orange). The data shown in this figure are obtained from simulations of hnRNPA1 LCD, LAF-1 RGG domain (WT and shuffled sequence) and Tau 2N4R performed at T = 310, 293, 323, and 277 K, respectively, and using rc = 4 nm. Error bars are standard deviations over trajectories of the systems at equilibrium.
With these considerations in mind, we use the CALVADOS 1 model with rc = 2 nm to run direct-coexistence simulations of IDPs for which csat has been measured experimentally (Table 4), i.e. variants of hnRNPA1 LCD, hnRNPA1 LCD∗ at various salt concentrations, and LAF-1 RGG domain variants. As we have shown that decreasing the range of the nonionic interactions disfavours PS, we perform these simulations at the experimental temperatures, which are lower by ∼30 K than those required to reproduce the experimental csat values when the model is simulated with rc = 4 nm (Figure 3E–G). The two-fold decrease in rc enables the model to quantitatively recapitulate the experimental csat data at the temperature at which the experiments were conducted. Notably, we show this for diverse sequences, across a wide range of ionic strengths, and for variants with different charge patterning and numbers of aromatic and charged residues. These results suggest that the range of interaction of the Lennard-Jones potential may be too large77. While the r−6 dependence is strictly correct for dispersion interactions between atoms, the nonionic potential of our model incorporates a variety of effective nonbonded interactions between residues, and hence the Lennard-Jones potential is not expected to capture the correct interaction range31.
Since CALVADOS 1 was developed using rc = 4 nm, we examined whether reoptimizing the model with the shorter cutoff could result in a comparably accurate model. As detailed in the Methods Section, we performed a Bayesian parameter-learning procedure32 using an improved algorithm, an expanded training set (Table 1), and rc = 2 nm. Figure S5 shows that the new model tends to underestimate the csat values of the most PS-prone sequences. We hypothesize that during the optimization the reduction of attractive forces due to the shorter cutoff is overcompensated by an overall increase in λ. We tested this hypothesis by performing the optimization with increasing values of rc, in the range between 2.0 and 2.5 nm, and found that the csat values predicted from simulations performed with rc = 2.0 nm tend to increase with the rc used for the optimization (Figure 5).
Figure 5. Saturation concentrations, csat, as a function of the cutoff used to optimize the model.
csat values are calculated from simulations performed using rc = 2.0 nm whereas the models are optimized using rc = 2.0, 2.2, 2.4, and 2.5 nm. Horizontal dotted lines represent experimental csat values from the references reported in Table 4.
Using rc = 2.4 nm for the optimization resulted in a model with improved accuracy compared to CALVADOS 1 (Figure 6), especially for the PS of LAF-1 RGG domain and the −23S+23T variant of A1 LCD. To test the robustness of the approach, the optimization was carried out starting from λ0 = 0.5 for all the amino acids (Figure 6) and from λ0 =M1 (Figure S6 and S7). The difference between the resulting sets of optimal λ values (Figure S7A) is lower than 0.08 for all the residues and exceeds 0.05 only for S, T and A. The model obtained starting from λ0 = 0.5 is more accurate at predicting PS propensities and will be referred to as CALVADOS 2 hereafter.
Figure 6.
(A) Comparison between λ sets of CALVADOS 1 (orange) and CALVADOS 2 (blue). (B) Distribution of the relative difference between experimental (Table 1) and predicted radii of gyration, 〈Rg〉, for CALVADOS 1 (orange) and CALVADOS 2 (blue). (C) Comparison between saturation concentrations, csat, at 293 K of variants of hnRNPA1 LCD measured by Bremer, Farag, Borcherds et al.58 (closed circles) and corresponding predictions of CALVADOS 1 (open orange circles) and CALVADOS 2 (open blue squares). (D–F) Correlation between csat from simulations and experiments for (D) A1 LCD variants, (E) A1 LCD* WT at [NaCl] = 0.15, 0.2, 0.3 and 0.5 M and (F) variants of LAF-1 RGG domain (Table 4).
The λ values of CALVADOS 1 and 2 differ mostly for K, T, A, M, and V, whereas the smaller deviations (|∆λ| < 0.09) observed for Q, L, I, and F (Figure 6A) are within the range of reproducibility of the method (Figure S7A). Although CALVADOS 1 was optimized using rc = 4 nm, predictions of single-chain compaction from simulations performed using rc = 2 nm are more accurate for CALVADOS 1 than for CALVADOS 2. This result can be explained by the opposing effects of decreasing the cutoff and calculating Rg values as 〈Rg〉 instead of . In fact, the values predicted by CALVADOS 1 are strikingly similar to the 〈Rg〉 values predicted by CALVADOS 2 (Figure 7A).
The correlation between experimental and predicted Rg values for the 67 proteins of Table 1 and Table 2 is excellent for both CALVADOS 1 and 2 (Figure 7B). On the other hand, CALVADOS 2 is more accurate than CALVADOS 1 at predicting PS propensities, as evidenced by Pearson’s correlation coefficients of 0.93 and 0.82, respectively, for the experimental and predicted csat values of the 26 sequences of Table 4 (Figure 7C).
Figure 7.
(A) Relative difference between experimental and predicted radii of gyration for CALVADOS 1 (orange) and CALVADOS 2 (blue). Full and hatched bars show ( – )/ where is calculated as the mean or the root mean square , respectively. The vertical dashed line splits the plot into the 51 and 16 sequences or solution conditions of the training set (Table 1) and test set (Table 2), respectively. Error bars represent the experimental error relative to . χ2 values in the legend are averages over 67 different sequences or solution conditions (Table 1 and Table 2). (B and C) Comparison between experimental and predicted (B) Rg (Table 1 and Table 2) and (C) csat values for CALVADOS 1 (orange) and CALVADOS 2 (blue). Pearson’s r coefficients are reported in the legend. Small squares in C show the same data as in Figure 6C–F whereas the large upward triangle, downward triangle, and circle show values for A2 LCD, FUS LCD, and Ddx4 LCD, respectively, at the conditions reported in Table 4.
Capturing the interplay between short-range nonionic and long-range ionic interactions is essential for accurately modeling the PS of IDPs58,78–80. Our results show that the decrease in the range of the nonionic potential reported in this work does not significantly perturb the balance between ionic and nonionic forces. In fact, CALVADOS 1 and 2 accurately predict the PS propensities of A1 LCD* at various salt concentrations, as well as the csat of variants of A1 LCD and LAF-1 RGG domain with different charge patterning (Figure 3E–G and Figure 6D–F). Moreover, CALVADOS 1 and 2 recapitulate the effect of salt concentration and charge patterning on the chain compaction of A1 LCD∗68 and p27-C constructs67, respectively (Figure 8A and 8B).
Figure 8.
Comparison between experimental Rg values and predictions of CALVADOS 1 (orange) and CALVADOS 2 (blue) for (A) A1 LCD* at different salt concentrations (50 mM < cs < 500 mM) and (B) p27-C constructs of different charge patterning (0.1 < κ < 0.8). Experimental conditions and references are reported in Table 2. (C) Predictions of CALVADOS 2 direct-coexistence simulations of the PS of constructs of the 1–80 N-terminal fragment of yeast Lge1 simulated at cs = 100 mM. Protein concentration profiles are shown as a function of the long side of the simulation cell for WT (blue), -11R+11K variant (orange), and -14Y+14A variant (green).
In the model, ionic interactions are also truncated and shifted. At the cutoff distance of 4 nm, the ionic energy decreases with increasing salt concentration and amounts to ±2.7 J mol–1 at cs = 150 mM and 20 °C. However, this energy is ~ 43 times larger at cs = 10 mM, which suggests that the model may considerably underestimate the strength and range of charge-charge interactions at low salt concentrations. To investigate this aspect, we performed single-chain and direct-coexistence simulations using a longer cutoff of 6 nm for the ionic interactions. The change in cutoff has a small effect on both the Rg (Figure S8A) and the csat values predicted for systems at cs = 150 mM (Figure S8C). Conversely, simulations at low salt concentration are considerably affected by the increase in cutoff. For the PRE data of A2 LCD at cs = 5 mM, we observe an improvement in the agreement with experiments (Figure S8B). Instead, the accuracy of the phase behaviour predicted for A2 LCD at cs = 10 mM decreases significantly as the csat value shows a ∼100-fold increase (Figure S8C). Since the vast majority of the available Rg and csat data in our training and test sets was measured at cs ≈ 150 mM, we are currently unable to further assess or improve the accuracy of the model at low salt concentrations.
As additional test systems, we considered constructs of the 1–80 N-terminal fragment of yeast Lge1, which have been recently investigated using turbidity measurements71. CALVADOS 2 correctly predicts that the WT Lge11–80 construct undergoes PS at the experimental conditions, albeit with a hundred times larger csat (50 ± 6 µM at cs = 100mM) compared to experiments (< 1 µM). In agreement with experiments, CALVADOS 2 predicts that mutating all the 11 R residues to K increases csat by over one order of magnitude whereas mutating the 14 Y residues of the 1–80 fragment to A abrogates PS (Figure 8C).
4 Conclusions
In the context of a previously developed Cα-based IDP model (CALVADOS), we show that neglecting the long range of attractive Lennard-Jones interactions has a small impact on the compaction of a single chain while strongly disfavouring PS. The effect can be explained by the smaller number of neglected pair interactions for a residues in an isolated chain compared to the dense environment of a condensate. Moreover, we find that the effect of reducing the range of interaction by a factor of two is relatively insensitive to sequence length and composition. Therefore, decreasing the cutoff of the Lennard-Jones potential of the Cα-based model engenders a similar generic effect on chain compaction and PS as a corresponding increase in temperature. We take advantage of this finding to solve the temperature mismatch of the CALVADOS model. Namely, we decrease the cutoff of the nonionic interactions from 4 to 2 nm and obtain accurate csat predictions at the experimental conditions, whereas simulations at temperatures higher by 30 °C were required in the original implementation. Finally, we used the shorter cutoff to reoptimize the stickiness parameters of the model against experimental data reporting on single-chain compaction. The small expansion of the chain conformations is overcompensated by an overall increase in stickiness so that the resulting model tends to underestimate the experimental csat values. By systematically increasing the cutoff used in the development of the stickiness scale, we find that performing the optimization using rc = 2.4 nm results in a model (CALVADOS 2) which yields accurate predictions from simulations run using rc = 2 nm at the experimental conditions. We present CALVADOS 2 as an improvement of our previous model by testing on sets of experimental Rg and csat data comprising 16 and 36 systems, respectively, which were not used in the parameterization of the model.
Ethics and consent
Ethical approval and consent were not required.
