# Active learning of the thermodynamics–dynamics tradeoff in protein
condensates

Yaxin An  Present address: Department of Chemical Engineering, Louisiana State
University, Baton Rouge, LA 70803 USA  Department of Chemical and Biological
Engineering, Princeton University, Princeton, NJ 08544 USA
Department of Chemistry, Princeton University, Princeton, NJ 08544 USA
Michael A. Webb<EMAIL_ADDRESS>Department of Chemical and Biological
Engineering, Princeton University, Princeton, NJ 08544 USA
William M. Jacobs<EMAIL_ADDRESS>Department of Chemistry, Princeton
University, Princeton, NJ 08544 USA

###### Abstract

Phase-separated biomolecular condensates exhibit a wide range of dynamical
properties, which depend on the sequences of the constituent proteins and
RNAs. However, it is unclear to what extent condensate dynamics can be tuned
without also changing the thermodynamic properties that govern phase
separation. Using coarse-grained simulations of intrinsically disordered
proteins, we show that the dynamics and thermodynamics of homopolymer
condensates are strongly correlated, with increased condensate stability being
coincident with low mobilities and high viscosities. We then apply an “active
learning” strategy to identify heteropolymer sequences that break this
correlation. This data-driven approach and accompanying analysis reveal how
heterogeneous amino-acid compositions and non-uniform sequence patterning map
to a range of independently tunable dynamical and thermodynamic properties of
biomolecular condensates. Our results highlight key molecular determinants
governing the physical properties of biomolecular condensates and establish
design rules for the development of stimuli-responsive biomaterials.

## Introduction

Biomolecular condensates—also known as “membrane-less organelles”—are self-
organized structures within the cytoplasm and nucleoplasm of living cells [1,
2, 3]. Condensates play diverse roles in a wide variety of biological
processes [4, 5, 6], in large part because of their ability to concentrate
proteins, RNAs, and small molecules in chemically specific environments [7, 8,
6, 5]. Nonetheless, the ability of condensates to tune the rates of
biochemical reactions is not only determined by their biochemical
compositions. Dynamical properties, such as viscosities and molecular
mobilities, also affect how condensates exchange biomolecules with the rest of
the cell [9, 10, 11] and interact with other intracellular structures, such as
chromatin [12] and cytoskeletal components [13]. Dynamical properties are also
directly relevant to condensate assembly and disassembly kinetics [14] and
irreversible aging processes [15]. Recent experiments have demonstrated that
dynamical properties vary widely across different condensates, both in vitro
and in vivo [11, 16].

Many condensates are believed to form spontaneously as a result of phase
separation taking place close to a local thermodynamic equilibrium [1, 17].
The thermodynamic stability of a condensate, which can be quantified by the
critical temperature of a biopolymer solution [18] or by the saturation
concentration required for phase separation [9, 19], is an equilibrium
property related to the multiplicity and strengths of the “multivalent”
interactions among the biomolecules comprising a condensate [9, 20, 21].
Intuitively, one might expect that the thermodynamic stability of a condensate
is correlated with its dynamical properties, which also emerge from the
collective effects of biomolecular interactions in the condensed phase. This
anticipated relationship represents a “tradeoff” between stability and
dynamics, since stronger thermodynamic driving forces are expected to
correlate with lower molecular mobilities [22]. Experimental evidence suggests
that such a correlation exists for certain ribonucleic condensates, such as
Arg/Gly-rich repeat protein/RNA condensates [11] and LAF-1/RNA condensates
[16]. However, understanding the general relationship between thermodynamic
stability and the internal dynamics of biomolecular condensates requires
additional and systematic exploration [23, 24].

Determining whether a stability–dynamics tradeoff extends to phase-separating
biopolymers is challenging due to the enormous diversity of biopolymers that
can form condensates. Here, we focus on condensates composed of intrinsically
disordered proteins (IDPs), a class of proteins that tend not to adopt stable
secondary and tertiary structures and are known to be essential components of
many naturally occurring condensates [25]. To study the properties of IDP
sequences, we employ a coarse-grained (CG) model that represents IDPs at
residue-level resolution [26]. Models of this type have been shown to
reproduce experimental measurements of both single-molecule and condensed-
phase properties [18, 27, 28, 29], suggesting that CG IDP models capture
essential physics of disordered polypeptides and retain sufficient chemical
specificity to predict differences among condensates composed of various IDP
sequences. Thus, CG IDP models seem well-suited to systematically explore the
relationship between IDP stability and condensate dynamics.

In this article, we combine molecular dynamics simulations and machine
learning to navigate the sequence space [30, 31] of IDPs and elucidate the
relationship between the phase behavior and internal dynamics of single-
component IDP condensates. We first establish a strong correlation between
proxies for condensate stability and dynamics in homomeric polypeptides, which
is understood using simple physical arguments. We then deploy a
computationally efficient machine-learning strategy based on “active
learning”[32, 33, 34] to identify IDP sequences that break this correlation,
exhibiting faster internal dynamics than homomeric polypeptides that assemble
into equally stable condensates. In this way, we identify “Pareto-optimal” IDP
sequences, meaning that sequence perturbations cannot enhance the dynamics
further without reducing the stability of the condensate. Finally, we examine
sequence features of Pareto-optimal sequences and perform a counterfactual
analysis[35, 36] to identify the sequence determinants of the limiting
thermodynamics–dynamics tradeoff. Taken together, our results demonstrate how
sequence design can be used to tune thermodynamic and dynamical properties
independently in the context of biomolecular condensates.

## Results

### Thermodynamic and dynamic properties of homomeric polypeptides are
strongly correlated

To establish a baseline expectation for the thermodynamics–dynamics tradeoff
in IDP condensates, we begin by studying a collection of homopolymeric
sequences at $300\text{\,}\mathrm{K}$. Specifically, we consider homopolymers
consisting of each amino-acid type with chain lengths of $N=20$, $30$, $40$,
or $50$ monomers. Molecular dynamics simulations are performed using the CG
model developed by Mittal and coworkers [26, 28], which treats IDPs at the
single-residue level in implicit aqueous solvent. Amino-acid residues interact
via a combination of bonded and nonbonded pair potentials, the latter of which
account for both electrostatic interactions and short-ranged hydrophobic
forces. We use a Debye screening length of $1\text{\,}\mathrm{n}\mathrm{m}$,
corresponding to an ionic strength of $0.1\text{\,}\mathrm{M}$. See “Model of
intrinsically disordered proteins” in Materials and Methods for further
details.

Throughout this work, we utilize the second-virial coefficient, $B_{2}$, as a
proxy for the thermodynamic stability of the condensed phase for IDPs that
undergo phase separation. $B_{2}$ quantifies the net attractive or net
repulsive intermolecular interactions of a pair of molecules in dilute
solution [37]. It is well established that $B_{2}$ can also be used to predict
the critical point of simple atomic and molecular fluids [38], colloidal
suspensions [39, 40], and homopolymer solutions [41]. A strong correlation
between $B_{2}$ and the critical temperature for IDP phase separation has also
been reported, albeit for a limited set of IDP sequences [18]. In all these
systems, a negative $B_{2}$ value—implying net attractive interactions—is
necessary for phase separation to occur. Although this condition turns out to
be insufficient to guarantee phase separation in heteropolymer solutions [42,
43], we find that $B_{2}$ nonetheless correlates strongly with the difference
between the coexisting phase densities for IDP sequences that do phase
separate, as we demonstrate using molecular dynamics simulations below. In
practice, we compute $B_{2}$ by calculating the potential of mean force,
$u(r)$, between the centers of mass of two polymer chains at
$300\text{\,}\mathrm{K}$ using adaptive biasing force simulations [44] (Fig.
1a; see also “Physical property calculations” in Materials and Methods).
$B_{2}$ is then obtained from the equation

$B_{2}=2\pi\int_{0}^{\infty}dr\,r^{2}\left[1-e^{-\beta u(r)}\right],$ (1)

where $r$ represents the distance between the center of mass of each chain and
$\beta\equiv 1/k_{\text{B}}T$. For convenience, we report $B_{2}$ relative to
a reference volume, $V_{0}=$5529\text{\,}\mathrm{\SIUnitSymbolAngstrom}^{3}$$,
equal to the pervaded volume of an ideal polymer chain,
${V_{0}=(4\pi/3)b_{0}^{3}(N/6)^{3/2}}$ [41], with the CG equilibrium bond
length, $b_{0}=$ $3.8\text{\,}\mathrm{\SIUnitSymbolAngstrom}$, and a chain
length of $N=50$.

Figure 1: Coarse-grained simulations of homomeric polypeptides predict a
strong relationship between condensate stability and internal dynamics. (a)
The potential of mean force, $u(r)$, between the centers of mass of two
polymer chains at $300\text{\,}\mathrm{K}$. Three example curves are shown for
net repulsive, weakly attractive, and strongly attractive 20-mers,
$\text{A}_{20}$, $\text{F}_{20}$, and $\text{Y}_{20}$, respectively. (b)
Empirical equation-of-state curves at $300\text{\,}\mathrm{K}$ for the three
homomeric polypeptides in panel a. Only two of these 20-mers, $\text{F}_{20}$
and $\text{Y}_{20}$, are predicted to phase separate to form condensed phases
with the indicated densities $\rho_{\text{c}}^{\text{F}_{20}}$ and
$\rho_{\text{c}}^{\text{Y}_{20}}\\!$, respectively. (c) The dimensionless
self-diffusion coefficients, $D/D_{0}$, of homopolymers in the condensed phase
are anticorrelated with their scaled dimensionless second-virial coefficients,
$-B_{2}N^{1/2}/V_{0}$. Points are shown on a log–log plot for the 29
homopolymers (out of 80 simulated) that phase separate at
$300\text{\,}\mathrm{K}$. Self-diffusion coefficients are obtained from
canonical-ensemble simulations conducted at the condensed-phase density,
$\rho_{\text{c}}$, determined from the equation-of-state analysis. Statistical
errors reflect the standard error of the mean and are comparable to the symbol
size.

Next, we develop a method to confirm whether an IDP sequence phase separates
and, if so, to determine dynamical properties in the condensed phase. We
perform an equation-of-state (EOS) analysis to approximate the condensed-phase
density, $\rho_{\text{c}}$, for polymers with negative second-virial
coefficients (Fig. 1b; see also “Physical property calculations” in Materials
and Methods). The EOS is a relationship between the osmotic pressure and
volume of the polymer solution at constant temperature. In macroscopic
systems, the pressure is a monotonically increasing function of the density at
constant temperature. However, in finite-size simulations, the EOS exhibits a
non-monotonic “van der Waals loop” [38, 45] at temperatures for which phase
separation occurs; fundamentally, this non-monotonicity arises due to the
inability of a small system to form a stable interface between coexisting bulk
phases. We therefore use short simulations of 100 polymer chains at
$300\text{\,}\mathrm{K}$ to compute the pressure as a function of density, and
we identify non-monotonic behavior as evidence of a phase-separated region.
Making the approximation that the coexistence pressure is near zero,
corresponding to a dilute low-density phase, we approximate the condensed-
phase density, $\rho_{\text{c}}$, as the highest-density root of a non-
monotonic EOS. By contrast, we consider a strictly non-negative EOS as
evidence that phase separation does not occur at $300\text{\,}\mathrm{K}$.
Following this approach, we find that 29 out of the 80 homomeric polypeptides
phase separate at 300 K. We later show that this method, which is amenable to
high-throughput simulations, is accurate to within a few percent of
$\rho_{\text{c}}$ values determined via direct-coexistence slab simulations
[28].

Finally, we use the self-diffusion coefficient, $D$, of a tagged chain within
the bulk condensed phase as a simple measure of the internal condensate
dynamics. We calculate $D$ from the mean-squared displacement of a tagged
chain in a canonical-ensemble simulation conducted at the condensed-phase
density, $\rho_{\text{c}}$ (see “Physical property calculations” in Materials
and Methods). This quantity is negatively correlated with the viscosity of the
condensed phase in phase-separated homopolymer solutions (Fig. S1 in the
Supplementary Materials). However, the greater statistical uncertainties in
the viscosity calculations make $D$ the more practical choice for quantifying
internal condensate dynamics in high-throughput simulations. We report self-
diffusion coefficients relative to a reference value,
$D_{0}=$0.42\text{\times}{10}^{-9}\text{\,}\mathrm{m}^{2}\mathrm{/}\mathrm{s}$$,
which corresponds to an ideal Rouse chain with chain length $N=50$,
$D_{0}=\tau k_{\text{B}}T/N\bar{M}$ [41], where
$\tau=$1\text{\,}\mathrm{p}\mathrm{s}$$ is the damping time in our Langevin
simulations and
$\bar{M}=$118\text{\,}\mathrm{g}\mathrm{/}\mathrm{m}\mathrm{o}\mathrm{l}$$ is
the average molecular weight of the 20 amino acid types.

Plotting $D/D_{0}$ versus $B_{2}/V_{0}$, we observe a strong negative
correlation between these quantities across all homopolymeric sequences that
are determined to phase separate (Fig. 1c). Surprisingly, this correlation is
stronger than the direct correlation between the condensed-phase density,
$\rho_{\text{c}}$, and either of these quantities alone (Fig. S2 in the
Supplementary Materials). This indicates that even though the condensed-phase
density at $300\text{\,}\mathrm{K}$ cannot be predicted solely on the basis of
the second-virial coefficient, this quantity nonetheless captures the key
physics necessary to predict the internal dynamics of a homopolymer condensate
at its equilibrium density.

To rationalize the empirical relationship between $D/D_{0}$ and $B_{2}/V_{0}$,
we consider a Rouse model [41] of an unentangled polymer melt. This model can
describe the condensed phases in these simulations, which have densities on
the order of $1\text{\,}\mathrm{g}\mathrm{/}\mathrm{m}\mathrm{l}$. In this
model, the self-diffusion coefficient is inversely proportional to the total
frictional force experienced by the chain, which typically follows an
Arrhenius law for a thermally activated process [24, 41]. It is thus
reasonable to relate the logarithm of the total friction to the reversible
work required to separate a pair of chains, which is proportional to
$\log(-B_{2})$ in the limit of strong attractive interactions. Accounting for
the number of interacting chains within the pervaded volume of a tagged
polymer [41], we therefore propose that the condensed-phase self-diffusion
coefficient for phase-separating homopolymers should scale as
$D\sim-N^{1/2}B_{2}$. Empirically, we find that including the prefactor
$N^{1/2}$ in the scaling relationship indeed improves the correlation across
the 29 phase-separating homopolymers that we simulated (Fig. 1c and Fig. S3 of
the Supplementary Materials). We thus conclude that this simple model captures
the essential relationship between the second-virial coefficients and
condensed-phase self-diffusion coefficients of homomeric polypeptide
sequences.

### Short polypeptides from known intrinsically disordered regions do not form
condensed phases

Having established a correlation between $D/D_{0}$ and $B_{2}/V_{0}$ for
homomeric polypeptides, we next examine the behavior of heteromeric
polypeptides based on extant IDP sequences. To do so, we perform CG
simulations of 1,266 short ($20\leq N\leq 50$) polypeptide sequences from
DisProt [46], a manually curated database of disordered regions of proteins.
This range of sequence lengths is used throughout this study to limit the
total computational expense of these simulations. Although these polypeptides
feature diverse sequence characteristics (Fig. S4 in the Supplementary
Materials), none of these polypeptide sequences form condensed phases
according to the criteria set forth in the prior section. This observation
could be attributed to potential bias within the DisProt database towards
soluble sequences [47], the imposition of an upper limit on sequence lengths
probed in this study, which excludes many longer extant IDP sequences that are
known to play important roles in intracellular phase separation [46, 48], and
our focus on single-component systems (e.g., exclusion of RNAs that may be
required to observe phase-separation [6]).

Of all analyzed DisProt sequences with $20\leq N\leq 50$, we find that the
vast majority (c.a. 80%) have positive $B_{2}$ (Fig. S5 in the Supplementary
Materials), while the remaining 255 sequences exhibit monotonically increasing
EOS curves, indicating that phase separation does not occur at
$300\text{\,}\mathrm{K}$ when simulated using this CG IDP model. Nonetheless,
we note that some of these $B_{2}$ values are more negative than some of the
phase-separating homomeric polypeptides considered in Fig. 1. This observation
highlights the importance of performing the EOS analysis to confirm whether a
polypeptide undergoes phase separation at $300\text{\,}\mathrm{K}$.

Figure 2: Active-learning approach to heteromeric sequence design. (a) The
workflow for designing new IDP sequences by integrating coarse-grained
modeling, machine learning (ML), and genetic-algorithm-based optimization.
Active learning is an iterative process consisting of (i) performing molecular
dynamics simulations to compute three physical properties of IDPs ($B_{2}$,
phase behavior, and $D$); (ii) learning the predictive features of IDP
sequences by training ML models to estimate these three physical properties;
(iii) designing new sequences by maximizing a fitness function; and (iv)
returning to step (i) to test the physical properties of the proposed IDP
sequences. The data set and ML models are subsequently updated via successive
iterations. (b) The relationship between the scaled condensed-phase self-
diffusion coefficients, $D/D_{0}$, and scaled second-virial coefficients,
$-B_{2}/V_{0}$, of all the designed IDP sequences that undergo phase
separation. The color of each grid element indicates the number of distinct
sequences within that region. The red highlighted area represents a one-
standard-deviation region spanned by the Pareto-optimal heteromeric sequences
(red stars). Homopolymer sequences (white circles) are shown for comparison.
(c) Convergence analysis of the active-learning approach. The hypervolume
under the Pareto front (blue symbols, left axis) increases with active-
learning iterations, while the fraction of newly updated data points at the
Pareto front (black symbols, right axis) tends to zero as the Pareto front
converges. Filled symbols/solid lines show the convergence of the entire
Pareto front, while open symbols/dashed lines describe the convergence of the
Pareto front in the regime where $-B_{2}/V_{0}\lesssim 50$. A final
“exploitation” round (shaded region) of sequence design confirms convergence
of the hypervolume, introducing only a few new sequences into the Pareto front
representation.

### Active learning identifies limits of tunability for condensate
thermodynamic-dynamic properties

Since the lack of phase separation precludes analysis of condensed-phase
dynamical properties using short DisProt sequences, we set out to generate
novel heteromeric sequences and explore the thermodynamics–dynamics tradeoff
beyond homopolymers. To this end, we adopt an active learning scheme [32] that
iteratively selects and simulates new polypeptides to probe this relationship
methodically and efficiently (Fig. 2a). The selection of new polypeptides is
guided by Bayesian optimization, which has been usefully deployed in materials
design, including the design of heteropolymer sequences, by informing the next
“best” sequence(s) to characterize [30, 49, 50, 51, 52, 53]. Within the
context of the thermodynamics–dynamics tradeoff, the “best” sequences are
those considered as Pareto optimal [54, 55], meaning that no other polypeptide
sequence can be found that simultaneously exhibits an increased $D$ and a
decreased $B_{2}$.

Our active learning scheme integrates Bayesian optimization with supervised
machine-learning models and a genetic algorithm for designing polypeptides.
Polypeptides are chosen for characterization based on a policy of expected
hypervolume improvement (EHVI) [56], which has previously been demonstrated to
converge towards a true set of Pareto-optimal sequences [57]. These Pareto-
optimal sequences are referred to as the Pareto front. Optimizing EHVI over
sequence space is facilitated by a genetic algorithm that leverages machine-
learning (ML) models, trained using the results of prior CG simulations, to
predict physical properties (i.e., phase behavior, $B_{2}$, and $D$) from
sequence features of prospective polypeptides. In particular, the ML models
employ a 30-dimensional feature vector that incorporates the amino-acid
composition [58, 59] and ten additional sequence-level descriptors, including
the sequence length, $N$; the average hydrophobicity per residue,
$\bar{\lambda}$; the sequence hydropathy decoration parameter, SHD [60]; the
fraction of positively and negatively charged residues, $\bar{q}_{+}$ and
$\bar{q}_{-}$; the net charge per residue, $|q|$; and the sequence charge
decoration, SCD; the mean-field second virial coefficient [60],
$B_{2}^{\text{MF}}$; the Shannon entropy, $S$; and the average molecular
weight, $\bar{M}$. All the ten sequence-level feature values are normalized
and standardized as described in “Design of novel intrinsically disordered
proteins” in Materials and Methods. Each iteration of active learning involves
training ML models with the current simulation data, generating 96 new
polypeptide sequences via the genetic algorithm, and subsequently simulating
these sequences according to the protocols introduced previously. Additional
technical details are reported in Materials and Methods, and complete
descriptions of ML model generation and validation are provided in the
Supplementary Materials (see Fig. S6).

We find that this active-learning approach rapidly identifies a representative
set of Pareto-optimal polypeptides (Fig. 2b and Fig. S7). First, using a ML
model for $B_{2}$ trained using DisProt sequence data, we perform an initial
round of simulations involving IDPs chosen to minimize $B_{2}$ only (Iteration
0). These simulations establish relevant data for training a ML model for $D$,
enabling subsequent multiobjective optimization via EHVI (see “Design of novel
intrinsically disordered proteins” in Materials and Methods). We
quantitatively assess convergence of the active-learning approach by
monitoring both the hypervolume (HV) under the current approximation of the
Pareto front and the fraction of newly added Pareto-optimal sequences relative
to the current number of Pareto-optimal points (Fig. 2c). Active learning
initially leads to a rapid expansion of the Pareto front (Iterations 1–3),
whereas subsequent iterations (Iterations 4–7) result in substantially
diminished expansion and mostly inconsequential additions to the Pareto front,
especially along the weakly-interaction portion ($-B_{2}/V_{0}\lesssim 50$) of
the front (dashed lines in Fig. 2c). In these latter rounds, the expansion of
the front occurs primarily via the addition of sequences with extremely
negative $B_{2}$ and low $D$. We therefore halt optimization based on EHVI at
this point and perform one last iteration of sequence design based on pure
exploitation of the ML models, which confirms convergence of the Pareto front
(Fig. 2c). The final data set includes 2,114 polypeptide sequences, including
35 Pareto-optimal sequences (denoted P1–P35, see Table S1) that provide a good
representation of the true Pareto front of heteromeric polypeptides.

By contrast with homomeric polypeptides, active learning identifies “designed”
heteropolymer sequences spanning a range of $D$ values at fixed $B_{2}$ and
vice versa (Fig. 2b). Although heteromeric polypeptides follow the same trend
of decreasing $D$ with decreasing $B_{2}$ as found for as homomeric
polypeptides, the data cannot be neatly collapsed by any obvious scaling. For
strongly attracting polypeptides ($-B_{2}/V_{0}\gtrsim 50)$, all sequences
possess limited condensed-phase mobility. Most of these, including the
homomeric sequences Y40, W40, and W50, are within the statistical error region
spanned by Pareto-optimal heteromeric sequences (Fig. 2b). However, more
weakly interacting polypeptides at the Pareto front ($-B_{2}/V_{0}\lesssim
50$) exhibit $D$ values that are roughly twice those of homopolymer sequences
with the same $B_{2}$. For example, at a fixed $B_{2}/V_{0}\approx 10$ in Fig.
2b, the reduced diffusion coefficient $D/D_{0}$ increases from $\sim 0.4$ to
$\sim 0.9$. (We note that, because the Bayesian optimization only targets
expansion of the Pareto front, it is likely that the maximum dynamic range of
$D$ values at fixed $B_{2}$ can exceed a factor of two.) Thus, Fig. 2b
illustrates the overall prospects and limitations of tuning the
thermodynamic–dynamic tradeoff via sequence design. Specifically, when IDPs
interact relatively weakly, sequence design might be used to fine-tune the
dynamical properties of single-component IDP condensates that possess similar
thermodynamic stabilities (here assessed by $-B_{2}$) or, alternatively, to
design IDP condensates that likely undergo phase separation below disparate
critical temperatures while exhibiting similar dynamical properties.

Figure 3: Validation of Pareto-optimal IDP sequences using direct-coexistence
simulations. (a) Snapshots of slab-geometry direct-coexistence simulations at
$300\text{\,}\mathrm{K}$ for five representative IDP sequences from the
thermodynamics–dynamics Pareto front (see Fig. 2b): P1 (blue), P8 (orange),
P14 (red), P19 (cyan), and P35 (purple). Slab simulations of 1000 chains are
performed in the canonical ensemble for
$4.5\text{\,}\mathrm{\SIUnitSymbolMicro s}$ (P1) and
$2\text{\,}\mathrm{\SIUnitSymbolMicro s}$ (P8, P14, P19, and P35). (b) Average
density profiles, measured orthogonally to the interface between the dilute
and condensed phases, obtained from the simulations in panel a. The final
$3\text{\,}\mathrm{\SIUnitSymbolMicro s}$ of the P1 trajectory and the final
$0.5\text{\,}\mathrm{\SIUnitSymbolMicro s}$ of the P8, P14, P19, and P35
trajectories are analyzed to compute these density profiles. Shaded areas
indicate each bulk condensed phase (except for P14 and P19 for clarity),
defined as the region where variations in the average density are less than
$0.02\text{\,}\mathrm{g}\mathrm{/}\mathrm{m}\mathrm{l}$. (c) The correlation
of the dimensionless second virial coefficient $-B_{2}/V_{0}$ and the density
difference between the condensed and dilute phases, $\rho_{C}-\rho_{D}$. (d)
Comparison of the dimensionless condensed-phase self-diffusion coefficients,
$D/D_{0}$, computed via the equation-of-state method within the active-
learning framework (AL) or via the direct-coexistence (DC) simulations shown
in panel a. The dashed line shows a linear fit to the data with a coefficient
of determination of $\text{R}^{2}=0.996$ and a slope of $0.94$.

### Coexistence simulations validate Pareto-optimal sequences discovered via
active learning

We validate the IDP sequences designed via active learning by performing
large-scale molecular dynamics simulations of phase coexistence. We carry out
direct coexistence simulations of five representative sequences from the
thermodynamics–dynamics Pareto front (P1, P8, P14, P19, and P35) using 1000
chains in a slab geometry [18, 27] at $300\text{\,}\mathrm{K}$. Simulation
durations of at least $2\text{\,}\mathrm{\SIUnitSymbolMicro s}$ are required
to reach equilibrium, resulting in the formation of an interface between
equilibrium condensed and dilute phases (Fig. 3a). The interfaces between the
coexisting dilute phases also become sharper as we move along the Pareto front
from P1 to P35 (Fig. 3b), which is consistent with a trend towards increasing
thermodynamic stability. Importantly, we also find that the differences
between the condensed and dilute-phase densities are anticorrelated with
$B_{2}$, which decreases from P1 to P35 (Fig. 3c). Both observations are
consistent with the previously reported correlation between critical
temperatures and second-virial coefficients in single-component IDP solutions
[18] and support our use of $B_{2}$ as a proxy for the condensed-phase
thermodynamic stability.

Importantly, these simulations recapitulate the condensed-phase dynamical
properties of Pareto-optimal sequences predicted by our active-learning
approach. From equilibrated simulations of phase coexistence, we determine the
mean condensed-phase density and compute the self-diffusion coefficient of
chains within the bulk of the condensed phase (i.e, away from the interface
with the dilute phase). The physical properties determined in this way compare
favorably with the values determined from the EOS analysis (Fig. 3c and Fig.
S8 in the Supplementary Materials) and single-phase simulations of condensed-
phase dynamical properties (Fig. 3d). In all cases, quantities computed via
these different methods differ by less than 10% and typically agree within the
statistical error. This agreement demonstrates that the high-throughput
simulations utilized by our active-learning approach are sufficiently accurate
compared to more computationally expensive direct-coexistence simulations.
Furthermore, these results suggest that our active-learning protocol (Fig. 2c)
indeed converges to a representative sequences on the thermodynamics–dynamics
Pareto front for this CG IDP model.

### Sequence features change drastically along the thermodynamics–dynamics
Pareto front

To better understand what governs the thermodynamics–dynamics tradeoff for
heteromeric polypeptides (Fig. 2b), we analyze the sequence characteristics of
IDPs on the Pareto front (Fig. 4a, stars). An aggregate comparison among the
Pareto-optimal polypeptides (P1–P35) is based on the pairwise cosine
similarity of their feature vectors, where higher scores indicate greater
similarity. For clarity, we decompose sequence similarity into the
contributions from _(i)_ the amino-acid composition (Fig. 4b, lower triangle)
and _(ii)_ the various sequence-derived descriptors (Fig. 4b, upper triangle).
The average cosine similarities among pairs of Pareto-optimal sequences are
_(i)_ 0.60 and _(ii)_ 0.56, respectively. For reference, the average cosine
similarities among pairs of DisProt sequences ($20\leq N\leq 50$) are _(i)_
0.57 and _(ii)_ 0.22; these reference values are indicated by the vanishing
points of the color bars in Fig. 4b.

Based on relative similarities, we propose that the Pareto-optimal sequences
roughly divide into three groups: P1–P8, P9–P20, and P21–P35, outlining a
high-$B_{2}$/high-$D$ regime, a transition regime, and a low-$B_{2}$/low-$D$
regime, respectively (top-left, middle, and bottom-right in Fig. 4a,b). While
P1–P8 generally exhibit stronger composition and sequence-descriptor
similarity than randomly selected pairs of disordered regions from DisProt,
they significantly differ from the other groups. Collectively, P9 to P35 tend
to show strong composition and sequence-descriptor similarities with one
another. However, the transition group (P9 to P20) has lower intra-group
similarity compared to the low-$B_{2}$/low-$D$ group (P21 to P35) (e.g., see
also dashed lines in Fig. S9 of the Supplementary Materials), which motivates
their distinction.

Figure 4: Physicochemical features of sequences at the thermodynamics–dynamics
Pareto front. (a) Representative IDP sequences along the Pareto front (stars)
divide into three regimes. Also shown are counterfactual sequences (circles)
of three representative Pareto-optimal sequences, P1 (blue), P8 (orange), and
P35 (purple); see text for the criteria used to select counterfactual
sequences. (b) Cosine similarities between the amino-acid composition
component of the feature vectors (lower triangle, left color bar) and between
the sequence-descriptor component of the feature vectors (upper triangle,
right color bar) of the 35 representative Pareto-optimal sequences. (c)
Variations of selected sequence descriptors (degree of polymerization, $N$;
average hydrophobicity, $\bar{\lambda}$; hydropathy distribution, SHD;
absolute net charge per residue, $|\bar{q}|$; the total charge,
$\bar{q}_{+}+\bar{q}_{-}$; and charge distribution, SCD) along the Pareto
front. Dashed lines indicate the mean value of each feature across all phase-
separating sequences studied via active learning. (d) Average feature
differences between each Pareto-optimal sequence and its counterfactuals, Eq.
(2). Error bars represent standard errors of the mean, and dashed lines
indicate zero feature difference.

However, we find that neither composition nor descriptor-based similarity
alone guarantees similar macroscopic properties among Pareto-optimal
sequences. On the one hand, Pareto-optimal sequences can exhibit high
variability in sequence characteristics despite occasionally possessing
similar physical properties. This is particularly true in the transition
regime, which is typified by high intra-group dissimilarity; for instance, P13
and P17 differ significantly in terms of sequence descriptors relative to
sequences that are directly adjacent on the Pareto front. This variability in
sequence characteristics also suggests a higher degeneracy of near-Pareto-
optimal polypeptides in the transition regime compared to the
high-$B_{2}$/high-$D$ and low-$B_{2}$/low-$D$ regimes, where the
characteristics of Pareto-optimal sequences tend to be more consistent. On the
other hand, sequences with high overall composition and sequence-descriptor
similarity can exhibit vastly different macroscopic properties. For example,
extremal sequences P1 and P35 have greater sequence-descriptor similarity than
P1 does with several other sequences in the high-$B_{2}$/high-$D$ regime. Even
within the higher-degeneracy transition regime, P9 and P20 maintain overall
high composition and sequence-descriptor similarity despite a 2.4-fold change
in $B_{2}$ and a 1.9-fold change in $D$ between these points on the Pareto
front. Taken together, these observations suggest that there is substantial
flexibility in choosing combinations of sequence features to achieve near-
Pareto-optimal condensate properties.

We next consider whether specific sequence characteristics, rather than
overall sequence similarities, might control property changes across the
Pareto front. In this way, we find that a subset of sequence descriptors are
most relevant and that the relative importance of each descriptor depends on
the regime. In Fig. 4c, we highlight the trends along the Pareto front of key
sequence descriptors including sequence length, $N$; average sequence
hydropathy, $\bar{\lambda}$; sequence hydropathy decoration, SHD; the absolute
net charge per residue, $|\bar{q}|$; the fraction of charged residues (or the
total charge), $\bar{q}_{+}+\bar{q}_{-}$; and sequence charge decoration, SCD.
Data for other features that have previously been cited as relevant to IDP
condensates [9, 60, 61, 62] are shown in Fig. S9 in the Supplementary
Materials. The highlighted descriptors display an array of behaviors when
moving across the Pareto front. $N$ is roughly sigmoidal with $-B_{2}$,
starting at small $N$ in the high-$B_{2}$/high-$D$ regime (P1–P8) before
gradually but noisily increasing to plateau at large $N$ in the
low-$B_{2}$/low-$D$ regime (P27–P35). By contrast, $\bar{\lambda}$ and SHD
exhibit U-shaped behavior, first generally decreasing in the
high-$B_{2}$/high-$D$ regime and then increasing once entering the transition
region. Interestingly, the net charge, $|\bar{q}|$, of Pareto-optimal
sequences is always close to zero, but there are qualitative differences in
total charge, $\bar{q}_{+}+\bar{q}_{-}$, and SCD in different regimes. In
particular, $\bar{q}_{+}+\bar{q}_{-}$ starts at finite fractional charge at
high-$B_{2}$/high-$D$ and initially increases before diminishing to zero in
the transition and low-$B_{2}$/low-$D$ regimes; SCD displays similar behavior,
but with the opposite sign.

The collective observations from Fig. 4c suggest two important guiding
principles regarding the thermodynamics–dynamics tradeoff of protein
condensates. First, the trends in charge-related features indicate that
Pareto-optimal polypeptides are polyampholytic only in the
high-$B_{2}$/high-$D$ regime. Polypeptides in this regime also tend to be
shorter and less hydrophobic than other phase-separating IDPs generated during
active learning. Therefore, polyampholicity in combination with nonuniform
charge distribution and weak hydrophobicity is key to achieving phase
separation with short (and therefore fast-diffusing) chains. Second, and by
contrast, Pareto-optimal polypeptides are largely composed of uncharged,
hydrophobic amino acids in the low-$B_{2}$/low-$D$ regime. In particular,
sequences P27–P35 are predominantly composed of tyrosine, tryptophan, and
phenylalanine, which have the highest hydrophobic character in the CG model
and are also among the more massive amino acids. Greater hydrophobicity and
longer chain lengths therefore result in both stronger inter-chain attraction
and more sluggish diffusion in the low-$B_{2}$/low-$D$ regime. It is notable
that these insights are consistent with prior computational and experimental
studies of protein-condensate rheology [9, 63, 64]. Therefore, these
principles may provide insights into general property trends resulting from
sequence modifications, even though they stem from analyzing a specific set of
Pareto-optimal polypeptides.

### Pareto-optimality is characterized by regime-dependent sequence
determinants

We perform a counterfactual analysis [65, 66, 36] of the Pareto-optimal IDPs
to gain further insight into the sequence determinants of the heteropolymer
thermodynamics–dynamics tradeoff. In this context, a “counterfactual” is an
IDP sequence that bears strong sequence similarity to a Pareto-optimal
sequence but exhibits distinctly different physical properties (Fig. 4a). Such
counterfactuals can therefore be used to attribute Pareto-optimality to small
variations in IDP sequences. Here, we consider a polypeptide to be a
counterfactual of a specific sequence on the Pareto front if _(i)_ the cosine
similarity of their feature vectors is greater than 0.9 and _(ii)_ the
counterfactual candidate is within a dimensionless distance of 0.15–0.3 away
from the Pareto front; for this calculation, we use a Euclidean distance in
the in the standard-normalized $B_{2}$–$D$ plane (see Supplementary Materials
for details). Requiring that the dimensionless distance is in the range
0.15–0.3 ensures that the properties of counterfactual sequences are
statistically resolvable from those on the Pareto front while remaining in a
similar regime of physical property behavior. Because each Pareto-optimal
polypeptide may have several counterfactuals satisfying the aforementioned
criteria, we examine the average feature differences,

$\langle\Delta
x_{k}\rangle\equiv\frac{1}{n_{\mathcal{C}_{ij}}}\sum_{j=1}^{n_{\mathcal{C}_{ij}}}(x^{(\mathcal{P}_{i})}_{k}-x^{(\mathcal{C}_{ij})}_{k}),$
(2)

where $k$ denotes a particular feature of a feature vector $\vec{x}$,
$\mathcal{P}_{i}$ is the $i$th Pareto-optimal sequence with its value of
$x_{k}$ indicated as $x_{k}^{(\mathcal{P}_{i})}$, $\mathcal{C}_{ij}$ is the
$j$th counterfactual of the $i$th Pareto-optimal sequence with its value of
$x_{k}$ indicated as $x_{k}^{(\mathcal{C}_{ij})}$, and $n_{\mathcal{C}_{ij}}$
is the total number of counterfactuals of the $i$th Pareto-optimal sequence.
Plotting $\langle\Delta x_{k}\rangle$ for all features versus $-B_{2}/V_{0}$
(Fig. 4d and Fig. S10 in the Supplementary Materials) allows us to ascertain
whether specific sequence descriptors distinguish Pareto-optimal versus non-
Pareto-optimal sequences.

In the low-$B_{2}$/low-$D$ regime, we find that Pareto-optimal polypeptides
possess several distinctive features relative to their counterfactual
sequences. In particular, Pareto-optimal sequences within this regime are
generally longer, feature more hydrophobic residues, and have larger values of
the sequence hydropathy decoration parameter, SHD, than their counterfactuals.
Pareto-optimal sequences in this regime also tend to have fewer charged
residues and a lower net charge than their counterfactuals. Since the Pareto-
optimal sequences in this regime are also neutral (Fig. 4c), this implies that
the counterfactual sequences possess a net charge despite having near zero
values of the sequence charge decoration parameter, SCD (Fig. 4c,d). Our
counterfactual analysis therefore indicates that hydrophobic interactions are
the dominant determinants of Pareto-optimal sequences in the
low-$B_{2}$/low-$D$ regime.

However, outside the low-$B_{2}$/low-$D$ regime, differences between
counterfactuals and Pareto-optimal sequences are not as easily resolved. In
the high-$B_{2}$/high-$D$ regime, no single sequence descriptor differs
consistently between Pareto-optimal and counterfactual sequences. Instead,
variations in $B_{2}$ and $D$ manifest via subtle manipulations, such as
exchanging residues with comparable hydrophobicity but different masses (Fig.
S10 in the Supplementary Materials). Within the transition regime, there are
statistically significant differences between features of Pareto-optimal
polypeptides and counterfactuals, but the directions of such differences are
not necessarily preserved across the Pareto front. In particular, the
hydrophobicity differences, $\langle\Delta\bar{\lambda}\rangle$, and total
charge differences, $\langle\Delta(\bar{q}_{+}+\bar{q}_{-})\rangle$, fluctuate
between positive and negative values within this regime (P9–P20 in Fig. 4d).
This consistently enhanced variability across multiple features (see also Fig.
S10 in the Supplementary Materials) can be attributed to the greater
degeneracy of sequence space near the Pareto front in the transition regime,
as discussed in the previous section. Overall, this counterfactual analysis
indicates that the relationship between sequence features and Pareto
optimality is nuanced, since relatively subtle changes from counterfactuals
dictate behavior in the high-$B_{2}$/high-$D$ regime and myriad feature
combinations exist for achieving near-optimality in the transition regime.

## Discussion

Using an active-learning approach, we have investigated the extent to which
thermodynamic and dynamical properties are coupled in IDP condensates. Our
approach uniquely combines high-throughput simulations of a coarse-grained
model of IDPs, Bayesian active learning strategies for predicting physical
properties from IDP sequences, and numerical optimization techniques to
efficiently explore sequence space. Overall, this approach confirms that a
tradeoff between thermodynamic stability and internal condensate dynamics is a
general, albeit tunable, feature of protein-based condensates. Importantly,
this work provide a robust assessment of condensate dynamics via a global
analysis of IDP sequence space, rather than relying on systematic sequence
mutations alone.

The essential insight of our active-learning study is that heteromeric
polypeptides, by virtue of their diverse physicochemical properties, tend to
decouple the thermodynamic and dynamical properties of IDP condensates. In our
study, the correlation between proxies for thermodynamic stability and
condensed-phase dynamics is strongest in the case of homomeric polypeptides.
In particular, both the thermodynamic stability and the condensed-phase
dynamics of homomeric polypeptide condensates can be predicted with high
accuracy on the basis of the second-virial coefficient and the chain length
alone. However, this correlation between thermodynamic and dynamic properties
weakens when we consider heteromeric IDP sequences, indicating that we can
design sequences to tune the thermodynamic and dynamic properties of IDP
condensates independently. Consistent with this notion, we observe that the
second-virial coefficient is not always sufficient for predicting phase
separation of heteromeric IDP sequences. Ultimately, by identifying a
representative set of Pareto-optimal polypeptides—heteromeric sequences that
outline a limiting tradeoff boundary with respect to the second-virial and
diffusion coefficients—we find that the condensed-phase self-diffusion
coefficient can be substantially increased relative to homomeric polypeptides
with the same second-virial coefficient.

Our systematic approach also provides insight into specific sequence
descriptors that govern Pareto-optimal sequences. However, we find that these
sequence determinants are surprisingly nuanced and differ depending on the
position of a sequence along the thermodynamics–dynamics Pareto front. Whereas
the extremal regions (high-$B_{2}$/high-$D$ and low-$B_{2}$/low-$D$) display
consistent trends with respect to key sequence features, the transition regime
is characterized by high sequence variability, which likely reflects greater
degeneracy with respect to sequence space within this regime. Chain length,
hydrophobicity, and hydrophobic patterning are consistently important
determinants of Pareto-optimality. By contrast, total charge and charge
patterning descriptors are primarily relevant in the high-$B_{2}$/high-$D$
regime that is dominated by polyampholytes, and enhancements to charge
patterning tend to suppress polypeptide diffusion. Finally, our counterfactual
analysis indicates that relatively minor sequence perturbations can result in
large changes to both $B_{2}$ and $D$, particularly in the transition regime.
These findings suggest that different design principles are relevant in
different regimes of physicochemical properties, and they highlight the
importance of considering multiple sequence features when optimizing
polypeptide sequences to achieve specific physical properties.

We emphasize that our numerical results apply to single-component protein
condensates and rely on the accuracy of the underlying coarse-grained IDP
model, which lacks physicochemical features such as directional sidechain
interactions, hydrogen bonds, and secondary structure. These limitations of
the underlying coarse-grained model influence the Pareto-optimal sequences,
many of which look dissimilar from known naturally occurring IDPs, despite
having overall low sequence complexity. In this sense, our designed sequences
should be interpreted as polypeptide-like heteropolymers with hydrophobic,
electrostatic, and excluded-volume interactions as opposed to literal amino-
acid sequences. Coarse-grained models are also known to predict dynamical
properties with reduced fidelity compared to all-atom models with many more
degrees of freedom [67, 68, 69]. It is likely that a greater range of
variability of dynamical properties would be predicted by a model that
incorporates such features, especially if such modifications were to allow for
more accurate descriptions of low-density condensed phases [16]. It is also
likely that the range of variability of dynamic properties would be increased
by examining IDP sequences with lengths greater than the upper limit of 50
used in this study. Nonetheless, despite these considerations, we expect that
our findings regarding the existence of a thermodynamics–dynamics tradeoff, as
well as the relative importance of chain length, charge patterning, and
hydrophobicity along the Pareto front, are robust with respect to alternative
IDP models [29, 70, 71, 72, 73, 74] and actual polypeptide-like polymers.
Experiments and higher-resolution simulations [67, 75, 76, 77] will be needed
to test such predictions.

A major strength of our approach is its ability to explore IDP sequence space
efficiently and characterize relationships among disparate biophysical
properties. Crucially, the strategy of combining high-throughput simulations
and Bayesian-guided optimization is independent of any limitations of the
underlying coarse-grained model. Moreover, we have demonstrated that this
approach can efficiently converge the Pareto front for a pair of biophysical
properties that depend on a complex sequence of computations. We anticipate
that this framework could be applied to other combinations of biophysical
traits and may also be useful for studying multicomponent and multiphasic
systems [78, 79, 80]. Our work therefore establishes a promising approach for
engineering custom biomolecular condensates with tunable thermodynamic and
dynamic properties.

## Materials and Methods

### Model of intrinsically disordered proteins

The coarse-grained (CG) model developed by Regy et al. [26] represents IDP
sequences, comprising 20 amino-acid types, in implicit solvent. The force
field consists of a harmonic bonded potential as well as short-range van der
Waals (vdW) and long-range electrostatic (el) nonbonded interactions,

$U_{\text{tot}}=\sum_{i}k_{\text{b}}(r_{i,i+1}-b_{0})^{2}+\sum_{i,j}\phi^{\text{vdw}}(r_{ij})+\sum_{i,j}\phi^{\text{el}}(r_{ij}),$
(3)

where $r_{i,i+1}$ is the distance between bonded residues $i$ and $i+1$, the
force constant is
$k_{\text{b}}=$10\text{\,}\mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{/}\mathrm{(}\mathrm{m}\mathrm{o}\mathrm{l}\cdot\AA^{2}\mathrm{)}$$,
and the equilibrium bond length is $b_{0}=3.82$ Å. The vdW interaction takes
the Ashbaugh–Hatch functional form,

$\phi^{\text{vdW}}(r_{ij})=\begin{cases}\phi^{\text{LJ}}(r_{ij})+(1-\lambda_{ij})\epsilon&\quad
r_{ij}\leq 2^{1/6}\sigma_{ij}\\\ \lambda_{ij}\phi^{\text{LJ}}(r_{ij})&\quad
r_{ij}>2^{1/6}\sigma_{ij}\\\ \end{cases}$ (4)

where $r_{ij}$ is the distance between nonbonded residues $i$ and $j$, and
$\phi^{\text{LJ}}(r_{ij})=4\epsilon[(\sigma_{ij}/r_{ij})^{12}-(\sigma_{ij}/r_{ij})^{6}]$
is the Lennard–Jones potential with $\epsilon=0.2$ kcal/mol. The interaction
strength, $\lambda_{ij}=(\lambda_{i}+\lambda_{j})/2$, and distance,
$\sigma_{ij}=(\sigma_{i}+\sigma_{j})/2$, parameters are determined from the
hydropathy scaling factors, $\lambda_{i}$, and vdW diameters, $\sigma_{i}$, of
the residues, respectively [26]. Lastly, screened electrostatic interactions
take the form

$\phi^{\text{el}}(r_{ij})=\frac{q_{i}q_{j}}{4\pi Dr_{ij}}e^{-\kappa r_{ij}},$
(5)

where $q_{i}$ is the charge of residue $i$, $D=80$ is the dielectric constant
of the solvent, and $\kappa=10$ Å is the Debye screening length.

### Physical property calculations

Second-virial coefficient. Second-virial coefficients are determined by
calculating the potential of mean force, $u(r)$, as a function of the center-
of-mass (COM) distance, $r$, between an isolated pair of chains. These
calculations are carried out using the adaptive biasing force (ABF) method
[44] implemented in the COLVARS package in LAMMPS [81]. A Langevin thermostat
is applied to maintain a constant temperature of $300\text{\,}\mathrm{K}$. Two
molecules are placed in a cubic simulation box with dimensions $300\times
300\times 300$ Å. ABF simulations are performed as a function of $r$, which
varies from $1\text{\,}~{}\mathrm{\SIUnitSymbolAngstrom}$ to
$100\text{\,}~{}\mathrm{\SIUnitSymbolAngstrom}$. Force statistics are stored
in bins of width $0.2\text{\,}\mathrm{\SIUnitSymbolAngstrom}$. The biasing
force is applied after 1000 samples are collected in each bin, after which a
final production run is performed for $5\text{\,}\mathrm{\SIUnitSymbolMicro
s}$ with a timestep of $10\text{\,}\mathrm{fs}$. The second-virial coefficient
is then obtained from Eq. (1). Thirty independent simulations are performed
for each free energy calculation to determine the statistical error.

Condensed-phase density. Considering the high computational costs of slab
simulations, we developed an alternative method to determine whether a CG
sequence undergoes phase separation at 300 K. To this end, we look for
evidence of a van der Waal’s loop in the equation of state (EOS) of a small
system. Specifically, we compute the pressure using canonical-ensemble
simulations of 100 chains with periodic boundary conditions at densities
ranging from
$0.2\text{\,}\mathrm{g}\mathrm{/}\mathrm{m}\mathrm{l}$–$1.2\text{\,}\mathrm{g}\mathrm{/}\mathrm{m}\mathrm{l}$.
A Langevin thermostat is applied to maintain a constant temperature of
$300\text{\,}\mathrm{K}$. Simulations are carried out for
$100\text{\,}\mathrm{ns}$ with a timestep of $10\text{\,}\mathrm{fs}$. We
assume a near-zero coexistence pressure to minimize statistical and fitting
errors, so that phase separation can be identified by the existence of a
negative pressure at a finite density. This approximation is justified by the
fact that the coexisting dilute phase is characterized by a low polymer
density and is nearly ideal, implying that the coexistence pressure is also
close to zero. In cases where negative pressures are observed, we fit a cubic
spline to the EOS to determine the condensed-phase density, $\rho_{\text{c}}$,
defined as the highest density at which the pressure crosses zero. The
pressure values are bootstrapped with replacement 50 times to determine the
average condensate densities and their standard errors.

Self-diffusion coefficient. To determine the self-diffusion coefficient of a
chain within the condensed phase, canonical-ensemble simulations of 100 chains
are performed at the condensed-phase density, $\rho_{\text{c}}$, with periodic
boundary conditions. The long-time behavior of the mean-squared displacement
is then used to compute the self-diffusion coefficient,

$D=\lim_{\Delta t\to\infty}{\frac{1}{6\Delta t}\langle|\vec{r}_{t+\Delta
t}-\vec{r}_{t}|^{2}\rangle},$ (6)

where $\vec{r}$ is the position of the COM of a tagged chain. The temperature
is maintained at $300\text{\,}\mathrm{K}$ using a Langevin thermostat with a
damping frequency of $1\text{\,}\mathrm{p}\mathrm{s}^{-1}$. Simulations are
carried out for $100\text{\,}\mathrm{ns}$ with a timestep of
$10\text{\,}\mathrm{fs}$. For each diffusion-coefficient calculation, the
statistical error is quantified using 30 independent simulations.

### Design of novel intrinsically disordered proteins

Overview of framework. We deploy an active learning [32, 55] framework based
on Bayesian optimization [34] to identify IDP sequences that define a Pareto
front of the properties $B_{2}$ and $D$. The active learning framework
utilizes coarse-grained simulations to generate data for $B_{2}$ as well as
pressure-density data and $D$, as appropriate. These data are used to train
three different machine-learning (ML) models: two Gaussian process regressors
[82] (GPR) and one random forest (RF) binary classifier [83]. These ML models
are used to make surrogate predictions during the optimization step of active
learning, during which new IDP sequences are selected for simulation. Sequence
optimization is carried out using a genetic algorithm, as described below.
Following optimization, simulations are performed to generate new data and
begin the next iteration of active learning.

Machine-learning models. GPR models are trained to predict $B_{2}$ and $D$,
while a RF model is trained to predict whether a given IDP sequence will
undergo phase separation. For a given sequence $i$, the input for all ML
models is $\vec{x}^{(i)}$, size-explicit augmented fingerprint as a feature
vector [58]. In particular, the feature vector is 30-dimensional, consisting
of 20 features reflecting the composition of the amino acids and ten sequence-
level chain descriptors. The elements of the feature vector are all normalized
via linear transformation methods to be unit-scale. Additional details
regarding the featurization, training, and hyperparameter optimization for all
ML models are provided in the Supplementary Materials. Although model accuracy
is not a primary objective, the three trained ML models show high accuracy in
predicting the aforementioned properties (see Fig. S6 in the Supplementary
Materials) at the conclusion of active learning. The coefficients of
determination $R^{2}$ of the two GPR models are as high as $0.975\pm 0.015$
and $0.970\pm 0.015$, and the accuracy of the RF classifier is 0.96.

Sequence optimization. Different fitness functions were used for sequence
design by a genetic algorithm depending on the active-learning iteration.
Denoting the iteration number by the index $k$, the fitness function
$f(\vec{x}^{(i)})$ is

$f(\vec{x}^{(i)})\\!=\\!\begin{cases}-B_{2}(\vec{x}^{(i)})&\text{if }k=0\\\
\text{EHVI}(\vec{x}^{(i)})\times\text{RF}(\vec{x}^{(i)})\times\alpha(\vec{x}^{(i)})\\!\\!\\!\\!&\text{if
}0<k<8\\\ -\tilde{B}_{2}(\vec{x}^{(i)})+\tilde{D}(\vec{x}^{(i)})&\text{if
}k=8,\end{cases}$ (7)

where EHVI($\vec{x}^{(i)}$) is the expected hypervolume improvement
acquisition function, RF($\vec{x}^{(i)}$) is the output label of the random
forest classifier (evaluating to 1 if phase separation is predicted and to 0
otherwise), $\alpha(\vec{x}^{(i)})$ is a scaling function that biases against
generation of sequences that have high similarity to previously simulated or
proposed sequences, and $\tilde{A}$ indicates a standard-normalized
transformation of the property $A$ based on the data acquired up to the given
iteration. For each iteration, 96 unique sequences are generated. In Eq. (7),
we compute the EHVI [56] for a given sequence $\vec{x}^{*}$,

$\text{EHVI}(\vec{x}^{*},\mathcal{P}^{\prime})=\sum_{i=1}^{n+1}\Psi(\tilde{D}^{(i)},\tilde{D}^{(i)},\tilde{D}^{*},\sigma^{*}_{\tilde{D}})\\\
\times\biggl{[}(\tilde{B}_{2}^{(i)}-\tilde{B}_{2}^{(i-1)})\Phi\Bigl{(}\frac{\tilde{B}_{2}^{*}-\tilde{B}_{2}^{(i)}}{\sigma^{*}_{\tilde{B}_{2}}}\Bigr{)}\\\
\qquad+\Bigl{(}\Psi(-\tilde{B}_{2}^{(i-1)},-\tilde{B}_{2}^{(i-1)}-\tilde{B}_{2}^{*},\sigma^{*}_{\tilde{B}_{2}})\\\
\qquad-\Psi(-\tilde{B}_{2}^{(i-1)},-\tilde{B}_{2}^{(i)},-\tilde{B}_{2}^{*},\sigma^{*}_{\tilde{B}_{2}})\Bigr{)}\biggr{]},$
(8)

where $\tilde{B}_{2}^{*}$ and $\tilde{D}^{*}$ are predicted values of the
second-virial and diffusion coefficients, respectively, obtained from the GPR
models; $\sigma^{*}_{\tilde{B}_{2}}$ and $\sigma^{*}_{\tilde{D}}$ are
corresponding uncertainty estimates from the GPR models;
$\mathcal{P}^{\prime}=\\{\vec{x}^{(0)},\vec{x}^{(1)},\ldots,\vec{x}^{(n)},\vec{x}^{(n+1)}\\}$
represents the current $n$-point approximation of the Pareto front, augmented
with additional reference points $\vec{x}^{(0)}=(X,-\infty)^{\top}$ and
$\vec{x}^{(0)}=(-\infty,X)^{\top}$ to facilitate the hypervolume calculation;
$\Phi(s)$ denotes the cumulative probability distribution function for a
standard normal distribution of $s$; and $\Psi$ is a function defined as

$\Psi(a,b,\mu,\sigma)=\int_{-\infty}^{b}\frac{(a-z)}{\sigma}\frac{1}{\sqrt{2\pi}}e^{-\frac{(z-\mu)^{2}}{2\sigma^{2}}}dz.$
(9)

The similarity penalty $\alpha(\vec{x^{(j)}})$ is defined to be

$\alpha(\vec{x}^{(j)})=\frac{j-1}{\sum_{k=1}^{j-1}\left\\{\frac{1}{2}[1-\vec{x}^{(j)}\cdot\vec{x}^{(k)}/(|\vec{x}^{(j)}||\vec{x}^{(k)}|)]\right\\}^{-1}},$
(10)

where $j=2,3,\ldots,96$ are indices of candidate sequences generated in a
given iteration. When a sequence is highly similar to one or more previously
generated sequences, $\alpha$ tends to zero. For $j=1$, we take $\alpha=1$
(i.e., no penalty is applied).

The active learning process is seeded with data obtained from simulations of
1,266 disordered sequences reported in DisProt [84]. We select and simulate
all sequences with sequence lengths between 20 and 50 residues (inclusive),
taking care to remove duplicated samples. Since none of these sequences are
determined to phase separate, we generate new sequences by maximizing a
fitness function that involves only $B_{2}$ in iteration 0. The resulting data
enable the training of ML models for $B_{2}$, $D$, and phase separation. Then
in subsequent iterations, we maximize a fitness function based principally on
the expected hypervolume improvement (EHVI) to identify polypeptides that
outline a Pareto front of $-B_{2}$ versus $D$. To confirm the convergence of
the Pareto front, a final exploitation-only round of optimization is performed
by maximizing a simple linear function of the standard-normalized properties
$-\tilde{B}_{2}$ and $\tilde{D}$.

For each iteration, 96 “child” sequences are successively generated that aim
to maximize the relevant fitness function in Eq. (7). To facilitate
convergence towards selecting each child sequence, 96 independent trials of
sequence optimization are executed in parallel, and only the sequence with the
best fitness score is chosen for explicit simulation; overall, this means that
$96\times 96$ sequence optimizations are performed in each iteration.

In iteration 0, the initial “parent” sequences are DisProt sequences with
negative $B_{2}$. In all subsequent iterations, the parent sequences are
Pareto-optimal sequences identified from the previous iteration. Candidate
child sequences are derived from the parent sequences via a combination of
crossover and mutation moves. In addition, “deletion” and “growth” moves are
performed, which either remove a portion of the current sequence or add a
portion of the current sequence within itself. These moves are executed
_independently_ in the order of crossover, mutation, deletion and growth with
probabilities of 0.5, 0.8, 0.2, and 0.5, respectively. These moves are
performed for 100 steps, which is sufficient for the genetic algorithm to
converge to a local maximum of the fitness function (Fig. S11 in the
Supplementary Materials). Additional details can be found in the Supplementary
Materials.

## Acknowledgments

The authors acknowledge assistance from Princeton Research Computing at
Princeton University, which is a consortium led by the Princeton Institute for
Computational Science and Engineering (PICSciE) and the Princeton University
Office of Information Technology’s Research Computing.

#### Funding:

This work was a grant from the Princeton Biomolecular Condensate Program to
M.A.W. and W.M.J.

#### Author contributions:

Conceptualization: M.A.W., W.M.J.; Methodology: Y.A., M.A.W., W.M.J.;
Investigation: Y.A., M.A.W., W.M.J.; Visualization: Y.A., M.A.W., W.M.J.;
Supervision: M.A.W., W.M.J.; Writing (original draft): Y.A., M.A.W., W.M.J.;
Writing (review and editing): Y.A., M.A.W., W.M.J.

#### Competing interests:

The authors declare that they have no competing interests.

#### Data and materials availability:

All data needed to evaluate the conclusions of the paper are available in the
main text or in the supplementary materials. Full datasets compiled from this
work are provided in Ref. [86]. Demonstrative Jupyter notebooks for essential
calculations are also provided at https://github.com/webbtheosim/condensate-
public.

## References

  * [1] C. P. Brangwynne, C. R. Eckmann, D. S. Courson, A. Rybarska, C. Hoege, J. Gharakhani, F. Jülicher, A. A. Hyman, Germline p granules are liquid droplets that localize by controlled dissolution/condensation, Science 324, 1729–1732 (2009).
  * [2] E. Gomes, J. Shorter, The molecular language of membraneless organelles, J. Biol. Chem. 294, 7115–7127 (2019).
  * [3] T. Mittag, R. V. Pappu, A conceptual framework for understanding phase separation and addressing open questions and challenges, Mol. Cell 82, 2201–2214 (2022).
  * [4] J.-Y. Kang, Z. Wen, D. Pan, Y. Zhang, Q. Li, A. Zhong, X. Yu, Y.-C. Wu, Y. Chen, X. Zhang, et al., LLPS of fxr1 drives spermiogenesis by activating translation of stored mrnas, Science 377, 1-13 (2022).
  * [5] G. He, T. GrandPre, H. Wilson, Y. Zhang, M. C. Jonikas, N. S. Wingreen, Q. Wang, Phase-separating pyrenoid proteins form complexes in the dilute phase, Commun. Biol. 6, 1-10 (2023).
  * [6] Y. Shin, C. P. Brangwynne, Liquid phase condensation in cell physiology and disease, Science 357, 1-11 (2017).
  * [7] N.-N. Deng, W. T. Huck, Microfluidic formation of monodisperse coacervate organelles in liposomes, Angew. Chem. 129, 9868–9872 (2017).
  * [8] T. Kojima, S. Takayama, Membraneless compartmentalization facilitates enzymatic cascade reactions and reduces substrate inhibition, ACS Appl. Mater. & Interfaces 10, 32782–32791 (2018).
  * [9] M. K. Hazra, Y. Levy, Biophysics of phase separation of disordered proteins is governed by balance between short-and long-range interactions, J. Phys. Chem. B 125, 2202–2211 (2021).
  * [10] Y. G. Zhao, H. Zhang, Phase separation in membrane biology: the interplay between membrane-bound organelles and membraneless condensates, Dev. Cell 55, 30–44 (2020).
  * [11] I. Alshareedah, M. M. Moosa, M. Pham, D. A. Potoyan, P. R. Banerjee, Programmable viscoelasticity in protein-rna condensates with disordered sticker-spacer polypeptides, Nat. Commun. 12, 1–14 (2021).
  * [12] L. Cai, G. G. Wang, Through the lens of phase separation: intrinsically unstructured protein and chromatin looping, Nucleus 14, 1-15 (2023).
  * [13] K. L. Price, M. Presler, C. M. Uyehara, D. C. Shakes, The intrinsically disordered protein spe-18 promotes localized assembly of msp in caenorhabditis elegans spermatocytes, Development 148, 1-13 (2021).
  * [14] L. Mediani, F. Antoniani, V. Galli, J. Vinet, A. D. Carra, I. Bigi, V. Tripathy, T. Tiago, M. Cimino, G. Leo, et al., Hsp90-mediated regulation of dyrk3 couples stress granule disassembly and growth via mtorc1 signaling, EMBO Rep. 22, 1-22 (2021).
  * [15] A. Agarwal, S. K. Rai, A. Avni, S. Mukhopadhyay, An intrinsically disordered pathological prion variant y145stop converts into self-seeding amyloids via liquid–liquid phase separation, Proc. Natl. Acad. Sci. 118, 1-10 (2021).
  * [16] M.-T. Wei, S. Elbaum-Garfinkle, A. S. Holehouse, C. C.-H. Chen, M. Feric, C. B. Arnold, R. D. Priestley, R. V. Pappu, C. P. Brangwynne, Phase behaviour of disordered proteins underlying low density and high permeability of liquid organelles, Nat. Chem. 9, 1118–1125 (2017).
  * [17] M. Feric, N. Vaidya, T. S. Harmon, D. M. Mitrea, L. Zhu, T. M. Richardson, R. W. Kriwacki, R. V. Pappu, C. P. Brangwynne, Coexisting liquid phases underlie nucleolar subcompartments, Cell 165, 1686–1697 (2016).
  * [18] G. L. Dignon, W. Zheng, R. B. Best, Y. C. Kim, J. Mittal, Relation between single-molecule properties and phase behavior of intrinsically disordered proteins, Proc. Natl. Acad. Sci. 115, 9929–9934 (2018).
  * [19] S. Das, A. N. Amin, Y.-H. Lin, H. S. Chan, Coarse-grained residue-based models of disordered protein condensates: utility and limitations of simple charge pattern parameters, Phys. Chem. Chem. Phys. 20, 28558–28574 (2018).
  * [20] S. Das, Y.-H. Lin, R. M. Vernon, J. D. Forman-Kay, H. S. Chan, Comparative roles of charge, $\pi$, and hydrophobic interactions in sequence-dependent phase separation of intrinsically disordered proteins, Proc. Natl. Acad. Sci. 117, 28795–28805 (2020).
  * [21] H. Y. J. Fung, M. Birol, E. Rhoades, Idps in macromolecular complexes: the roles of multivalent interactions in diverse assemblies, Curr. Opin. Struct. Biol. 49, 36–43 (2018).
  * [22] S. E. Harding, P. Johnson, The concentration-dependence of macromolecular parameters, Biochem. J. 231, 543–547 (1985).
  * [23] R. Laghmach, I. Alshareedah, M. Pham, M. Raju, P. R. Banerjee, D. A. Potoyan, RNA chain length and stoichiometry govern surface tension and stability of protein-rna condensates, Iscience 25, 1-15 (2022).
  * [24] G. M. Wadsworth, W. J. Zahurancik, X. Zeng, P. Pullara, L. B. Lai, V. Sidharthan, R. V. Pappu, V. Gopalan, P. R. Banerjee, RNAs undergo phase transitions with lower critical solution temperatures, bioRxiv (2022).
  * [25] C. J. Oldfield, A. K. Dunker, Intrinsically disordered proteins and intrinsically disordered protein regions, Annu. Rev. Biochem. 83, 553–584 (2014).
  * [26] R. M. Regy, J. Thompson, Y. C. Kim, J. Mittal, Improved coarse-grained model for studying sequence dependent phase separation of disordered proteins, Protein Sci. 30, 1371–1379 (2021).
  * [27] G. L. Dignon, W. Zheng, J. Mittal, Simulation methods for liquid–liquid phase separation of disordered proteins, Curr Opin Chem 23, 92–98 (2019).
  * [28] G. L. Dignon, W. Zheng, Y. C. Kim, R. B. Best, J. Mittal, Sequence determinants of protein phase behavior from a coarse-grained model, PLoS Comput. Biol. 14, 1-23 (2018).
  * [29] J. A. Joseph, A. Reinhardt, A. Aguirre, P. Y. Chew, K. O. Russell, J. R. Espinosa, A. Garaizar, R. Collepardo-Guevara, Physics-driven coarse-grained model for biomolecular phase separation with near-quantitative accuracy, Nat. Comput. Sci 1, 732–743 (2021).
  * [30] M. A. Webb, N. E. Jackson, P. S. Gil, J. J. de Pablo, Targeted sequence design within the coarse-grained polymer genome, Sci. Adv. 6, 1-10 (2020).
  * [31] A. J. Gormley, M. A. Webb, Machine learning in combinatorial polym. chem., Nat. Rev. Mater. 6, 642-644 (2021).
  * [32] J. S. Smith, B. Nebgen, N. Lubbers, O. Isayev, A. E. Roitberg, Less is more: Sampling chemical space with active learning, J. Chem. Phys. 148, 1-10 (2018).
  * [33] C. Kim, A. Chandrasekaran, A. Jha, R. Ramprasad, Active-learning and materials design: the example of high glass transition temperature polymers, MRS Communications 9, 860–866 (2019).
  * [34] B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, N. De Freitas, Taking the human out of the loop: A review of bayesian optimization, Proceedings of the IEEE 104, 148–175 (2015).
  * [35] G. P. Wellawatte, A. Seshadri, A. D. White, Model agnostic generation of counterfactual explanations for molecules, Chem. Sci. 13, 3697–3705 (2022).
  * [36] G. P. Wellawatte, H. A. Gandhi, A. Seshadri, A. D. White, A perspective on explanations of molecular prediction models, J. Chem. Theory Comput. 19, 2149–2160 (2023).
  * [37] M. K. Hazra, Y. Levy, Affinity of disordered protein complexes is modulated by entropy–energy reinforcement, Proc. Natl. Acad. Sci. 119, 1-11 (2022).
  * [38] D. A. McQuarrie, Statistical Mechanics (Harper Collins, New York, 1976).
  * [39] G. Vliegenthart, H. N. Lekkerkerker, Predicting the gas–liquid critical point from the second virial coefficient, J. Chem. Phys. 112, 5364–5369 (2000).
  * [40] R. Tuinier, G. A. Vliegenthart, H. N. Lekkerkerker, Depletion interaction between spheres immersed in a solution of ideal polymer chains, J. Chem. Phys. 113, 10768–10775 (2000).
  * [41] M. Rubinstein, R. H. Colby, et al., Polymer physics, vol. 23 (Oxford university press New York, 2003).
  * [42] S. Rekhi, D. S. Devarajan, M. P. Howard, Y. C. Kim, A. Nikoubashman, J. Mittal, Role of strong localized vs weak distributed interactions in disordered protein phase separation, J. Phys. Chem. B. 127, 3829–3838 (2023).
  * [43] A. Z. Panagiotopoulos, Phase separation and aggregation in multiblock chains, The Journal of Chemical Physics 158, 1–9 (2023).
  * [44] J. Comer, J. C. Gumbart, J. Hénin, T. Lelièvre, A. Pohorille, C. Chipot, The adaptive biasing force method: Everything you always wanted to know but were afraid to ask, J. Phys. Chem. B 119, 1129–1151 (2015).
  * [45] K. Binder, B. J. Block, P. Virnau, A. Tröster, Beyond the van der waals loop: What can be learned from simulating lennard-jones fluids inside the region of phase coexistence, Am. J. Phys 80, 1099–1109 (2012).
  * [46] A. Hatos, B. Hajdu-Soltész, A. M. Monzon, N. Palopoli, L. Álvarez, B. Aykac-Fas, C. Bassot, G. I. Benítez, M. Bevilacqua, A. Chasapi, L. Chemes, N. E. Davey, R. Davidović, A. K. Dunker, A. Elofsson, J. Gobeill, N. S. G. Foutel, G. Sudha, M. Guharoy, T. Horvath, V. Iglesias, A. V. Kajava, O. P. Kovacs, J. Lamb, M. Lambrughi, T. Lazar, J. Y. Leclercq, E. Leonardi, S. Macedo-Ribeiro, M. Macossay-Castillo, E. Maiani, J. A. Manso, C. Marino-Buslje, E. Martínez-Pérez, B. Mészáros, I. Mičetić, G. Minervini, N. Murvai, M. Necci, C. A. Ouzounis, M. Pajkos, L. Paladin, R. Pancsa, E. Papaleo, G. Parisi, E. Pasche, P. J. B. Pereira, V. J. Promponas, J. Pujols, F. Quaglia, P. Ruch, M. Salvatore, E. Schad, B. Szabo, T. Szaniszló, S. Tamana, A. Tantos, N. Veljkovic, S. Ventura, W. Vranken, Z. Dosztányi, P. Tompa, S. C. E. Tosatto, D. Piovesan, DisProt: intrinsic protein disorder annotation in 2020, Nucleic Acids Res. 48, D269–D276 (2019).
  * [47] J.-F. Yu, Z. Cao, Y. Yang, C.-L. Wang, Z.-D. Su, Y.-W. Zhao, J.-H. Wang, Y. Zhou, Natural protein sequences are more intrinsically disordered than random sequences, Cellular and Molecular Life Sciences 73, 2949–2957 (2016).
  * [48] Q. Li, X. Peng, Y. Li, W. Tang, J. Zhu, J. Huang, Y. Qi, Z. Zhang, Llpsdb: a database of proteins undergoing liquid–liquid phase separation in vitro, Nucleic Acids Res. 48, D320–D327 (2020).
  * [49] K. Shmilovich, R. A. Mansbach, H. Sidky, O. E. Dunne, S. S. Panda, J. D. Tovar, A. L. Ferguson, Discovery of self-assembling pi-conjugated peptides by active learning-directed coarse-grained mol. simul., J. Phys. Chem. B 124, 3873–3891 (2020).
  * [50] M. J. Tamasi, R. A. Patel, C. H. Borca, S. Kosuri, H. Mugnier, R. Upadhya, N. S. Murthy, M. A. Webb, A. J. Gormley, Machine learning on a robotic platform for the design of polymer–protein hybrids, Adv. Mater. 34, 1-12 (2022).
  * [51] S. Kosuri, C. H. Borca, H. Mugnier, M. Tamasi, R. A. Patel, I. Perez, S. Kumar, Z. Finkel, R. Schloss, L. Cai, M. L. Yarmush, M. A. Webb, A. J. Gormley, Machine-assisted discovery of chondroitinase ABC complexes toward sustained neural regeneration, Adv. Healthc. Mater. 11, 1-12 (2022).
  * [52] K. M. Jablonka, G. M. Jothiappan, S. Wang, B. Smit, B. Yoo, Bias free multiobjective active learning for materials design and discovery, Nat. Commun. 12, 1-10 (2021).
  * [53] A. Tran, J. Tranchida, T. Wildey, A. P. Thompson, Multi-fidelity machine-learning with uncertainty quantification and bayesian optimization for materials design: Application to ternary random alloys, J. Chem. Phys. 153, 1-8 (2020).
  * [54] M. Emmerich, J.-w. Klinkenberg, The computation of the expected improvement in dominated hypervolume of pareto front approximations, Rapport technique, Leiden University 34, 7–3 (2008).
  * [55] M. Zuluaga, A. Krause, M. Püschel, $\epsilon$-pal: An active learning approach to the multi-objective optimization problem, J. Mach. Learn. Res. 17, 3619–3650 (2016).
  * [56] K. Yang, M. Emmerich, A. Deutz, T. Bäck, Multi-objective bayesian global optimization using expected hypervolume improvement gradient, Swarm Evol. Comput. 44, 945–956 (2019).
  * [57] B. P. MacLeod, F. G. Parlane, C. C. Rupnow, K. E. Dettelbach, M. S. Elliott, T. D. Morrissey, T. H. Haley, O. Proskurin, M. B. Rooney, N. Taherimakhsousi, et al., A self-driving laboratory advances the pareto front for material properties, Nat. Commun. 13, 1-10 (2022).
  * [58] R. A. Patel, C. H. Borca, M. A. Webb, Featurization strategies for polymer sequence or composition design by machine learning, Mol. Syst. Des. Eng. 7, 661-676 (2022).
  * [59] R. A. Patel, M. A. Webb, Data-driven design of polymer-based biomaterials: High-throughput simulation, experimentation, and machine learning, ACS Appl. Bio Mater. (2023).
  * [60] W. Zheng, G. Dignon, M. Brown, Y. C. Kim, J. Mittal, Hydropathy patterning complements charge patterning to describe conformational preferences of disordered proteins, J. Phys. Chem. Lett. 11, 3408–3415 (2020).
  * [61] A. Statt, H. Casademunt, C. P. Brangwynne, A. Z. Panagiotopoulos, Model for disordered proteins with strongly sequence-dependent liquid phase behavior, J. Chem. Phys. 152, 1-10 (2020).
  * [62] A. H. Mao, S. L. Crick, A. Vitalis, C. L. Chicoine, R. V. Pappu, Net charge per residue modulates conformational ensembles of intrinsically disordered proteins, Proc. Natl. Acad. Sci. 107, 8183–8188 (2010).
  * [63] I. Alshareedah, et al., Sequence-encoded grammars determine material properties and physical aging of protein condensates, bioRxiv (2023).
  * [64] S. Rekhi, C. G. Garcia, M. Barai, A. Rizuan, B. S. Schuster, K. L. Kiick, J. Mittal, Expanding the molecular language of protein liquid-liquid phase separation, bioRxiv (2023).
  * [65] S. Wachter, B. Mittelstadt, C. Russell, Counterfactual explanations without opening the black box: Automated decisions and the gdpr, Harv. JL & Tech. 31, 841-887 (2017).
  * [66] S. Dandl, C. Molnar, M. Binder, B. Bischl, Multi-objective counterfactual explanations, International Conference on Parallel Problem Solving from Nature 1, 448–469 (2020).
  * [67] N. Salvi, A. Abyzov, M. Blackledge, Solvent-dependent segmental dynamics in intrinsically disordered proteins, Sci. Adv. 5, 1-10 (2019).
  * [68] J. F. Rudzinski, Recent progress towards chemically-specific coarse-grained simulation models with consistent dynamical properties, Computation 7, 1-28 (2019).
  * [69] S. Dhamankar, M. A. Webb, Chemically specific coarse-graining of polymers: Methods and prospects, J. of Poly. Sci. 59, 2613–2643 (2021).
  * [70] A. Garaizar, J. R. Espinosa, Salt dependent phase behavior of intrinsically disordered proteins from a coarse-grained model with explicit water and ions, J. Chem. Phys. 155, 1-13 (2021).
  * [71] A. P. Latham, B. Zhang, Maximum entropy optimized force field for intrinsically disordered proteins, J. Chem. Theory Comput. 16, 773–781 (2019).
  * [72] Z. Benayad, S. von Bülow, L. S. Stelzl, G. Hummer, Simulation of FUS protein condensates with an adapted coarse-grained model, J. Chem. Theory Comput. 17, 525–537 (2020).
  * [73] G. Tesei, T. K. Schulze, R. Crehuet, K. Lindorff-Larsen, Accurate model of liquid–liquid phase behavior of intrinsically disordered proteins from optimization of single-chain properties, Proc. Natl. Acad. Sci. U.S.A. 118, 1-10 (2021).
  * [74] G. Tesei, K. Lindorff-Larsen, Improved predictions of phase behaviour of intrinsically disordered proteins by tuning the interaction range, Open. Res. Europe 2, 94 (2022).
  * [75] W. Zheng, G. L. Dignon, N. Jovic, X. Xu, R. M. Regy, N. L. Fawzi, Y. C. Kim, R. B. Best, J. Mittal, Molecular details of protein condensates probed by microsecond long atomistic simulations, J. Phys. Chem. B 124, 11671–11679 (2020).
  * [76] N. Galvanetto, M. T. Ivanović, A. Chowdhury, A. Sottini, M. Nüesch, D. Nettels, R. Best, B. Schuler, Ultrafast molecular dynamics observed within a dense protein condensate, bioRxiv (2022).
  * [77] L. A. Gruijs da Silva, F. Simonetti, S. Hutten, H. Riemenschneider, E. L. Sternburg, L. M. Pietrek, J. Gebel, V. Dötsch, D. Edbauer, G. Hummer, et al., Disease-linked TDP-43 hyperphosphorylation suppresses TDP-43 condensation and aggregation, EMBO J. 41, 1-23 (2022).
  * [78] W. M. Jacobs, Theory and simulation of multiphase coexistence in biomolecular mixtures, J. Chem. Theory Comput. 19, 3429–3445 (2023).
  * [79] S. Do, C. Lee, T. Lee, D.-N. Kim, Y. Shin, Engineering dna-based synthetic condensates with programmable material properties, compositions, and functionalities, Sci. Adv. 8, 1-14 (2022).
  * [80] P. Y. Chew, J. A. Joseph, R. Collepardo-Guevara, A. Reinhardt, Physical determinants of multiphase organisation in multi-component protein/rna condensates, Biophys. J. 122, 295a–296a (2023).
  * [81] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, et al., LAMMPS-a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comput. Phys. Commun. 271, 1-34 (2022).
  * [82] C. E. Rasmussen, C. K. Williams, et al., Gaussian processes for machine learning, vol. 1 (Springer, 2006).
  * [83] G. Biau, E. Scornet, A random forest guided tour, Test 25, 197–227 (2016).
  * [84] A. Hatos, B. Hajdu-Soltész, A. M. Monzon, N. Palopoli, L. Álvarez, B. Aykac-Fas, C. Bassot, G. I. Benítez, M. Bevilacqua, A. Chasapi, et al., Disprot: intrinsic protein disorder annotation in 2020, Nucleic Acids Res. 48, D269–D276 (2020).
  * [85] Y. Zhang, A. Otani, E. J. Maginn, Reliable viscosity calculation from equilibrium molecular dynamics simulations: A time decomposition method, J. Chem. Theory Comput. 11, 3537–3546 (2015).
  * [86] M. Webb, W. Jacobs, Y. An, W. Oliver, Thermodynamic and dynamics data for coarse-grained intrinsically disordered proteins generated by active learning. DataSpace. https://doi.org/10.34770/6tnm-7b56 (2023).