## Abstract

We analyse a generic motility model, with the motility mechanism arising by contractile stress due to the interaction of myosin and actin. A hydrodynamic active polar gel theory is used to model the cytoplasm of a cell and is combined with a Helfrich-type model to account for membrane properties. The overall model allows consideration of the motility without the necessity for local adhesion. Besides a detailed numerical approach together with convergence studies for the highly nonlinear free boundary problem, we also compare the induced flow field of the motile cell with that of classical squirmer models and identify the motile cell as a puller or pusher, depending on the strength of the myosin–actin interactions.

## 1. Introduction

Living cells move themselves around using different strategies, well adapted to their environment. A full understanding of the mechanisms behind cell motility is still missing but remains central for many biological and biomedical processes. Various generic mechanisms have been proposed to describe motility in different situations. Many eukaryotic cells for example move using a crawling motion. Here, motility results mainly from polymerization and depolymerization of actin filaments. The underlying treadmilling process, if combined with local adhesion of the cell on a substrate, leads to macroscopic motion. The treadmilling process and the associated crawling motion have been studied from a microscopic point of view, see [1–3] for a review on existing mathematical models. Continuum models, which allow for spatial and temporal resolution, have been considered for such a crawling motility mechanism in [4–8]. All these approaches use a reaction–diffusion system along the cell membrane and/or within the cytoplasm to effectively account for actin polymerization and combine it with a mechanical or hydrodynamic model for cell dynamics. This allows the description of the morphology and evolution of eukaryotic cells and for it to be linked to realistic signalling networks, as, for example, considered in [7,9].

Other motility mechanisms are less explored, but necessary in situations in which local adhesion is less evident, such as for cells moving in martigels [10,11] or freely swimming microorganisms. We here consider a motility mechanism arising by contractile stress due to the interaction of myosin and actin. Microscopically, myosin motor complexes use the energy from ATP hydrolysis to grab onto neighbouring actin filaments and exert stress. This process is also known for eukaryotic cells, where it shapes the rear of the cell, but it can also lead to motility itself. Here, the exerted stress is contractile and leads to a microscopic quadrupole flow around the myosin–actin complexes. A hydrodynamic active polar gel theory is developed to model these phenomena on a continuum level [12–14]. If considered in a confinement, a splayed polarization of the filaments can occur and has already been used as a route to motility [15–17]. All these studies consider a droplet. In the first case, with a surface tension using a numerical approach based on hybrid lattice Boltzmann simulations, in the second, the same setting is considered using a stream-function finite difference scheme and in the third a droplet of fixed shape is considered using an analytic description.

We will here extend the approach to include also bending properties of a cell membrane, which, however, turns out to be of less relevance for the motility mode within the considered parameter regime. The focus of the paper is a detailed computational study of the motility mechanism due to myosin–actin interactions. We explain the model used, which is here formulated in a phase-field description; demonstrate thermodynamic consistency of the overall model (without the active components); consider an adaptive finite-element discretization in space and a semi-implicit time discretization for the system of equations; and show convergence studies for critical parameters. As the considered motility mode results from a physical instability, a stable numerical discretization is essential for a detailed analysis. The simulation code is used to demonstrate the robustness of the motility mechanisms and detailed parameter studies are provided to contribute to a better understanding of the mechanisms behind cell motility for environments without local adhesion. We also analyse the flow field induced by the motile cell and compare it with a squirmer model, which allows the identification of the motile cell as a puller or pusher, depending on the strength of the myosin–actin interactions. We further discuss possible extensions of the model, e.g. combinations of myosin–actin interactions with actin polymerization. All simulations are restricted to two dimensions. The described model can also be used for three-dimensional cell motility, where the myosin–actin interactions are assumed to dominate and treadmilling only plays a minor role. However, computational studies require an adequate preconditioner/solver for the system and its development is still current research. As already exemplarily shown in [15], the motility mode remains persistent in three dimensions and we expect a similar robustness of the instability. However, a quantitative comparison of critical parameters, as well as comparisons with fluid flow measurements of moving cells will require computational intensive three-dimensional simulations.

## 2. Mathematical model

The used model is an extension of the considered approach in [15] and provides a generic route to study individual processes leading to cell motility. We will focus here on myosin–actin interactions as a source for cell motility. We review the equations and highlight the modifications.

### 2.1. Energy

We consider the free energy of the system
2.1which consists of the energy of the filament network *E*** _{P}** in the cytoplasm of the cell

*Ω*

_{cp}(

*t*), described by an orientation field

**P**, which is the mesoscopic average orientation of the actin filaments, the surface energy

*E*

_{S}of the cell membrane

*Γ*(

*t*), described by a phase-field variable

*ϕ*and the kinetic energy

*E*

_{kin}inside and outside of the cell, characterized by the velocity

**u**. For the sake of simplicity, we consider in the derivation equal density

*ρ*and viscosity

*η*for the cytoplasm and the fluid outside

*Ω*

_{out}(

*t*), which is considered as an isotropic Newtonian fluid, so that 2.2with

*Ω*=

*Ω*

_{cp}(

*t*) ∪

*Γ*(

*t*) ∪

*Ω*

_{out}(

*t*).

The phase-field variable is chosen, such that *ϕ* ≈ 1 in the cytoplasm and *ϕ* ≈ −1 in the fluid outside. The cell membrane is implicitly defined by the zero level set of *ϕ*. In [15], the cell has been considered as a droplet for which the surface energy reads
2.3where *W*(*ϕ*) = (1/4)(*ϕ*^{2} − 1)^{2} denotes the double-well potential, *ɛ* determines the interface thickness and *σ* is the surface tension. We here also take bending energy of the cell membrane into account and use the Helfrich [18], or modified Willmore energy in a phase-field approximation [19,20]
2.4where *b*_{N} denotes the bending rigidity and *H*_{0} the spontaneous curvature. We will set *H*_{0} = 0 for simplicity. If *ɛ* tends to zero *E*_{S,CH} → *σ*∫_{Γ} d*s* [21] and *E*_{S,W} → *b*_{N}∫_{Γ}(*H* − *H*_{0})^{2} d*s* [22] with *H* the mean curvature. We will consider the combination of both surface energies
2.5for which *Γ*-convergence for *ɛ* → 0 was shown in [23].

The energy of the filament network is given by [15]
2.6The gradient term with the positive Frank constant *k* is a simplification of a general distortion energy formulation from the theory of liquid crystals, with the assumption of the same value of the stiffness associated with splay and bend deformations (e.g. [24]). Linking *ϕ* to the second term allows restriction of **P** to the cytoplasm: if *ϕ* < 0, the minimum is obtained for |**P**| = 0 and thus the term does not contribute to the energy, and for *ϕ* > 0 the term forms a double-well with two minima with |**P**| = 1 and the form specified by the parameter *c*_{0}. The last term in equation (2.6) guarantees for *β*_{0} > 0 that **P** points outwards in normal direction to the cell boundary. This is expected to be of relevance for polymerization and depolymerization of actin filaments and used in [1,8,25], but for the motility mode considered here a strong preference of the orientation of **P** at the cell boundary cannot be seen. In [15], it is argued that small *β*_{0} values can resemble the effect of a weak external field. We will therefore consider both cases *β*_{0} = 0 and 0 < *β*_{0} ≪ 1 as in a more general approach with a combination of myosin–actin interactions and treadmilling *β*_{0} > 0 will be required anyhow. Figure 1 provides a schematic picture of the used variables.

Before we introduce the governing equations, we consider the energies in a non-dimensional form. We consider the characteristic values for space velocity and energy with characteristic length *L*, characteristic velocity *U* and fluid viscosity *η*. This yields a timescale and a pressure We further define the constants *c*_{1} = (*c*_{0}*L*^{2}/*k*) and *β* = (*β*_{0}*L*/*k*) and the dimensionless quantities:

— Reynolds number

— capillary number

— bending capillary number

— a polarity number Pa = (

*ηUL*/*k*)— an active force number Fa = (

*ηU*/*ζL*),

where *ζ* > 0 describes a contractile and *ζ* < 0 an extensile stress. Dropping the notation, we obtain the energies in a non-dimensional form
which are used in the following.

### 2.2. Governing equations

The equations are based on [15]. We denote the variational derivative or chemical potential of the orientation field and the phase field by and respectively.

#### 2.2.1. Orientation field equation

The orientation field equation considers a polar liquid crystal theory combined with generalized hydrodynamics, see [26,27] and [28–30] for a review, and is given by
2.7where the left-hand side is the co-moving and co-rotational derivative where the vorticity tensor defined as takes rotational effects from the flow field into account, where The deformation tensor and the non-dimensional constant *ξ* relates the coupling between the orientation field and the flow field and describes the alignment on **P** with the flow, where *ξ* > 0 for rod-like and *ξ* < 0 for oblate cells. Furthermore, *κ* = *η*_{rot}/*η* is a scaling factor between rotational and dynamic viscosity. The non-dimensional chemical potential reads
2.8

#### 2.2.2. Phase-field equation

We consider the phase field as an implicit representation of the cell surface and consider a regularized advection equation for the phase-field variable *ϕ* with the advected velocity given by the fluid velocity **u**. The introduced diffusion term is scaled with a small mobility coefficient *γ* > 0. The evolution equation reads
2.9with non-dimensional chemical potential
2.10with
2.11which describes the influence of the orientation field and
2.12which accounts for the bending and surface tension effects with
2.13and
2.14introduced to write the higher-order equation for *ϕ* as a system of second-order equations for *ϕ*, *μ*, *ψ*.

#### 2.2.3. Flow equations

The physics of the flow are described by the Navier–Stokes equations
2.15with hydrodynamic stress tensor The viscous stress is
2.16The active stress is
2.17which describes the phenomenologically introduced activity [28,31], with denoting the rescaled phase-field function, which serves as an approximation of a characteristic function for *Ω*_{cp}(*t*), with in *Ω*_{cp}(*t*) and in *Ω*_{out}(*t*). The third term, which describes the stress coming from the distortions of the filaments, reads
2.18For the Ericksen stress, we consider the divergence to be defined through
2.19which describes the stress coming from the cell surface as well as from the filaments as a result of their energy minimizing behaviour [32,33]. This term also follows for the considered case *E*_{S} = *E*_{S,CH} from the explicit form used in [15].

#### 2.2.4. Initial and boundary conditions

We consider a cell in a canal and take a rectangular domain *Ω*. We assume periodic boundary conditions on the left and right boundary for all variables. At the upper and lower boundary, we use homogeneous Neumann boundary conditions: and as well as Dirichlet boundary conditions **u** = 0 and *ϕ* = −1. The initial condition for *ϕ* is the implicitly described initial cell shape with *r* the signed distance function to the membrane *Γ*(0) and for **P** we apply an equal aligned filament network where *δ* is a vector-valued random number generated following an uniform distribution on the interval [−0.05, 0.05] in order to break the symmetry. For all simulations, we start with a circular cell with the radius *R* = 5 which is placed in the centre of *Ω* = [0, 160] × [0, 40]. The initial condition for the orientation field is

#### 2.2.5. Material parameters

We consider the following material parameters, see table 1, which are adapted from [7,15] and the references therein. The low Reynolds number allows the flow equation to be restricted to a Stokes system.

#### 2.2.6. Analytical results and numerical treatment

Neglecting all active terms, the proposed system of equations fulfil thermodynamic consistency. This is shown in appendix A. If we further neglect the orientation field (**P** = 0), the model reduces to a phase-field approximation used for vesicle–fluid interactions (e.g. [34–36]). Further neglecting the bending forces by considering only *E*_{S} = *E*_{S,CH}, we obtain ‘Model H’ in the classification of Hohenberg & Halperin [37]. If *ε* tends to zero, this special case converges to a two-phase flow problem with a jump condition for the fluid stress tensor and a continuity condition for the fluid velocity **u** at the interface (e.g. [38]). Even if this analysis cannot easily be carried over to the full system, the last condition is expected to hold and thus guarantees that fluid cannot flow through the membrane.

The system of partial differential equations is discretized using the parallel adaptive finite-element toolbox AMDiS [39,40]. We further explore an operator splitting approach, allowing the subproblems of the flow field, the orientation field and the phase-field evolution to be solved separately in an iterative process. In time, a semi-implicit discretization is used, which, together with an appropriate linearization of the involved nonlinear terms, leads to a set of linear systems in each time step. Details are described in appendix B.

## 3. Simulations

### 3.1. Motility due to contractile and extensile stress

As in [15], motility can be achieved by means of a spontaneous splay deformation. It is a two-stage process, with an elongation of the cell as a consequence of a quadrupolar straining flow resulting from the active stress tensor The elongation stops, if the surface forces characterized by Ca and Be balance the active stress. The orientation field **P**, which remains rather uniform during the elongation, starts to fluctuate, which induces a shear flow parallel to the orientation field and a spontaneously splay instability. The splayed configuration breaks the axial symmetry of the system and transforms the quadrupolar flow in a dipolar flow with two large vortices running across the cell, which has an influence on the cell shape and causes the cell to move with constant shape and at constant velocity along the symmetry axis (figure 2).

For completeness, we also demonstrate an example for cell motility due to extensile stress. Here, the vortices are reverse and the cell is stretched in the *x*_{1}-direction. Together with the active stress, which now generates a flow normal to the filaments, a bend instability occurs, describing an alignment of the filaments along the curved shape of the cell. This results in a downward motion (figure 3). The only modification needed to achieve this is 1/Fa = −3/2.

The shape and the direction of both instabilities depend on the initial conditions as well as small disturbances due to external influences. Figure 4 shows the opposite splay instability (first row) and the opposite bend instability (second row). Although the cell moves in the contrary direction, the velocity profile has a similar shape as before.

All these results qualitatively agree with Tjhung *et al*. [15]. We now turn to more quantitative comparisons and test the robustness of the instabilities.

### 3.2. Onset of motility

In any case, motility is only possible if the strength of the myosin–actin interactions exceeds a critical value. We obtain a critical activity parameter 1/Fa_{crit} ≈ 0.75. Below 1/Fa_{crit}, no instability occurs and the cell does not move. This is at least the case for *β* = 0 and in qualitative agreement with [15]. The bending capillary number Be does not influence the behaviour within the considered parameter regime. However, a quantitative comparison with the results in [15], where 1/FA_{crit} ≈ 0.5 is measured, cannot be achieved as not all parameters used in [15] are known and the critical value turns out to be highly sensitive to various parameters, which will be analysed below. Figure 5 shows the upper branch of the bifurcation diagram separating a stationary state from a splayed and moving state by plotting the constant velocity of the cell. For *β* > 0, the transition to a immotile cell is smoothed out. We no longer have a sharp transition and observe motility also below 1/Fa_{crit}, again in agreement with [15].

The onset of the instability and the time required to reach a constant shape moving with constant velocity depends on the used parameters. The stronger the myosin–actin interactions, the faster this shape is reached. This effect is most pronounced for *β* = 0 and decreases for *β* > 0. The time to reach a constant shape moving with constant velocity also depends on membrane properties of the cell. While the bending capillary number Be only plays a minor role in the considered parameter regime, the influence of the capillary number Ca is significant. The smaller the surface tension, the longer it takes to reach the desired shape. Again, this effect is less pronounced for *β* > 0.

### 3.3. Convergence tests

All obtained results are very sensitive to various parameters. The motility results from a splay or bend instability, which, for example, is heavily influenced by the elasticity of the filament network, related to the Frank constant *k*, which is here carefully chosen together with other physical parameters to observe the instability. Due to this sensitivity on the physical parameters, we would like to consider the influence of numerical parameters on the described phenomena.

We consider convergence tests. As we are primarily interested in cell motility, we first consider a parameter regime for which our cell becomes motile and moves with a constant shape and constant velocity. We consider the case of contractile stress and thus, movement in horizontal direction. We use shape and velocity for validation and measure the following quantities:

— the

*x*_{1}-coordinate of the centre of mass,

— the mean velocity of the cell

— as an average of the

*x*_{1}-component of the velocity in*Ω*_{cp}, where and— the circularity of the cell, which is defined as the quotient of the perimeter of an area-equivalent circle and the perimeter of the cell where

*B*(*ϕ*) is the perimeter of the cell.

We used absolute values for all quantities and the following error norm: ‖*e*‖_{2} = ((∑_{I}|*q*_{t,ref} − *q*_{t}|^{2})/(∑_{I}|*q*_{t,ref}|^{2}))^{1/2}, where *q*_{t} is the temporal evolution of quantity *q*. The solution on the finest grid serves as reference solution *q*_{t,ref}. Table 2 shows the relative error norms as well as the relative order of convergence (ROC) for the desired quantities if *ɛ* is reduced. We consider two cases *β* = 0 and *β* = 0.05. Together with *ɛ*, we also refine the mesh size to guarantee the same number of grid points within the diffuse interface layer for all simulations and the time step to ensure the same relation between mesh size and time step. The time interval is *I* = [0, 500]. Other parameters are obtained from table 1. We see essentially first-order convergence, the higher numbers in ROC are probably due to fortunate circumstances. Figure 6 show the shape and position for various *ɛ*, visualizing the convergence and confirming the choice of *ɛ* = 0.21 for the previous and further studies.

The second test considers the onset of motility. How sensitive is the obtained critical parameter 1/Fa_{crit} on The relation is shown in figure 7. A deeper analysis of the interface profile, as shown for a one-dimensional cut of a cell in figure 8 explains this dependency as |**P**| is slightly more smeared out than *ϕ*. This has an influence on the active stress Its divergence is reduced at the interface for increasing *ɛ* and therefore a larger activity is needed to initiate the instability.

### 3.4. Influence of different viscosities

Up to now we have considered equal density and viscosity for the cytoplasm and the fluid outside. The model can easily be extended to relax this restriction. We hereby follow a typical extension of ‘Model H’, taking *ρ* = *ρ* (*ϕ*) and *η* = *η* (*ϕ*). As shown in [41], the results for this approach are comparable to other more advanced approaches. In the following, we only consider variations in *η* and define with an appropriate function interpolating between *η*_{out} and *η*_{cp}, which are rescaled dimensionless numbers corresponding to the viscosity in the fluid outside and the cytoplasm, respectively. Figure 9 shows the dependency of 1/Fa_{crit} on the values of *η*_{out} and *η*_{cp}. Decreasing the viscosity, but keeping both values equal, leads to a reduction of the required activity for motility, but increasing the viscosity in the cytoplasm and keeping the viscosity in the outside fluid constant, in all cases, leads to an increase of the required activity. This can be explained by the necessity to induce a characteristic flow pattern in *Ω*_{cp} to induce the instability, which becomes harder to achieve for larger viscosities.

The viscosity also has an influence on the cell velocity. The reached stationary velocity *v*_{cell} increases if *η*_{out} is reduced. For more realistic parameters, with an even larger ratio of *η*_{out}/*η*_{cp}, we thus expect faster moving cells. The slope of the corresponding bifurcuation branch, as in figure 5, above 1/Fa_{crit} is reduced if *η*_{cp} is increased. The sharp transition to motility for *β* = 0 and the smoothed out transition for *β* > 0 remain.

## 4. Discussion

We already emphasized that this model describes cell motility without adhesion. Can we relate the motility mode to any freely swimming microorganism? In order to answer this question, we first compare the induced flow field with theoretical predictions for a squirmer model ([42,43] and e.g. [44]). The surface tangential velocity for a circular squirmer in a co-moving frame in polar coordinates is given by
4.1where *n* determines the velocity of the cell, whereas *m* = *n*_{2}/*n*_{1} defines whether the swimmer is a pusher (*m* < 0), a puller (*m* > 0) or a neutral (stealth) swimmer (*m* = 0), and *α* is the angle between the swimmers fixed swimming axis and the vector pointing to the surface. Figure 10 shows the surface tangential velocity for different swimmers, where we choose *n*_{1} = 0.15 as well as *m* = 0 (stealth), *m* = 0.5 (puller) and *m* = −0.5 (pusher). The profiles significantly differ with the extrema in that part of the swimmer, which is responsible for the motion. In case of a puller, it is the cell front (0 < *α* < *π*/2) and (3*π*/2 < *α* < 2*π*), whereas as the pusher is driven by the rear, so the extrema appear for (*π*/2 < *α* < 3*π*/2). For a neutral swimmer, the extrema are at *π*/2 and 3/2*π*.

We now compare these results with our simulations. We therefore extract the surface tangential velocity in the co-moving frame from our simulations. We use a contractile stress and consider for such that the stationary profile and velocity is reached. Figure 11 shows the profile for various parameters 1/Fa and *β* = 0.05. In comparison with the analytical results, we find puller dynamics for 1/Fa ≤ 0.5, similarities to neutral swimmers for 1/Fa = 0.75 and pusher dynamics for 1/Fa ≥ 1. For *β* = 0, we qualitatively obtain the same results for 1/Fa ≥ 1/Fa_{crit} and thus only pusher dynamics. The corresponding velocity profiles from the squirmer model are obtained from a data fit (figure 11): *n*_{1} = 0.086, *m* = 0.357 (puller), *n*_{1} = 0.172, *m* = 0.059 (neutral) and *n*_{1} = 0.291, *m* = −0.139 (pusher), respectively. Although we are comparing results for nearly circular shapes, see figure 12 for the corresponding stationary profiles, with that from analytic results for circular shapes, we observe a reasonable agreement.

The analytical flow field of a circular squirmer particle can be described by a superposition of a uniform background velocity, in our case, the constant velocity of the moving cell *v*_{cell}, a Stokeslet, a stresslet and a source doublet. In [45], this is used to identify typical experimental flow fields. We here consider the same approach and use the velocity field of a circular cell with centre of mass in a co-moving frame, given by
4.2where is the polar axis, scaled with the distance , the unity vector in *x*_{1}-direction and **I** the identity matrix. We prepared our numerical solution: , and claim outside the circular cell shape with radius *R* = 5 to determine *v*_{cell}, *A*_{st}, *A*_{str} and *A*_{sd}. Table 3 shows the parameters obtained from the data fit. For 1/Fa = 0.5, the stresslet parameter *A*_{str} is negative which indicates a puller-like velocity profile and for 1/Fa = 1, *A*_{str} is positive, indicating a pusher-like velocity profile. For 1/Fa = 0.75, the data fit suggests a low puller-like velocity profile. However, we should keep in mind that we compare velocity profiles of a circular and a non-circular shape. This discrepancy can be seen by analysing the relative error between the numerical results and the fitted analytical solution (figure 13). The maximum of the error appears at the part of the cell where it is compressed and does not overlap with the circular shape.

Even if a transition from puller-like to pusher-like dynamics can be observed for increasing actin–myosin interactions, the flow characteristics are much less developed than in typical squirmer models [44] and are dominated by the Stokeslet contribution. Within the analytical treatment of a circular droplet in [17], it was found that the droplet behaves like a puller. However, for the small splay considered, the corresponding flow field is not sufficient for motility and it is the quadrupole moment that characterizes the motility mechanism, resembling the motility mechanism of a squirmer. This is consistent with our findings for low 1/Fa.

In [45], the same fitting approach is used to analyse the flow topology for swimming microorganisms, such as *Cloamydomonas reinhardtii* and *Volvox carteri*. Here, the flow is also strongly dominated by the Stokeslet contribution and puller-like dynamics are only mildly developed. However, for a quantitative comparison of our results with the flow fields of such microorganisms, or that of bacteria, which typically show pusher-like dynamics, more experimental data are required. It would be interesting to investigate how predictions of the considered model in three dimensions compare with such measured flow fields in the future.

## 5. Conclusion

We here review and extend a proposed generic model for cell motility [15], which is based on spontaneous symmetry breaking in active polar gels. It models the interaction of myosin and actin as the driving mechanism for motility and does not require adhesion. The model is extended to include further membrane properties, in particular bending properties, which, however, turn out to be of minor relevance for motility in the considered parameter regime. Detailed numerical studies are performed and convergence studies considered to demonstrate the stability of the used algorithm, which is based on a phase-field description. The results clearly indicate the independence of the physical instabilities, the splay or bend instability, which are responsible for cell motility in the considered model, and possible numerical instabilities and show the robustness of the motility mode. With this confidence in the model and the developed numerical algorithm, the results are compared with model and experimental data for swimming microorganisms. Within certain parameter regimes a transition from puller-like to pusher-like dynamics can be found for increasing myosin–actin interactions, demonstrating the generic properties of the model. A quantitative comparison with swimming microorganisms is not yet possible and besides the lack of available experimental data, requires three-dimensional simulations and probably further model extension. One possible way to extend the model is a combination of the myosin–actin interactions with the treadmilling process of actin polymerization and depolymerization, described in the introduction. However, qualitative similarities with generated flow fields of microorganisms, such as *V. carteri* could already be found. The simulated flow field as well as the measured flow field is dominated by the Stokeslet contribution. In [45], it is argued that this behaviour will have an effect on the rheology of suspensions of such microorganisms. With these properties, suspensions of our modelled cells would probably behave more like suspensions of sedimenting particles, as higher-order moments are negligible in flow fields dominated by the Stokeslet contribution. However, if this assumption holds, or the weakly developed puller- or pusher-like dynamics in the considered model are already sufficient to observe typical phenomena in active fluids, as, for example, phase separation, have to be tested.

## Funding statement

W.M. and A.V. acknowledge support from the German Science Foundation through grant no. Vo-899/11. We further acknowledge support from the European Commission within FP7-PEOPLE-2009-IRSES PHASEFIELD and computing resources at JSC through grant no. HDR06.

## Acknowledgement

Simulations were carried out at ZIH at TU Dresden and JSC at FZ Julich.

## Appendix A

**A.1. Thermodynamic consistency**

Without the active terms, the proposed system of equations is thermodynamically consistent. To show this, we consider A 1with A 2 A 3 A 4which yieldswhere we have used the definition for and , which show that the integrals involving these terms vanish, and the identity from which follows that .

## Appendix B. Numerics

The system of partial differential equations is discretized using the parallel adaptive finite-element toolbox AMDiS [39,40].

**B.1. Time discretization**

We split the time-interval *I* = [0, *T*] into equidistant time instants 0 = *t*_{0} < *t*_{1} < … and define the time steps *τ* : = *t*_{n+1} − *t*_{n}. Of course, adaptive time steps may also be used. We define the discrete time derivative *d*_{t} · ^{n+1}: = ( · ^{n+1} − · ^{n})/*τ*, where the upper index denotes the time-step number and, for example, *ϕ*^{n}: = *ϕ* (*t*_{n}) is the value of *ϕ* at time *t*_{n}. In each time step, we solve:

(1) the flow problem for and

*p*^{n+1}:(2) The orientation field for : where we linearize

(3) The phase-field evolution for

*ϕ*^{n+1},*μ*^{n+1},*ψ*^{n+1}: where we again linearize the nonlinear terms by a Taylor expansion of order one, e.g. ((*ϕ*^{n+1})^{2}− 1)*ϕ*^{n+1}= ((*ϕ*^{n})^{2}− 1)*ϕ*^{n}+ (3(*ϕ*^{n})^{2}− 1)(*ϕ*^{n+1}−*ϕ*^{n}).

**B.2. Fully discrete finite-element scheme**

The fully discrete finite-element scheme follows in a straightforward manner. A *P*^{2}/*P*^{1} Taylor–Hood element is used for the Stokes problem, all other quantities are discretized in space using *P*^{2} elements. The obtained linear system, for which the direct unsymmetric multifrontal method UMFPACK is used, is solved in each time step. We use an adaptively refined triangular mesh with a high resolution along the cell membrane to guarantee at least five grid points across the diffuse interface as well as a high resolution within the cytoplasm to appropriately resolve the orientation field. The criteria to refine or coarsen the mesh is purely geometric and related to the phase-field variable *ϕ*. Due to the use of adaptivity, we need to interpolate the old solution defined on onto the new mesh To do this without violating the conservation of cell volume, we solve in every adaption step, with and defined on and *ϕ*^{n,old} on We use a multi-mesh strategy [46] to virtually integrate the first term on the finest common mesh which guarantees a constant cell volume as long as time steps are appropriately chosen. We require the interface not to propagate over a whole element within one time step. With this restriction, all numerical tests show that is conserved.

- Received February 22, 2015.
- Accepted April 1, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.