Introduction

Biomolecular condensates formed through Liquid-liquid phase separation (LLPS) are increasingly recognized as fundamental regulators of diverse cellular processes, including stress responses Van Der Lee et al. (2014); Feric et al. (2016); Shin and Brangwynne (2017); Boeynaems et al. (2018); Burke et al. (2015); Biamonti and Vourc’h (2010); Riback et al. (2017), cellular signaling Wippich et al. (2013); Su et al. (2016), gene expression regulation Morimoto and Boerkoel (2013); Hnisz et al. (2017), and chromatin organization Strom et al. (2017). Dysregulation of condensate formation, in contrast, has been implicated in a wide range of human diseases, most notably cancer and numerous neurodegenerative disorders Molliex et al. (2015); Patel et al. (2015); Boeynaems et al. (2018); Conicella et al. (2016); Osterburg and Dötsch (2022). Intrinsically disordered proteins (IDPs), which lack stable tertiary structures and often contain multivalent interaction motifs, are key molecular drivers of LLPS Molliex et al. (2015); Uversky et al. (2015); Zhou et al. (2018); Dignon et al. (2020); Elbaum-Garfinkle et al. (2015). Understanding the molecular principles governing IDP-mediated LLPS is therefore essential for elucidating the mechanisms underlying both physiological condensate function and pathological phase transitions. A central yet incompletely understood aspect of LLPS is the multi-step desolvation of biomolecular chains as they transition from a homogeneously dispersed solution to a condensed phase. How this desolvation energetics contributes to the thermodynamics and kinetics of LLPS remains a major open question.

Because experimentally characterizing the microscopic molecular events of LLPS at sub-molecule level is still challenging for most experimental assays, a broad spectrum of computational approaches has been developed to characterize biomolecular phase separation. These range from mean-field Flory (1942); Huggins (1941) and random-phase-approximation polymer theories Lin et al. (2017); Lin and Chan (2017); McCarty et al. (2019) to diverse molecular simulation frameworks based on atomistic or coarse-grained (CG) force fields Dignon et al. (2018b,a); Davtyan et al. (2012); De Jong et al. (2013); Baul et al. (2019). Although all-atom molecular dynamics (MD) simulations can capture the microscopic processes of biomolecules with unprecedented temporal and spatial resolutions, they become computationally prohibitive when applied to biomolecular phase separation involving large number of protein chains. This challenge has motivated the development of CG models that reduce computational cost while retaining the essential physicochemical features governing phase behavior. These CG frameworks include the sticker-and-spacer paradigm and its lattice Monte-Carlo implementation LASSI Harmon et al. (2017); Ruff et al. (2021); Choi et al. (2019), the hydrophobicity-scale (HPS) family and its refinements Dignon et al. (2018b,a); Regy et al. (2021); Dannenhoffer-Lafage and Best (2021), and Mpipi, which explicitly incorporates cation-π interactions Joseph et al. (2021). Advanced multiscale models have also been pioneered to predict the coupled phase behavior of IDPs and chromatin with remarkable accuracy Espinosa et al. (2020). Higher-resolution or experimentally fine-tuned CG force fields (e.g., AWSEM-IDP Wu et al. (2018), Martini3 Zhang et al. (2023); Souza et al. (2021), COCOMO Valdes-Garcia et al. (2023); Jussupow et al. (2025), AICG2+ Li et al. (2014); Bian et al. (2024), and SIRAH Darré et al. (2015)) further capture hydrogen bonding, secondary-structure propensity, and sequence-dependent packing. Notably, the maximum entropy optimized force field (MOFF) integrates maximum entropy biasing and energy landscape theory, successfully predicting the experimental dimensions and macroscopic phase separation of both ordered and disordered proteins Latham and Zhang (2019, 2021). More recently, Tesei et al. introduced the CALVADOS model, which optimizes residue-specific interaction parameters within the HPS framework using SAXS and paramagnetic relaxation enhancement (PRE) data across a large IDP dataset Tesei et al. (2021); Tesei and Lindorff-Larsen (2022); Tesei et al. (2024), thereby improving the accuracy and transferability of CG models for predicting IDP conformations and LLPS behavior.

Despite these advances, most residue-level CG models rely on implicit solvent representations. This simplification neglects the contribution of water molecules in mediating inter-residue interactions, which can be crucial for the thermodynamics and kinetics of biomolecular folding and phase transitions Cheung et al. (2002); Liu and Chan (2005); Karanicolas and Brooks III (2002); Tarus et al. (2006); Wu et al. (2011b). In particular, the desolvation penalty associated with electrostatic and hydrophobic interactions exhibits marked differences between implicit and explicit solvent representations Salari and Chong (2012); Wu et al. (2011a). More importantly, conventional implicit-solvent CG models cannot capture the multi-step desolvation events that accompany the transition from a dilute solution to a dense condensate in LLPS. Recent theoretical and experimental studies highlight that water reorganization and solvent-mediated entropy changes are not merely passive consequences of LLPS but active driving forces that shape condensate stability, internal structure, and dynamics Mukherjee and Schäfer (2023); Mukherjee et al. (2024); Li and Hou (2023). Nevertheless, implicit-solvent CG models remain widely used because they strike an effective balance between computational efficiency and residue-specific resolution. Incorporating explicit desolvation physics into such models thus represents a particularly promising strategy for enhancing their physical realism without sacrificing tractability in LLPS simulations. Notably, residue-level desolvation terms have already been successfully employed in protein folding studies Cheung et al. (2002); Liu and Chan (2005); Ferguson et al. (2009); Chen and Chan (2014); Zhang and Chan (2010); Karanicolas and Brooks III (2002); Gasic and Cheung (2020); Camilloni et al. (2016). Within a structure-based modeling framework, Cheung et al. introduced a desolvation barrier separating the direct-contact minimum from the solvent-separated minimum in pairwise residue-residue potentials Cheung et al. (2002). Simulations demonstrated that this minimalist desolvation model captures water-expulsion events during the final stages of folding. Chan and colleagues further showed that incorporating a desolvation barrier markedly enhances folding cooperativity Liu and Chan (2005); Ferguson et al. (2009); Chen and Chan (2014); Zhang and Chan (2010). These previous successes provide a strong foundation for the explicit implementation of desolvation effect in coarse-grained models to more realistically investigate the critical contribution of desolvation in LLPS.

In this study, we advance the coarse-grained modeling of biomolecular LLPS by explicitly incorporating desolvation terms into the energy function of CG model. By performing molecular simulations, we systematically investigate how desolvation influences phase behavior, thermodynamic properties, and the microscopic dynamics underlying condensate formation as well as the internal dynamics of the equilibrated dense phase. We further optimize the desolvation parameters using all-atom MD simulations and experimental R_g measurements, yielding a refined framework for realistically modeling sequence-dependent phase behavior of biomolecules at microscopic resolution.

Results

Explicit modeling of desolvation effects in coarse-grained models

In order to elucidate the role of water molecules in mediating inter-residue interactions, we first performed all-atom MD simulations with explicit solvent for a set of amino acid analogues and computed the potential of mean force (PMF) along inter-residue distances (r). Methane, methanol, acetate ion, ammonium ion and methanamide were chosen to represent non-polar, polar, negatively charged, positively charged, and backbone-like residues, respectively. These amino acid analogues cover the key features of most amino acid types. The simulations were performed under ambient conditions (300 K and 1 atm), and the sampled snapshots were then analyzed to determine the inter-molecule distances and the corresponding PMF profiles. Additional simulation details are provided in the Methods section.

Figure 1A shows the PMF profile for the methane system. In addition to the deep potential well corresponding to the direct contact between two methane molecules (r ∼ 4 Å), there exists a shallow potential well at larger distance (r ∼ 7 Å). This shallow minimum arises from transient, solvent-separated pseudo contacts between the solute residues, hereafter referred to as solvent-separated contact. To form a direct contact between the solutes, the water molecules mediating the solvent-separated contacts need to be excluded, which leads to a desolvation barrier between the two potential wells. Similar features have been reported in previous all-atom MD simulations Liu and Chan (2005). The PMF profiles extracted based on the simulations for other amino acid analogues exhibit the same overall shape (Figure 1B and Figure 1—figure supplement 1), but the positions of the direct-contact minimum (r_dc), the solvent-separated minimum (r_ss), and the desolvation barrier (r_b) vary depending on the analogue molecules. The depths of the two potential wells (ϵ and ϵ_ss) and the height of the desolvation barrier (ϵ_b) likewise differ across analogue pairs. The desolvation effect is also intrinsically modulated by solvent polarity, which can be realized experimentally by changing ionic strength or adding cosolvents. To elucidate this relationship, we performed additional simulations of the methane system in water with systematically varied polarity. For simplicity, the solvent polarity was tuned by uniformly scaling the partial charges of the hydrogen and oxygen atoms in the water molecules. This approach directly modulates the strength of solvent-solvent electrostatic interactions, which in turn dictates the stability of the hydration shell surrounding the solute. We compared the standard water model (1.0× charges) with reduced-polarity models (0.9× and 0.8× charges) (Figure 1A). The results reveal that a reduction in solvent polarity leads to a simultaneous decrease in the height of the desolvation barrier and the depth of the solvent-separated minimum. These observations underscore the importance of explicitly modeling the desolvation energetics and exploring the effects of different desolvation strengths on the thermodynamics and kinetics of protein LLPS.

Figure 1.

In this work, we introduce desolvation terms into the HPS model developed by Dignon et al. (2018b). In this framework, each residue is simplified to a single spherical bead. The HPS energy function includes the bonded term and nonbonded term (see Materials and methods). The nonbonded term is described by the short-range pairwise interactions and long-range electrostatic interactions. To account for the desolvation effect, we added two Gaussian functions centered at r_b and r_ss to represent the desolvation barrier and the solvent-separated contact, respectively. The resulted nonbonded energy function with desolvation effect is given by

Here Φ(r, ϵ, λ, σ) denotes the nonbonded term in the original HPS model energy function Dignon et al. (2018b). The parameters λ and σ correspond to the hydrophobicity scale values and effective bead sizes, respectively, for the 20 amino acids.

We first employ a homopolymer protein model with a chain length of 50 residues to illustrate the desolvation effect. For each residue, we use the parameter set extracted by averaging among the intrinsically disordered proteins in the Disprot dataset Aspromonte et al. (2024), which gives an amino acid mass of m_residue = 110 amu (atomic mass unit), a hydrophobicity parameter λ = 0.640, and an effective bead size σ = 0.536 nm. The parameter ϵ_b governs the height of desolvation barrier associated with expelling water molecules between residues, whereas ϵ_ss determines the depth of the local minimum corresponding to a solvent-separated contact. Together, these parameters shape the desolvation potential and modulate the energetic cost of excluding hydration water. By tuning the parameter values of ϵ_ss and ϵ_b, we generate energy functions with varying degree of desolvation contributions (Figure 1C), enabling systematic investigation of how desolvation influences the phase behavior and thermodynamics of IDPs. Unless explicitly stated otherwise, all thermodynamic quantities are reported in terms of a reduced temperature, defined as T ^∗ = k_BT/ϵ. Here, ϵ denotes the characteristic interaction energy scale of the system, which is set to ϵ = 0.2 kcal/mol following the convention of the HPS model Dignon et al. (2018b).

Thermodynamic impact of desolvation on the phase diagram

To capture the essential physics of IDP phase separation, we employed a coarse-grained model of 125 homopolymer chains using a slab simulation protocol (Methods). As a control, we first performed molecular dynamics simulations within the HPS model framework without considering the desolvation terms across a range of temperatures (Figure 2A-C). Representative snapshots clearly illustrate the transition from a homogeneous solution to a stable condensate as the temperature decreases (Figure 2B). We extracted time-averaged density profiles along z-axis to identify coexisting dense and dilute phases and thereby constructed the binodal curve (Figure 2A,C), which reproduces phase behaviors reported in previous molecular simulation studies Dignon et al. (2018b). We then repeated the slab simulations while explicitly including desolvation energetics with varying ϵ_b and ϵ_ss to investigate how desolvation reshapes the phase diagram and condensate structure.

Figure 2.

To elucidate how desolvation regulates thermodynamics of LLPS, we systematically mapped phase diagrams under varying desolvation parameters (Figure 2D and E). Increasing the desolvation barrier ϵ_b monotonically lowers the critical temperature (Figure 2D), suggesting a reduced phase separation propensity. Analysis of residue-residue radial distribution functions showed that higher ϵ_b suppresses density distribution in the barrier region, suggesting an entropic penalty that restricts the conformational volume available for residue-residue association (Figure 2—figure supplement 1C). By contrast, deepening the solvent-separated well ϵ_ss elevates (Figure 2E), which is attributed to an enthalpic stabilization driven by water-mediated bridging as shown by the enhanced population of solvent-separated configurations (Figure 2—figure supplement 1D). These opposing effects, entropic restriction versus enthalpic gain, suggest that the desolvation potential emerges as a key thermodynamic regulator that tunes macroscopic phase behavior via the balance of conformational entropy and solvent-mediated enthalpy (lower panels of Figure 2D and E).

Because the critical temperature is an intrinsic, sequence-specific property of a given system, it is also important to examine how desolvation strength alters the overall shape of the phase diagram when is held constant. To do so, we rescaled the temperature axis by , which effectively normalizes overall energetic strength and isolates desolvation-specific effects on the phase boundary. For the desolvation models with different barrier height ϵ_b, the normalized binodal curves nearly overlap (Figure 2F). By contrast, varying the solvent-separated well depth ϵ_ss results in distinct deviations in the density of the condensed phase (Figure 2G). Specifically, stronger water-mediated interactions lower equilibrium density by biasing the ensemble toward configurations with larger inter-residue separations (Figure 2—figure supplement 1E, F). Thus, solvent-separated interactions directly regulate molecular packing and, by introducing hydrated configurations, prevent artificial over-compaction as typically encountered in coarse-grained models with implicit solvent. These results suggest that explicitly including desolvation enables coarse-grained models to more faithfully capture the physicochemical structure of biological condensates.

Conformational modulation by desolvation in dense and dilute phases

We next investigate how desolvation modulates the conformational distribution of IDPs across dilute and condensed phases. To quantify the global chain dimensions, we calculated the radius of gyration (R_g) as a primary structural metric. The resulting R_g distributions for protein chains in both phases are presented in Figure 3B and C. Notably, the protein chains adopt more compact conformations in the dilute phase (smaller ⟨R_g⟩), whereas those in the condensed phase exhibit more extended conformations (larger ⟨R_g⟩). This conformational expansion upon condensation reflects a transition from the intra-chain dominated interactions to the inter-chain dominated interactions as typically encountered during phase of separation of IDPs Dignon et al. (2018b); Wei et al. (2017); Hazra and Levy (2021). Within the condensed phase, protein chains preferentially form inter-chain contacts that lower the free energy of the system, which compensate for the entropic cost associated with reduced translational and conformational freedom. The R_g distributions also suggest that chain conformations in the dilute phase are much more sensitive to variations in desolvation parameters than in the condensed phase (Figure 3B and C). Consequently, the variation in the conformational difference arises predominantly from the conformational changes of isolated chains in the dilute phase. Physically, ΔR_g captures the transition from intramolecular-dominated collapse in the dilute phase to intermolecular-mediated expansion in the dense phase, thus serving as a metric for the relative strength of inter-chain versus intra-chain interactions.

Figure 3.

To further elucidate the thermodynamic implications of desolvation-modulated conformations, we analyzed the relationship between the conformational shift ΔR_g and the thermal distance to the critical point, defined as . Strikingly, data from all simulated systems nearly collapse onto a single curve, revealing a strong universal correlation between the magnitude of conformational change and the thermal distance to the phase transition point (R² = 0.942, Figure 3D). This result demonstrates that the conformational response to phase separation is governed fundamentally by how far the system resides thermally from the critical point, irrespective of the specific microscopic desolvation energetics.

We rationalize this observation using the Flory-Huggins theory Huggins (1942); Flory (1942) as a bridge between the macroscopic thermal distance and the microscopic interaction strength. In this theoretical framework, the thermodynamic state is defined by the interaction parameter χ = ϵ_eff/(k_BT), where ϵ_eff represents the effective pairwise energy scale specific to each desolvation model Flory (1953). By invoking the critical condition ϵ_eff = k_BT_cχ_c, where χ_c is a universal constant determined solely by the fixed chain length N, we eliminate the system-specific energy scale ϵ_eff.

This substitution reveals that the deviation from the critical interaction parameter scales strictly with the normalized thermal distance, [χ(T_sim) − χ_c] ∝ (T_c − T_sim)/T_sim. The thermodynamic driving force χ can then be related to the structural observable ΔR_g. Since ΔR_g captures the structural transition from a intrachain-interaction dominated state in the dilute phase to an interchain-interaction dominated state in the dense phase, we posit that this conformational shift scales linearly with the excess interaction strength, expressed as ΔR_g ∝ [χ(T_sim) − χ_c]. Combining these relations leads to a theoretical prediction where the thermal distance is related to the conformational change scaled by the simulation temperature

This derivation explains the approximate linearity observed in Figure 3D while also accounting for the systematic deviation caused by the temperature coefficient, which is evident from the varying slopes between low-temperature and high-temperature data points. To validate our theoretical derivation and the underlying linear response assumption, we removed this temperature factor and replotted the data against the normalized thermal distance . As shown in the inset of Figure 3D, this operation yields a nearly perfect linear relationship (R² = 0.986). This statistical enhancement confirms the above argument that ΔR_g scales linearly with the interaction distance χ −χ_c, thereby establishing that the conformational response of IDPs is intrinsically dictated by the thermal distance from the critical point.

We further show the relationship between the phase density difference (Δρ = ρ_dense − ρ_dilute) and the conformational difference (ΔR_g) between coexisting phases in Figure 3E. From the previously established proportionality T_c − T_sim ∝ ΔR_g, one can relate the density contrast Δρ to ΔR_g through the critical-scaling relation by Δρ = A(T_c − T)^β (Equation (6)), where β is the critical exponent. Substituting the linear relation T_c − T_sim ∝ ΔR_g into this expression yields Δρ ∝ (ΔR_g)^β. On logarithmic scales, a clear power-law dependence emerges (Figure 3E), with , which is consistent with our theoretical expectation. The fitted exponent (≈ 0.4) is in reasonable agreement with the theoretical critical exponent β = 0.325 based on 3D Ising model Rowlinson and Widom (2013). Collectively, these results reveal a fundamental coupling between microscopic conformational reorganization and macroscopic packing density during biomolecular condensation, underlying conformational reorganization as a structural signature of the phase transition.

Dynamic consequences of desolvation in phase separation

Beyond governing thermodynamic phase boundaries and conformational distribution, the desolvation potential fundamentally reshapes the dynamical landscape of IDP condensates. The interplay between inter-chain interactions and solvent exclusion not only dictates the stability of the dense phase but also modulates the transport properties within the condensate and the macroscopic kinetics of phase separation. To probe the chain mobility within the dense phase, we employed a strategy analogous to fluorescence recovery after photobleaching (FRAP) experiments following prior molecular dynamics simulation work Yamada and Takada (2023), tracking the spatiotemporal evolution of specific chains initially localized at the slab center (Figure 4A and B). For the standard HPS model without desolvation, the tagged chains achieve a uniform distribution that matches the equilibrium density profile within 50 ns (Figure 4B). This rapid relaxation confirms a liquid-like nature of the condensate. Incorporating desolvation effects significantly slows down the relaxation process. The spatial distribution of tracked chains remains biased toward the initial center position after 50 ns and requires a prolonged timescale (∼ 150 ns) to fully equilibrate. This slowing down suggests that desolvation introduces additional roughness to the energy landscape that govern molecular transport.

Figure 4.

To more rigorously quantify this mobility, we computed the time-averaged mean squared displacement (TAMSD) Metzler et al. (2014) for chains within the slab. The chain motion exhibits sub-diffusive scaling, TAMSD = K ⋅ Δt^α, with a characteristic exponent α ≈ 0.91 (Figure 4A, bottom right panel). This deviation from normal diffusion (α = 1) primarily reflects the viscoelastic nature of the crowded condensed phase with multi-valent dynamic interaction network, where chain motion is hindered by molecular crowding and chain constraints. At extended lag times, the TAMSD curves naturally plateau as they approach the geometric limit imposed by the slab thickness. Notably, our observed exponent is comparable to simulation results reported for protein diffusion in three-dimensional condensates (α ≈ 0.85) Watanabe et al. (2025), consistently highlighting the restricted nature of chain dynamics within dense phase. Given the robustness of α across different interaction parameters, we utilized the generalized diffusion coefficient K to characterize the transport rates.

The modulation of the diffusion coefficient reveals a clear dependence on the interaction parameters (Figure 4C). At a constant temperature, K increases with the desolvation barrier ϵ_b and decreases with the well depth ϵ_ss (Figure 4C, left panels). These seemingly divergent trends are effectively unified by the packing state of the condensate. As illustrated in Figure 4C (right panel), the diffusion coefficient is predominantly correlated with the bulk density (ρ_dense), accompanied by a secondary temperature modulation, where higher temperatures provide a modest increase in mobility beyond the baseline packing constraints. Consistent with the density trends established in Figure 2D and E, a higher ϵ_ss promotes tighter packing, increasing the effective friction, whereas a higher ϵ_b prevents close contacts, creating greater free volume that facilitates diffusion. These results indicate that, at a fixed absolute temperature, macroscopic density serves as the primary determinant of chain mobility.

However, this density-dominated view is incomplete because it convolutes intrinsic energy landscape features with thermodynamic state shifts. Since variations in ϵ_b and ϵ_ss significantly alter the critical temperature, comparisons made at a constant absolute temperature inherently involve systems positioned at different thermodynamic quench depths, defined as the normalized temperature deviation from the critical point. This highlights the distinct thermodynamic driving forces for different systems. To disentangle the intrinsic influence of the desolvation potential from these thermal effects, we renormalized our analysis. First, we scaled the simulation temperature by the critical temperature of each system and utilized the ratio as the control parameter for quench depth, noting that a lower value corresponds to a deeper quench. This ensures the thermodynamic equivalence across different interaction parameters. Second, to factor out the inherent linear dependence of Brownian motion on temperature where the diffusion rate scales with k_BT, we defined a reduced diffusion coefficient as detailed in the Supplementary Information. In this renormalized frame, the intrinsic signature of the energy landscape emerges clearly (Figure 4D). At the same quench depth, where thermodynamic driving forces are comparable, a higher desolvation barrier ϵ_b (darker blue points) yields a lower reduced diffusion rate . Crucially, this kinetic divergence becomes increasingly pronounced at deeper quenches (lower ) but effectively vanishes near the critical point, indicating that the regulatory role of the intrinsic energy landscape amplifies as the system moves further from the critical point. Notably, the standard HPS model exhibits much faster intrinsic mobility (Figure 4 red triangles), confirming that neglecting specific desolvation energetics artificially smooths the energy landscape and leads to an overestimation of chain dynamics. The trend shown in Figure 4C and D uncovers the physical reality masked by density effects. Specifically, it reveals that while a higher ϵ_b loosens the macroscopic packing, it simultaneously roughens the microscopic energy landscape. The desolvation barriers increase the activation energy required for local chain rearrangement, thereby suppressing the intrinsic mobility. Thus, the chain dynamics are governed by a hierarchy of mechanisms: macroscopic density sets the baseline mobility through free volume, while the roughness of the desolvation-mediated energy landscape imposes a fine-tuned kinetic constraint.

Desolvation modulated kinetic arrest of coarsening dynamics

Elucidating the molecular events during the early stage of phase separation is challenging experimentally. We therefore examine droplet-growth dynamics at the early stage of phase separation by performing molecular simulations with explicit consideration of desolvation effect. Initially, the system was equilibrated at a supercritical temperature (T ^∗ = 2.98) for 10 ns to ensure a homogeneous distribution. The temperature was then quenched to T ^∗ = 1.79, located deep within the two-phase coexisting region, to initiate spontaneous phase separation (Figure 4A). Visually, the system evolves from an initially uniform distribution into interconnected structures within 1 ns ∼ 2 ns, which subsequently breaks up into discrete high-density domains along z-axis. Over time, these domains gradually coalesce into larger slabs, exhibiting the characteristic morphological evolution of spinodal decomposition Cahn and Hilliard (1958).

To quantify the phase separation kinetics, we adopted a two-stage metric strategy. In the early stage, where interfaces are diffusive, we monitored the normalized density variance , see Supplementary Information) derived from the instantaneous density profile ρ(z) to track the growth rate of density fluctuations. As clear interfaces emerged, we switched to tracking the average slab size (⟨ξ⟩), defining dense domains as continuous regions where ρ(z) exceeds a threshold value (here the average system density was used as the threshold). To ensure statistical robustness, all kinetic metrics were averaged over six independent slab simulation replicas. The resulting kinetic profiles (Figure 4F and G) reveal two distinct regimes. First, immediately after the temperature quench, the system undergoes a rapid exponential growth in density variance (top panels). This behavior is characteristic of spinodal decomposition Cahn (1961). The results show that stronger effective attraction (higher ϵ_ss or lower ϵ_b) tends to accelerate these initial fluctuations. This is quantitatively evidenced by the monotonic decrease in the characteristic time scale τ_1/2, defined as the time required to reach 50% of the maximum normalized variance. Second, in the late stage (t > 10 ns), the domain growth transitions into a coarsening regime following a power-law scaling ⟨ξ⟩ ∼ t^α with α ≈ 0.26 (Bottom panels). Theoretically, the Brownian motion coalescence (BMC) mechanism predicts an asymptotic scaling of α = 1/3 Majumdar et al. (1994). Our observed exponent falls slightly below this theoretical value, likely due to the viscoelastic nature of the polymer solution. Specifically, the sub-diffusive dynamics of individual chains (TAMSD ∼ Δt^0.91) inherently constrain the macroscopic domain diffusion Doi et al. (1988), thereby retarding the coarsening rate below the standard Brownian prediction.

Intriguingly, between the initial spinodal decomposition and the late-stage Brownian coalescence, we observed a distinct intermediate plateau regime where the domain growth is temporarily arrested. As shown in the slab size evolution (Figure 4F and G, bottom), the domain growth stalls for a distinct period before entering the late-stage scaling regime. This kinetic arrest is visually corroborated by the density profiles at 5 ns in Figure 4E, where distinct domains are observed to form but fail to fuse immediately. The quantitative analysis in Figure 4F and G further reveals that these domains persist as distinct entities for durations that vary significantly depending on the desolvation parameters. This phenomenon is consistent with the characteristics of viscoelastic phase separation Tanaka (2000). It suggests that strong inter-chain interactions within the dense phase form a transient network structure. While surface tension drives the fusion of domains, the internal structural connectivity generates a resistance to the necessary shape deformation, thereby delaying coalescence. Crucially, the plateau duration extends significantly in systems with higher ϵ_ss or lower ϵ_b. This trend corroborates our viscoelastic interpretation, suggesting that stronger inter-chain attractions reinforce the transient network, thereby effectively prolonging the kinetic arrest. This kinetic trap is particularly pronounced in the standard HPS model (red trajectories in Figure 4F and G), due to its relatively stronger pairwise attractions. The sensitivity of kinetic arrest of coarsening dynamics to desolvation parameters underscores the explicit modeling of desolvation in coarse-grained model in order to achieve a more physically plausible molecular simulations of LLPS. Synthesizing these observations, Figure 4E provides a schematic illustration of the three-stage evolution. This diagram maps the progression of driving forces, highlighting the transition from the initial regime dominated by inter-chain interactions, to a final stage governed primarily by surface tension-driven coarsening.

Collectively, we presented a comprehensive analysis of the phase separation dynamics, bridging the microscopic diffusion of single chains to the macroscopic growth of condensed domains. Our analysis underlies the essential role of desolvation potential. Beyond modulating thermodynamic stability and conformational ensembles as discussed earlier, it acts as a critical regulator of phase separation dynamics.

Parameterization of desolvation model for IDPs

The results from the homopolymer model discussed above underscore the critical role of desolvation in regulating both the thermodynamics and dynamics of LLPS. Consequently, explicitly incorporating water-mediated interactions into residue-level coarse-grained models is essential for more realistically capturing the physics of protein phase separation. To this end, we developed a transferable parameterization strategy derived from all-atom simulations and validated it against experimental observables of IDPs.

Experimentally determining the desolvation parameters, specifically, the barrier height ϵ_b and the solvent-separated minimum depth ϵ_ss, remains challenging. To bridge this gap, we leveraged all-atom molecular dynamics simulations of amino acid analogues to extract the corresponding energy scales. Taking the depth of the direct-contact potential well (ϵ) as a reference scale, we expressed the desolvation parameters as ϵ_b = α_bϵ and ϵ_ss = α_ssϵ. By fitting the PMF profiles obtained from all-atom simulations (Figure 1—figure supplement 1), we determined the coefficients α_b and α_ss. Notably, these coefficients were highly consistent across different pairs of amino acid analogues. For simplicity and transferability, we adopted the averaged values α_b = 0.33 and α_ss = 0.06 as the initial baseline for our coarse-grained model.

To evaluate the effectiveness of these physics-derived parameters, we first integrated them into the framework of the standard HPS model. While the HPS model performs remarkably well in capturing sequence-dependent properties of LLPS, previous studies have noted that parameters optimized solely against single-chain R_g can face challenges in transferability to the condensed phase, often leading to an overestimation of phase separation propensity (i.e., higher critical temperatures, T_c) Dignon et al. (2018b). We assessed the model performance on a dataset of 47 IDPs taken from Tesei et al. Tesei et al. (2021). The deviation of simulated single-chain R_g from experimental SAXS data was quantified by , where σ_exp represents the experimental uncertainty. When the desolvation parameters derived from all-atom simulations (α_b = 0.33, α_ss = 0.06) were added to the HPS model without modifying its original ϵ, we observed a substantial improvement in conformational accuracy. The average discrepancy decreased drastically from 1,082.27 (original HPS) to 85.88 (HPS-desolvation) (Figure 5B). This improvement indicates that introducing the desolvation terms effectively mitigates the over-compaction tendency of the original HPS model likely by the competition between direct and water-separated contacts.

Figure 5.

We also incorporated the desolvation terms into CALVADOS2 Tesei and Lindorff-Larsen (2022), a force field whose hydrophobicity parameters (λ) were reparameterized using Bayesian optimization to reproduce experimental R_g and PRE data. Directly adding desolvation terms to CALVADOS2 led to substantial chain over-expansion , reflecting the fact that CALVADOS2 is already finely tuned to reproduce experimental data through a top-down optimization procedure and the inclusion of additional desolvation terms is incompatible with its parameterization. To address this, we implemented a systematic optimization strategy in which the desolvation barrier height (ϵ_b), solvent-separated minimum depth (ϵ_ss), and the overall energy factor ϵ were scanned over broad parameter ranges to identify the parameter set that minimizes deviations between simulated and experimental R_g values. The full optimization workflow is summarized in the flowchart in Figure 5A, and the resulting error landscapes are shown in Figure 5—figure supplement 1 and 2. The optimized parameter set (Table 1) achieved excellent agreement with experiments, reducing to 9.54 (Figure 5C), which is comparable to the accuracy of the original CALVADOS2 model in describing R_g values.

Table 1.

We further probed the phase behavior using the FUS low-complexity domain (FUS LC) as a test case. Experimental observations indicate that the critical temperature is around 25 ^◦C Burke et al. (2015). While the original HPS model overestimates this value with a critical temperature T_c ∼ 340 K, integrating the desolvation terms into the HPS framework (HPS-desolvation) results in a T_c of ∼ 290 K (Figure 5D) which aligns well with experimental observations. In addition to improving the prediction of the critical temperature, we observed that the desolvation model significantly mitigates the overestimation of condensate density common in implicit solvent coarse-grained simulations Regy et al. (2021); Benayad et al. (2020). In the CALVADOS2 framework, the inclusion of desolvation reduces the dense phase density while maintaining a critical temperature comparable to the original model (Figure 5D). Similarly, for the HPS model, when we rescaled the interaction strength ϵ to match the single-chain R_g and thereby align the critical temperature with the HPS-desolvation system (Figure 5—figure supplement 3), we observed the same behavior. Quantitatively, across both frameworks, the desolvation-modulated models consistently yield a packing density approximately 5% lower than the original models. This distinction underscores a limitation of standard Lennard-Jones potentials where a single parameter ϵ simultaneously controls chain dimension, critical temperature, and condensate density. In contrast, our model introduces tunable parameters ϵ_b and ϵ_ss that allow for accurate reproduction of experimental densities without compromising the prediction of phase boundaries. The observed density reduction is physically attributed to the increased population of solvent-mediated contacts relative to direct contacts. This mechanism introduces water-mediated porosity into the droplet and mimics the realistic hydration properties of protein condensates, thereby offering a more physically realistic description of IDP phase separation.

Discussions and Conclusion

Coarse-grained molecular simulations have become indispensable for elucidating the molecular principles underlying LLPS, complementing both theoretical frameworks and experimental observations. While residue-level models with implicit solvent have achieved remarkable success in capturing sequence-dependent phase behaviors of biomolecules, they often face a fundamental tradeoff due to the neglect of explicit solvent excluded volume and hydration effects. Specifically, interaction parameters calibrated to match single-chain dimensions or critical temperatures frequently result in an overestimation of the condensed phase density. This limitation suggests that standard Lennard-Jones potentials, which couple phase stability and density through a single energetic parameter, may lack the physical dimensionality required to fully describe the thermodynamics of hydrated protein condensates.

To address this challenge, we developed a coarse-grained framework that explicitly incorporates desolvation effects into the pairwise energy function. Application of this model to homopolymer systems demonstrates that desolvation substantially reshapes phase behavior. Incorporating desolvation modulates the critical temperature and alters the phase diagram by regulating intermolecular packing within the dense phase. Moreover, the model reveals a tight coupling between desolvation-controlled chain compaction during the formation of condensate phase and the thermal distance from criticality (temperature gap), providing a theoretical basis for linking single-chain conformational changes (ΔR_g) to macroscopic phase behavior. Dynamic analysis further captures the three characteristic stages of the early stage of LLPS, including spinodal decomposition, kinetic arrest, and coarsening. Crucially, it reveals how desolvation effect regulates the kinetics of these stages by modulating the effective inter-chain interactions and the roughness of the energy landscape. These results underscore the importance of accounting for desolvation energetics in governing both the thermodynamics and kinetics of LLPS.

Motivated by these insights, we further developed a transferable parameterization strategy for residue-level models. The resulting desolvation-aware CG model improves the accuracy of IDP conformational ensembles and LLPS thermodynamics while preserving the computational efficiency inherent to CG methodologies. By bridging the gap between the computational efficiency of implicit-solvent representations and the physicochemical realism of water-mediated interactions, this framework provides a mechanistic means to decouple the control of phase stability (T_c) from the regulation of condensed-phase density (ρ_dense). This decoupling effectively mitigates the over-compaction of dense phase commonly observed in implicit solvent simulations with CG models, enabling the molecular simulations qualitatively more consistent with the physical nature of biological condensates.

Despite these advances, the current model employs a single set of desolvation parameters (α_b, α_ss) for all amino acid types. While this simplification is effective for capturing global phase behavior, it overlooks residue-specific desolvation energetics that may influence selective interactions in heterogeneous or compositionally complex sequences. Future extensions of the model could incorporate residue-specific desolvation parameters derived from bottom-up parameterization or expanded experimental datasets, thereby enhancing predictive accuracy for sequence-dependent LLPS.

In summary, this work introduces a versatile and transferable extension to existing coarsegrained force fields that explicitly incorporates the physics of water-mediated interactions. By enabling efficient and realistic treatment of desolvation, the model provides a robust framework for simulating biomolecular condensates with improved thermodynamic and structural fidelity, as well as more accurate kinetic characterization. Future developments, such as residue-specific desolvation energetics and extensions to multicomponent systems involving nucleic acids, hold strong potential for deepening our mechanistic understanding of the complex phase behaviors that organize cellular biochemistry.

Methods

All-atom MD simulations

The all-atom MD simulations were performed using GROMACS. Ten amino acid analogues were solvated in a cubic water box under periodic boundary conditions. The OPLS-AA force field was used to model the solute molecules, and the TIP4P water model was employed to represent solvent. Na⁺ and Cl⁻ ions were added to neutralize the system, yielding a total system size of 1,010 atoms. After an initial energy minimization of 5 × 10⁴ steps, the system was equilibrated in the NVT ensemble at 298 K for 100 ps, followed by an NPT equilibration at 298 K and 1 atm for 1 ns. Temperature and pressure were controlled using the Nosé–Hoover thermostat and Parrinello–Rahman barostat Parrinello and Rahman (1981), respectively. A time step of 2 fs was used for all MD steps. The simulation details are provided in Table 2. Subsequently, a production simulation of 100 ns was performed under the same NPT conditions for structural and thermodynamic analysis. The intermolecular separation r between two amino acid analogues was defined as the distance between their C atoms. The potential of mean force (PMF) was computed as PMF(r) = −k_BT ln P (r), where P (r) is the probability distribution of intermolecular distances obtained from the production trajectory.

Table 2.

Coarse-grained simulations

We utilized the slab method Dignon et al. (2018b) to extract the densities of the coexisting dilute and dense phases. In these simulations, 100 protein chains were placed in a rectangular box of dimensions 15 nm × 15 nm × 280 nm. The density profile along the z-axis was computed to identify the condensed and dilute regions, from which thermodynamic quantities such as the critical temperature were subsequently determined. Simulations were performed over a range of different temperatures. To implement the pairwise interaction potentials of the HPS model, we employed the azplugins package (https://github.com/mphowardlab/azplugins) together with HOOMD-blue v2.9.7 Anderson et al. (2020), running with GPU acceleration. All simulations used a time step of 0.01 ps. Each trajectory was first equilibrated for 3 × 10⁷ steps, followed by a production run of 1 × 10⁸ steps. Additional details of the simulation framework are provided in Table 3.

Table 3.

Hydrophobicity Scale Model

Our desolvation model is developed based on the hydrophobicity scale (HPS) coarse-grained model by Dignon et al. (2018b), in which each residue is represented by a single bead. The potential energy includes bonded term, electrostatic term and short-range pairwise term. The bonded interaction term is given by a harmonic potential V ^b = k_b(r_i,i+1 − r⁰)² with the spring constant k_b = 10 kJ/Ų and the bonded length r⁰ = 0.38 Å. r_i,i+1 represents the distance between the consecutive residues i and i + 1 along the protein chain. The electrostatic interaction term is modeled by a Coulombic potential with Debye-Hückel Debye and Hückel (1923) electrostatic screening caused by salt ion concentration. The function is given by

where κ is the screening length and D = 80, which is the relative dielectric constant of water. For the short-range pairwise interaction, the HPS model uses a 6-12 Lennard-Jones potential with the interaction strength given by hydrophobicity scale Kapcha and Rossky (2014), which can describe the effective interaction between residues. The hydrophobicity value λ, ranging from 0 to 1, varies for different residues. The arithmetic average is used for the interaction parameter between two residues, i.e. hydrophobicity value λ_i,j = (λ_i + λ_j)/2 and amino acid size σ_i,j = (σ_i + σ_j)/2. The energy function of the short-range pairwise interaction is given by

where Φ_LJ is standard 6-12 Lennard-Jones potential

The free parameter ϵ is set as 0.2 kcal/mol, with which the molecular simulations can best reproduce the experimental results of R_g values for a set of IDPs.

Simulation Framework

We use slab simulation scheme in the coexistence simulations. Initially, 100 chains are randomly distributed in a cubic box with periodic boundary condition at 200 K. The box length is first scaled to 15 nm, which can effectively prevent protein chains from interacting with their periodic images generated by periodic boundary conditions. Then, the x- and y-dimensions are set fixed and the z-dimension is extended to 280 nm. The temperature is then gradually increased to the target temperature and the simulations are conducted at different temperatures for 1 μs in the NVT ensemble maintained by a Langevin thermostat, with a friction coefficient of 0.01 ps⁻¹. A time step of 0.01 ps is used for all the simulations. We use HOOMD-Blue v2.9.7 packages Anderson et al. (2020) to conduct the simulations, which utilizes both CPU and GPU resources to accelerate calculation.

The slab density profile along z-axis is then determined after the system gets equilibrated. In order to avoid the influence of periodic boundary conditions, we first move the center of mass of the system in the z-axis direction to the position of z = 0. We then determined the distribution of dense phase and dilute phase of the phase separated system using the defined thresholds. The area with density higher than 0.95 of the maximum density (ρ_max) is treated as dense phase, and their average density is regarded as the phase density of the dense phase (ρ_dense). For the dilute phase, the threshold is taken as ρ_min + 50 mg/mL, where ρ_min is the minimum density. The density of dilute phase (ρ_dilute) is then determined by the same method. If ρ_max − ρ_min < 50 mg/mL, the system is considered not to undergo phase separation.

Critical Temperature

The critical temperature was determined based on the vapor-liquid interfacial properties of Lennard-Jones chains Blas et al. (2008). The critical temperature T_c is obtained from ρ_dense and ρ_dilute at different temperatures by

where β is the critical exponent, which is set as 0.325 Rowlinson and Widom (2013) and A is the fitting parameter. Data with temperature below the critical temperature is used to fit this equation. The phase density at critical temperature is determined by the law of rectilinear diameters Davies (1912)

with B and C being the free parameters determined from the simulation data. At the critical temperature, the distinction between dense phase and dilute phase disappears, which means ρ_dense = ρ_dilute = ρ_c. The phase density is then calculated by ρ_c = B + C ⋅ T_c. A phase diagram which shows temperature change with the density of dense phase and dilute phase is then obtained.

Derivation of the Reduced Diffusion Coefficient

To investigate the intrinsic influence of the interaction potential on chain dynamics, we must isolate the structural contributions from the trivial thermal acceleration. While keeping the quench depth (T/T_c) constant ensures thermodynamic consistency, it necessitates varying the absolute temperature T, which inherently alters the base diffusion rate via thermal fluctuations. To resolve this, we adopt the constant quench depth approach to ensure thermodynamic consistency and introduce a normalized metric to eliminate the trivial thermal velocity contribution.

According to the fundamental Einstein relation, which connects the diffusion coefficient K to the thermal energy k_BT and the friction coefficient ζ:

Here, the friction coefficient ζ encapsulates all resistance to motion. Following the hybrid theoretical framework proposed by Macedo and Litovitz Macedo and Litovitz (1965), which unifies free-volume theory with activated rate theory, the friction term can be factorized into a packing-dependent contribution and an energy-barrier-dependent contribution:

where the first term represents the hindrance due to limited free volume (v_f), governed by the packing density ϕ, and the second term is the Arrhenius factor describing the activated hopping over an effective energy barrier E_eff, which corresponds to the roughness of the energy landscape induced by the desolvation potential.

In our analysis, by fixing the quench depth (T/T_c), we ensure that the macroscopic packing density ϕ remains comparable across different systems (as verified by the phase diagram in main text). Consequently, the free-volume contribution to friction is effectively invariant. Substituting this back into Equation (8), the diffusion coefficient scales as:

The pre-factor T represents the trivial thermal acceleration inherent to Brownian motion. To isolate the intrinsic energy landscape signature, we define the reduced diffusion coefficient by normalizing out this linear thermal contribution:

Thus, serves as a physically rigorous metric to isolate the intrinsic influence of the energy function on chain dynamics. By normalizing out the trivial linear thermal velocity contribution inherent to Brownian motion, this definition is necessary to decouple the specific kinetic modulation imposed by the desolvation potential from the background thermal acceleration and thermodynamic driving forces. Consequently, allows for a direct assessment of how the interaction potential fundamentally alters transport properties, independent of the trivial temperature effects required to maintain thermodynamic consistency.

Calculation of Normalized Density Variance

To quantify the rate of density fluctuation growth during the early stage of phase separation, we first computed the instantaneous spatial variance of the density profile along the z-axis, , defined as

where ρ(z, t) is the local density at position z and time t, and is the mean density of the system. To enable a fair comparison of kinetics across systems with distinct interaction parameters (ϵ_ss and ϵ_b), it is crucial to decouple the rate of phase separation from the difference in equilibrium densephase densities. Therefore, we normalized the instantaneous variance by its equilibrium baseline:

Here, the normalization factor is the time-averaged density variance calculated from the final 30% of the simulation process. This period corresponds to the steady state characterized by the formation of a well-defined slab with stabilized phase densities. By using this normalized metric , we effectively factor out the variations in final density contrast, ensuring that the observed differences in the time-evolution curves purely reflect the acceleration or deceleration of the kinetic process itself.

Desolvation Parameters

Due to the lack of experimental data, the empirical parameters α_b and α_ss were estimated using all-atom molecular dynamics simulations of amino acid analogues with explicit solvent. The potential of mean force (PMF) obtained from the simulation results is shown in Figure 1—figure supplement 1. Methane (CH₄), methanol (CH₃OH), acetate ion , ammonium ion , and methanamide (CH₃NO) were chosen as representative amino acid analogues for non-polar, polar, negatively charged, positively charged, and backbone-like residues, respectively. The α_b, α_ss values were then obtained by fitting the nonbonded energy function with desolvation effect (Equation (1)) to the experimental PMF results. The λ parameter values for these amino acid analogous were assigned based on the λ values of the corresponding amino acids. Specifically, we used λ_GLY = 0.649 for methane, λ_SER = 0.595 for methanol and λ_CYS = 0.595 for methanamide. The results are shown in Table 1.

Data availability

All simulation codes, input files, and scripts for plotting figures of the manuscript are openly available on GitHub at https://github.com/kaizhangnju/desolvation-CG-model. Additional data supporting the findings of this work are available from the corresponding author upon reasonable request due to large file size constraints.

Figure supplements

Figure 1—figure supplement 1.

Figure 2—figure supplement 1.

Figure 5—figure supplement 1.

Figure 5—figure supplement 2.

Figure 5—figure supplement 3.