Phase Separation Driven by Active Flows

We extend the continuum theories of active nematohydrodynamics to model a two-fluid mixture with separate velocity fields for each fluid component, coupled through a viscous drag. The model is used to study an active nematic fluid, mixed with an isotropic fluid. We find micro-phase separation, and argue that this results from an interplay between active anchoring and active flows driven by concentration gradients. The results may be relevant to cell-sorting and the formation of lipid rafts in cell membranes.

Introduction.-There is increasing evidence that phase ordering plays a fundamental role in biological processes such as morphogenesis and colony growth [1][2][3][4]. Phase ordering of cell types (Fig. 1) [5], membrane-less organelles [6] and RNA-protein mixtures [7] inside cells have been studied using free-energy prescriptions. However, biological matter is intrinsically out of thermodynamic equilibrium, and it is unlikely that biological phase ordering relies entirely on the framework of equilibrium thermodynamics.  [8]). (d) Phase ordering of extensile (magenta) and contractile (green) epithelial MDCK cells (courtesy of L. Balasubramaniam [9]).
Self-motile particles form an important class of nonequilibrium systems called active matter [10]. Active matter exhibits non-equilibrium phase segregation mechanisms, such as motility-induced phase separation (MIPS) [11][12][13] where active particles get trapped in regions of high particle density. Phase separation has been found in continuum models [14][15][16], in the presence of hydrodynamic interactions [17], or driven by aligning torques [18,19]. Active nematics consist of active rod-like particles with orientational order. The elongated particles generate active dipolar stresses which destabilise the aligned state and lead to active turbulence characterised by a chaotic and highly vortical velocity field [20][21][22]. Recent work has shown that bacterial biofilms [23,24], epithelial tissues [25,26], and microtubule-motor mixtures [27] can be modelled as active nematics. There is little work about phase ordering in active nematics [28,29], but recently microphase separation was predicted to occur in extensile gels near the nematic-isotropic transition point [30].
It is still not clear whether and how active stresses can lead to phase ordering in a binary mixture of an active nematic and a passive fluid in the absence of any thermodynamic driving force. To answer this question, we extend the continuum theories of active nematohydrodynamics to a two-fluid model with separate velocity fields, coupled through viscous drag, for each fluid component [31,32]. This allows us to model relative motion between the two fluids which is a requirement for activity-driven phase ordering. Two-fluid models that have been used in biological contexts include studies of cytoskeleton dynamics and growing biofilms [33][34][35].
Model.-We model a binary mixture of an active nematic fluid (component 1), and a passive isotropic fluid (component 2), in d dimensions. Each fluid component, denoted by i ∈ 1, 2, has a density field ρ i , velocity field u i and a chemical potential µ i . The active nematic is associated with a second rank tensor order parameter field [36], where S denotes the degree of nematic ordering, and n denotes the orientation of the local nematic director field.
Each component fluid i obeys the mass continuity equation and the momentum balance equations where ϕ = ρ 1 /ρ c is the concentration of the active fluid, and ρ c = ρ 1 + ρ 2 is the total fluid density. The forces acting on the fluids, on the rhs of Eq. (2), are thermodynamic forces, viscous drag between the component fluids, internal viscous dissipation σ visc,i = η i (∇u i )+(∇u i ) T − 2 d (∇ · u i )I , and, for the active nematic fluid, elastic and active forces [20], respectively. The stresses in the nematic are given by where F is the free energy, H is the molecular field [20] defined by H = − ∂F ∂Q + I d Tr ∂F ∂Q , λ is the flow aligning parameter, and ζ is the magnitude of the activity. Eq.
The nematic tensor field evolves according to [37] ∂ t Q + u 1 · ∇Q − S = ΓH (6) where S = (λE 1 + Ω 1 ) · Q + I/d + Q + I/d · (λE 1 − Ω 1 ) − 2λ Q + I/d (Q : ∇u 1 ) is the co-rotation term describing the response of the orientation field to strain and vorticity in the flow, while the right hand side describes the relaxation to a state of minimum free energy. Here, E 1 the rate of strain tensor and Ω 1 the vorticity tensor of the active fluid. Equilibrium is described by a free energy [20,21,28,38] where the nematic elastic constant K LC , C, a and κ are material parameters. The first term represents the pressure contribution to the free energy from the isothermal equation of state for the fluid [39]. The rest of the first line describes a Landau-Ginzburg free energy with a > 0 favouring a homogeneous mixed state (ϕ = 1/2) and κ > 0 penalizing interfaces. The second line is the Landau-de Gennes free energy describing an ordered nematic.
We change variables to give one equation for the combined fluid, which is to a good approximation incompressible, and one for the active fluid. To do this, we define the centre of mass velocity of the combined fluid u c = ϕu 1 +(1−ϕ)u 2 , and the relative slipping velocity between the fluids δu = u 1 −u 2 . Assuming that the viscous drag between the components is the fastest relaxation process in the system, δu << u 1 , u 2 , we neglect all terms of order (δu) 2 , and assume that the drag between the fluids is linear in δu. Moreover, recalling that µ i = ∂F/∂ρ i , the thermodynamic forces can be written as a stress tensor [38][39][40]. Assuming that both fluids have the same viscosity η, the internal viscous dissipation of the combined fluid is F visc,c = η c ∇ 2 u c where η c = 2η.
Adding Eqs. (1), (2) for each component then shows that the combined fluid satisfies The equations for the first (active) component are where G = γ(u c − u 1 ) + (1 − ϕ)∇δµ is the force applied by the passive component, and δµ = −δF/δϕ. Eqs. (8)(9)(10)(11) are solved using a lattice Boltzmann (LB) algorithm [38,40]. For the equation of the Q field (6) we use a finite difference (FD) approach. 5 FD steps are taken for each LB step to optimize simulation speed. Simulations are run in a periodic box of size 200 x 200 for 20,000 LB time-steps, with a random initial director configuration. Parameter values are ρ 1 = 20, Phase separation driven by activity.-We solve the equations of motion for a mixture which comprises equal concentrations of an active nematic fluid and a passive fluid. The mixture is incompressible so that the total fluid density remains constant, but the relative concentrations of individual components can vary in space and time. The two components are initially homogeneously mixed so that the concentration of the active component ϕ = 0.5 everywhere. However, above a threshold activity the system quickly self-organises into dynamicallyevolving regions which are, respectively, rich or poor in the active fluid component. A typical snapshot of the phase-ordered state is shown in Fig. 2(a), and Movie 1 in the Supplemental Material (SM) [41] shows its time evolution. The domains continually break up and re-form. They are elongated by extensile active flows [28] and then pulled apart by the chaotic active turbulent flows .
We measure the standard deviation of ϕ, denoted by ∆, to quantify the degree of phase separation. Fig. 2(b) shows how ∆ evolves with time towards a dynamical steady state. tion is entirely due to the activity, and requires no passive forces. The free energy in Eq. (7) is chosen to be a single well potential, a > 0, which favours mixing. To explain the mechanism causing the phase ordering consider a small fluctuation in concentration. From Eq. (5) this will lead to an active force tangential to the interface between the higher and lower activity regions F tangential = 2ζ|∇(Sϕ)| (m · n)(l · n) l where m and l are unit vectors normal (pointing away from the more active region) and tangential to the interface respectively, and n is a unit vector along the nematic director (Fig 3(a)) [28], see also SM, Sec. 1 [41]. This force sets up flows parallel to the interface. The nematic director rotates in the flow leading to a net alignment parallel (perpendicular) to the interface in extensile (contractile) systems. This is termed active anchoring [28,42].
Similarly the active force normal to the interface is F normal = −ζ|∇(Sϕ)| (2(m · n) 2 − 1) m. In the absence of anchoring, F normal averages to zero. However, when there is an active anchoring of the director, the force drives the active component up the concentration gradient towards the active nematic region. Hence the relative flow between the two components enhances concentration fluctuations to form active regions of concentration ϕ > 1/2. The same mechanism leads to domain formation in contractile systems, as the change in the direction of the active anchoring is compensated by the change in sign of the activity. This is illustrated schematically in Figs.
Varying the Model Parameters.-We now discuss more detailed numerical results investigating how the size of the phase separated domains and the difference in concentration between the two phases, ∆, depend on the model parameters. Fig. 5(a) shows that ∆ increases, but does not scale linearly, with the activity coefficient ζ. This is because at higher activities the flow is more turbulent with many defects which reduce the strength of the active anchoring. The formation of concentration gradients is opposed by the bulk free energy which scales as a, the relative drag γ, and the surface tension κ. We consider the balance between the activity and each of these in turn, choosing parameters where the given relaxation mechanism is dominant. Results are presented for extensile activity but they also hold for contractile systems since the same mechanisms are at play (SM, Fig. A2 [41]).
In Fig. 5(b) ∆ is plotted as a function of ζ/a 0.5 showing collapse onto a single curve. The collapse is a result of the balance between the driving force normal to the interface and the thermodynamic force which restores mixing. The former scales as ζ 2 because it depends on the strength of the anchoring and on the magnitude of the active force, both of which scale with ζ. The latter scales as a.
The data also lies on a single curve when ∆ is plotted against the ratio ζ/γ (Fig. 5(c)). This shows the expected inverse relation between the activity, which creates a slipping velocity between the active and passive fluids, and the viscous drag that dampens the slip.
A similar plot, obtained by changing the surface tension coefficient κ at constant a and γ, shows the best data collapse when ∆ is plotted against ζ/κ 0.8 (Fig. 5(d)). We do not have a simple argument for this dependence as several factors contribute. In addition to the thermodynamic penalty, an increase in surface tension weakens the activity gradient by making the interface more diffuse, and hence will alter the details of the active interfacial forces and active anchoring. The length-scale of the domains can be calculated numerically from the first zero of the spatial correlation function of the fluctuations in the concentration field. For fixed surface tension, Fig. 5e shows that the domain size is controlled by the usual active nematic length scale L act = K LC /ζ. This is, however, not the only length scale controlling the physics. The domain size can be changed by varying a second length scale, related to the surface tension, L surf = κη/ζ, as shown in Fig. 5f. Active droplets break up by forming bend or splay instabilities with typical lengthscale L act , or by self-shearing instabilities [28,43] with typical lengthscale L surf . Snapshots showing formation and dissolution of a droplet are shown in SM, Fig. A3 [41].
Threshold for phase ordering.-It is possible to obtain an exact expression for the activity threshold above which phase ordering occurs, in the limit of perfect surface anchoring, and ignoring interface curvature. The effects of curvature are discussed in Sec. 2 of the SM [41].
Consider a concentration gradient around ϕ = 1/2 in the x-direction, given by ∇ϕ =| ∇ϕ |x. Let the activity coefficient be ±|ζ|, where the + (−) sign denotes extensile (contractile) activity. Assuming strong active anchoring at the interface, this leads to a director field Q xx = ∓S nem √ ϕ, Q xy = 0. The normal active force at the interface is then F act = (3/2)S nem √ ϕ |ζ | |∇ϕ|x, and the restoring force driven by the free energy is F rest = −2aρϕ(1 − ϕ)|∇ϕ|x. This leads to an expression for the slip velocity δu x in the Stokes limit: The concentration gradient will increase for activities greater than |ζ| crit = 4aρϕ(1−ϕ) Fig. A4 [41].) Summary & outlook.-By introducing an active two-fluid model with viscous drag between the fluids we have found that a homogeneous mixture of an active nematic fluid and a passive fluid spontaneously phase orders into regions with different concentrations of each fluid component. This is driven dynamically by active flows set up at concentration interfaces, and is an example of liquid-liquid phase separation [44] out of thermodynamic equilibrium.
The two-fluid model reduces to the lyotropic formulation of nematohydrodynamics that has been successful in defining active and passive regions within a simulation, in the limit of infinite viscous drag between the fluids. However, because the lyotropic approach uses a single velocity field for both phases [28,45], it cannot account for relative flows between, or compressibility of, the two components.
Active nematics have been observed to preferentially adhere to a substrate or container wall, a phenomenon known as active wetting [46][47][48][49]. A straightforward extension of our argument shows that for an initially homogeneous active-isotropic mixture, the normal force generated at a boundary would cause wetting of the boundary by the active nematic for both extensile and contractile activity, even in the absence of surface thermodynamic forces. The wetting would then be enhanced by any thermodynamic planar (normal) anchoring for extensile (contractile) activity.
In future work it will be interesting to study mixtures where both the fluids are active. In addition to cell sorting this is relevant to the growth and dynamics of bacterial and algal colonies which contain more than one species [50]. Moreover, the relaxation of the incompressibility constraint will allow the study of how compressible flows affect defect behaviour, and whether they can stabilize structures such as asters or the cellular lumen.