Immersed Boundary Simulations of Active Fluid Droplets

We present numerical simulations of active fluid droplets immersed in an external fluid in 2-dimensions using an Immersed Boundary method to simulate the fluid droplet interface as a Lagrangian mesh. We present results from two example systems, firstly an active isotropic fluid boundary consisting of particles that can bind and unbind from the interface and generate surface tension gradients through active contractility. Secondly, a droplet filled with an active polar fluid with homeotropic anchoring at the droplet interface. These two systems demonstrate spontaneous symmetry breaking and steady state dynamics resembling cell motility and division and show complex feedback mechanisms with minimal degrees of freedom. The simulations outlined here will be useful for quantifying the wide range of dynamics observable in these active systems and modelling the effects of confinement in a consistent and adaptable way.


Introduction
Active fluids are ubiquitous in biology and include the cell cytoskeleton [1,2], bacteria films and even schools of fish [3,4]. The common trait of these systems is that they are driven out of equilibrium at the level of their constituents, which interact hydrodynamically with their neighbours. In this paper we use a continuum model of active fluids developed in [5][6][7] where the constituent fluid particles generate local force dipoles. These systems demonstrate aspects of cytoskeletal dynamics in cells or collections of microswimmers such as bacteria, which generate dipolar stresses in their surrounding medium. Recent experimental work has shown that it is possible to confine these systems to droplets or vesicles in vitro to isolate their dynamics and observe emergent behaviour [8][9][10][11]. With these systems as motivation we consider two limits of the active fluid model, firstly the case of an isotropically ordered active fluid confined to a droplet interface and secondly the case of a droplet of active fluid with strong local polar ordering. We are able to capture the full dynamics of these two systems in a single computational model that predicts interesting phase behaviour in both cases.
We have developed numerical simulations to model active fluids confined within and on the interface of deformable droplets in 2-dimensions. The simulations are based on the Immersed Boundary method introduced by Peskin [12] to model interactions between fluids and elastic interfaces in biological systems. This method models the droplet interface as a 1D Lagrangian mesh of points which interacts with a fixed fluid mesh via a numerical analogue of the Dirac delta function (as outlined in the Methods section). The advantage of this method over phase field methods (such as [13][14][15]) is that we can track and numerically conserve quantities defined on the interface naturally in the simulation and potentially model a wide ranges of interface properties. Moreover, we can easily incorporate external boundaries in a self-consistent manner. However, a limitation of this method is that it is significantly more difficult to consider droplet separation or combination, phenomena which are captured automatically by phase field methods.
Applying the immersed boundary method to active fluids is a natural step as the properties and geometries of biological interfaces can be complex. We simulate the coupled dynamics of active particles on the interface as well as in the droplet bulk in a self consistent way. Firstly, this allows us to exactly conserve the mass of active particles within the drop (as shown in Methods). This leads to interesting droplet dynamics including transitions from stationary to motile and back to stationary as a function of activity (see Results: Active Polar Fluid). Secondly, we show that this method can also consistently model an active polar fluid droplet and compare the results of these to previous studies in the literature as well as characterising previously unpredicted states. Furthermore, we have made improvements to previous fluid immersed boundary models by ensuring conservation of droplet area (to be consistent with incompressibility) and to reduce numerical instabilities in the boundary due to the finite sized mesh (as shown in Methods). These are necessary improvements in order to model these active systems where flows are generated locally.
The structure of the paper is as follows. We first introduce the analytical equations describing the system in the following section, and then the numerical implementation of these is presented in the Methods section. Following this we present tests of the simulations and compare with other simulation data. We then present some interesting results from these simulations for the two cases of active droplets discussed. These include spontaneous symmetry breaking leading to steady swimming states, droplet deformation, and oscillating states.
In the Discussion section we describe potential uses for these simulations and extensions to these that could be performed in the future to investigate even more biologically relevant systems. In particular we discuss various ways that confinement can be added to this model. Furthermore, we outline how to extend these simulations to 3 dimensions and how we expect the results presented here to generalise to the 3D case. Finally, we also discuss potential experimental realisation of these simple systems in vitro and how these simulations can be used to simulate multicomponent active droplets or droplets with active and passive components.

Fluid Droplet Model
We construct numerical simulations of active fluid droplets using an Immersed Boundary (IB) method for modelling fluid-fluid interactions. The IB method was originally formulated for modelling fluid-structure interactions in biological problems, where the structure can be represented by elastic fibres [12,16]. In our simulations, we follow the work of [17] which adapts this to fluid-fluid interactions in a 2D domain. We begin by considering an infinitely thin, closed interface C represented by the position vector X = (X(s, t), Y(s, t)) where s is a Lagrangian coordinate along the curve C and t is time. The force density on this interface due to the surface tension is given in 2D by: F b ðs; tÞ ¼ r s gðs; tÞτðs; tÞ ½ ð 1Þ where γ is the surface tension, τ is the tangent vector to the surface and r s = @/@s. This acts as an external force on the fluids, which can be translated to a 2D force density via the 2D Dirac Delta function as follows where x is the position vector in the medium. The 2D Dirac Delta is defined by the identity Z O gðxÞdðx À XÞd 2 x ¼ gðXÞ ; where O = O 0 + O 1 is the total fluid domain, consisting of fluids labelled 0 and 1. In the case of incompressible Stokes' flow, which we assume here, force balance in the fluid dictates: where η is the fluid viscosity, v is the flow velocity and P the pressure. At the length and velocity scales of the cell cytoskeleton, inertia is negligible. An estimate of the Reynolds number using the parameter tables (Tables 1, 2 and 3 in Methods) and the density of water ρ w = 1000kg m −3 gives Re % r w zDmR 2 0 =Z 2 0 % 10 À11 (assuming the maximum velocity scale is v % zΔμR 0 /η 0 ). Finally, if we assume that the fluids are completely immiscible, then the interface velocity will be the same as the velocity in the medium at that point, such that: These simple equations model the evolution of a Newtonian fluid droplet embedded in another Newtonian fluid, with a varying surface tension γ(s, t) at the interface. In [17] surface tension decreases with surfactant concentration [18][19][20], however here we consider the case where it depends on a concentration of active matter on the interface generating contractile stresses. Thus the model is similar to models of anti-surfactants on an interface that locally increase surface tension [21,22]. The difference in an active system is that energy is consumed locally by the particles and converted to work, and it is assumed that the chemical fuel concentration can be treated as a constant everywhere.

Active Isotropic Interface
We consider an active liquid crystal dispersed on the interface such that the fluid particles are ordered along the tangent line of the interface. This is the 2D limit of the case of a 3D fluid droplet with an isotropically ordered liquid crystal confined to the plane of the interface (the 3D case is discussed analytically in [23]). As in the model of the actin cortex in [24] we consider that this active fluid is contractile enough to generate a negative pressure in the (1D) boundary layer for small concentrations, and that a passive pressure will become stronger at larger concentration. This is in qualitative agreement with experiments on reconstituted crosslinked actomyosin in vitro which, within a range of myosin and crosslinker density, initially contracts and increases in density until a steady state size is reached [25]. In this case, we consider the lowest order contribution to be linear in c as we are using a one-component description. These contractile forces, when confined along the tangent line of the droplet interface, give rise to a change in surface tension equivalent to the stress in the active fluid: where γ 0 is the bare surface tension of the fluid drop. The activity z c Δμ is a phenomenological parameter quantifying the strength of the active contractile z c Δμ < 0 or extensile z c Δμ > 0 forces with passive repulsion force proportional to B ! 0. It will be useful to consider results with respect the effective activityzDm ¼ z c Dm þ Bc 0 where c 0 is the stationary steady state concentration. Note also that there is potentially a passive pressure term linear in c(s, t) absorbed into z c Δμ, but as all of these terms are phenomenological this does not change the model. In general the order of the active term in powers of c requires details of the microscopic interaction responsible, but for the scope of this paper we consider this linear case only. As in [17] conservation of mass on the moving interface C(t) gives: d dt Z CðtÞ cðs; tÞdlðs; tÞ ¼ Z Cð0Þ @ @t cðs; tÞ @Xðs; tÞ @s where |@X(s, t)/@s| is the stretching factor, relating the length of a boundary element dl to its initial length ds such that dl(s, t) = |@X(s, t)/@s|ds. Eq (8) implies that: @cðs; tÞ @t þ cðs; tÞ j@Xðs; tÞ=@sj @ @t @Xðs; tÞ @s ¼ 0 : The second term in Eq (9) is the Lagrangian advection term. Diffusion of the mass on the closed interface is incorporated with the addition of the term @c/@t = D@ 2 c/@l 2 where D is the diffusion coefficient on the interface. Including this term, the evolution of the concentration c becomes: @ @t cðs; tÞ @Xðs; tÞ @s ¼ D @ @s @cðs; tÞ=@s j@Xðs; tÞ=@sj þ qðs; tÞ : The final term q(s, t) captures any binding to or unbinding from the boundary, which we define in Eq (13).
We explicitly consider the coupling of the interface concentration to a bulk concentration ρ in the internal fluid. We only consider this concentration as active when bound to the interface, because an isotropic stress generates no flow in an incompressible medium. In the case of the cytoskeleton active fluid, our system could model a two component system of motor proteins and actin filaments, where the isotropic filament network is concentrated on the interface and remains fixed while a population of motor proteins in the droplet can bind and unbind from this network. Then the equations describing an internal bulk concentration can be incorporated into this model by defining a level-set function H such that: Note that length element dl = |@X/@s|ds changes as the boundary deforms. Then, following [26] we define the concentration field of the unbound active particles as Hρ and hence the advection-diffusion equation describing this quantity is: The unbound concentration of active material at the boundary ρ b therefore is defined as: We assume for this model that the binding rates obey a linear scaling with local concentration. However it is likely that this is not a realistic assumption given the large gradients in concentration field that the simulations predict. Investigating and incorporating more realistic binding behaviour goes beyond the scope of this work and would demand a model informed by the microscopic details of the system in question.
Eq (12) ensures that the total mass of active material is conserved such that: Thus these equations describe advection-diffusion inside the droplet with containing boundary conditions.

Active polar fluid droplet
We independently incorporate an active fluid in the bulk of the drop. As mentioned in the previous section, an isotropic active fluid cannot generate flow if it is incompressible, so we consider the case of an active liquid crystal with strong polar ordering. Thus, we define a unit polarisation vector p which describes the direction of local polar ordering. As shown in [5] we can model the passive dynamics of this system by the Frank free energy functional with appropriate boundary terms. In this context we can write this as: The first term in Eq (16) is the one constant approximation of the Frank free energy density. The second term is used to replace the Lagrange multiplier which ensures |p| = 1 as used in [15]. This term is minimised when p is a unit vector and so for large c b this is realised everywhere inside the droplet except near to any topological defects. This can be seen by writing the polarisation as p ¼ jpjðrÞ½ cos ðbðr; yÞÞr þ sin ðbðr; yÞÞθ so that in the large anchoring limit the free energy is minimised for β = 0 everywhere, and |p| must satisfy the Ginzburg-Landau equation [27]: This ordinary differential equation (ODE) is non-linear due to the |p| 3 term and hence we are only able to solve it numerically. Assuming that |p| = 1 at the droplet boundary (strong anchoring) and |p| = 0 at the origin we find the solutions plotted in Fig 1. In the limit c b ! 1 this term fulfils the role of the Lagrange multiplier h 0 k and strictly enforces |p| = 1 on O 0 and |p| = 0 on O 1 . In the simulations, we have to use finite c b , hence there is some finite gradient of polarisation across the interface and around the defect, which sets the effective core size of the defect in the system. (16) is minimised for filament polarisation p aligned with the outward normal direction near the interface. We consider this polar boundary condition as actin networks are found to order in such a way near the cell membrane in experiments [28]. This term is the source of polar ordering in the simulations as all the other terms are symmetric for p ! −p.

The final term in Eq
This model now looks similar to many phase field models of (active and passive) liquid crystal emulsions in the literature [13][14][15]29]. The key difference in the immersed boundary formulation is that the interfacial width does not depend on the surface tension and is a numerical approximation to a sharp interface, and the level-set function H is enslaved to the position of this interface. This means we can usefully change the properties of the boundary without requiring a change to any other part of the calculation.
Thus, if we assume that activity is confined to the drop we can write the active stress as [5] s act ij ¼ ÀzDmH p i p j À The passive distortion stress in the liquid crystal is defined in terms of the molecular field h [30] The molecular field h = −δF H /δp acts to minimise the free energy functional F H with respect to p, and ν is the coupling constant between the flow gradient and polarisation. The final term is the Ericksen stress, a generalisation of the hydrostatic pressure for complex fluids, which in this case is given by [7]: In addition, there is a stress term due to the polarisation across the interface given by (as used in [14]): The stress then is incorporated into the force balance so that Eq (4) becomes: Finally the evolution of the polarisation field is governed everywhere by [5]: where the vorticity tensor ω ij = (@ i v j − @ j v i )/2 and strain rate u ij = (@ i v j + @ j v i )/2 act to reorient the polar particles.

Results and Analysis
In this section we look at some preliminary tests performed to check the accuracy of the simulations and results that are pertinent to the research questions posed in the introduction. Namely, we quantify symmetry breaking transitions and steady state behaviour including droplet motility and deformation which resemble aspects of cytoskeletal dynamics driven by active stresses. We also look at the role of certain parameters on the phase transitions we observe and relate this to similar systems discussed in the literature. We normalise the length and time scales of the simulations with respect to the initial droplet radius R 0 and the viscosity η and following some checks of accuracy choose step sizes of h = 0.075, Δs % h/2 and Δt = 5 × 10 −4 . We choose a grid length of L x = L y = 9, which gives M x = M y = 120, which has prime factors of 2, 3, and 5, resulting in efficiency savings in the FFT algorithm. A full list of the parameter value ranges used in these simulations can be found in the methods section. We relate these values to SI units by setting the length, time and energy scales to cytoskeleton relevant values such that R 0 = 10 μm and η = 10kPa s [31], and zΔμ = −10. . .10kPa [32].

Fluid Droplet Deformation
We first test the fluid simulation code alone without any active terms to check that it agrees with the expected behaviour of a simple fluid drop in 2D. We consider a droplet with an initially deformed shape given by: where y 0 kþ1=2 ¼ ðk þ 1=2ÞDs (see Fig 2). We can then measure the evolution of this shape at time n using the Fourier coefficients: whereX is the point position in the coordinate system with the droplet centre of mass at the origin and y n k ¼ tan À1 ðỸ n k =X n k Þ. Generally, we compute the first few orders of Fourier coefficients l = 2, 3, 4, but in this particular test we just present a n 2 as the shape index as this is initially set to 1/2 by Eq (24). The evolution of this shape index is presented in Fig 2(b) for several values of surface tension γ. As one would expect the droplet returns to a spherical shape, and the time it takes to do so is linearly proportional to the surface tension. Reducing the grid spacing from h = 0.075 to h = 0.025 makes no discernible difference to the droplet evolution and nor does increasing the grid size L x, y . Thus our choice of L x, y = 9R 0 appears sufficient to render interactions with the droplet's periodic image negligible.
We also note that once the droplet reaches steady state, there is still some flow in the system even though analytically the flow should be zero. This arises due to inaccuracies in the mapping to the Cartesian mesh. We find that these flows are very small and do not affect the droplet shape. However, these flow are on average only an order of magnitude smaller than those produced by the random perturbation that we apply, so they could affect the symmetry breaking direction.

Active Boundary Simulations
In this section we present the steady states observed in the active isotropic interface model. First, we consider the case where the active fluid is bound to the interface, setting k off = 0 and (a) Black line shows the boundary conformation at time t = 0 given by Eq (24). Blue arrows show the initial flow field, scaled for visibility. (b) Evolution of the shape factor a n 2 (see Eq (25)) with respect to the simulation time t (log scale) for values of the surface tension γ (in simulation units) spanning two orders of magnitude.
doi:10.1371/journal.pone.0162474.g002 quantify the steady states observed. Following this we increase this binding rate and observe the differences this coupling can make to the observed steady state.
In all of these cases, in order to enable the droplet to break symmetry we apply a perturbation to the system at time t = 2.5. We do this by altering the concentration at the interface at the beginning of the time-step by a small amount. In order to conserve mass, we perturb the concentration using a Fourier series: where θ k is as defined in Eqs (25) and (26) and the coefficients a kl and b kl are assigned randomly generated values between −10 −3 and 10 −3 . Activity and diffusion only. We find that for all values of contractile activity tested above a threshold value, a swimming steady state is observed with the droplet remaining circular (see Fig 3) as we predict analytically in [23]. We find it necessary to use a non-zero value of B on the interface in order to ensure numerical stability and physically reasonable results. In Fig 3  we show how the dependence of steady state droplet swimming speed against effective activitỹ zDm only depends weakly on the parameter B. The effect on droplet speed of the parameter B is weak because while increasing B reduces the magnitude of the gradients in c (and hence γ) along the interface, these gradients are spread over more of the droplet interface and still generate flow towards the concentration peak.
We also see from Fig 3 that the activity dependence of the speed is typical of active steady state phenomena [14] initially showing a linear dependence on activity near the threshold and a weaker dependence at higher activity values. This suggests an optimum swimming speed efficiency, above which extra swimming speed begins to cost more active energy.
The phase diagram in Fig 4(a) shows how the activity threshold depends on the diffusion of c on the interface. This dependence is approximately linear and is well characterised by the Linear Stability Analysis presented in [23]. Interestingly, Fig 4(b) shows that even though diffusion increases the active energy required to break symmetry, faster diffusion rates result in faster swimming droplets at the same activity values. Further, we see that at large enough activity, we observe the formation of two peaks in the concentration profile, which quickly coalesce to form one.
Coupling to a bulk concentration. Setting the binding rates k on,off to be non-zero can completely alter the dynamics of the system. In the case of slow bulk diffusion or binding rate compared to the active flow timescale, we see that this coupling simply dampens the instability discussed in the previous section and reduces the final droplet speed. This is because the active concentration is recycled by the internal fluid, unbinding more at peak concentrations and diffusing slowly through the drop (see Fig 5). We fix the bulk diffusion D b and increase the binding rate k b = k on = k off and observe that the feedback from the internal concentration changes the observed behaviour. The flow field in the swimming drop reference is rearward near the interface and forward through the middle. This means that the active particles unbinding at the droplet rear are advected through the middle of the drop and this feedback can lead to the  Immersed Boundary Simulations of Active Fluid Droplets nucleation of a second peak in the concentration. In the simulations this results in a continuous transition with respect tozDm from a steady motile state to 'wandering' motile states (characterised by the periodic nucleation and coalescence of a second concentration peak, snapshot shown in Fig 5(b)) and stationary two-peak solutions (Fig 5(c)). Movies of these simulations can be found in the Supplementary Information.
We see that the wandering state is really a transition phase from steady motility to a stable two-peak configuration. This result is very interesting as it suggests a simple feedback mechanism such as advection-diffusion of material through the droplet is sufficient to stabilise the two-peak solution. Physical intuition suggests that two peaks of this nature should always be unstable due to the hydrodynamic interaction between them. Thus these stable solutions highlight a simple physical mechanism that is generic to such droplet systems. However, as outlined in [23], this transition is not simply a function of the effective Peclét number in the drop, since faster diffusion in the drop decreases the activity required to spontaneously form multiple peaks in concentration on the interface. However, it is clear that the binding rate and advection terms need to be large to stabilise the two-peak solution.

Active Polar Fluid
Comparison with published results. We can test that the simulations of a droplet of active polar fluid are working by comparing to the hybrid lattice Boltzmann (LB) simulations published in [14]. By adjusting our simulation parameters to be the same as in the hybrid LB simulations of an active polar fluid we can compare the steady states directly. As in reference [14] we Immersed Boundary Simulations of Active Fluid Droplets set the following parameters as: K = 0.04, c b = 2.5, R 0 = 17, Δt = 1, h = 1, W 1 = 0, η = 5/3, Γ = 1 and ν = −1.1. We initialise both sets of simulations with a small polarisation in the positive x direction inside the droplet, and perturb the systems by applying a small normal anchoring W 1 = 0.0001 for the first 5000 time-steps. We set the surface tension to γ = 0.15 which is approximately the effective surface tension in the phase field model in [14]. We also note that, the active stress in the hybrid LB simulations is proportional to zΔμφ where φ = 2 inside the droplet, and φ = 0 outside. Therefore, when we compare activity values in Fig 6(a), the activity values used in the hybrid LB simulations are multiplied by 2 for consistency.
We can see from Fig 6 that we have good agreement between the two sets of simulations, with similar steady state velocity trends and droplet shape over a range of contractile activity values. We see that the symmetry breaking threshold agrees well between the simulations, however the non-linear dependence of droplet speed on activity is slightly different. This may in part be due to the estimates of surface tension which have no analytical analogue in the LB simulations, as well as the only approximate incompressibility achieved by LB (which is exact in our model).
Symmetry Breaking: Motion and Deformation. In this section we consider specifically the case with a large polar anchoring term W 1 = 50h in order to investigate the effect of this coupling between internal polarisation and shape. We observe that with this boundary condition we only see an active motile state for extensile activityzDm > 0. This is in contrast to simulations by Giomi et al. [15] which show that, with normal anchoring an active nematic droplet shows swimming states for contractile activity. This suggests that it is the difference between nematic and polar order that is responsible for the difference. The stationary state of the system with polar order consists of an aster with defect q = +1 at the droplet centre, whereas in the nematic case it consists of an elongated droplet with q = ±1/2 defects along the long-axis of the droplet. This difference in stationary state therefore changes the geometry such that the stability threshold for motion changes sign. This change appears to be robust for the values of ν we have tested (±1.1 corresponding to 'rod-like' and one to 'disc-like' flow aligning particles [30]).
As predicted in [23] the activity threshold for the onset of motion does not depend on the surface tension, however the droplet steady state speed does. As we can see in Fig 7 the motile steady state is elongated and the shape deformation depends on the surface tension. In Fig 7 we plot the k = 2 mode shape parameter from Eqs (25) and (26) (essentially measuring droplet elongation) and final droplet velocity against the surface tension. Clearly, due to the strong anchoring at the interface, the polarisation profile is coupled to the boundary shape, Immersed Boundary Simulations of Active Fluid Droplets so when the surface tension is small the polarisation profile is able to deform the droplet more and reach a faster steady state. As surface tension is increased the deformation reduces and goes towards a circular droplet, and in this limit the droplet velocity begins to plateau at its minimum value. This shows a positive feedback between the shape and the motion in these droplets, but in this case it does not affect the symmetry breaking threshold. This defect driven motion is similar in appearance to those predicted for crawling active droplets in [33,34], despite the difference in geometry and flow boundary conditions.
For contractile activity we see that the lowest order symmetry breaking mode corresponds to a non-motile state where the droplet appears to initiate division. The rate of deformation of the droplet interface is exponential and is modelled well by the linear stability analysis at small times. In some cases we see that this deformation is halted at some steady state 'dumbbell' shape, where the defect is elongated into a line defect along the long axis of the droplet and can separate into 3 distinct defects (of charge +1, -1 and +1 respectively) shown in Fig 8. At larger activity values or lower γ we see that the opposite sides of the immersed boundary are almost brought into contact. However, these simulations do not yet have the capability of modelling division of the droplet as this is not a part of the IB model, which is a shortcoming in comparison to phase field models. Thus in these cases the simulation method breaks down. In most cases we test the deformation does stop and we observe further symmetry breaking of the internal polarisation field such that the droplet front enlarges and the droplet swims with the 3-defect configuration. However we do not characterise these dynamics here as they are complex and depend strongly on the value and sign of the liquid crystal coupling parameter ν.
In principle one could impose conditions that would allow division of the boundary, such that there were two separate immersed boundaries in the simulation. To do this, one would need to calculate the separation between non-neighbouring boundary points k 1 and k 2 , and once they were within a certain distance (say ds) separate the boundary into 2 boundaries one from k 1 + 1. . .k 2 the other from k 2 . . .k 1 .
Onset of Active Turbulence. It is well known that active systems can display turbulent like behaviour at zero Reynolds' number when the the equilibrium relaxation time is long compared to the timescale of the active dynamics. In terms of activity, the threshold for this behaviour is proportional to ηK/ΓL 2 where K is the one Frank elastic constant and L is the system size [3,4,35]. These turbulent-like states have been observed experimentally in bacterial suspensions and reconstituted cytoskeletal networks [9,36] where the active flows can spontaneously create and annihilate polarisation defects in pairs and generate vortices in the flow. We observe the onset of such dynamics in an active polar fluid droplet by reducing the elastic constant in our simulations to K = 10 −3 . We find that due to the strong anchoring condition at the  interface, and relatively high surface tension the polarisation field remains approximately fixed at the boundary. However, reducing the constant K reduces the equilibrium interaction between these filaments and those near the droplet centre which show dynamics that are only weakly coupled to the interface.
This transition to turbulent motion is easily observed if we consider the case of contractile activity which for low activity values stabilises the defect at the droplet centre. Snapshots of the dynamics are displayed in Fig 9. We see that due to the low K value the stationary configuration becomes unstable but the flows are not large enough to significantly deform the drop. The result is unsteady dynamics of the polarisation field in the centre of the drop that are mostly decoupled from the interface dynamics. This results in stretching and distortion of the original +1 defect into a line and then further pairwise creation of ±1 defects (shown in Fig 9 by the points where four black stripes intersect). As we scale the droplet size up, or decrease K, one would expect that these turbulent dynamics would become faster as the effective energy cost of deformation becomes smaller relative to the activity. Note that due to the polar nature of the system modelled here, only defects of integer topological charge can be generated, so these have fundamentally different dynamics to the half integer defects observed in active nematics.

Extension to 3 dimensions
This method can be extended to a full 3D simulation, but the steps required to do so go beyond the scope of this paper. The equations related to the fixed fluid grid, such as the fluid flow Eq (22) and polarisation dynamics Eq (23) remain unchanged in 3D and the same combination of Fast Fourier Transform and Runge-Kutta methods (see Methods) can be utilised to solve for the dynamics of these quantities. However, the fluid-fluid interface is a 2D surface in the generic 3D case, and so this requires a different representation from the 1D Lagrangian mesh used here. A direct extension would be to represent the fluid interface as a 2D triangular mesh. Then the boundary points will evolve in exactly the same way as before. However, rather than removing single points from the interface (as outlined in Methods) one would have to remove individual triangles below a threshold area and replace the mass contained within that triangular area. Similarly, triangles with excessive eccentricity or area should be split along their short axis into two triangles to ensure stability. In general one could use any standard re-meshing algorithm ensuring that the concentration on the interface is numerically conserved and the spatial distribution is unchanged. Such a mesh has been used to model surfactant coated drops in [37]. An alternative method is to only model the Immersed Boundary points implicitly via the Cartesian fluid mesh (known as the Immersed Interface Method [38]). Then one has to explicitly consider the stress jump in the normal direction at the interface in order to impose the boundary force [39]. This normal force is then imparted on the Cartesian mesh within the interfacial region (defined by the numerical delta function range). Similarly, the boundary region and associated level-set function are advected due to flow normal to the interfacial region. Thus the approach is closely related except that the boundary mesh is implicitly defined rather than an explicit Lagrangian mesh.
We expect that the dynamics reported in this paper will be representative of what would be seen in 3D. In [23] a linear stability analysis of these two systems predicts the same qualitative symmetry breaking behaviour despite the flow being quantitatively changed by the geometry. The main qualitative difference one would expect to observe is a deformation of the droplets with an active interface in 3-dimensions, and an investigation of how this deformation couples to concentration dynamics would be an interesting investigation for future work.

Summary and Conclusions
In this paper we have detailed the numerical method that we have used to simulate the evolution of active droplets in 2D. We have used an Immersed Boundary (IB) method to model the fluid boundary as a Lagrangian mesh of points. This allows us to explicitly simulate a concentration of isotropic active contractile particles on the interface that alter the surface tension. Further, we use the immersed boundary to define a level-set function which defines the inside and outside of the drop. Using this we can define an active polar liquid crystal phase in the droplet interior with anchoring at the droplet interface.
Our results reveal steady state migrating states in both types of active droplet. In all cases we investigate the steady state droplet velocity, which scales non-linearly with activity, levelling off at larger activity values. This means there is a generic reduction in swimming efficiency as activity is increased in these systems as observed for the flow rate in other active droplet and bulk systems [40].
The first system we consider is a contractile isotropic active fluid on the droplet interface and observe that without any binding and unbinding the only observed active dynamic steady state is a swimming droplet with one peak in concentration. We see that for small binding rates the steady state remains unchanged but reduces the droplet velocity. At higher binding rates feedback from the fluid flow stabilises a second concentration peak and the droplet can reach a stationary two-peak configuration for large enough activity. This shows how a simple feedback mechanism can stabilise more complex behaviour, since at low binding rates if two peaks are formed in the concentration they always attract and coalesce to form a single peak at the droplet rear. These higher order states show that this active interface shares much of its internal dynamics with models of thin layers of active fluids [41,42], and so provides a useful insight into how such a system evolves immersed in an external fluid. Thus, it is a useful step in modelling the dynamics of the actin cortex in cells which many active thin film models are based on.
Secondly, we simulate an active polar fluid droplet with strong anchoring at the interface. We see that if when we take the Frank elastic constant to be significant (or equivalently a small droplet) we see transitions to highly ordered dynamic states, including motion, rotation and symmetric deformation. Due to the geometry of the droplets and the polar anchoring condition, we find strikingly different behaviour for the contractile and extensile regimes, with extensile activity promoting translational symmetry breaking and contractile promoting symmetric modes of deformation. Using a smaller Frank constant (equivalent to a larger droplet) we observe the onset of turbulent dynamics, discussed in depth for active nematics in [43]. Due to the confinement inside the droplet the dynamics are still constrained and we see that the number of defects in the polarisation field rarely exceeds 3 for the range of activities simulated here. However, the dynamics are unsteady due to the slow relaxation of the molecular field compared to the active flow timescale.
The results discussed here are consistent with previous studies of other active fluid droplet systems [14,15]. Additionally, given the increase in experimental progress in producing stable active gels and fluids in vitro [25], and confining such systems to droplets and vesicles [8][9][10] it is possible that these simplified models can be tested experimentally. In particular the active fluid interface model may be achievable by coating a fluid droplet or the inside of a vesicle with a stable isotropic crosslinked actin layer. Active contraction of the layer could be introduced by trapping myosin motor complexes inside the droplet/vesicle using existing microfluidic techniques. Such motor complexes may have different binding kinetics than the linear case assumed here, however one would expect to observe qualitatively similar behaviour.
A useful extension of these simulations that is beyond the scope of the work we present here would be to simulate these systems in 3 dimensions. From our analysis in [23] we would expect much of the observed dynamics to remain qualitatively the same in 3D, with the exception that we would expect deformation of the droplets under surface tension gradients. In particular it would be interesting to simulate the droplet deformation when we predict steady two peak solutions in the concentration of an active isotropic fluid on the droplet interface (see Fig 5) which could equally manifest as a ring or two poles of higher concentration in this geometry.
The simulations presented here can also be used to model interactions with external boundaries or elastic solids in a self-consistent way through the IB method. Similarly, one could investigate the dynamics of multiple active droplets simply by adding more Lagrangian boundary meshes. These boundaries interact via their coupling to the fluid flow, and hence one would be able to simulate confinement of these droplets without necessarily sacrificing the periodic solution method outlined in the Methods section. A preliminary investigation into such external boundaries is discussed in [44].
Finally, the modular nature of the code means that these simulations are capable of simulating a coupled active interface and an internal active polar fluid. The dynamics of such a system will be even more complex and so the results presented here will be relevant to understanding those systems in future. Similarly, this will be relevant also to systems of passive liquid crystal drops propelled by active surfactants as studied in [45].
In conclusion, the results we present here demonstrate complex dynamics in confined active fluids with minimal degrees of freedom. These simple one-component models demonstrate steady state dynamics akin to those observed in cells, cell fragments and reconstituted cell components that are driven by the active cell cytoskeleton. The numerical method we outline in this paper, based on the Immersed Boundary method enables continuum level simulation of these systems and has a wide range of potential applications in this field.

Methods
In this section we outline explicitly the numerical method used in the simulations and the parameter value ranges studied.

Coordinate System
The bulk fluid domain O is defined on a staggered periodic Cartesian grid of points with spacing h in each direction. The 3 overlaid grids correspond to positions (x i , y j ), (x i−1/2 , y j ), and (x i , y j−1/2 ) where x i−1/2 = x i − h/2 and the same for y. The points x i and y j run from x 0 = −(L x − h)/2 where L x, y are the lengths of the simulation box in the x and y directions and M x, y = L x, y /h is the total number of grid points in each direction). The interface C is described by a staggered grid of Lagrangian points s k where s k = kΔs with k = 0. . .N(t) − 1. These points have corresponding Cartesian coordinates (X k , Y k ) and (X k+1/2 , Y k+1/2 ) = (X k+1 + X k , Y k+1 + Y k )/2. The step size is defined as Δs = l 0 /N(0) where l 0 is the initial boundary length and N(0) the initial number of boundary points. Points can be added or removed from this mesh during the simulation, which will be outlined later in this section. Fig 10 gives an outline of which quantities are defined on which grid. We integrate over time with step Δt such that t = nΔt for n = 0. . .T. We will denote the time-step at which a quantity is defined by a superscript where necessary (e.g. X nþ1 k ), if none is given then all quantities in that equation will be at the same arbitrary time-step n.
Staggered grids are commonly used for fluid simulations of this nature in order to achieve consistency between first and second order gradients without the need to use a 5-point stencil in 1D (or 13-point in 2D). We define the gradient operator on the bulk fluid as r h g i, j (@ x , @ y )g i, j = (g i+1, j − g i, j , g i, j+1 − g i, j )/h so that if G i, j = r h g i, j then the x-component of the function G i, j is defined on the grid (x i+1/2 , y j ) and the y-component on (x i , y j+1/2 ). Values of G i, j can then be interpolated to other lattices where required. The Laplacian operator can then be defined as such that it is on the same grid as the function g. Similarly, the boundary derivative @ s g k = (g k+1 − g k )/Δs is defined on the k + 1/2 lattice, such that @ 2 s g k ¼ @ s ð@ s g kÀ1 Þ ¼ ðg kþ1 À 2g k þ g kÀ1 ÞDs 2 .

Immersed Boundary Equations
We use the 2D numerical delta function δ h (r − X) described in [12] to couple boundary quantities to the bulk fluid and vice versa. This delta function is given by the product of two 1-dimensional functions F: where the function F acts as a numerical analogue of the 1D analytical Dirac delta function. It is desirable to use a function that is only non-zero for a few points, as this aids computational efficiency, hence a restriction regularly chosen is F(r) = 0 if r < −2 or r > 2. In [12] several constraints are imposed that one would expect from a delta function, namely that it be a continuous function in r and: X i¼even where λ is a constant. Note that the condition in Eq (30) gives the identity ∑ i F(r − i) = 1 but also ensures a symmetric distribution of F. The derivation of F is outlined in reference [12], here we just give the resulting function: From this it is straightforward to check that λ = 3/8. It should be noted that this distribution maintains its properties for all real values of r and is not restricted to integer values.
Using this delta function one can distribute the force density on the boundary F ðbÞ k to the fluid. First, we discretise Eq (1) to give: where the tangent vector (defined on the X k+1/2 grid) is given by The corresponding force density components in the fluid then can be defined as: f ðy;bÞ where the (x) and (y) superscripts are used to denote components of the vector. Note that the x component is defined on the x i−1/2, j grid and the y component on the x i, j−1/2 , this is the same for the components of the velocity vector. We compute the velocity components v ðxÞ i;j and v ðyÞ i;j , and pressure P i,j by solving the incompressible Stokes equation (i.e. Eq (4) such that rÁv = 0). To do this we employ a projection method where we first define the intermediate solutions v 0 and P 0 such that where f tot is the sum of all the force densities acting on the fluid. First, we solve for v 0 ; utilising the periodicity of the lattice we can write any function g i, j in terms of a 2D Discrete Fourier Transform (DFT): Here we have used the notation I ¼ ffiffiffiffiffiffi ffi À1 p to avoid confusion with the index notation. Using Eqs (40) and (28) we can re-write Eqs (38) and (39) so that for each value of α and β: Using standard NAG library functions [46] to compute the Fast Fourier Transforms (FFTs), we find the DFTs of the components of f tot i;j . Looping over all values of α and β we can use Eqs (42) and (43) to assign values to the DFT of v 0 i;j . We then compute the inverse DFT ofv 0 to acquire values for v 0 i;j . Following this we calculate r h Á v 0 i;j and compute its DFT. Then we can calculate P 0 i;j similarly. Finally, we relate these quantities back to the velocity and pressure as follows: This is a computationally efficient method for solving the Stokes equation as each set of assignments only require M x M y operations. The time taken for the FFT computations scales approximately as M x M y log 10 (M x M y ) so this is the limiting complexity in this process.
Once the bulk velocity is computed we can calculate the velocity of the interface using Eq (6): V ðyÞ k ¼ where i 0 is the value of (X k + (L x − h)/2)/h rounded down to the nearest integer and j 0 is the same but for Y k . Generally in IB methods this interface is updated using a simple forward first order Euler scheme, i.e. X nþ1 k ¼ X n k þ DtV k . Instead we use an explicit fourth order Runge-Kutta method as outlined later in this section. We find that this method improves the numerical stability of the algorithm compared to a forward Euler update, particularly in the case where the boundary position is coupled to an internal polarisation field.

Active boundary equations
We define a concentration of active particles c k on the boundary points X k+1/2 which alters the surface tension γ k according to Eq (7). The concentration dynamics are then governed by Eq (10) which we discretise using the Crank-Nicholson scheme (as in [17,26]): The quantity q contains the binding and unbinding rates for the concentration and is defined further below in Eqs (59) and (60). We update the concentration c after updating the boundary position X so that X n+1 is defined. Then the set of Eqs in Eq (50) can be represented by a periodic tridiagonal matrix: Immersed Boundary Simulations of Active Fluid Droplets Rearranging Eq (50) we can read off the expressions for l,d,u, and b: We solve this matrix directly using standard NAG linear algebra functions [46]. We test that the mass is conserved numerically by outputting the change in total mass in the simulation between time-steps and find that these values are randomly distributed around zero with maximum deviations on the order of machine precision.

Coupling to a passive bulk fluid
In order to numerically model the equations in the Model section we need a numerical analogue for the level-set function H. Taking the gradient of both sides of Eq (11) we see that by definition: Thus, as in [26,47], to find the numerical equivalent of the function H we solve the following Poisson equation: We again solve this equation using the same FFT method used in to solve Eqs (42), (43) and (44). However, as we set the magnitude of the zero-wavenumber mode of the Fourier Transform to zero, this only defines H up to a constant (more exactly it defines H so that its average value is zero). This can be corrected simply in the simulation by recording the value of H at any point where rH = 0 and H < 0 (a point strictly outside the drop) and then subtracting the value of H at this point from the value of H at all points in the simulation. This then results in a profile where all points outside the droplet are approximately equal to zero. However, as discussed in [48] this can lead to numerical issues when dividing by H, and so we redefine the level set function at all points as H ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi H 2 þ 2 p setting = 10 −6 as a sufficiently small discrepancy to avoid noticeable 'leaking' of the internal concentration.
We can then define a concentration field that is contained within the droplet by the quantity Hρ. Thus, in order to conserve mass and ensure that the concentration is contained within the droplet we can use the same equations as used in [26] to describe a soluble surfactant. We update ρ from Eq (12) using the following Crank-Nicholson scheme (as used in [26]): where D h is a central difference operator on the staggered grid (as defined in [26]). The exchange terms q n k and q nþ1 k are given by: P M y j¼0 Hi; j n r n i;j d h ðx i;j À X n kþ1=2 Þh 2 À k off c n k ð59Þ q nþ1 As we update the boundary positions and concentration before this step, X nþ1 k , c nþ1 k and H nþ1 i;j are already defined. Analytically we know that the system should reach a stable steady state where r ¼ k off c 0 =k on inside the drop and c = c 0 on the interface. In these simulations however, the bulk concentration evaluated at the interface is i¼0 P M y À1 j¼0 H n i;j r n i;j d h ðx i;j À X n kþ1=2 Þh 2 , and so the equivalent steady state for this numerical system is given by: Where the droplet shape and position is assumed fixed in this passive case and the angled brackets indicate the average value over the boundary such that: Thus, we initialise the simulation with a concentration c k = c 0 = 1 on the interface and ρ = 0 everywhere except at x i = y j = 0 where ρ = ρ 0 everywhere so that the concentration in the system is directly proportional to H and is in the stable steady state configuration (in the absence of activity or external flows).

Active Polar Fluid
We define the polarisation field p i, j , h i, j and stress components s ðact;dis;intÞ i;j from Eqs (18)-(21) on the centred lattice (x i , y j ) and interpolate any gradients that are evaluated on the other lattices therein. First, we calculate the molecular field and the free energy density using Eq (16) from the identity h = −δF H /δp: where the central gradient operator r hc is defined by r hc g i,j = (g i+1,j − g i−1,j , g i,j+1 − g i,j−1 )/2h. From these we calculate the components of the stress a i;j ¼ act i;j þ dis i;j þ int i;j , and then finally the force density contribution from these stresses: As with the velocity field, the x-component of the force density is defined on the x i−1/2, j lattice and the y-component on the x i, j−1/2 lattice. This force density is then added to the boundary force density in Eq (38) before calculation of the velocity. Once the velocity is computed, the polarisation update equation can be computed from Eq (23) where iÀ1=2;j þ @ y v ðxÞ iþ1=2;j þ @ y v ðxÞ iÀ1=2;jÀ1 þ @ y v ðxÞ iþ1=2;jÀ1 i u ðyyÞ i;j ¼ @ y v ðyÞ i;jÀ1=2 : The polarisation is updated using the coupled fourth order Runge-Kutta method at the same time as the boundary is updated using Eqs (48) and (49), detailed in the "Algorithm Overview" section.

Adding and Removing Boundary Points
As the immersed boundary in these simulations is the interface between two fluids, there is no resistance to neighbouring points getting too close or too far apart, as there would be in an elastic boundary. Therefore, it is necessary to add and remove points from the simulation to keep the boundary well defined and maintain accuracy. We remove the point X k whenever |@ s X k−1 | < 0.65 and add a point between X k and X k+1 when |@ s X k | > 1. 35. See Fig 11 for a visual guide to this process.
To remove a point X k we reassign X 0 m ¼ X mþ1 for all k m < M − 1, then reduce the total number of boundary points to N 0 = N − 1. Conservation of mass then dictates: and again c 0 m ¼ c mþ1 for all k m < N − 1. The dashed quantities all denote the new values after removal of the boundary point.
To add a point between X k and X k+1 we reassign X 0 mþ1 ¼ X m for all k < m < N then add the point X 0 kþ1 ¼ ðX kþ1 þ X k Þ=2 with concentration c 0 kþ1 ¼ c k to conserve mass (c k remains unchanged).
We notice that when large flow fields are generated in the simulations, inaccuracies in the boundary positions can accumulate that need to be addressed. Firstly, we observe 'buckling' of neighbouring boundary points in areas where the boundary flow converges, which results in points not being removed from the simulation even though they are overlapping. Thus, we also add a condition for the removal of a boundary point if jr 2 s R k j > d max where, after some testing of different boundary shapes, we arrived at d max ¼ 0:05R 0 =Ds 2 , where we have defined R k to be the distance of the boundary point from the droplet centre of mass. Therefore any excessive second order gradients in the boundary position are effectively smoothed out. However, we choose the limit so that smooth boundary deformations are still permitted, such as those observed in our results for an active polar fluid drop.

Droplet Area Conservation
We observe that the droplets can lose significant area over the whole simulation time, despite the fluid incompressibility. This effect is observed in the fluid immersed boundary simulations presented in [17] and [26] but is very slight and gradual (on the order of 10 −4 % area loss over a whole simulation) as the size of the flows is imposed externally. It is possible to observe much more rapid area loss in our simulations due to the continual addition and removal of interface points, and so to prevent this, and be consistent with the incompressibility condition, we n + 3/4 values of the polarisation and boundary condition. In between each fractional update we repeat steps 1-4. Then finally the updated polarisation and boundary are given by: We run the simulation, iterating over steps 1-9 usually 10 6 times, unless we deem that longer is required. A full list of the parameter value ranges is given below (Tables 1, 2

and 3):
Supporting Information S1 Video. Steady state swimming droplet from Fig 5(a). Colour gradients show concentration on the interface from black (low) to yellow (high). (MPEG) S2 Video. 'Wandering' droplet from Fig 5(b). Colour gradients show concentration on the interface and in the bulk, from black (low) to yellow (high). (MPEG)