Capturing chemical reactions inside biomolecular condensates with reactive Martini simulations

The study examines how biomolecular condensates formed via liquid–liquid phase separation (LLPS) modulate chemical reactions. Condensates are implicated in organizing cellular biochemistry and may have played roles in protocell evolution by recruiting molecules and serving as reaction crucibles that alter reaction rates and specificity. While experiments have shown condensates affect enzymatic and redox processes, their interior dynamics are challenging to probe, motivating molecular dynamics (MD) approaches. Atomistic simulations are costly for large, slow condensate systems, so coarse-grained (CG) methods like the Martini force field are widely used for LLPS. An extension, reactive Martini, enables modeling chemical reactions via virtual sites and dummy beads and has been validated for benzene-1,3-dithiol oligomerization via disulfide formation. This work aims to demonstrate reactive Martini simulations of chemical reactions within condensates: (1) reproduce pH-dependent LLPS of a synthetic peptide (LFssFL), (2) quantify partitioning of benzene-1,3-dithiol and peptide-functionalized variants into condensates via PMFs and experiments, and (3) assess how condensates alter macrocycle ring size distributions and reaction rates.

Prior studies established condensates as dynamic organizers influencing enzyme activity and redox reactions, and implicated in early-life chemistry. Computationally, CG models (including optimized IDP models) have provided insights into LLPS. The Martini CG force field is broadly applied across biomolecular and materials contexts, including condensates. Reactive Martini extends Martini to model chemical reactions using tabulated potentials, reproducing benzene-1,3-dithiol disulfide oligomerization. Work on peptide-based condensates (e.g., minimalistic dipeptide stickers and spacers) informs how sticker hydrophobicity and spacer solvation govern phase behavior. Studies also highlight challenges in hydration strength in CG models for IDPs and the need for parameter tuning to better match experimental hydration and partitioning.

Computational: The LFssFL synthetic peptide (two leucine–phenylalanine dipeptides linked by a disulfide spacer) was modeled with Martini 3. At neutral pH, the two N-termini were modeled as protonated (Q5 beads, +1 charge each); at high pH, deprotonated (uncharged P6 beads) to mimic pH-induced LLPS. The disulfide spacer mapped to SP1 beads. For selected simulations, peptide–water interactions were increased by adding a backbone virtual site interacting with water (ε = 0.465 nm, σ = 0.1 kJ/mol) to enhance hydration. Systems were prepared with Polyply and Gromacs; a pre-partitioned geometry for a condensate droplet used Polyply’s gen_coords. Unreactive simulations (LLPS characterization) used standard Martini parameters: LJ and Coulomb cutoffs 1.1 nm, T = 300 K with velocity-rescaling thermostat (τ = 1 ps), P = 1 bar with Berendsen (τ = 4 ps) for equilibration then Parrinello–Rahman (τ = 12 ps) for production, compressibility 3×10^-4 bar^-1, time steps 10 fs (equilibration) and 20 fs (production). Diffusion was analyzed via incoherent scattering functions; fits compared by AIC to determine single vs double exponential decays. Partitioning free energies: PMFs were computed with the accelerated weight histogram (AWH) method along the distance between the center of mass of a pre-equilibrated condensate slab and either a benzene-1,3-dithiol molecule or peptide-functionalized XGLKFK. Each PMF was averaged over 3 independent 1 μs simulations. A variant with increased condensate hydration assessed sensitivity to water content. Reactive simulations: Reactive Martini parameters (tabulated potentials) modeled disulfide bond formation between benzene-1,3-dithiol monomers to form macrocycles. Gromacs 2018.8 with group cutoff scheme enabled tabulated potentials; neighbor list 1.2 nm updated every 10 steps; nonbonded cutoffs 1.1 nm; time step 10 fs; LINCS order 8 with 2 iterations. Production runs were in triplicate. A reduced reactive potential well depth (from 60 to 20 kJ/mol) was also tested to allow reversible bond breaking (thiol–disulfide exchange-like behavior), generating new tabulated potentials via the reactive Martini notebook. Simulation setups: (i) Pre-partitioned systems with reactants enriched inside an LFssFL condensate droplet (~94 mM benzene-1,3-dithiol) assessed product distributions; (ii) Aqueous systems at equal total peptide count (and a higher-concentration small water box) served as controls to disentangle concentration effects; (iii) Out-of-equilibrium co-assembly simulations introduced varying peptide loadings (e.g., 2%, 6%, 12% weight; examples with 500, 1500, 2500 peptides) with 216 reactive molecules to probe concurrent LLPS and reaction kinetics over up to 500 ns–1 μs. Experimental partitioning assay: LFssFL condensates (3 mM, pH 8.2, 200 mM borate) were prepared; XGLKFK 1mer, 3mer, and 4mer were spiked to 0.2 mM each. Samples incubated 30 min at 4 °C, then 60 μL aliquots analyzed by UPLC to determine total concentrations. After centrifugation (5000 rcf, 30 min, 4 °C), supernatant concentrations were measured by UPLC. Condensate phase volume was computed by difference. Calibration curves (0.05–0.2 mM) enabled concentration determination; three independent repeats per component. Partition coefficients were computed as the ratio of condensate to supernatant concentrations. Preoxidation of 1mer to 3mer/4mer minimized hybrid product formation with LFssFL while leaving minimal 1mer to permit thiol–disulfide exchange.

- LLPS reproduction: Martini 3 simulations recapitulated the pH-dependent phase separation of LFssFL: no LLPS at neutral pH; spontaneous LLPS at high pH. Diffusion analysis via ISFs required two exponential decays for LLPS (heterogeneous) vs one for homogeneous neutral pH systems. The high-pH condensate slab contained ~17% water by weight, lower than experimental reports (>60%), consistent with underestimation of peptide hydration in standard Martini. - Partitioning thermodynamics: PMF for benzene-1,3-dithiol showed favorable transfer into the condensate by ΔG ≈ −9 ± 1 kJ/mol, reduced to −7 ± 1 kJ/mol with increased condensate hydration. Peptide-functionalized XGLKFK exhibited ΔG ≈ −10 ± 2 kJ/mol with a distinct interfacial minimum, suggesting interfacial accumulation. Experimentally, XGLKFK 1mer had a partition coefficient Kp = 5.6 ± 0.9 (ΔG ≈ −4.3 ± 0.3 kJ/mol). Larger macrocycles also partitioned: 3mer Kp = 26.4 ± 1.3; 4mer Kp = 7.1 ± 1.2, indicating size-dependent, nontrivial partitioning. - Macrocycle product distributions: Reactive Martini simulations starting with reactants inside a condensate produced a ring size distribution shifted toward larger macrocycles versus purely aqueous systems; while 3mers remained modal, larger rings (including 10–14mers) formed within condensates. Control aqueous simulations at matched higher reactant concentration reproduced the shift toward larger rings, indicating that the primary driver is increased local concentration inside condensates. - Reversibility test: Reducing the reactive potential well depth (60 → 20 kJ/mol) enabled observable bond breakage, yielding formation and breakup of large rings; 2mers and 3mers remained dominant, approaching equilibrium populations in the last ~20 ns. - Out-of-equilibrium co-assembly: When phase-separating peptides were introduced into reacting systems, maximum ring sizes increased (up to 7mers versus 5–6mers in references) after ~100 ns. Increasing peptide weight fractions (2%, 6%, 12%) modestly influenced distributions and stabilized ring formation counts. Reacted rings partitioned into nascent droplets and coalesced into single condensates by ~1 μs in examples. Monomer consumption completed by ~25 ns in a system with 216 reactive molecules; consumption rates were higher at greater peptide concentrations, indicating that condensate formation chaperones monomer partitioning, raises effective concentration inside droplets, and accelerates reaction rates.

The simulations and experiments collectively show that condensates preferentially recruit benzene-1,3-dithiol and peptide-functionalized reactants, increasing their local concentration and thus altering both kinetics and thermodynamics of macrocyclization. The shift to larger ring sizes inside condensates is primarily explained by concentration effects, as confirmed by aqueous controls at matched concentrations. Reaction rates correlate with peptide abundance during co-assembly, supporting a chaperone-like role of condensates in recruiting reactants to dense phases where reactions proceed faster. Environmental factors intrinsic to condensates, notably water content and interfacial properties, modulate partitioning free energies and likely influence product distributions. Allowing bond reversibility by tuning reactive well depth demonstrates the framework’s flexibility to mimic thiol–disulfide exchange conditions and highlights how dynamic exchange can impact steady-state ring populations. Together, these results address the hypothesis that condensates can act as reaction crucibles: they concentrate reactants, speed reactions, and bias product distributions toward larger macrocycles.

This work establishes reactive Martini as a practical tool to simulate chemical reactions inside biomolecular condensates. The study (i) reproduces pH-dependent LLPS of a synthetic peptide (LFssFL), (ii) quantifies favorable partitioning of benzene-1,3-dithiol and XGLKFK into condensates with agreement between PMFs and experimental partitioning, and (iii) demonstrates condensate-driven changes in macrocycle ring size distributions and accelerated reaction rates via enhanced local concentrations. Sensitivity to condensate environment (e.g., hydration) and the ability to model reversible disulfide chemistry by tuning reactive potentials were shown. Future directions include extending to other reactions and condensate systems (e.g., validated with WGR-1), improving reactive potentials to better match benzenedithiol macrocycle formation and stack stabilization (e.g., incorporating peptide β-sheet interactions), leveraging titratable Martini for pH-dependent dynamics, and overcoming current performance limits (e.g., adopting Verlet scheme in Gromacs) to access longer timescales for concurrent phase separation and reactions.

- Hydration underestimation: The Martini peptide model underestimates condensate water content (~17% vs >60% experimentally), impacting partitioning free energies and potentially product distributions. - Timescale and performance: Reactive Martini currently relies on tabulated potentials with the group cutoff scheme in Gromacs 2018.8, limiting accessible timescales; LLPS completion can be slower than reaction completion, complicating full co-assembly equilibration. - Reaction model: The primary model does not include disulfide bond breaking; reversibility was mimicked only by reducing the reactive potential well depth. Accurate pH-dependent thiol–disulfide exchange kinetics are not explicitly modeled. - Potential optimization: Observed very large macrocycles (10–14mers) may reflect model parameters; further refinement of reactive potentials and inclusion of peptide side-chain β-sheet interactions could improve realism. - Simplified pH treatment: Protonation states were fixed (homogeneous) to mimic pH changes, not dynamically titrated.