Fig1.A: mixed binary mixture, rounds with different colors denote different molecules. Fig1.B: demixed binary mixture, one component forms a dense liquid droplet

Things might get a little different when interaction among molecules are considered. Under the conditions of constant temperature, volume and particle numbers, the system is at equilibrium when the Helmholtz free energy F is the lowest. Based on regular solution model, the free energy density f takes the following form in the unit of k_BT[1]:

f =ϕ ln⁡ϕ +(1-ϕ )ln⁡(1-ϕ )++χϕ (1-ϕ )+λ/2 |∇ϕ|², (1)

where ϕ is the volume fraction of one component (let us say component A), χ is a parameter characterizing the strength of intermolecular interactions and λ is related to the surface tension between interfaces. The volume fraction of A is defined as the volume of A molecules divided by the total volume of the system. In a binary system, the volume fraction of the other component, let us say component B, becomes naturally 1-ϕ.
First let’s focus on the symmetric part of f , i.e f₀ =ϕ ln⁡ϕ+(1-ϕ)ln⁡(1-ϕ)+χ ϕ(1-ϕ), and see how its shape changes as we vary χ. When A and B attract, χ is less than 0; on the contrary, when A and B repulse, χ is greater than 0. Easily seen from fig2, when χ<2, f₀ only has one minimum; when χ>2, f₀ has two minima and one maximum. A bifurcation takes place when χ=2 and essentially alters the feature of free energy density. Fig2: plot of symmetric part of free energy density under various χ

For better illustration, we scrutinize two typical cases, χ=0 and χ=4. Fig3 shows the free energy density for mixed state in blue solid line and demixed state in green dotted line. Obviously seen from fig3, when χ=0, for any initial concentration represented by ϕ₀, it always requires extra free energy for the system to demix into any two separate states ϕ₁ and ϕ₂, where the green dotted line is higher than any point on the blue solid line between ϕ₁ and ϕ₂; when χ=4, there exists a range of ϕ₀ to separate into two demixed compositions ϕ₁ and ϕ₂, where the green dotted line is lower than any point on the blue solid line between ϕ₁ and ϕ₂. This is the situation where phase separation could happen spontaneously.
Fig3: The blue solid line and the green dotted line represent the free energy density for mixed state and demixed state respectively. (A)χ=0. The green dotted line is always higher than the blue solid line, indicating an extra energy for separation; (B)χ=4. The green dotted line is below the blue solid line, making spontaneous phase separation possible.

To be more precious, we can specify the conditions where separation could happen. According to fundamental work on liquid-liquid phase separation, when d²f/dϕ²<0, any local perturbation will result in spontaneous separation, such a formation is named spinodal decomposition; when d²f/dϕ²>0 and between the two minima, only sufficiently large global perturbation can make phase separation happen, such an approach is called nucleation; and their boundary is named spinodal line. Now, if the free energy density function is symmetric, when ϕ lies outside the two minima, phase separation cannot happen. The boundary between whether phase separation can take place or not is called binodal line.

Based on the criteria above, we plot the phase diagram of a binary mixture. In fig4, we represent the initial concentration by ϕ in the x-axis and vary χ in the y-axis. The region confined by spinodal line is the unstable region, for separation can take place under any local perturbation. In contrast, the region between the binodal line and the spinodal line is the metastable region, where only sufficiently large global perturbation can initiate separation.
Fig4: Binary phase diagram. Binodal line shows the boundary between whether phase can separate or not. Spinidal line shows the boundary between two different formations: spinodal decomposition and nucleation. The area confined by spinodal line is the unstable region while the area between bimodal line and spinodal line is the metastable region.

As aforementioned, our system is a ternary mixture system consisting of two multivalent proteins and water, which is a bit more complicated. To capture the basic feature of three-component phase separation, we use a similar theoretical model for simulation. The free energy density f is now written in the unit of k_B T as:

f=ϕ₁ ln⁡ϕ₁ +ϕ₂ lnϕ₂ +ϕ₃ ln⁡ϕ₃ +χ_(1-2) ϕ₁ ϕ₂+χ_(1-3) ϕ₁ ϕ₃ +χ_(2-3) ϕ₂ ϕ₃ +λ_1/2 |∇ϕ₁ |^2+λ_2/2 |∇ϕ₂ |^2+λ_3/2 |∇ϕ₃ |^2 where ϕ₁ denotes the first multivalent protein (FKBP for example), ϕ_2 denotes the second multivalent protein (Frb for example) and ϕ₃ denotes water. An intrinsic relation of the three is given by: ϕ_1+ϕ_2+ϕ₃ =1 Phenomenologically speaking, ϕ₁ and ϕ_2 condensate together and separate from ϕ₃ . For computational convenience, χ_(1-2) and χ_(1-3) are taken value above 2 and χ_(2-3) below 2. A ternary phase diagram is calculated in a similar way by determining whetherϕ₁ and ϕ₃ separate and whether ϕ_2 and ϕ₃ separate. The results are shown in fig5.
Fig5: Ternary phase diagrams. (A)χ_(1-2)=0; χ_(1-3)=8; χ_(2-3)=3; (B) χ_(1-2)=0; χ_(1-3)=8; χ_(2-3)=4.

The phase diagram only provides a rough approximation on where phase separation can happen, it is insufficient to predict what happens after the separation. Hence, we further recur to the continuum model first proposed by Cahn and Hilliard to simulate a dynamic process.
∂c/∂t=∇⋅(M∇μ) μ=df/dc-λ∇^2 c

We use finite element method in a 100×100 mesh and select Neumann boundary condition to solve the partial differential equations above. Crank-Nicolson method is used for time-stepping with a footstep of 2.0×〖10〗^(-3). The initial composition is given by adding a perturbation of strength 〖10〗^(-2) to a homogenous sate. A typical result is given as followed:
Fig6: Simulation for dynamic evolution of phase separation under χ_(1-2)=0, χ_(1-3)=3, χ_(2-3)=4. (A) The concentration distribution of ϕ_1; (B) The centration of ϕ_1 on the sampled red dotted line in (A).

We further adjust the interaction strength between two proteins, thus affecting both χ_(1-3) and χ_(2-3). As χ increases, indicating a stronger interaction, the time for phase separation to occur is decreased, which is in accordance with our experimental results.
Fig7: Simulation for phase separation under different interaction strengths. The stronger the interaction, the faster separation emerges.