Phase separation in soft repulsive polymer mixtures: foundation and implication for chromatin organization

Given the wide range of length scales, the analysis of polymer systems often requires coarse-graining, for which various levels of description may be possible depending on the phenomenon under consideration. Here, we provide a super-coarse grained description, where polymers are represented as a succession of mesosopic soft beads which are allowed to overlap with others. We then investigate the phase separation behaviors in a mixture of such homopolymers based on mean-field theory, and discuss universal aspects of the miscibility phase diagram in comparison with the numerical simulation. We also discuss an extension of our analysis to mixtures involving random copolymers, which might be interesting in the context of chromatin organization in a cell nucleus.


Introduction
Phase separations in polymer solution and blend have a long history of research due to its importance in fundamental science as well as industrial applications [1][2][3][4][5][6][7] .Recently, its pivotal role in the field of biophysics has been recognized as a basic mechanism to organize various cellular and nuclear bodies [8][9][10][11][12][13][14][15][16][17][18][19] .Here, given the complexity in biological systems, standard approaches such as the Flory-Huggins theory to analyze the phase separation do not always suffice, and various extensions or modifications may be called for depending on the phenomena under consideration.
In this paper, we provide one such example, where we investigate the phase separation behavior of polymer mixtures made of mesoscopic segments.Our work has been motivated by the recent attempts to simulate chromatin organization in cellular nucleus.In ref. 20 , Fujishiro and Sasai constructed a polymer model of the whole genome of human cells, where each chromatin is modeled as a succession of soft-core monomers.Here, individual monomers (beads) represent ∼ 10 2 kbp of DNA, which is much larger than the conventional monomers defined in standard theory or simulation of polymer systems.They argued that the interaction between such mesoscopic segments is soft and repulsive, and the imbalance in such repulsion in systems with e.g., eu-and hetero-chromatic monomers would trigger the phase separation.Similar modelings of large scale behavior of chromatin with soft-core potential naturally arises after the coarse-graining, hence have been employed in other works as well [21][22][23][24][25] , where the soft potential incorporates the entropic effect relevant to the mesoscipic segments.
How can we describe such a phase separation phenomena in chromatin theoretically?The immediate complication lies in the copolymer nature of the chromatin model 21,26 .However, even if we let aside the sequence effect and restrict our attention to a binary mixture of homopolymers, the application of a Flory-Huggins theory is hampered because of the allowed overlap between monomers due to the soft-core nature of the inter-monomer potentials.A key insight would thus be obtained by the phase behavior of the mixture of soft particles.This problem has been extensively studied by groups of Likos, Löwen and Kahl [27][28][29] .Very recently, Staňo, Likos and Egorov have extended the framework to the system of chains of soft beads 30 .Although their primal target is a mixture of linear polymers and ring polymers (or polycatenanes), we expect that the same physics applies to chromatin system, too.
Our first aim is thus to recapitulate and to numerically validate the theoretical framework for the binary mixture of polymers made of soft monomers, which allows one to analyze the phase behavior.In Sec. 2, we introduce the mean-field free energy for our system.From the analysis of the free energy, we present in Sec. 3, the miscibility phase diagram and compare it with the result from molecular dynamic (MD) simulations.Sec. 4 is devoted to discussions on universal aspects of the phase diagram, comparison with a conventional Flory-Huggins theory, and connection to the Gaussian core model.Building on the framework, we also discuss its extension to a system containing copolymers in Sec. 5.

Free energy of the soft repulsive polymer mixture
Following Staňo 30 , we adopt the following free energy density for the mixture of polymers modeled as the succession of soft beads which represent mesoscopic segments where k B T is thermal energy, c x and N x are, respectively, the number density of beads and the chain length (number of beads per chain) of component x(= a or b).Parameters χ xy (> 0) represent the strength of the repulsive interaction between beads x and y.Note that in this representation, the parameters χ xy have a unit of volume, and we measure them in unit of the volume of individual beads.In other words, we assume, for simplicity, the characteristic size (σ ) of beads a and b are equal, which is taken to be the unit of length.Although there is no attraction, the phase separation may be induced by the asymmetry in the repulsion, i.e., χ aa ̸ = χ bb .At first sight, Eq. ( 1) looks as a free energy in the second virial approximation valid for low concentrations.As we shall show below, however, the free energy (1) is capable of describing the phase separation in concentrated regime (c a + c b )σ 3 > 1 (see Seq. 4.3 for discussion).
Let us first clarify a mathematical aspect relevant to the phase equilibria condition in the system described by the free energy Eq.(1).If the homogeneous mixture of polymer A and polymer B separates into phase 1 and phase 2, the demixed state is specified by the concentrations of both components in respective phases, i.e., (c b ).The number of unknowns is thus n u = 4.
On the other hand, the phase equilibria between two phases indicates the equalities of chemical potentials µ between two phases (α = 1 or 2) for both components (x = a or b), i.e., µ b and also the mechanical balance ensured by the equality of pressure P (1) = P (2) , where To check the validity of the free energy prediction, we have performed numerical simulations of the polymer mixture.Briefly, the system is a mixture of two types of linear homopolymers A and B, where A (B) polymer is made of a succession of N a (N b ) beads of type a (b) (see Appendix for details of the simulation model).To represent the soft repulsion between monomers, we employ the Gaussian potential, see Eq. (15) in Appendix, where the strength of the repulsion between x-bead and y-bead is ε xy in unit of k B T .
The numerically determined phase boundary shown in Fig. 2  (a b ) = (0.25, 0.25), and further increase in concentration leads to a well-defined phase separated structure (Fig. 2 (b)).Remarkably, the numerically determined phase diagram resembles that predicted by the free energy analysis (Fig. 1).More specifically, we find that the numerical and analytical phase diagrams almost overlap under the correspondence χ xy /ε xy ≃ 2.5 between interaction parameters in free energy and interaction strengths in simulation.

General aspects of phase diagram
In Sec. 3, we have shown one example how the phase boundary alters with the change in chain length or the interaction strength (Fig. 1 (b)).To clarify the dependence of the shape of phase diagram on system parameters in a more systematic way, it is desir- able to find out universal aspects inherent to the model described by the free energy (1).
To this end, we introduce the rescaled concentrations ca = c a N a χ aa and cb = c b N b χ bb , which enables us to rewrite Eq. ( 1) as where irrelevant linear terms in concentrations are dropped, and coefficients are The above free energy density is invariant under the parameter changes which keep k 1 and k 3 = χ ab / √ χ aa χ bb constant.These conditions are satisfied by the following transformations where p, q, k are positive real numbers.Therefore, with the change in parameters according to Eqs (4) and ( 5), the phase diagram drawn in terms of rescaled concentrations remains the same.A similar analysis has been done by Staňo et.al 30 .Indeed, from the stability analysis of the uniform state, one finds the spinodal curve and the critical point is determined by Eq. ( 6) together with Eq. ( 6) indicates that a necessary condition for the phase separation k 3 > 1 ⇔ χ ab > √ χ aa χ bb , and Eq. ( 7) determines the location of critical point on the spinodal curve 30 .In Fig. 3, we demonstrate a collapse of the phase diagram upon rescaling.

Comparison with Flory-Huggins theory
It is instructive to compare the present theory with the standard Flory-Huggins theory for polymer blends.The Flory-Huggins free energy (per lattice site) for a blend of polymer A and B is written as 1 2 f where φ x is the volume fraction of component x, and a nondimensional parameter χ measures the nature and the strength of interaction.The incompressiblity condition enforces φ a + φ b = 1.Since χ is usually positive, corresponding to attraction among the like-species, such an interaction acts as a driving force to the phase separation.The free energy (8) reduces to that of polymer solution in the limit N b = 1 where the component b represents a solvent.Figure 4 shows the phase diagram calculated from Eq. ( 8).When we fix the interaction parameter at the value χ > χ c , where χ c = ( 1/Na + 1/Nb) 2 /2 is the critical value for the phase separation, the phase diagram as a function of φ a is one-dimensional line with two points φ a representing the phase boundaries.If the overall concentration falls in between these two points, the uniform state is metastable (outside spinodal region) or unstable (inside spinodal region) and the system phase separates into the A-poor (dilute) and A-rich (concentrated) phases with the volume fractions φ (1) a and φ (2) a , respectively.Note that the dimensionality of the phase boundary at fixed χ is d pb = 0, i.e., points, which is a consequence of the equality of number of conditions (n c = 2, i.e., µ a and µ a and φ (2) a ).We note that a set of conditions µ a and µ a and Π (1) = Π (2) , where is the osmotic pressure, the use of which may be more common in polymer solution, where the component b is regarded solvent.These two methods are equivalent due to the relation which follows from the incompressible condition 1 , where v 0 is the volume of the monomers and solvents.In contrast, in our description of polymer mixture with soft potential, the solvent degrees of freedom is already integrated out, and c a and c b are independent variables without constraint, i.e., free from the incompressible condition.One can conceive that our system under consideration is a three component system (two solute A and B plus solvent), and the free energy density (1) represents a mesoscopic description after coarse-graining.
It is also known that the interaction part in the Flory-Huggins free energy initially takes the form (χ aa /2)φ 2 a + (χ bb /2)φ 2 b + χ ab φ a φ b .Rewriting it into the form of Eq. ( 8) with a single parameter χ = χ aa /2 + χ bb /2 − χ ab is made using the incompressibility condition.Again, it does not apply to our soft polymer description.Since our system originally possesses three components, we naturally need three interaction parameters to characterize the system.We also note that the critical χ parameter in blends of long polymers N A , N B ≫ 1 is χ c → 0 in Flory-Huggins theory.A necessary condition for the phase separation in this limit is thus χ > 0 ⇔ (χ aa + χ bb )/2 > χ ab .In contrast, as discussed in Sec. 4, the corresponding condition in our soft polymer mixtures is χ ab > √ χ aa χ bb independent of the chain length 30 .

Relation with Gaussian core model
In our model, polymers are described as N successive soft beads, where these beads are already mesoscopic entities with their internal degrees of freedom integrated out.We note here that, starting from a microscopic model, there is a freedom to choose N, i.e., the degree of the coarse-graining.Although the extreme limit of the choice is N = 1, in which individual polymers are described as single soft particles, the validity of such a description may break down at high concentration 31 .
It is known that the effective pair potential between two isolated polymer coils in dilute solution can be well approximated by a simple Gaussian potential where ε ≃ 2 and the width R is of order of the gyration radius the coil 32 .The fact that the energy scale of the interaction is of (a) order of thermal energy indicates entropic nature of the interaction.It has been shown that the above potential also provides a reasonable description for the effective interaction in semidilute solution, where polymers are overlapped.Thermodynamic properties of fluid composed of soft particles interacting through Eq. ( 10), i.e., the Gaussian core model, has been analyzed in detail by Louis et.al. 32 .Phase separation in binary mixture of such fluids has been also extensively studied [27][28][29] .As discussed in Sec.One may then conclude that the introduction of the "polymerization index" N a and N b might be auxiliary for the description of homopolymer mixtures.However, there are, at least, two reasons we need the polymeric description.First, N = 1 description is known to suffer from a significant concentration dependence of the effective repulsive interactions once polymer coils start to overlap deep in semidilute concentration regime.This motivates the multisegment description with N > 1 31 , where a suitable choice for N would be guided by the overlapping condition for mesoscopic segments.Second, once there arises a characteristic length scale in the problem, we need the polymeric description as a succession of beads with appropriate degree of coarse-graining.In Sec. 5, we provide one such example, where we analyze the effect of modulation in local physical properties along polymers, (i.e., due to post translational modification) on the phase separation.
Another point deserving comment is the relation between the strength of the interaction potential ε xy and the interaction parameter χ xy in our free energy (1).We have shown in Sec.3 that the simulation results quantitatively match with the free energy prediction under the relation χ xy /ε xy ≃ 2.5.Since the free energy (1) takes apparently the same form as the virial expansion up to second order, one may expect that the ε xy − χ xy relation would be obtained from χ xy = − (e −U xy (r)/k B T − 1)d⃗ r.As emphasized in ref. 32 , however, the free energy ( 1) is based on the random phase approximation (RPA) 33 .As such, it becomes more and more accurate in higher concentration regimes in contrast to the second virial approximation 33 .In fact, unlike the virial expansion, the quadratic form of the free energy in concentrations is a consequence of the RPA closure, where the direct correlation functions, which appears in the Ornstein-Zernike relation, is independent of the concentrations.The analysis of equation of state with RPA leads to the identification χ xy = (1/k B T ) U xy (r)d⃗ r = π 3/2 σ 3 xy ε xy .Given the resultant ratio χ xy /ε xy = π 3/2 is considered to be an upper bound compared to a more accurate estimate e.g., obtained from hypernetted chain closure 32 , we find our result χ xy /ε xy ≃ 2.5 reasonable and providing an overall consistency of the soft core model description of the phase separation based on the free energy Eq.( 1).

Mixtures with random copolymers
So far discussed is a foundation for the coarse-grained description of polymer mixture, where polymers are represented as succession of soft mesoscopic beads.In particular, we have focused on the phase behavior of mixture of homopolymers.As stated in Introduction, however, one of the main motivations to necessitate such a description is its relevance to describe the large scale behavior of chromatin in cellular nucleus.In this section, we would like to discuss a simple extension of our theory, which may be linked to a certain aspect of chromatin organization in living cells.
It is known that interphase chromatin in early embryo is quite homogeneous inside nucleus, which is, in certain sense, reminiscent to a uniform solution of homopolymers 34 .With the progress of the development stage, however, several characteristic structures, such as heterochromatin foci and transcriptional factories, start to appear 35 .Responsible for such structure formations would be a phase separation, which is driven by local alternation of chromatin monomers caused by e.g., post-translational modification.
The change in the chemical state in chromatin monomers likely induces the modulation of physical properties along chromatin polymer, which could be represented by a copolymer model.Since the variation in repulsive forces primarily reflects the difference in density of core-bearing monomers within the coarse-grained segments 20,36 , the segment "a" represents regions where chromatin exists in a relaxed, less condensed state, reminiscent of euchromatin, while the segment "b" corresponds to more condensed regions akin to heterochromatin.The structure formation under consideration could thus be treated by the appearance of copolymers in the matrix of homopolymers.With this in mind, let us consider a mixture of homopolymers H (with length N h ), which composes of type a beads only, and copolymer C (with length N c ) , which composes of types a and b beads.The monomer concentration of H and C polymers are c h and c c , respectively.For the analytical tractability in a simple mean-field description, we assume the latter to be a random copolymer, which is characterized by the fraction α of b beads, i.e., the number of b beads in a copolymer C is αN c .
The free energy of the mixture is written as which takes the same form as Eq. ( 1) except for the appearance of new interaction parameters.While χ hh = χ aa trivially from the definition of the homopolymer H, the others χ cc and χ hc are nontrivial, which appear instead of χ bb and χ ab , respectively.Given the randomness in the sequence of the copolymer, we can evaluate these interaction parameters as mean values of the inter-beads interactions χ aa , χ bb , χ ab ;  In Fig. 5, we show phase diagrams obtained from the free energy (11) with Eqs. ( 12) and ( 13) for a fixed α.Note that α = 0 reduces to a homopolymer solution (with only type a bead), and α = 1 corresponds to a blend of homopolymers A and B analyzed in earlier sections.Here we show the cases for α = 0.3, 0.5.As expected, the region for phase separation enlarges with the fraction α.In addition, the results depend on a relative stiffness between beads a and b.As shown, the system is more prone to phase separation when the matrix polymer is softer χ aa < χ bb reflecting the asymmetry in the phase diagram of homopolymer mixtures (Sec.3).
To check the validity of the free energy prediction, we again performed numerical simulations using Gaussian potentials to represent bead-bead soft repulsions.In Figs. 6 and 7, we compare the theoretical phase diagram in Fig. 5 with simulation results, where the repulsion strengths for Gaussian potentials are set to be ε xy = χ xy /2.5 between beads x and y as determined from the result of homopolymer mixtures.As shown, the agreement is rather sat-isfactory, demonstrating that the overall trend of phase separation is well captured by the proposed free energy.In Figs. 6 and 7, we also show the spatial profiles of the monomer concentration c h of homopolymer together with the corresponding typical snapshots.

Discussions and Summary
Numerical simulations of large scale chromatin organization in cellular nucleus often adopt highly coarse-grained models, in which chromatin polymer is represented as a succession of soft beads [20][21][22][23][24][25] .Unlike a model which employs nucleosomes as monomers of chromatin polymer, each of beads here represents a substantial amount of nucleosomes, thus regarded as a mesoscopic entity, allowing mutual overlaps with entropic penalty.It has been shown that the effective interaction between such mesoscipic segments is soft and repulsive, qualitative feature of which is well approximated by the Gaussian potential 20 .
We have considered binary mixtures of such soft repulsive polymers and investigated how the imbalance in repulsive interactions between different species leads to the phase separation.After summarizing universal aspects of the phase diagram based on invariant property of the free energy upon changes in parameter values, we have extended the theory to mixtures including random copolymers, which may have some implications to chromatin phase separation.
As discussed in Sec. 5, random copolymer model is inspired by epigenetic modification of chromatin.This modification is performed and maintained by enzymes, thus includes energy consuming nonequilibrium process.In this sense, our description based on equilibrium framework should be considered a useful effective description to elucidate the impact of phase separation on chromatin organization.The same remark applies to many of current biophysical modelings of chromatin, not only for its structural organization but also for its dynamics.Yet, there are several other works, which emphasize possible impacts of various nonequilibrium effects on chromatin 18,[38][39][40][41][42][43][44][45] .Perhaps, some of these effects associated with nonequilibrium activities could be described by effective equilibrium models.Such a strategy may well work to understand some aspect of the problem, but may fail to capture other aspects.In our opinion, it remains to be seen how and when nonequilibrium factors are critically important in chromatin biophysics.The same comment would apply to topological constraints, another factor presumably important in chromatin, but not explicitly included in our description 46,47 .In this regard, it is interesting to note that as discussed in 42 , the similar physics as described in present paper may be important in blend of non-concatenated ring polymers, where the soft repulsion arises from the so-called topological volume due to topological constraints [48][49][50] .
As possible extensions of our work, we first note that our theory deals with the macro-phase separation, hence does not capture the possible appearance of mesophases.However, the occurrence of the "micro-phase separation" is naturally expected in copolymer systems, the elucidation of which should provide further insight into the problem of chromatin organization in nuclei.Secondly, although we have only analyzed bulk property based on meanfield theory, we expect that the effect of correlations and interfacial properties at phase boundary and near a confining wall could be analyzed by following an approach outlined in ref. 30 .It would be interesting to see how such an analysis can be compared to the chromatin spatial profile, e.g.near the nuclear membrane.
Finally, we point out that there are several studies on compressibility effects in polymer solutions, which become evident, for instance, in pressure-induced phase separation 51 .Here, interesting phenomena such as the acousto-spinodal decomposition have been predicted 52 .Although comparison with these studies may be interesting, we note that the loss of incompressible condition in our description results after coarse-graining, i.e., integrating out solvent degrees of freedom.Therefore, to address the kinetic effect, we need to take properly solvent effects into account.initial spatial concentration profile to prepare the phase separated initial state.
To perform various statistical analysis, we sampled microstates of the system every 1000 steps after the system reaches equilibrium, i.e., at 2 × 10 6 steps (simulation runs in Sec. 4) and 3 × 10 6 steps (Sec.5) after setting the interaction strengths to appropriate values.However, for the case of simulation runs given in Sec. 4 at (c (0) a , c (0) b ) = (0.3, 0.3), the sampling starts at 9 × 10 6 steps, and for the case of Sec. 5 at α = 0.3, (ε aa , ε bb , ε ab ) = (2.5, 0.5, 1.5), (c (0) h , c (0) c ) = (0.5, 0.5), the sampling starts at 7 × 10 6 steps.In all simulations, we collect 1001 independent samples of particle configurations in equilibrium.For the simulation at (c (0) a , c (0) b ) = (0.3, 0.3), as only one exception, production runs end at 18 × 10 6 steps, and particle coordinates are sampled every 1000 steps from 9 × 10 6 to 18 × 10 6 steps, during which 9001 independent samples of particle configurations in equilibrium are collected for statistical accuracy improvement in the vicinity of the critical point.We have confirmed that the physical properties of the simulation system are not significantly changed when we start the simulation from the homogeneously mixed initial states.

Conflicts of interest
'There are no conflicts to declare.
c b ), leading to the number of condition to determine the phase equilibria n c = 3. Comparing the number of unknowns and that of conditions, we expect that the dimensionality of the phase boundary, i.e., binodal is d pb = n u − n c = 1. 3 Phase diagram of the soft repulsive polymer mixture In Fig. 1 (a), we show an example of the miscibility phase diagram obtained from the free energy Eq.(1) under the fixed interaction parameters.As we have discussed, the phase diagram is twodimensional spanned by c a and c b , in which the uniform state (bottom left) and the demixed state (upper right) are separated by onedimensional phase boundary.Tie lines, which connect (c b ) in demixed state, are negatively sloped, indicating the phase separation is a segregative type.As expected, the region for demixing widens with the increase in either chain length or the repulsion strength, see Fig. 1 (b).Note that despite the symmetry in the chain length N a = N b in the examples shown here, the phase diagram exhibits the asymmetry about the diagonal c a = c b .The asymmetry is caused by the difference in physical properties of type a and b beads, which leads to the phase rich in softer beads b being more concentrated than the other.Such a feature can be made more evident by re-plotting the phase diagram in total concentration c = c a + c b and the composition ψ = c a /c plane (Fig. 1 (c)).
) is obtained with the interaction strength ε aa = 2, ε bb = 1, ε ab = 1.5 and the chain length N a = N b = 20, where the overall concentrations are varied from (c

Figure 1
Figure 1 Miscibility phase diagram obtained from the free energy (1).(a) Binodal (solid curve) and spinodal (dashed curve) in the case of repulsion parameters χ aa = 2, χ bb = 1, χ ab = 1.5 and the chain length N a = N b = 20.Two curves meet at the critical point marked by a circle.Some examples of tie lines are also shown.(b) Dependence of phase diagram on chain length and repulsion parameters, where binodals obtained from different set of parameters are shown.(c)Replot of the phase diagram (a) in the plane spanned by total concentration c = c a + c b and the composition ψ = c a /c.

Figure 2
Figure 2 Phase boundary (solid circle) of miscibility phase diagram obtained from numerical simulation with interaction parameters ε aa = 2, ε bb = 1, ε ab = 1.5 and the chain length N a = N b = 20.Open symbols (square, circle and triangle) indicate the overall concentrations of the system simulated; square or triangle indicates that the uniform state is stable or unstable, while circle represents the vicinity of the critical point.Overlapped solid curve is the binodal curve obtained from mean field theory with χ aa = 5, χ bb = 2.5, χ ab = 3.75.(b) Spatial profile of concentration of a-bead under various overall concentrations.(c)Histogram of number density of a-bead, where horizontal axis represents the deviation from the uniform cooncentration c a − c (o) a .
b ) and the number of unknowns (n u = 2, i.e., φ

Figure 4
Figure 4 (a) Miscibility phase diagram for a polymer blend obtained from Flory-Huggins free energy Eq.(8).(b) Under a fixed χ-parameter (χ > χ c ), the phase boundaries are points (d pb = 0) on a line.

4 . 1 ,
our free energy (1) can formally be mapped to that case by (N a , N b , c a , c b ) → (1, 1, c a /N a , c b /N b ).The analysis in Sec.4.1 may then indicate that the miscibility phase diagram is intact if we simultaneously transform the interaction parameters as (χ aa , χ bb , χ ab ) → (χ aa N 2 a , χ bb N 2 b , χ ab N a N b ).

Figure 7
Figure 7 Quantitative comparison of theoretical phase diagrams against numerical simulation.Parameters are the same as those in Fig. 6 except for the fraction of b beads in copolymer, which is α = 0.5 here.These conditions correspond to those in Fig. 5 (c) and (d), recpectively.(Top) Binodal (solid curve) with the critical point marked by solid circle.Open symbols indicate the overall concentration (c (0) h , c (0) c ) adopted in numerical simulations performed under the parameter correspondence ε xy = χ xy /2.5, where open squares and triangles, respectively, indicate the one-phase and two-phase regions.In the latter, concentrations after phase separation are shown by solid triangles, which are connected by tie lines.(Middle) Spatial profile of monomer concentration of H-polymer.(Bottom) Typical snapshots obtained in simulations with (c (0) h , c (0) c ) = (0.2, 0.2), where red beads represent type-a beads contained in H-polymer, while yellow or blue beads represent type-a or type-b beads contained in C-polymer, respectively.The snapshots were rendered using OVITO 37 .