Theory of droplet ripening in stiffness gradients

Liquid-liquid phase separation is an important mechanism for compartmentalizing the cell's cytoplasm, allowing the dynamic organization of the components necessary for survival. However, it is not clear how phase separation is affected by the complex viscoelastic environment inside the cell. Here, we study theoretically how stiffness gradients influence droplet growth and arrangement. We show that stiffness gradients imply concentration gradients in the dilute phase, which transport droplet material from stiff to soft regions. Consequently, droplets dissolve in the stiff region, creating a dissolution front. Using a mean-field theory, we predict that the front emerges where the curvature of the elasticity profile is large and that it propagates diffusively. This elastic ripening can occur at much faster rates than classical Ostwald ripening, thus driving the dynamics. Our work shows how gradients in elastic properties control the size and arrangement of droplets, which has potential applications in soft matter physics and plays a role inside biological cells.


Introduction
Phase separation has recently been established as a crucial mechanism for organizing membrane-less organelles in biological cells. [1][2][3][4] These membrane-less organelles, or biomolecular condensates, often posses properties of liquid-like droplets, like internal rearrangement, spherical shapes, and monomer exchange with their surroundings. However, these droplets exist in complex environments that cannot be described as simple liquids. For example, the cytoskeleton in the cytosol 5,6 as well as the chromatin in the nucleoplasm 7 have solid-like properties, which could affect droplets. Indeed, recent experiments showed that droplets typically form in chromatin-poor regions in the nucleus. 8 This suggests that the elastic properties of chromatin suppress droplet formation. However, it is difficult to disentangle the effect of the elastic surroundings from other potential processes that could affect droplet formation in the complex case of a living cell.
Physical experiments with synthetic materials can help to understand how an elastic environment affects droplet growth. For instance, oil droplets growing in a homogeneous PDMS gel form mono-disperse emulsions and the stiffness of the gel controls droplet size. 9 There are several advantages of these experiments: first, the system is not driven, implying it relaxes toward equilibrium after preparation. Second, the gel is strongly cross-linked, so it does not spontaneously rearrange on the time scale of the experiments. 10 Viscous relaxation is thus negligible. Third, droplets are large compared to the mesh size, implying that the gel can be described by a continuum theory. Taken together, these properties allow to isolate the effects of an elastic environment on droplets.
The aim of the present paper is to understand how spatially varying stiffnesses affect droplet dynamics. This is motivated by recent experiments showing that stiffness gradients lead to elastic ripening, where droplets dissolve in stiffer regions. 11 Moreover, these experiments revealed a dissolution front invading stiffer regions, while the material of the dissolving droplets accumulated in softer regions. We already developed a theoretical description of the situation, which is based on the assumption that the gel exerts a pressure onto droplets that is proportional to the local stiffness. 11 Numerical simulations of this theory showed excellent agreement with the measured data. In the present paper, we analyze this model numerically and analytically to understand the details of elastic ripening. In particular, we use a simplified, coarse-grained version to derive scaling laws for where and when the dissolution front starts and how it evolves in time. These simple equations allow us to predict how the model parameters affect elastic ripening in more complex situations, e.g., when the elasticity profile varies in all spatial dimensions. interact by exchanging material via the dilute phase. This exchange is driven by the difference between the concentration in the dilute phase and the concentration right outside a droplet's surface. We determine the latter by considering a single droplet in an homogeneous elastic environment.

Thermodynamics of isolated droplets
Motivated by recent experimental results, 11 we describe a threecomponent system of two liquids and a network of cross-linked polymers. At high temperatures, the two liquids are miscible and form a homogeneous gel with the polymer network. When the temperature is lowered, one liquid forms droplets by segregating from the other two components. Therefore, this phase separation can effectively be described as a binary system; droplets with a high volume fraction f in of the segregating liquid coexist with a dilute gel phase of lower volume fraction f. In the case of thermodynamically large phases, these equilibrium volume fractions can be determined by minimizing the total free energy of the system using the Maxwell construction; 12 see ESI. † In particular, the equilibrium conditions imply that the pressures and chemical potentials are identical in both phases.
In the case of a single spherical droplet of radius R embedded in an isotropic elastic matrix, the droplet exhibits an additional pressure P due to both surface tension and elastic effects. In the simple case of a dense droplet phase and an ideal dilute phase, the resulting equilibrium volume fraction in the dilute phase can be approximated as (1) see ESI. † Here, c in is the concentration of the segregating fluid inside the droplet, k B is Boltzmann's constant, T is absolute temperature, and f 0 denotes the equilibrium volume fraction of the dilute phase in a thermodynamically large system for P = 0. Consequently, eqn (1) implies that the additional pressure P increases the equilibrium volume fraction in the dilute phase. The pressure P on a single spherical droplet of radius R embedded in an elastic gel is where g is the surface tension and P E is the pressure exerted by the isotropic elastic gel, which is related to its stress-strain curve. We here consider droplets that are much larger than the gel's mesh size, implying that the gel experienced large strains during droplet growth. Such situations typically arise when the gel has a maximal stress P C that it can sustain. For example, when growing droplets fracture the gel, P C is the critical stress that is required for fracturing. In the simplest case, P C is proportional to the Young's modulus E of the gel, P C = zE. The proportionality constant z can be determined experimentally or from theory. For example, Neo-Hookean theory 13,14 predicts z = 5/6, while the silicon gels used in the elastic ripening experiments exhibit z E 0.5. 15 Taken together, we here focus on the case where the pressure exerted by the gel is a simple function of the stiffness, P E = zE. For simplicity, we assume that this pressure is also exerted on shrinking droplets, ignoring small hysteresis effects due to gel fracturing. 15 To see whether surface tension g is important in the ripening experiments, we next estimate the relevant pressure differences between droplets. The pressure difference generated by surface tension is roughly gDR/R 2 , where R is the mean droplet radius and DR denotes the difference in the radii. Conversely, the pressure difference created by the external gel can be estimated by the typical stiffness difference DE. Surface tension is thus negligible if the ratio is small. In the ripening experiments, 11 we have g = 4 mN m À1 , DR E 5 mm, R E 10 mm, and DE = 740 kPa, implying K E 3 Â 10 À4 , so that Ostwald ripening is indeed negligible. Taken together, the equilibrium volume fraction f eq can thus be approximated by where Ê = c in k B T/z is the relevant stiffness scale. This expression allows us to determine the volume fraction f eq outside a droplet embedded in a gel described by a local stiffness E.

Dynamical equations of emulsions
We now consider an emulsion of droplets embedded in a gel with spatially varying stiffness E( x). We describe the system by the droplet radii R i and their positions x i , as well as the volume fraction f( x) in the dilute phase. The thermodynamics discussed in the previous section imply that the equilibrium volume fraction f eq right outside each droplet depends on its position. The difference between f eq and f drives a diffusive flux between the droplet and the dilute phase. Since we are only interested in dynamics on length scales larger than the droplet radii, we evaluate all involved quantities at the droplet position x i . Consequently, the material efflux J i integrated over the droplet surface can be expressed as where D is the molecular diffusivity. 12 This flux drives changes in droplet radius, where we used that the volume fraction f in inside the droplet is much larger than the fraction f eq outside. This equation describes how a droplet grows by taking up material from its immediate surrounding. On large length scales, material diffuses in the dilute phase, implying where V i = (4p/3)R i 3 are the individual droplet volumes. Here, the last term accounts for the material exchange with droplets. Note that we consider immobilized droplets whose positions x i are constant, since the elastic gel restricts the motion of the droplets.
In principle, our model allows overlapping droplets, but since they do not move and their initial positions are far apart, overlapping does not occur in practice. Consequently, eqn (5) and (6) fully describe how an emulsion of droplets evolves over time.
The dynamics of the system is governed by two diffusive fluxes that act on different length scales. Locally, material is exchanged between the droplets and the dilute phase by the flux J i . Conversely, transport on longer length scales only happens in the dilute phase by the redistribution flux ÀDrf.

Numerical simulations show dissolution fronts
We simulated eqn (5) and (6) to understand the effects of an elasticity gradient on the emulsion dynamics. Motivated by elastic ripening experiments, 11 we consider a system consisting of two regions of different stiffness E soft and E stiff . When the two regions are put side-by-side, a transition region emerges.
To capture this, we model the stiffness profile in the entire system by where x denotes the coordinate perpendicular to the interface and w is the width of the transition region. Generally, we chose parameter values motivated by experiments. 11 For instance, we initialized the simulations with tiny droplets distributed uniformly without overlap in the whole system and imposed a uniform volume fraction field f( x) in the dilute phase. Fig. 1 shows the time course of a typical simulation. Starting from the homogeneous initial condition, the system quickly forms two separate regions aligned with the stiffness profile (upper panel). Here, the stiff side exhibits smaller droplets and larger volume fractions in the dilute phase compared to the soft side. Droplets then start dissolving in the transition region and a dissolution front moves into the stiff side. Simultaneously, droplets grow on the soft side of the transition region while droplets far into the soft side remain unchanged. The exact same dynamics have been observed in the elastic ripening experiments; 11 see Fig. 2. These dynamics can be understood qualitatively by considering the diffusive fluxes in the system.
In the initial stage, the system is supersaturated everywhere, f 4 f eq . Consequently, material is transferred from the dilute phase to the droplets until a local equilibrium is reached, f = f eq . Eqn (4) implies that f eq is smaller for softer regions, so more material is absorbed by the droplets. We thus observe larger droplets in softer regions (see Fig. 1), consistent with experimental observations. 9 After the initial, local equilibration, material redistribution on longer length scales sets in. Since the stiffer side exhibits larger volume fractions f in the dilute phase, material is transported to the soft side. Consequently, on the stiff side, f drops below the local equilibrium volume fraction f eq , droplets shrink, and eventually dissolve. This process starts close to the transition region, since the redistribution flux is driven by gradients in f, which do not exist further away. Once droplets start disappearing in the transition region, droplets further away begin to be affected and a dissolution front forms that invades the stiff side. All the material redistributed from the stiff side is absorbed by the droplets on the soft side close to the transition region, which effectively shield all the other droplets on the soft side.

A coarse-grained model explains the dissolution dynamics
To understand the front's dynamics quantitatively, we next consider a simplified version of our model. Here, we do not describe the dynamics of individual droplets, but rather focus on the fractions of material in droplets and in the dilute phase. We thus introduce the coarse-grained volume fraction c of material contained in droplets, where the integrals are performed over a small discretization volume centered at x. Note that the discretization volume needs to be large enough to contain multiple droplets, but also small compared to the characteristic length scales of the elasticity gradient. Introducing the local mean droplet volume V( x,t), we can express c as c = f in nV, where n is the local droplet  The dynamics of the coarse-grained system follow from eqn (5) and (6). We show in the ESI, † that the dynamical equations are The first equation describes the local exchange of material between droplets and their surroundings, while the second equation accounts for the redistribution of material over long length-scales. Fig. 3 shows that the results of the numerical simulation of the coarse-grained model are virtually indistinguishable from that of the detailed model. Therefore, the coarse-grained model captures the essential features of the elastic ripening process. In particular, the dynamics of the dissolution front are governed by the material distribution, while individual droplets are irrelevant.

Dissolution starts at high curvatures of the elasticity profile
We now use the coarse-grained model to understand where and when the dissolution front appears, i.e., when droplets first disappear. In the numerical simulations, we observe that f quickly approaches f eq by local equilibration between the droplets and the dilute phase. Assuming that the system is initialized with f = f init and c = c init , the volume fractions approach f E f eq and c = c 0 with c 0 E c init + f init À f eq after this initial stage. Consequently, the volume fraction f in the dilute phase is controlled by the stiffness profile E( x), while the remaining material is absorbed in the droplets.
After the initial equilibration stage, the inhomogeneities in E, and thus in f, drive diffusive fluxes toward the soft side. However, we observe that these fluxes mostly affect c and hardly change f before the first droplets disappear. To understand the dynamics in this stage, we approximate eqn (9b) by q t c E Dr 2 f eq . Consequently, c evolves as (10) and the curvature of the equilibrium field f eq , set by the elasticity profile, controls droplet dynamics. In particular, droplets grow in convex regions (r 2 f eq 4 0), while they shrink in concave ones. Note that eqn (10) only holds when c 4 0, since otherwise droplets would be absent and the flux in the dilute phase changes f; see eqn (9b). We can use eqn (10) to estimate the time and position of the start of the dissolution front. In particular, droplets dissolve after a time t * ( x) E Àc 0 ( x)/(Dr 2 f eq ), when all material is removed from the droplet phase. The dissolution front starts at the earliest of these time points, t start = min x -(t * |t * Z 0), which is given by Fig. 2 Our model quantitatively captures the dynamics of dissolution fronts in the elastic ripening experiments. 11 Shown are the experimental and numerical mean droplet radii as a function of the distance from the interface on the stiff side for five different time points. The model parameters are E stiff = 0.0341Ê, E soft = 3.2 Â 10 À4 Ê, w = 0.37c, and the mean droplet radius on the stiff side is 0.09c; see ESI, † for further details on the comparison. Fig. 3 The coarse-grained model (red lines, eqn (9)) captures the detailed dynamics of the full model (blue lines, eqn (5) and (6)). Shown are the profiles of the volume fractions f in the dilute phase (upper panels) and the fraction c contained in droplets (lower panels) for three different time points t. The model parameters are given in Fig. 1.
where the minimum is constrained to regions where t * Z 0, i.e., where droplets shrink (r 2 f eq o 0). The location corresponding to the minimum denotes the starting position of the front. Eqn (11) highlights that the front appears where droplets are small and sparse (small c 0 ) as well as dissolve quickly (strongly negative curvature r 2 f eq ). The starting time t start can be estimated for the simple stiffness profile given by eqn (7). In particular, the droplet volume fraction c 0 will be close to the value c stiff deep into the stiff side; the curvature is approximately r 2 f eq B Df/w 2 , where Df = f eq (E stiff ) À f eq (E soft ) denotes the difference in the equilibrium volume fractions between the two sides and w is the width of the transition region. Using these estimates, eqn (11) suggests a time scalê which should govern the starting time of the front. In contrast, the starting position s start is dominated by the location of the largest negative curvature of f eq . For the simple stiffness profile given in eqn (7), this position should scale with the width w of the transition region. We test the prediction of eqn (11) and the scaling discussed above by comparing to numerical simulations of the full model; see Fig. 4. The collapse of the starting times shown in the left panel indicates thatt is the relevant time scale for this process. Moreover, the actual prediction t start given in eqn (11) is within a factor of two of the measured data. This analysis shows that the front appears earlier for larger stiffness differences (larger Df) between the two sides, for narrower transitions regions (smaller w), as well as when there is less material in the droplet phase (small c stiff ). Fig. 4b shows the corresponding starting positions, which clearly are determined by the width w of the transition region. The data collapse indicates that neither the absolute stiffnesses nor the droplet size affects where the front appears.

Two fronts move in opposite direction initially
Shortly after the first droplets dissolved, the surrounding droplets continue to shrink and disappear. Consequently, the region devoid of droplets expands in all directions. The material of the shrinking droplets is transported toward the soft side by diffusive fluxes. Initially, the material will accumulate where the equilibrium field f eq has the largest positive curvature (maximal r 2 f eq ); see eqn (10). The accumulating material is absorbed by the droplets in this region, which thus grow; see Fig. 1. The fact that droplets grow on the soft side implies that the front moving from s start toward this side slows down. Conversely, the front moving in the opposite direction can continue invading the stiff side.

Dissolution front moves diffusively at late times
We next analyze the late time dynamics of the front invading the stiff side. Specifically, we consider the case where the width L of the region devoid of droplets is large compared to the width w of the transition region (L c w). At this stage, the front moving toward the soft side is virtually stationary. We can thus understand the dynamics of the front invading the stiff side by analyzing the width L of the region devoid of droplets (where c = 0). The dynamics in this region are governed by a simple diffusion equation of the volume fraction f in the dilute phase; see eqn (9b). At the stiff side of this region, the dissolving droplets impose the equilibrium fraction f eq (E stiff ) as a boundary condition for the diffusion equation. The corresponding boundary condition at the soft side can be approximated by f eq (E soft ). For simplicity, we focus on slow fronts where the diffusion equation is in a stationary state, so the fraction f in the region devoid of droplets is governed by where Df = f eq (E stiff ) À f eq (E soft ) and x denotes the distance from the boundary on the soft side. The dynamics of L can be obtained by considering the change of the total amount of material on the stiff side (x Z 0). Because of material conservation, this change is equal to diffusive flux at x = 0, which can be determined from eqn (13). We show in the ESI, † that this implies q t L = a/(2L) with where c stiff is the droplet volume fraction deep into the stiff side. Consequently, the region devoid of droplets expands as where t 0 is such that L(t 0 ) = 0. This equation implies a diffusive motion of the front with diffusivity a. Fig. 4 The starting time t start and position s start of the front scale with the predicted time and length scales, respectively. Results from numerical simulations of the full model (dots) for various parameters are compared to the theoretical predictions from the coarse-grained model (gray line). The time point t start of the first dissolving droplet is shown in panel (a) as a function of the predicted associated time scalet given in eqn (12). The theoretical prediction given by eqn (11) is shown for c = 0.09f 0 and w = 1.45c. Panel (b) shows the associated starting position s start together with the equivalent prediction following from eqn (11), which is shown for c = 0.09f 0 and E stiff = 0.15Ê. The remaining parameters are E soft = 10 À4 Ê and f in = 1.
We test the theoretical prediction given in eqn (15) by comparing to numerical simulations of the full model. Fig. 5 shows the recorded times and positions when droplets dissolved (gray symbols), thus marking the dissolution fronts. The fronts start in the transition region on the stiff side and then move in opposite directions. The front moving toward the soft side slows down and comes to a halt on the soft side of the transition region, as predicted in the previous section. Conversely, the front invading the stiff side is quicker and does not stop. We measure its speed by fitting eqn (15) to the front positions deep into the stiff side to extract a and t 0 ; see Fig. 5a. Since the model explains the measured data at late times, we conclude that the front moves diffusively.
The fitted front diffusivity a is presented in Fig. 5b as a function of the relevant non-dimensional parameter Df/c stiff . This parameter compares the strength Df of the elastic ripening to the fraction c stiff of material that needs to be removed from the droplets. Consequently, the front is faster for larger Df/c stiff . The theoretical prediction for a, given in eqn (14), matches the data well for a o 2D. The fact that the front diffusivity a needs to be smaller than or comparable to the molecular diffusivity D is not surprising since we assumed that the front is slow enough for the region devoid of droplets to be in a stationary state; see eqn (13). Consequently, our theory predicts a maximal front diffusivity of 4D while the simulations indicate that faster fronts are possible.

Elasticity profiles spatially control droplets
So far, we considered elasticity profiles that only vary in one direction. However, our full model (eqn (5) and (6)) and the coarse-grained model (eqn (9)) are more general. In particular, eqn (11) implies that droplets first disappear in regions of strongly negative curvature r 2 f eq . Then, a front of dissolving droplets moves in all directions. The material of the dissolving droplets accumulates in regions of large positive curvature r 2 f eq , where droplets grow and remain for a long time. Consequently, where droplets grow or shrink is governed by f eq ( x) and thus the elasticity profile E( x). Fig. 6 shows a numerical simulation of the full model for a two-dimensional elasticity profile (left panel). A detailed simulation of a similar pattern has already identified that droplets accumulate in the soft valleys 8 and the time course shown in the right panels of Fig. 6 confirms that droplets follow the dynamics described above. A movie of this simulation, as well as one for a more complex elasticity profile, can be found in the ESI. † Taken together, this shows that we can engineer elasticity profiles to locate droplets in precise arrangements.

Conclusions
We presented a theoretical description of elastic ripening, which agrees quantitatively with experimental data. 11 Therefore, our theory identifies the driving mechanism of elastic ripening: the elastic matrix exerts a pressure onto droplets that raises the equilibrium concentration in their surroundings; gradients in this concentration then lead to diffusive material transport in the system. Surprisingly, the droplets do not start to dissolve in regions where the stiffness is maximal nor where its gradient is largest. In fact, our coarse-grained model reveals that droplets initially shrink faster where the curvature of stiffness is larger.  (15). Model parameters are E stiff = 0.09Ê and c = 0.09f 0 . Remaining parameters are given in Fig. 1. (b) The front diffusivity a (dots) determined from fitting to numerical simulations is compared to the prediction (line) given by eqn (14). Simulations were done for c stiff /f 0 = 0.03, 0.09, 0.30 for various E stiff , while the remaining parameters are the same as in panel (a). However, at late times droplets only remain in the softest regions. Taken together, we find complex dissolution dynamics, where for instance two fronts move in opposite directions, as in the experiments. 11 In particular, we derived scaling laws for the time when the first droplets disappear (eqn (12)) and for the late-stage dynamics of the dissolution front (eqn (15)). Taken together, these laws can be used to predict the effect of elastic ripening for various setups.
Elastic ripening competes with Ostwald ripening since both effects are driven by pressure differences between droplets. Their relative importance is quantified by K, given by eqn (3), when DE is the relevant stiffness difference in the droplets' surroundings. For instance, although elastic ripening dominates initially in the simulation shown in Fig. 6, Ostwald ripening will eventually occur between the remaining islands because they have the same elastic properties (DE E 0). A similar situation occurs in heterogeneous elastic materials, where DE corresponds to local stiffness variations. Thus, both elastic ripening and Ostwald ripening can happen in the same system, on different length-and time-scales.
The elastic ripening in stiffness gradients is similar to other droplet coarsening dynamics in gradient systems. For instance, concentration gradients, e.g., of regulating species that compete for mRNA binding, 16 have been shown to bias droplet locations in experimental 1 and theoretical studies. [17][18][19] Similarly, other external fields, like temperature gradients created by local heating 20,21 or even gravity 22 could be used to control droplet arrangements. Such systems can be analyzed using approaches that are similar to the ones presented here.
We showed that elastic ripening allows to control droplet arrangements, which could for instance be used in technical applications for micropatterning or for creating structural color. Moreover, our theory can help to understand the localization of biomolecular condensates in biological cells. For instance, elastic ripening explains experiments where droplets have been induced in the stiff regions of heterochromatin, but moved into softer regions immediately. 8 We expect that similar processes happen in the cytosol, where biomolecular condensates should be less likely where the cytoskeleton is dense. Interestingly, there are counterexamples, like centrosomes that localize to regions of high microtubule density [23][24][25] or ZO-1 clusters that concentrate in the acto-myosin cortex. 26,27 This seems to contradict elastic ripening, but in both examples the condensates interact with the elastic matrix: centrosomes bind the tubulin of microtubules 28 and the ZO-1 protein interacts with the F-actin of the cortex. 29 Consequently, there are two competing gradients in this situation: droplets are repelled by the stiffness of the surrounding matrix but are attracted by its molecular components. Indeed, when the actin-binding domain of ZO-1 is removed, the clusters do not accumulate in the cortex anymore, but are more broadly distributed, 26 as predicted by elastic ripening.
Our theoretical description of elastic ripening can be naturally extended to include other effects. In fact, the dynamics described by eqn (5) and (6) already include surface tension effects when the pressure given by eqn (2) is used. Although we did not analyze the impact of surface tension since it is negligible in the experiments, it will become important on longer timescales or when surface tensions are large. Moreover, the elastic properties of the matrix might be more complex than considered here. For instance, the cytoskeleton can shown strain-stiffening, which might arrest droplet growth, as well as visco-elastic effects, which allow to relax elastic stresses. 6 The latter effect can be captured by the theory of viscoelastic phase separation, which affects the coarsening behavior. 30,31 This stress-relaxation, as well as the reduced stress due to hysteresis effects, 15 could slow down elastic ripening. Finally, the droplets themselves can possess elastic properties. This is particularly true in biological condensates, 32 which sometimes even form solid-like aggregates 33 that potentially cause diseases. 4 All these effects might be crucial for understanding the behavior of biomolecular condensates in cells.

Conflicts of interest
The authors declare no conflicts of interest.