Self-organization and shape change by active polarization in nematic droplets

Active forces occurring within cells can drive crucial biological processes that involve spontaneous organization and shape change, such as cell division. Motivated by recent in vitro experiments of nematic droplets of cytoskeletal filaments and motors that self-organize and divide, we present a minimal hydrodynamic model that combines the nonequilibrium kinetics of motor-filament interactions with equilibrium nematic phase separation. The motors organize within droplets and structure filaments into polarized aster defects. At large motor activity, they can even deform or divide the droplet, or form multi-aster chains of droplets. Our predicted phase diagram recapitulates these experimentally observed shapes.

Introduction. Active mechanical forces enable living systems, particularly animal cells, to move, change shape, organize components, and divide. Subcellular cytoskeletal assemblies, comprising polar filaments and molecular motors that transduce biochemical reactions to generate active mechanical forces, drive these processes [1]. Understanding the general physical principles of living matter provides insight into cell biology, as well as guides the engineering of artificial cells that exhibit spatiotemporal organization of components and spontaneous shape change characteristic of cell division [2]. Model in vitro systems of purified cytoskeletal proteins, which capture elements of cell biological phenomena with only a fraction of the biochemical complexity occurring in vivo, exhibit a rich array of collective phenomena [3,4] that motivates bio-inspired active matter theory [5].
Recently, phase segregated macromolecular droplets have emerged as model systems to investigate spatiotemporal organization in biological cells and phase separation has been proposed as a possible primitive means of subcellular organization in protocells [6]. These macromolecular liquids are typically composed of disordered proteins and nucleic acids, and consequently, internal droplet order as well as motor activity are absent. In contrast, recent experiments indicate that biopolymers with high aspect ratio, such as cytoskeletal filaments, can form phase separated droplets with orientational order, as in liquid crystals, because of the alignment of filaments in the dense phase [7]. This nematic order confers an equilibrium spindle shape to these droplets [8,9], known as tactoids, which arises from a competition of droplet surface tension, the tendency of the filaments to align with the interface, and elasticity arising from bulk nematic order [10,11].
The ordered structure in biomolecular fluids influences the emergence of collective phenomena in active systems [5]. When confined to a droplet, active forces lead to non-equilibrium phenomena such as droplet shape change [12][13][14][15][16][17][18][19], motility [20,21] and dynamics governed by the geometry of the confining droplet [22]. In addition to these active nematic fluids, the directed "walking" of motors on filaments, along with motor-based filament crosslinking, can lead to polar order, where filaments prefer to point in the same direction at the mesoscopic scale, and form self-organized defect structures, such as asters and vortices [23]. Such polar active states are fruitfully described by hydrodynamic theories [24][25][26][27][28][29][30]. In recent experiments, myosin motors were shown to self-organize at the midplane of the aforementioned actin-filament based nematic droplets [31]. When sufficiently active myosin motors are present, they deform the droplet, even splitting it into two. A simple free energybased model of a nematic droplet, considering only the mutual alignment of the filaments and motors and adhesion of the droplet to the motor complex, was invoked to capture these key behaviors, but relied on arguments specific to the shape and structure of both the motor complexes and the droplet [8]. On the other hand, active mechanical forces, specifically the directed sliding of filaments by motors leading to their sorting by polarity [1], are expected to play a role in the dynamics of the organization. The effects of such active emergence of local polar order and associated defects, within a nematic droplet with a preferred orientational axis, have not yet been theoretically explored.
In this Letter, we combine continuum modeling that describes the structure of equilibrium nematic droplets with an active mechanical model for how motors move and slide filaments according to polarity, while selforganizing into localized asters. Using numerical simulations, complemented by theoretical analysis, we show that the resulting motor-filament self-organization destabilize droplets, giving rise to a rich array of experimentally observed structures including deformed, divided and multi-lobed droplets that can be generated by tuning one of two motor activity-dependent parameters.
1D polarity sorting model. To build up intuition on the role of active forces in motor self-organization, we first 2 examine a simplified one dimensional setup of motors interacting with filaments, as sketched in Fig. 1(i). Here, red lines depict actin filaments, blue circles show myosin II motors, while the black arrows indicate the motors respective direction of motion. This scenario arises, for instance, in contractile actin bundles [32,33] and when actin filaments are locally oriented along a nematic director [34]. Since motors walk towards a specific end of the polar filament (the barbed or plus end), we consider a density of left and right pointing filaments, denoted by n + and n − respectively. Momentum balance for the motor-filament system requires that a right (left) pointing filament is pushed to the right (left) by the motor as the motor moves in a single direction along the polar filament. The active motion of motors of density m thus gives rise to mass fluxes of both motors and filaments. In the dilute regime, the fluxes of right and left pointing filaments can be expressed as J + = ζn + m and J − = −ζn − m respectively, where ζ is a parameter related to the active motor force. The net flux of motors is written as −v 0 (n + −n − )m, where v 0 is the self-propulsion velocity of motors.
Including the diffusion of filaments and motors with coefficients D and D m respectively, the flux conservation equations can be written in terms of filament density ρ = n + + n − and the 1D polarization p = (n + − n − )/ρ 0 of filaments, where ρ 0 = ρ is the average density, giving FIG. 1. Illustration of steps in the active polarity-sorting model. (i) Motors "walk" towards barbed ends of filaments which they slide in the opposite direction as a result of momentum conservation. (ii) Polarization and (iii) motor density at steady state as a function of position in a one-dimensional model (Eqs. S7-S6) that shows motor localization to center of bundle by polarity sorting. (iv) In the 2D nematic droplet model, the filaments are initially oriented along the long axis without any polarity preference. Motors bind to filaments and advect them according to their polarity. (v) At steady state, motors gather at the droplet midplane and sort filaments into an aster that pinches the droplet.
The steady-state solution of Eqs. (2)-(3), (see SI I-II), is the one dimensional equivalent of an aster, where the populations of right and left pointing filaments are completely sorted to the right and the left of the aster, and is shown in Fig. 1(ii)-(iii). Note that in this minimal model, we have not considered the usual role of myosin motors as active crosslinkers which create pairwise forces on anti-parallel filaments, leading to self-straining flows proportional to local polarization [35]. This simple model reproduces the key experimental observation that motors migrate to the center of the droplet, and predicts that they induce strong polarity sorting. 2D nematic droplet with motor-induced polarization. To explore this prediction of active polarity sorting in nematic droplets and its implications for droplet shape, we build a more realistic hydrodynamic description for a suspension of filaments and motors on a frictional substrate that damps out large-scale fluid flows. In contrast with a thin fluid film, the filaments we seek to describe aggregate into droplets with free interfaces that separate the high density nematic from the low density isotropic phases, depicted in Fig. 1(iv). This is conveniently described by a nondimensional "phase field" corresponding to filament density, ψ, where ψ > 0 ( ψ < 0) describes the interior (exterior) of the droplet. Filaments in the high density droplet interior align in orientation, described by the 2D nematic order parameter, Q ij . These two ingredients result in nematic droplets with tactoid shape at equilibrium [36]. The active motion of filaments and polarity sorting induced by motors results in a net polarization within the droplet. Unlike in the 1D model, polar order in 2D is non-conserved and can be induced by motor-driven torques or relaxed by rotational diffusion. Generalizing Eq. (1) -(3) to 2D, and observing usual principles of conservation and symmetry [5], we obtain the dynamical equations, Eq. (4)-(5) describe the conservation of motors and filaments, respectively, and include both active and passive fluxes on the right side. The motor flux includes diffusion in the first term, and active motion of motors with a propulsion velocity v 0 in the second, Additionally, we include binding-unbinding kinetics for the motors (third and fourth term in Eq. (4)), where motors bind with rate k on wherever filaments exist (expressed by θ(ψ), the Heaviside step function) from a "reservoir" of free motors in solution, and unbind with rate k off , but are not restricted to diffuse within the droplet. Eq. (5) includes a flux created by a free energy of inter-filament interactions (full form given below) as well as active flux induced by motors advecting filaments in their polarity direction. Eq. (6) describes the induction of polarization by torques caused by gradients in motor and filament density, and its relaxation through rotational diffusion. Seen in the framework of the Toner-Tu hydrodynamic theory that describes the flocking of polar active particles [37], where p has the status of both an orientational order parameter as well as a fluid velocity, the −ζ 0 ∇ψ gives the gradient of pressure, and −ζ p ∇m is the gradient of an active stress created by motors [38]. Terms similar to this latter also arise in the theory of chemotactic colloids [39] and equilibrium polar liquid crystals [40],in both of which cases a chemical concentration can guide polarization. . Microscopically, the ζ p term captures the preferred orientation of filaments towards regions with higher motor density by the binding and crosslinking by motors, while ζ originates from the sliding forces exerted by motors. Note that the ζ specifies both polarization and active flux in the one dimensional model, whereas the ζ p term arises in the two dimensional model because of the additional orientational degree of freedom. The parameters v 0 , ζ and ζ p then all depend on motor activity but can in principle be varied independently by tuning motor properties such as size, shape, processivity and crosslinking. We assume in Eq. (7) that the nematic order is strong and arises from equilibrium forces. The timescales for relaxation towards equilibrium are specified by the "mobility" coefficients, Γ ψ , Γ p and Γ Q .
The equilibrium dynamics assume a coupled phase transition in ψ and nematic order Q, and a relaxation of p included in the total free energy, The density free energy Eq. (8b) models phase separating droplets according to standard Cahn-Hilliard dynamics. Ignoring corrections for curved interfaces [41], we use droplet surface tension ("line tension" in 2D) and its interfacial width, ξ = 2B/ν 2 [42]. The free energy for the polarization, Eq. (8c), includes two relaxation terms, α p , β p > 0 (corresponding to lack of spontaneous polar order) and an elastic term, κ p . Eq. (8d) is the Landau-de Gennes free energy for the nematic order, Q ij , with elasticity κ Q . All equilibrium couplings between fields are written in Eq. (8e), where the first term controls the density-driven isotropic-nematic transition and induces nematic order within the droplet. The second term is a "weak anchoring" that aligns the nematic parallel to the droplet interface for A > 0. The third term aligns the polarization with the nematic or-  der. Overall, this model recapitulates the key elements of nematic phase separation of the filaments and their coupling with motor activity. Results. We now employ numerical simulations to explore the consequences of motor activity on the dynamics and morphology of phase-separating nematic droplets. Starting from an initially circular droplet of radius R 0 , we integrate Eqs. 4-7 until a non-equilibrium steady-state is reached (see SI III for simulation details and parameters). Figure 2 presents typical simulation results. We obtain an elongated droplet with the nematic aligned along its long axis, as shown by the density and nematic director plots for a typical case in Fig. 2(i). The motor density accumulates at the core of the aster it induces with an outward polarization as shown in Fig. 2(ii). We first explore the interplay between motor-generated active forces, which tend to distort the equilibrium structure, and surface tension which resists such deformation. To this aim, we perform simulations with varying mo- surface tension ∞ 20 40 60 ≥ p theoretical analysis single tactoid, one aster two connecting tactoids, one aster three connecting tactoids, two asters full division (iv) three tactoids tor induced polarization, ζ p , relative to surface tension, γ. The resulting time sequences of observed droplet shapes, starting from an unpolarized droplet, are shown in Fig. 2(iii). When ζ p is low compared to surface tension, Fig. 2(iii a), motors localize towards the center to form an aster, which only slightly polarizes and deforms the droplet. At intermediate ζ p (iii b), motors localize more strongly, resulting in a stronger elongation and the appearance of a constriction of the droplet around the midplane between strongly polarized lobes. Increasing the ζ p further can lead to two distinct scenarios: In (iii c), the aster divides and a third, central lobe with strong polarity gradient emerges between two constrictions. In (iii d), finally, ζ p is sufficiently strong to induce full division of the initial droplet into two polarized daughter droplets. The short-time dynamics of this model thus reproduces the motor centering, aster formation and polarity sorting features of the 1D model (Fig. 1), while at longer times a complex diversity of morphologies emerges from the interplay between surface tension and motor activity.
To rationalize this rich phenomenology, we build a morphological phase diagram in Fig. 3(i) in ζ p -γ while keeping other parameters constant. Each distinct steady state structure is shown in Fig. 3(ii)-(v) (top). Interestingly, each of these morphologies corresponds to experimentally observed shapes, as shown in Fig. 3(ii)-(v) (bottom). Confirming the qualitative findings described in Fig. 2(iii), at medium to high surface tension and low ζ p , we find motors form asters in the midplane of the undeformed droplet (blue stars), whereas higher ζ p increases the influence of the centered aster on droplet shape, pinching it into two lobes (orange circles). This is consistent with the experimental observation that motors always localize to droplet midplane, but only deform the droplet when there are more active motors [31]. Qualitatively, the motors at the aster core splay the filaments, which are anchored to the interface, and which therefore results in an inward pinching of the interface. The motor induced polarization (second term in Eq. 6) is derivable from a functional, −(ζ p /Γ p ) d 2 x m∂ i p i , corresponding to an effective free energy for motor-induced spontaneous splay. This corresponds to a negative contribution to surface energy [43] (see SI IV for the derivation), that drives the droplet shape instability when ζ p m 0 /Γ p > γ . We thus expect the transition between undeformed and pinched droplets to occur across a dividing line, ζ p ∼ γ, which is confirmed by our phase diagram data in Fig. 3(i).
By further increasing ζ p and lowering surface tension, we find droplets with two asters which pinch them into three-lobed structures (Fig. 3(iv), green upward triangles). Importantly, we also find experimental realizations of such multiply pinched, equi-lobed structures, shown in Fig. 3(iv) bottom panel. Such linear "strings of tactoids" connected by multiple motor clusters have not been previously reported and are evocative of fibers with periodic contractile units such as in muscle or anomalous, multipolar biological spindles [44].
Multi-aster states such as aster lattices are in fact a generic feature of bulk active polar fluids [26,30], but here, we report and analyze their occurrence within droplets. To explore the aster-forming instability arising from the feedback between motor flux and motor induced polarization in Eq. 4 and Eq. 6, we perform a linear stability analysis for an incompressible, polar fluid featuring only p and m (see SI V). The analysis yields a characteristic spacing between asters in bulk, l c , such that l c ∼ 1/ m 0 ζ p , that decreases with motor density and motor induced polarization. Based on this, we expect multiple asters can be accommodated within the droplet when its major axis, l d , is long enough compared to l c . We confirm this numerically in (Fig. 3(vi)) by showing that the critical motor concentration for transition from one to two asters for varying droplet size (but all other parameters constant), does indeed decrease with initial droplet size, with an expected scaling of m 0 ∼ 1/l 2 d . Since the aspect ratio of a droplet at equilibrium increases with decreasing surface tension, the droplets get longer and can accommodate multiple asters when l d > l c ∼ ζ −1/2 p is satisfied. This explains why lower surface tension favors the two-aster state (green triangles) in Fig. 3(i). At even higher ζ p and lower γ, we find a region of the phase diagram in Fig. 3(i) where the droplet is fully divided (purple diamonds).
To see how droplet division can be enhanced, we now explore an alternative shape instability in our model Eq. 5, whereby an active density flux pushes filaments away from the aster. By varying the strength of this active flux, ζ, and the surface tension γ, while keeping ζ p (= 5.0) fixed, we obtain the phase diagram in Fig. 3(vii), which shows undeformed single tactoids (blue stars), pinched droplets (orange circles) and fully divided droplets (purple diamonds). Here, we find that droplet division occurs over a wide range in parameter space, showing that ζ is a good control parameter to trigger droplet division, while ζ p can be used to obtain multiasters states.
To analyze how the active flux parameter, ζ, for full division scales with γ, we note that if the motor and polarization are "fast variables" that relax quickly to their steady state aster solution (shown in the SI VI), the active density flux term in Eq. (5), −ζ∇ · (mp) can be obtained from the variation of an effective "free energy" functional, −ζ d 2 x ωψ . Here, ω is a scalar potential that gives the steady state polarization corresponding to an aster through mp = −∇ω. It can be shown that this effective free energy term makes a negative contribution to the droplet surface energy (detailed in the SI VI) that scales with ζ and which can then destabilize the droplet when sufficiently strong compared to surface tension. We thus expect the transition to fully divided droplets to scale as ζ ∼ γ, which is confirmed by our simulations in Fig. 3(vii).
Conclusion. We build and explore a minimal theoretical model for motor-filament droplets which captures four different shapes and show them experimentally in the actomyosin droplet system. We predict a phase diagram of expected droplet shapes based on motor activity and filament interactions, which can be tested in future experiments by systematically varying these properties. Unlike other theoretically proposed active droplet division mechanisms, the droplet shape changes we predict and show in experiment, do not require fluxes of chemical reactants [13] or large-scale fluid flows and defect dynamics [12], that arise in active nematic morphodynamics [45]. Instead, the emergence of local polar order within a nematic and the consequent spatial localization of activity, is crucial in our model, which was not explored in a previous equilibrium nematic model for this system [31]. In addition to in vitro cytoskeletal spindles, our results also have implications for the physical principles behind cell division and the organization of the biological spindle, where motor-induced active organization of filaments [9], shape change [46], co-existence of polar order with nematic order [47] and even formation of multi-lobed structures [48] can occur. We expect our work to inform strategies to use localization of active agents to achieve self-actuated shape morphing in materials [49] and synthetic cells.

I. DERIVATION AND SOLUTION OF ONE DIMENSIONAL MODEL
The 1D model for the dynamics of filament density, filament polarization and motor density given in the main text follow from the conservation equations: Our one dimensional model (Eqs. (1)- (3)) is simplified by assuming incompressibility, i.e. the density is ρ = ρ 0 . In the steady state Eqs. (1)-(3) then reduce to The Eqs. (S4)-(S5) are solved with "natural boundary"conditions that the motor density and polarization as well as their corresponding fluxes decay to zero at x = ±∞ to yield: where ξ = 2Dm v0 is the typical aster length scale and m 0 = Dξ ζρ0 is the aster strength. Figure 1(ii-iii) shows the solutions Eqs. (S6)-(S7), which show the one dimensional equivalent of an out pointing aster. The motors gather at the center of the aster as expected. Note that in 1D, the flux at steady state is a constant throughout the system, and that natural boundary conditions require this flux to vanish. In other words the diffusive flux of the motors and filaments is balanced by the corresponding active flux that advects them, at steady state. This competition sets the width of the distribution of motors at the "aster": ξ ∼ Dm v0 . In a confined droplet geometry, and at higher dimensions, solutions with nonzero fluxes are possible and seen in our numeric solutions.

II. 1D NUMERICS
To test the minimal 1D model in a finite domain, we numerically integrated the following equations Here, we included the following free energy which has a phase field for the density ρ such that a one dimensional droplet is formed and F 1d has decaying terms for the polarization p. We integrate Eq.(S8)-(S10) numerically using a a finite difference method with a forth order Runge-Kutta time discretization and periodic boundary conditions. The simulation parameters can be found in table I. Figure S2 show the motor concentration ( Fig. S2(a)), polarization ( Fig. S2(b)) and density (Fig. S2(c)) for a simulation at low activity ζ = 1.5. We find the emergence of a single droplet with a centering aster.
Going to higher activity ζ = 3.5, we find that the droplet can divide by means of the active motion of filaments. Figure S1 shows a time series of the division of a droplet. Here, the blue solid line show the motor concentration, the dashed orange line the polarization, and the dashed dotted green line the density.

III. SIMULATION METHOD
We solve Eq. (4)-(7) numerically using a finite difference scheme with a fourth order Runge-Kutta time integration on a 74 × 74 grid. We use a timestep of δt = 0.01 and the various parameters are given in Table III, while we vary surface tension, γ, and the strength of motor-induced density flux and molecular torque: ζ and ζ p . We initialize the system with a circular droplet of radius R 0 .
The parameters used in the simulations in the main text are given in Table II and Table III .

IV. ASTER CONTRIBUTION TO SURFACE ENERGY
The formation of an aster is actively induced by the motor gradient term in Eq. 6, which can be shown to arise from an effective free energy contribution, ζ p /Γ p d 2 x p·∇m, to the polarization equation. Integrating by parts and setting the surface term to zero, this term can be seen to be equivalent to a spontaneous splay energy, −ζ p /Γ p d 2 x m∇ · p.
For a uniform motor concentration, m(x) = m 0 , this spontaneous splay term can be transformed into a surface term using the Green's theorem (which is also the 2D version of the divergence theorem) to give, −ζ p m 0 /Γ p dln · p, where the line element dl is along the boundary of the droplet andn is the unit normal to the droplet boundary. Since the p points outwards in the asters, we getn · p > 0, which therefore contributes a negative line tension to the droplet interfacial energy. Comparing with the droplet line tension energy, γ dl, and assuming that polarization has a value p 0 at the droplet boundary, we arrive at the condition, ζ p p 0 m 0 > Γ p γ for destabilization of droplet shape by Fig.3(i) Fig.3(i) and Fig.3(vii)) and to compute the critical motor concentration for three aster splitting as function of initial droplet size (Fig.3(vi)).
the spontaneous splay induced by motor activity.

V. ASTER FORMATION IN BULK ACTIVE POLAR MODEL
The gain more insight into the formation of multiple asters we study a simplified model for polar active fluids in bulk. We consider the following equations for motor concentration and polarization field, ∂ t p = −ζ p ∇m + Γ p (κ p ∆p + α p p − β p p|p| 2 ). (S14) Here, we use ∆ ≡ ∇ 2 is the Laplace-Beltrami operator In this simplified model, we consider a spontaneous long range order in the polarization instead of nematic order. We also neglect motor binding kinetics, and instead set the total number of motors to a constant average value, m 0 , which corresponds to the steady state of motor binding kinetics: m 0 = k on ψ 0 /k of f . Further, we suppress the density kinetics by going to the incompressible limit and setting the compressibility ζ 0 = 0. Equations (S13)-(S14) are integrated using a finite difference method with periodic boundary condition in a square box with side length l. For the time integration, we use a fourth order Runge-Kutta method with a timestep, δt = 0.01. We initialize the system with an aster in the polarization field p = (cosφ, sinφ) T and a uniform motor concentration, m 0 . All simulation parameters are given in Table IV.
First, we study the system size dependence by varying the box size l, and initial motor concentration m 0 = 0.5. In Fig. S3 the number of asters as a function of system size l is shown. The single aster state from the initial condition is destabilized above a critical system size l c = 54 and we find a monotonic increase in number of asters. Linearzing Eq.(S13)-(S14) around the polarized state p = e y + δp, m = m 0 + δm and transforming into Fourier space with wavevector k = (k x , k y ) T gives ∂ t δm = −D m |k| 2 δm + iv 0 m 0 k · δp + iv 0 ik y δm, (S15) where we neglected the Landau terms. This system has a maximal eigenvalue whose wavevector is given by The corresponding wavelength L max = 2π/k max determines the grid spacing of asters in our simulation. Hence, we expect to find more than one aster when the system is larger then 2L max = 52.2 which agrees reasonably well with the system size we find in our simulations (see Fig. S3).
Second, we investigate the number of motors in the system while keeping the system size constant at l = 32. We initialized the system with an aster in the polarization field, and varied the initial motor concentration m 0 . Figure S4 shows the number of asters as a function of motor concentration. We find a critical motor concentration m 0 = 1.5 above which there are multiple asters emerging. Using the maximal wavevector determined Eq.(S17) by our linear stability analysis we can compute the motor concentration at which the corresponding wavelength L max is half the system size, such that it can support two asters. The resulting critical motor concentration, has a value of m 0 = 1.3, which agrees reasonably well with the critical motor concentration we find from our simulations.
In Fig S5, we show the critical motor concentration at which one aster becomes unstable as a function of systems size. We show critical concentrations from both our simulations and our linear stability analysis, which agree well. The critical motor concentration becomes smaller with larger system size.

VI. STEADY STATE ASTER SOLUTION
In order to find a steady state solution corresponding to an aster, we focus on the motor and polarization equations, which at steady state are given by, 0 = −ζ p ∇m + Γ p κ p ∆p, (S19) 0 = ∇ · [D m ∇m + v 0 (mp)]. (S20) Here, we neglected binding kinetics in the motor equation as well as the decaying terms for the polarization p and set compressibility ζ 0 → 0. A solution to Equation (S20) can be obtained by satisfying, ∇ln(m) = −v 0 /D m p. We insert this into the Laplacian of Eq.(S19) and find a solution, 0 = λ m m + ∆ln(m), where λ m = ζpv0 ΓpκpDm . This is known as the Liouville-Bratu-Gelfand equation [50][51][52] and has been solved for a number a physically relevant cases [53]. A single aster solution is given by