A mechanism for cell motility by active polar gels

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 to consider 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.


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 e.g. [1,2] and [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,5,6,7]. 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 to describe the morphology and evolution of eukaryotic cells and link it to realistic signaling networks, as e.g. considered in [8,7].
Other motility mechanisms are less explored, but necessary in situations in which local adhesion is less evident, such as for cells moving in martigels [9,10] 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 on neighboring 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 myosinactin complexes. A hydrodynamic active polar gel theory is developed to model these phenomena on a continuum level, see [11,12,13]. If considered in a confinement, a splayed polarization of the filaments can occur and has already been used as a route to motility [14,15,16]. 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 extent the approach to include also bending properties of a cell membrane and consider an adaptive finite element approach.
The focus of this paper is a detailed computational study of the motility mechanism due to myosin-actin interactions in 2D cells. We explain the used model, 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 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 analyze the flow field induced by the motile cell and compare it with a squirmer model, which allows to identify 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. Basically, the described model can also be used for 3D 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.
2. Mathematical model. The used model is an extension of the considered approach in [14] 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 which 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 equal density ρ and viscosity η for the cytoplasm and the fluid outside Ω out (t), which is considered as an isotropic Newtonian fluid, so that 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 φ. For convenience, we will also introduceφ = 0.5(φ + 1) as a rescaled phase-field function, which serves as an approximation of a characteristic function for Ω cp (t). In [14] the cell has been considered as a droplet for which the surface energy reads where 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 [17], or modified Willmore energy in a phase-field approximation [18,19] where 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 → σ Γ ds [20] and E S, [21] with H the mean curvature. We will consider the combination of both surface energies for which Γ-convergence for ε → 0 was shown in [22]. The energy of the filament network is given by [14] The 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, see e.g. [23]. Linking φ to the second term allows to restrict 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 eq. (2.6) guarantees for β 0 > 0 that P points outwards in normal direction to the cell boundary. Fig. 2.1 provides a schematic picture of the used variables. Before we introduce the governing equations, we consider the energies in a nondimensional form. We consider the characteristic values for space x = Lx, velocity u = Uû and energy E = ηU L 2Ê , with characteristic length L, characteristic velocity U and fluid viscosity η. This yields a time scale t = L Ut and a pressure p = ηU Lp . We further define the constants c 1 = c0L 2 k and β = β0L k and the dimensionless quantities:  notation we obtain the energies in a nondimensional form which are used in the following.

Governing equations.
The equations are based on [14]. We denote the variational derivative or chemical potential of the orientation field and the phase field by P = δE δP and φ = δE δφ , respectively. 2.2.1. Orientation field equation. The orientation field equation results from a classical Toner-Tu field theory [24,25] which is combined with flow-orientation couplings, see e.g. [26,27] for a review, and is given by where the left hand side is the co-moving and co-rotational derivative where the vorticity tensor defined as Ω = 1 2 (∇u − ∇u) takes rotational effects from the flow field into account, where ∇u = (∂ j u i ) (i,j) . The deformation tensor D = 1 2 (∇u+∇u ) and the nondimensional constant ξ relates the coupling between the orientation field and the flow field and describes the alignment of the cell with the flow, where ξ > 0 for rod-like and ξ < 0 for oblate cells. Furthermore, κ = η rot /η is a scaling factor between rotational and dynamic viscosity, to be taken as 5/3. The nondimensional chemical potential reads

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 with nondimensional chemical potential which describes the influence of the orientation field and which accounts for the bending and surface tension effects with introduced to write the higher order equation for φ as a system of 2nd order equations for φ, µ, ψ.

Flow equations.
The physics of the flow are described by the Navier-Stokes equations with hydrodynamic stress tensor σ = σ viscous +σ active +σ dist +σ ericksen . The viscous stress is The active stress is which describes the phenomenologically introduced activity [26,28], withφ = 0.5(φ + 1) denoting the rescaled phase-field function, which serves as an approximation of a characteristic function for Ω cp (t), withφ ≈ 1 in Ω cp (t) andφ ≈ 0 in Ω out (t). The third term which describes the stress coming from the distortions of the filaments, reads For the Ericksen stress we consider the divergence to be defined through and describes the stress coming from the cell surface as well as from the filaments as a result of their energy minimizing behavior, see [29,30]. This term also follows for the considered case E S = E S,CH from the explicit form used in [14].
2.2.4. Thermodynamic consistency. Neglecting all active terms, the proposed system of equations fulfill thermodynamic consistency. To show this, we con-siderĖ where we have used the definition for ∇ · σ ericksen and σ dist , which show that the integrals involving these terms vanish, and the identity u×(∇×u) = ∇(|u| 2 )−(u·∇)u from which follows that u · (−Re(u · ∇)u) = 0.
2.2.5. Initial and boundary conditions. We consider a rectangular domain Ω and split the boundary into two parts. 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: ∇P · n = ∇P · n = 0, and ∇µ · n = ∇ψ · n = 0 as well as Dirichlet boundary conditions u = 0 and φ = −1. The initial condition for φ is the implicitly described initial cell shape φ = tanh(r/( √ 2ε)), with r the signed distance function to the membrane Γ(0) and for P we apply an equal aligned filament network P = (P 1 , P 2 ) + δ, where δ is a vector-valued random number generated following a 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 center of Ω = [0, 160] × [0, 40]. The initial condition for the orientation field is P = (1, 0) + δ.
2.2.6. Material parameters. We consider the following material parameters, see Tab. 2.1, which are adapted from [14,7]  3. Numerics. The system of partial differential equations is discretized using the parallel adaptive finite element toolbox AMDiS [31]. We restrict our simulations to the two-dimensional case. We further explore an operator splitting approach, allowing to solve the subproblems of the flow field, the orientation field and the phase-field evolution separately in an iterative process. In time, a semi-implicit discretization is used, which, together with an appropriate linearization of the involved non-linear terms, leads to a set of linear systems in each time step.
3.2. Fully discrete finite element scheme. The fully discrete finite element scheme follows in a straight forward 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 T h 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 T n h onto the new mesh T n+1 h . To do this without violating the conservation of cell volume, we solve φ n,old , θ = φ n,new , θ in every adaption step, with θ and φ n,new defined on T n+1 h and φ n,old on T n h . We use a multi-mesh strategy [32] to virtually integrate the first term on the finest common mesh T n h ∪ T n+1 h , which guarantees a constant cell volume as long as time steps are appropriately chosen.

Simulations.
4.1. Motility due to contractile and extensile stress. As in [14] 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 σ active . 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, see   Table  2.1 and we changed 1/Fa=1.125 and take β = 0 (no explicit forcing for P to point outwards at the cell boundary).
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, see Fig. 4.2. 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. Fig 4.3 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.

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.76. Below 1/Fa crit no instability occurs and the cell does not move. This is at least the case for β = 0 and in agreement with [14]. Fig. 4.4 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 .
The onset of the instability and the time required to reach a constant shape moving with constant velocity depends on the used parameters. As stronger the myosinactin interactions, as 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, 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.

Convergence tests.
All obtained results are very sensitive to various parameters. The motility results from a splay or bend instability, which e.g. is heavily influenced by the elasticity of the filament network. 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 center of mass, • the mean velocity of the cell as an average of the x 1 -component of the velocity in Ω cp , where u = (u 1 , u 2 ) , 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 . Tab. 4.1 show 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.1. 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 Tab. 2.1. We see essentially first order convergence, the higher numbers in ROC are probably due to fortunate circumstances. Fig. 4.5 show the shape and position for various ε, visualizing the convergence and confirming the choice of ε = 0.21 for the previous and further studies.  Table 4.1: Relative error norms and convergence orders for critical parameters, upper part β = 0 and lower part β = 0.1.
The second test considers the onset of motility. How sensitive is the obtained critical parameter 1/Fa crit on ε? The relation is shown in Fig. 4.6. A deeper analysis of the interface profile, as shown for a 1D cut of a cell in Fig. 4.7 explains this dependency as |P| is slightly more smeared out than φ which will have an influence also on the onset of instability.  5. 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 [33,34] and e.g. [35]. The surface tangential velocity for a circular squirmer in a co-moving frame in polar coordinates is given by where n 1 determines the velocity of the cell, whereas m = n2 n1 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 5.1 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 u T = (u 1 − v cell , u 2 ) φ(x,t)=0 fort > 0 such that the stationary profile and velocity is reached. Figure 5   Polarization field (first row), velocity field in co-moving domain (second row). With increasing activity, the splay instability is enhanced, which moves the maximum of the velocity field along the interface (indicated by red angular arrow) from the front to the rear, visible also through the position of the vortices in the cell (indicated by blue arrow), which are located more towards the front for puller dynamics and more towards the rear for pusher dynamics.
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 [36] 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 center of mass x cm = (0, 0) in a co-moving frame, given by where r = x/r is the polar axis, scaled with the distance r = x 2 1 + x 2 2 , e 1 the unity vector in x 1 -direction and I the identity matrix. We prepared our numerical solution: u = (u 1 − v cell , u 2 ) T , x = (x 1 − x cm , x 2 − y cm ) T and claim |u − v| → min outside the circular cell shape with radius R = 5 to determine v cell , A st , A str and A sd . Table 5.1 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 analyzing the relative error |u − v|/v cell between the numerical results and the fitted analytical solution, see  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 [35] and are dominated by the Stokeslet contribution. In [36] the same fitting approach is used to analyze 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.
6. Conclusion. We here review and extend a proposed generic model for cell motility [14], 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. 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 fitted analytical solution |v| (middle) and the relative error (right) for 1/Fa=0.5 (first row), 1/Fa=0.75 (second row) and 1/Fa=1 (third row). For the analytical solution as well as the error analysis we approximated the cell shape by a circle, with radius R = 5 obtained from the initial condition. The data fit indicates 1/Fa=0.5 as puller and 1/Fa=1 as pusher. For 1/Fa=0.75 we get a puller like velocity profile, where we expected a neutral swimmer, which of course can be a result of the approximated circular shape of the cell. the physical instabilities, the splay or bend instability, which are responsible for cell motility in the considered model, and possible numerical instabilities. 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 microorganism is not yet possible and besides the lack of available experimental data, requires probably further model extension. One possible way to extend the model is a combination of the myosin-acting interactions with the treadmilling process of acting polymerization and depolymerization, described in the introduction. However, qualitative similarities with generated flow fields of microorganisms, such as Volvox carteri could already be found. In both cases, the flow field is dominated by the Stokeslet contribution. In [36] it is argued that this behavior is going to have an effect on the rheology of suspensions of such microorganisms. With these properties, suspensions of our modeled cells would probably behave more like suspensions of sedimenting particles.