Liquid–liquid phase separation driven compartmentalization of reactive nucleoplasm

Phase separation is a fundamental mechanism for the emergent order in an ordinary and biological matter [1, 2]. Recently, phase separation of bio-molecules has also become a cornerstone physical mechanism for understanding the intracellular organization [2–4]. A wide range of membrane-less compartments are found to form through biomolecular phase separation, including nucleoli [5–7], stress granules [8–10], chromatin domains [11–13] and transcriptional centers [14, 15]. The primary components driving the intra-cellular phase separation are proteins and nucleic acids with multi-valent interaction centers [16, 17].

The stickers-and-spacers framework has emerged as a viable model explaining the existence of a broad class of sequence encoded driving forces of disordered proteins which serve as nucleating centers for biomolecular condensates [3, 18]. These newly appreciated abilities of proteins and nucleic acids for forming large-scale liquid bodies is offering fresh avenues for understanding mechanisms of the coordinated action of biomolecules in gene regulation and cellular organization that go beyond single-molecule action. It is well known that eukaryotic nuclei are rich in disordered proteins linked with transcription and chromatin architecture reorganization activities [19, 20]. There have been several proposals of the functional roles that may include catalysis of biochemical reactions, noise buffering and inducing ultra-sensitive signals [21–23]. However, understanding the mechanistic picture that links phase separation to the functional gene regulatory processes inside the nucleus has remained elusive [22] due to heterogeneous, multi-scale, and non-equilibrium nature of the nuclear environment [24].

In this work, we propose a minimal model of a reactive nucleoplasm (figure 1) with an objective to illustrate, in a proof of principle manner; (i) how spatially resolved non-equilibrium reactive en-events generate qualitatively distinct from equilibrium phase behavior and (ii) how the interplay of various kinetic timescales in the system impacts gene expression patterns. The minimal reactive nucleoplasm model consists of a ternary solution filled with incompressible fluid consisting of 'active' protein, RNA components, and 'passive' nucleoplasmic buffer. The mathematical formulation of the model is based on a generic ternary diffuse interface model of one-step transcription/translation and phase-separation inside the nucleoplasm. The model resolves the formation and dissolution of protein–RNA droplets as well as reactive events of influx/creation and outflux/degradation of proteins and RNA. By an extensive set of simulations exploring the interplay of timescales in the system, we found a broad kinetic regime dominated by length-scale heterogeneity of phase-separated droplets which correspond to heterogeneous gene expression patterns.

Zoom In Zoom Out Reset image size

Various simple mathematical models have been used to explore generic aspects of equilibrium and non-equilibrium phase-separation as well as the physical properties of resulting condensates [25–33]. The most relevant to the present contribution are the work by Berry et al [29] , Tang and Yang [33], and Glotzer [32], where authors have considered a binary and ternary fluid model of nucleoplasm, which couples phase separation with first-order exchange among soluble and insoluble components. Authors have simulated different stages of spinodal decomposition and explored its impact on Ostwald ripening of droplets. It was shown that kinetics of the RNA flux accelerates the ripening of droplets, thereby showing a link between the thermodynamics of phase separation and kinetics of droplet formation. Other recent notable studies worth highlighting here include the work by Yamamoto et al [30], which have studied the non-spatial model of architectural RNA phase separation coupled with non-equilibrium production. Finally, an important theoretical framework by Ilker and Joanny [34] establishes and equivalence of phase-separation kinetics with the Cahn–Hilliard effective temperature.

Some of the key distinctions of the present model from previous studies include (i) spatially resolved formulation of gene expression and phase separation (ii) explicit connection of phase-separated patterns with transcription and translation models of gene expression (iii) exploring the impact of dynamic turnover in RNA–protein–nucleoplasm ternary diagram on global patterning. Thus, the present work clarifies as a first step the dual nature of nuclear order and gene expression, and provides a useful theoretical framework for understanding equilibrium and non-equilibrium origins of intra-nuclear patterning.

In order to explore the kinetics of phase separation of RNA–proteins–nucleoplasm reacting mixture, we use a local thermodynamic approach that describes the phase separation of a ternary mixture based on the Flory–Huggins model [35, 36] for polymer solutions coupled with chemical reaction-diffusion equations (figure 1). In this section, we describe the mathematical formulation of the model of a 'minimal reactive nucleoplasm': where the liquid–liquid phase separation of single RNA and protein components is coupled with reactive events of transcription, translation, and degradation.

A minimal one-step model of independent, unregulated transcription and translation [37] is used to set the lifetimes of protein and RNA components comprising the nucleoplasmic milieu defined by the following the chemical reactions:

$$\begin{align}\hfill & \varnothing {k}_{\mathrm{R}}{\to }\mathrm{R}{k}_{\mathrm{d}}{\to }\varnothing ,\hfill \\ \hfill & \varnothing {k}_{\mathrm{P}}{\to }\mathrm{P}{k}_{\mathrm{d}}{\to }\varnothing ,\hfill \end{align} \tag{ 1 }$$

where the k_R, k_P and k_d are the rate coefficients for transcription, translation and degradation, respectively; The ∅ symbol is used to denote the degradation of RNA and protein. The kinetics of phase separation of ternary reactive protein–RNA–nucleoplasm mixture is given by the following set of reaction–diffusion equations [32]:

$$\begin{align}\hfill \frac{\partial {\varphi }_{1}}{\partial t}& =\nabla \left({M}_{1}\left(\left\{{\varphi }_{i}\right\}\right)\nabla \left(\frac{\delta F\left[\left\{{\varphi }_{i}\right\}\right]}{\delta {\varphi }_{1}}\right)\right)\hfill \\ \hfill & \quad +{f}_{1}\left(r\right)\left[{k}_{\mathrm{R}}\left(1-{\varphi }_{1}-{\varphi }_{2}\right)-{k}_{\mathrm{d}}{\varphi }_{1}\right]\hfill \end{align} \tag{ 2 }$$

$$\begin{align}\hfill \frac{\partial {\varphi }_{2}}{\partial t}& =\nabla \left({M}_{2}\left(\left\{{\varphi }_{i}\right\}\right)\nabla \left(\frac{\delta F\left[\left\{{\varphi }_{i}\right\}\right]}{\delta {\varphi }_{2}}\right)\right)\hfill \\ \hfill & \quad +{f}_{2}\left(r\right)\left[{k}_{\mathrm{P}}\left(1-{\varphi }_{1}-{\varphi }_{2}\right)-{k}_{\mathrm{d}}{\varphi }_{2}\right],\hfill \end{align} \tag{ 3 }$$

where φ_i (i = 1, 2, 3) are the order parameter variables associated with the local concentration of ith-component forming the mixture, M_i are their corresponding mobility coefficients, and F is the free energy functional describing the thermodynamics of ternary mixture. The total local density in the system is considered constant, respecting the incompressibility condition expressed as: ${\sum }_{i=1}^{3}{\varphi }_{i}=1$. The indices 1, 2, and 3 denote RNA, protein, and nucleoplasm, respectively. The function f₁(r) and f₂(r) are used to define the spatially heterogeneous distribution of components mimicking RNA expression inside nucleus and protein translation and flow into nucleus. The RNA is created in the center of nucleoplasm thereby mimicking a transcriptional process [38], herein we used the function defined as ${f}_{1}\left(r\right)=\mathrm{exp}\left(-{\left(r-{r}_{\mathrm{c}}\right)}^{2}/a\right)$ where r is the distance from the center of the domain located at r_c. The coefficient a is a localization lengthscale which is set to a = 1. The protein component is flown into the nucleus from the nuclear boundaries thereby mimicking translation [38], herein we used the function defined as f₂(r) = 1 on the domain boundary ∂Ω and 0 otherwise. Heterogeneous mobility is an interesting aspect of nucleoplasm [39], however, and is undoubtedly deserving of a separate investigation. For simplicity, we assume the mobility coefficients for all components to be constant M₁({φ_i }) = M₂({φ_i }) = M.

The free energy functional taken in this model is based on the Flory–Huggins energy formulation of the ternary mixture system [40, 41], which is given by:

$$\begin{equation}F\left[\left\{{\varphi }_{i}\right\}\right]=\int \mathrm{d}{\Omega}\left[{f}_{\text{bulk}}\left({\varphi }_{1},{\varphi }_{2},{\varphi }_{3}\right)+\sum _{i}\frac{{\kappa }_{{\varphi }_{i}}}{2}{\left\vert \nabla {\varphi }_{i}\right\vert }^{2}\right]\end{equation} \tag{ 4 }$$

with

$$\begin{equation*}{f}_{\text{bulk}}\left({\varphi }_{1},{\varphi }_{2},{\varphi }_{3}\right)=\sum _{i}\frac{{\varphi }_{i}}{{N}_{i}}\enspace \mathrm{ln}\left({\varphi }_{i}\right)+\sum _{i\ne j}{\chi }_{ij}\enspace {\varphi }_{i}\enspace {\varphi }_{j},\end{equation*}$$

where f_bulk(φ₁, φ₂, φ₃) is the local free energy density, and the square-gradient coefficients ${\kappa }_{{\varphi }_{i}}$ are positive constants controlling the interfacial free energies. χ_ij are the interaction parameters between species i–j, and N_i is the degree of polymerization.

We note here that there also exist other mathematical formulations of the local free energy for describing the phase separation of ternary polymeric mixtures. For instance, the generalized multiwell potential used by Yang et al [42] to study multiphase systems where the local free energy density is the summation of all the double-well potential for each phase-field variable complemented by a polynomial term that connects all the variables altogether.

The dynamic equations are brought to a dimensionless form, which reveals the essential time- and length-scales of the problem. Let us introduce l, the characteristic length of the system, and τ is the characteristic time. The dimensionless form of the kinetic equations yields four dimensionless parameters that control the dynamics of phase-separation where components are undergoing chemical reactions. The τ_D/τ = l²/(M × τ) = 1 is the ratio between diffusion timescale and characteristic time which is fixed to be a unit. The τ/τ_T = τ × k_T is the ratio between characteristic time and time-scale associated with transcription. The τ/τ_P = τ × k_P is the ratio between characteristic time and time-scale associated with protein formation. The τ/τ_d = τ × k_d is ratio between characteristic time and time-scale associated with degradation. To solve the dynamic equations (2) and (3) in the dimensionless form, we use a fully implicit finite element C++ library from the multiphysics object-oriented simulation environment (MOOSE) [43]. The simulations were performed on a rectangular domain of dimensions L_x × L_y : (50 × 50), with periodic boundary conditions. A quadrilateral element QUAD4 with four nodes was used for domain meshing with the refinement. The total number of elements used for the fine mesh is 10 000. The time step of integration Δt is fixed at 0.05 (a.u.). In this work, we only consider the case of a symmetric ternary mixture undergoing a chemical reaction with χ₁₂ = χ₂₃ = χ₁₃ = 3 that ensures phase-separation of protein–RNA droplets in the absence of chemical reactions χ > χ_c, where χ_c denote the critical point. We also assumed that the chemical reaction takes place at the comparable scales with the phase separation process. Near the critical point χ_c, the mean-field Flory–Huggins free energy can be approximated through Taylor expansion via Landau form [32]. It is important to construct representative free energy functional that accounts for fluctuations of the phase-field variables that are dominant near the critical point. To express the critical behavior a new formulation of the free energy of mixing renormalized by the spatial variation has been proposed by Yamamoto et al [44, 45]. The chain lengths N_i , or degree of polymerization, are assumed to be the same: N₁ = N₂ = N₃ = 1. It is noted here that the entropic term of the local free energy is proportional to the inverse polymeric lengths N_i . As a consequence of the reduction of entropy for long-chain lengths, the phase diagrams and the spinodal curve will be modified in the way to increase the separating-phase region of the phase diagram. In this case, a small variation of the interaction parameter between species will lead easily to phase separation. The impact of phase-separation of small and large chains such as genomic regions and proteins/RNA could be an interesting study to address in future investigations. The other parameters of simulations were set to ${\kappa }_{{\phi }_{1}}={\kappa }_{{\phi }_{2}}=0.031\enspace 25$. The initial configuration of phase-field variables is generated by φ_i (r) = ⟨φ_i ⟩ + δφ_i , where ⟨φ_i ⟩ is an initial average concentration associated with φ_i , and δφ_i is a small random perturbation amplitude. For all simulations presented here, the initial average of species i are set to ⟨φ₁⟩ = 0.3 and ⟨φ₂⟩ = 0.09, with δφ₁ ∈ [−0.05, 0.05] and δφ₂ ∈ [−0.01, 0.01].

Here we report the findings obtained by analyzing the results of simulations with a minimal reactive nucleoplasm model. We have organized the results in three broad kinetic regimes, which are described by three orders of magnitude in degradation time-scale τ_d = τ, 10τ, 100τ with respect to the diffusion time-scale τ. These time-scales correspond to different rates of turnover of nucleoplasmic components from fast to slow. For each regime, we have carried out extensive grid-based kinetic parameter sweeps exploring the coupling of transcriptional and transnational time-scales on the backdrop of fixed phase-separating free energy landscape of RNA–protein–nucleoplasm components. Despite the incredible simplicity of the minimal reactive nucleoplasm model, the time course of simulations (figures 2–4) has revealed a non-trivial patterning of nucleoplasm which are a dramatic departure from equilibrium thermodynamics of ternary phase separation in the absence of spatially non-uniform reaction–diffusion.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We start by investigating the fast dynamical regime corresponding to τ_d = τ. Figure 2 shows the spatial and temporal evolution of RNA–protein components in a reactive nucleoplasm environment. Four interesting cases are highlighted in figure 2, corresponding to different time-scale separation between translation τ_P and transcription τ_R processes. When τ_P = τ_R = τ we observe rapid coarsening dynamics (droplet growth kinetics) and emergence of RNA droplets with stable 'protein front' by which we refer to the formation of protein layer surrounding nucleoplasm domain edges [figure 2(A)]. Coarsening dynamics for either protein front or RNA droplets are also observed for fast translation/slow transcription or slow translation/fast transcription, respectively [figures 2(B) and (C)]. Naturally, it is expected that the more rapid degradation timescales for RNA and protein lead to the disappearance of both RNA droplets and protein fonts [figure 2(D)]. Thus, we can conclude that in the fast dynamical regime, nucleoplasm patterns are entirely set by the kinetic parameters with thermodynamic free energy landscape taking the back seat.

The formation of RNA–protein patterns is more intricate in the intermediate dynamical regime corresponding to τ_d = 10τ (figure 3). In the figure 3, we show the patterns arising from competing transcription, translation, and degradation with different transcription and translation timescales. When both transcription and translation have comparable time scales [figure 3(A)] to that of diffusion, we find a rapid progression of protein front on the one hand and a rapid RNA seed droplet growth at on the other. In this regime, once the nucleoplasm reaches a non-equilibrium steady state, droplet patterning is now becomes dictated by the thermodynamics of Flory–Huggins interaction parameters, which set the extend of mixing between protein and RNA components. When both the transcription and translation have the same time scale but are now significantly slower than diffusion [figure 3(D)], then the system displays heterogeneity in protein RNA droplet distribution. The situation with faster transcription and the faster translation is predictable: [figures 3(B) and (C)] the system tends to increase in protein front and RNA nucleating center, respectively. We have only highlighted the most interesting cases for each dynamical regime. The simulation results for the full set of time-scales are presented in the supporting information (https://stacks.iop.org/PB/18/015001/mmedia) figure S1–S10. The temporal evolution of average concentration for different kinetic regimes with τ_R = τ_P = 100τ is shown in figure S14.

We now turn to the analysis of the slow dynamical regime corresponding to a degradation timescale τ_d = 100τ. In this regime, we once again highlight four interesting cases [figures 4(A)–(D)]. When the system has comparable timescales for transcription, translation, and diffusion τ_P = τ_R = τ, the nucleoplasmic component gets eliminated, and the system rapidly reaches a balance between protein front and RNA seed. Slowing down the translation τ_P = 100τ_R = 100τ by two orders of magnitude leads to a non-trivial patterning with the RNA seed, nucleoplasm and RNA droplets all con-existing at a non-equilibrium steady-state. On the other hand, slowing down the transcription [figure 4(C)] leads to a steady-state with significantly downsized RNA-seed at the expense of the peripheral protein front. Slowing down both transcription and translation [figure 4(D)] leads to a non-trivial patterning with protein front, RNA droplets, and nucleoplasmic environment, but this time with no dominant RNA-seed.

To quantify the emergence of different droplet patterns, we have summarized the phase behavior of a reactive nucleoplasm model via. Global kinetic phase diagram, figure 5. Here one can see a global picture of how the interplay of kinetic timescales is favoring uniform vs binary vs ternary phases.

Zoom In Zoom Out Reset image size

In order to quantify the length-scales of emergent patterns, we have analyzed the azimuthally averaged dynamic structure factors S(k, t) = ∫dk_Ω S(k, t) associated with each dynamical regime (see figure 6 and supporting information, figure S11–S13). Analyzing dynamic structure factors reveals characteristic length-scales of protein/RNA droplets as well as emergent dynamical heterogeneity manifesting in the fast or slow coarsening of the droplets. The computed structure factors show a broad range of kinetically controlled states where one has bi-modality of RNA (or protein) components. This bimodality or more broadly heterogeneity of distribution is directly linked to the heterogeneity of droplet sizes and shapes that one can quantify from simulation images (figures 3 and 4). Multi-modal gene expression is a feature often linked with phenotypic heterogeneity and which has been mostly explained by citing the underlying non-linear dynamics of dichotomous switching noise in a spatially uniform master equation formalism [46–48]. The results of figure 6 clearly show that transcriptional heterogeneity can originate purely from the spatially non-uniform nature of gene expression, which is being modulated by timescales of phase separation, transcription, and translation.

Zoom In Zoom Out Reset image size

Finally, we have computed the length scale of patterns summarizing transitional heterogeneity in three kinetic regimes of the minimal reactive nucleoplasm model (figure 7). The length-scales patterns are quantified by using pre-computed azimuthally averaged dynamic structure factors $R\left(t\right)=\frac{\int \mathrm{d}kS\left(\mathbf{\text{k}},t\right){k}^{-2}}{\int \mathrm{d}kS\left(\mathbf{\text{k}},t\right){k}^{-1}}$ [49]. The length-scale patterns R(t) quantify the dominant length-scales which emerge during the time-evolution of a reactive nucleoplasm. Analyzing the evolution of length-scale patterns for three dynamical regimes, we see clearly that for the fast dynamical regime, there is only one dominating length-scale, which is dictated by the rapid degradation kinetic timescale. For the intermediate and slow dynamical regimes, however, we find heterogeneity of length-scales, which emerges from the disparity in transcription and translational timescales. We note that this heterogeneity has both structural and dynamical manifestations, as one can see by analyzing the power low of length-scale patterns R(t) ∼ At^α . We find two exponents, one which is characteristic for the early coarsening stage (α ∼ 1/3) and second (α ∼ 3/8) for the later accelerated phase separating evolution toward a steady state. We note that power laws in the dynamical variables of the nucleoplasmic environment have been detected and characterized in a large number of experiments [39, 50]. These power-law dependencies, however, are hard to disentangle in terms of distinct contributions since their origin may come from any combination of phase-separation, polymeric effects, confinement, and non-equilibrium motorized activities in the nucleus. In this work, we have only managed to scratch at the surface of the fascinating dynamical patterning potential of the active nucleoplasmic environment. For future studies, it would be interesting to investigate the impact of differential mobility, hydro-dynamical coupling, and as well as investigate gene expression beyond a simple one-step model of unregulated reactions that have not been done in the present contribution.

Zoom In Zoom Out Reset image size

Recent in vitro experiments with binary and ternary protein and RNA mixtures [51–53] have shown the rich complexity of phases which can emerge through liquid–liquid phase separation via modulation of stoichiometry of components and point mutations. These experiments show potentially novel regulatory strategies which when combined with non-equilibrium cellular processes such as transcription and translation could produce coordinated gene regulatory and signaling actions.

To this end, in the present contribution we introduce a minimal reactive nucleoplasm model combining spatially dependent transcription and translation to liquid–liquid phase separation of RNA and protein components embedded in the backdrop of passive nucleoplasm buffer. We use the minimal reactive-nucleoplasm model to cleanly dissect how the interplay of transcription, translation, and degradation time scales couples with the liquid–liquid phase separation of RNA and protein components in a model ternary solution. By carrying out extensive grid-based sweeps of kinetic parameter space, we uncover various non-trivial patterning and length-scale heterogeneity compared to the classic Flory–Huggins thermodynamic picture for ternary polymeric solutions. Our central finding is the existence of a broad kinetic regime characterized by a slow turnover of components and timescale disparity between transcription and translation under which a phase separating system can display bi-modal distribution. The significance of this finding is that the observed heterogeneity of gene expression is linked directly with the heterogeneity of length-scales in phase-separated condensates. The main findings and the minimal reactive nucleoplasm model thus establishes a useful framework with which one can further elucidate the emergence of nucleoplasm patterns and phenotype heterogeneity from first principles modeling of phase-separation and reaction–diffusion processes.

Research reported in this publication was supported by the National Institute of General Medical Sciences of the National Institutes of Health under Award No. R35GM138243. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation Grant No. ACI-1548562 [54] on the Stampede2 machine at the Texas Advanced Computing Center (TACC) through allocation CTS190023. The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper. The authors also acknowledge financial support from Iowa State University.