Liquid–liquid phase separation driven compartmentalization of reactive nucleoplasm

The nucleus of eukaryotic cells harbors active and out of equilibrium environments conducive to diverse gene regulatory processes. On a molecular scale, gene regulatory processes take place within hierarchically compartmentalized sub-nuclear bodies. While the impact of nuclear structure on gene regulation is widely appreciated, it has remained much less clear whether and how gene regulation is impacting nuclear order itself. Recently, the liquid–liquid phase separation emerged as a fundamental mechanism driving the formation of biomolecular condensates, including membrane-less organelles, chromatin territories, and transcriptional domains. The transience and environmental sensitivity of biomolecular condensation are strongly suggestive of kinetic gene-regulatory control of phase separation. To better understand kinetic aspects controlling biomolecular phase-separation, we have constructed a minimalist model of the reactive nucleoplasm. The model is based on the Cahn–Hilliard formulation of ternary protein–RNA–nucleoplasm components coupled to non-equilibrium and spatially dependent gene expression. We find a broad range of kinetic regimes through an extensive set of simulations where the interplay of phase separation and reactive timescales can generate heterogeneous multi-modal gene expression patterns. Furthermore, the significance of this finding is that heterogeneity of gene expression is linked directly with the heterogeneity of length-scales in phase-separated condensates.

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 ultrasensitive signals [21][22][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, multiscale, 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 nonequilibrium reactive enevents 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 phaseseparation inside the nucleoplasm. The model resolves the formation and dissolution of protein-RNA droplets as well as reactive events of influx/creation and out-flux/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.
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][26][27][28][29][30][31][32][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 firstorder 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.

Minimal reactive nucleoplasm model
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: 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]: 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: ∑ i = 1 3 φ i = 1. The indices 1, 2, and 3 denote RNA, protein, and nucleoplasm, respectively. The function f 1 (r) and f 2 (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 (r) = exp(−(r − r c ) 2 /a) 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 2 (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 1 ({φ i }) = M 2 ({φ 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: where f bulk (φ 1 , φ 2 , φ 3 ) is the local free energy density, and the square-gradient coefficients κ φ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 2 /(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 objectoriented 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 χ 12 = χ 23 = χ 13 = 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 1 = N 2 = N 3 = 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 κ ∅1 = κ ∅2 = 0.031 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 ⟨φ 1 ⟩ = 0.3 and ⟨φ 2 ⟩ = 0.09, with δφ 1 ∈ [−0.05, 0.05] and δφ 2 ∈ [−0.01, 0.01].

Results
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.
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 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 tempo ral 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 high-light 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. 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.
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][47][48]. The results of figure 6 clearly show that transcriptional heterogeneity can originate purely from the spatially nonuniform nature of gene expression, which is being modulated by timescales of phase separation, transcription, and translation.
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 t = ∫ dkS k, t k −2 ∫ dkS k, t 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.

Conclusion
Recent in vitro experiments with binary and ternary protein and RNA mixtures [51][52][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 nonequilibrium 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 back-drop 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.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.  Evolution of phase-field variables φ i for three phase mixture undergoing the chemical reaction with a degradation time-scale fixed at the same of diffusion (τ d = τ ). From up to bottom: snapshots corresponding to simulation results with (A) τ R = τ P = τ; (B) τ R = 100τ P = 100τ; (C) τ P = 100τ R = 100τ ; and (D) τ R = τ P = 100τ . The color code in blue, green, and red indicates the RNA, protein, and nucleoplasm regions, respectively.  Evolution of phase-field variables φ i for three phase mixture undergoing the chemical reaction with a degradation time-scale fixed at the same of diffusion (τ d = 10τ). From up to bottom: snapshots corresponding to simulation results with (A) τ R = τ P = τ ; (B) τ R = 50τ P = 50τ; (C) τ P = 100τ R = 100τ ; and (D) τ R = τ P = 50τ . The color code in blue, green, and red indicates the RNA, protein, and nucleoplasm regions, respectively.  Evolution of phase-field variables φ i for three phase mixture undergoing the chemical reaction with a degradation time-scale fixed at the same of diffusion (τ d = 100τ ). From up to bottom: snapshots corresponding to simulation results with (A) τ P = 10τ R = 10τ; (B) τ R = 10τ P = 50τ; (C) τ R = 50τ P = 10τ; and (D) τ R = τ P = 100τ. The color code in blue, green, and red indicates the RNA, protein, and nucleoplasm regions, respectively.  The blue star indicates that the RNA domain becomes the dominant phase in the long timescale limit. The green star means that the proteins-droplets become the dominant phase in the long timescale limit. The pink star indicates three coexisting phases of the ternary mixture. In contrast, the red star indicates that RNA-proteins turn totally into the nucleoplasm.