Three-dimensional phase field model for actin-based cell membrane dynamics

The interface dynamics of a 3D cell immersed in a 3D extracellular matrix is investigated. We suggest a 3D generalization of a known 2D minimal phase field model suggested in [1] for the description of keratocyte motility. Our model consists of two coupled evolution equations for the order parameter and a three-dimensional vector field describing the actin network polarization (orientation). We derive a closed evolutionary integro-differential equation governing the interface dynamics of a 3D cell. The equation includes the normal velocity of the membrane, its curvature, cell volume relaxation, and a parameter that is determined by the non-equilibrium effects in the cytoskeleton. This equation can be considered as a 3D generalization of the 2D case that was derived in [2].


I. INTRODUCTION
Over the past decade, increasing attention has been paid to the formulation of computational models that describe the motility of 3D cells that crawl on flat substrate [3], [4] or invade 3D extracellular matrices (ECM) [5], [6], [7], which is the subject of this paper. This challenge has been mentioned about a decade ago in the review paper [8] as barely started field. Since that time, the keratocyte motility, including the lamellipodium waves dynamics [9], has been rather well explored using experimental and theoretical approaches.
Actin polymerization is a basic mechanism in 2D cell motility [10], [9]. When investigating the dynamics of 3D cell invasion in 3D ECM, one finds that there are several mechanisms or modes of migration that control 3D cell motility, and a cell may switch between them depending on cell intrinsic and extrinsic factors [11]. Among these mechanisms is the generation of protrusive force through hydrostatic pressure. This mechanism is called amoeboid and it does not need actin polymerization at the leading edge to generate protrusions. Another mechanism, which is the subject of this paper, is the Lamellipodium-based protrusion in 3D-ECMs which is called mesenchymal. This mechanism is based on actin polarization to form protrusion that enables cell migration in 3D ECMs, and it is analogous to the 2D actin-based cell motility. This migratory mode has been observed in some metastatic cancer cells while moving in complex 3D environments, see [11] and references therein.
Recently a minimal computational phase field model of 3D cell crawling on general substrate topography, not only on flat surface, has been formulated [12]. In their model they assume that the actin exists only nearby the substrate surface and vanishes far away. Therefore this model is a nontrivial generalization of their 2D model [1].
In the present paper we suggest a 3D generalization of the original 2D model developed by Aranson and coworkers [1], [13] (see also [2]). It is assumed that the actin exists over all the interface of the 3D cell, hence that model can and be used as a model of 3D cell surrounded by 3D ECM.
In the framework of our model, the cell geometry is de-scribed by a phase field coupled with a three-dimensional vector field of the actin network polarization. Thus, we suggest a model of actin-based 3D cell motility surrounded by 3D ECM. Relative to other phase field approaches, this model can be considered as a simple minimal model describing the 3D cell motility [14], [12].
In the present paper, we consider the case where the cell does not move as a whole but can change its shape. The structure of the paper is as follows: In Sec. II we present the minimal 3D phase field model. In Sec. III we investigate the dynamics of the spherical shape interface. In Sec. IV we consider the general shape interface. We derive a closed evolutionary nonlocal equation that describes the interface dynamics. Also, we investigate the stability of the spherical shape membrane. Finally, in Sec. V we present the conclusions.

II. FORMULATION OF THE NONLOCAL PROBLEM
We consider the following minimal phase-field model of self-polarization and motility of spherical cell shape. This model is a 3D generalization of the 2D model developed in the context of keratocyte motility in [1] and [13] (see also [2]).
where u(r, θ, ϕ, t) is the order parameter that is close to 1 inside the cell and 0 outside. The interface is defined by the relation u(r = ρ(θ, ϕ, t)) = 1/2. The three-dimensional polarization vector field P(r, θ, ϕ, t) = pr + qθ + wφ represents the actin orientations. It is assumed that P(r = 0) is close to zero, and the cell does not move as a whole. See [2] for more details about the formulation of this simplified version of the full model that was developed in [1].
The model contains several constant parameters: D u is the stiffness of diffuse interface, D p is the diffusion coefficient for P, α is the coefficient characterizing advection of u by P, β determines the creation of P at the interface, τ −1 is the inverse time of the degradation of P inside the cell, V 0 = 4πρ 3 0 /3 is the overall volume of the cell, µ is the stiffness of the volume constraint, and σ is the contractility of actin filament bundles. All the parameters listed above are positive. Notice that the model (1a)-(1e) is nonlocal due to the definition of δ.
Because of the spherical symmetry of the problem, we employ the spherical coordinate system with the corresponding differential operators (see Fig. 1), x = r cos θ cos ϕ, y = r cos θ sin ϕ, z = r sin θ, (2b) In the next sections we use arguments similar to those we used in the 2D case [2]. We begin our analysis with consideration of the spherical cell dynamics.

III. ROTATIONALLY SYMMETRIC CASE
Assume that our fields have a rotational symmetry, i.e., u = u(r, t), and P = p(r, t)r.
In the present paper, the basic assumption is that the ratio of the thickness of the cell wall (i.e., the width of the transition zone, where u(r, t) is changed from nearly 1 to nearly 0) to the size of the cell is small, see Fig. 1. In that case, the nonlocal term in (1b) can be estimated as The model (1a)-(1e) takes the form, We introduce the scaling that describes slow dynamics of large enough cell radius, and define the transition zone variable, Also we define, Consequently the chain rule yields It holds that The scaled equations (3a)-(3e) are, Because the motion of the front is influenced by its curvature 1/R, the term R t u ζ in (5a) should balance the curvature term 2u ζ /R(t) in the same equation. That justifies the choice of the time and radius scaling (4). Assume that the parameters have the scaling and introduce the expansions Substituting (6) and (7) into (5) and collecting terms of the same order, we obtain at the leading order system (8), Following the Ginzburg-Landau theory and Fourier transform method, we find the solutions of the system, See Fig. 2 for the plot of the fields in (9). The equation for u 1 at the order O( ) is The solvability condition of the latter equation yields the following closed equation for the interface dynamics R(t), see [2] for more details. The expression (9a) yields, (11) Hence equation (10) may be written in the form , where, It is more convenient to consider the following presentation of equation (12) of the front dynamics, The expression in the right-hand side of (15a) has the meaning of a "force" acting on the cell surface. Therefore, equation (15a) can be considered as a "force-velocity relation" of the motionless 3D cell immersed in 3D ECM (see [10]).
Recall that Ω 1 is positive and Ω 2 is negative, therefore, Ω(β) is a monotonically growing function of β, see (13), (14). The nonlinearity of the dependence of Ω on β is caused by the nonlinear term −σp 2 in the expression (5b) for δ. The physical origin of that term is the contraction of actin filament bundles.
In Fig. 3 we plot the function f (R) for two values of β; the graph shows the existence of two stationary radii, stable, R + (β) and unstable, R − (β), for β = 1, while no stationary states for β = 0.4. One can conclude that there exists a critical value β c such that there are no stationary solutions when β < β c . Below we find that β c = 0.637 for the chosen set of parameters.
The critical value β c has to satisfy three constraints: Fig. 4(b)). As a result we find that β c is the positive solution of the quadratic equation that can be found explicitly: notice that R * does not depend on R 0 ; for values of parameters indicated in Fig. 2, R * = 1.041 . Because Ω(β) is a monotonically growing function of β, R − (β) decreases with the growth of β, and therefore R − (β) < R * for any β > β c . On the contrary, R + (β) increases with the growth of β. For values of parameters indicated in Fig. 2 and R 0 = 1, we find β c = 0.637 . If β < β c , then f (R) < 0 for any R, therefore the cell radius decreases with time and tends to zero during a finite time (see Fig. 4(a)). The temporal evolution of R(t) in the case of β > β c depends on the relation between the initial radius R 0 and R − (β). If R 0 > R − (β), then R(t) → R + (β) as t → ∞; if R 0 < R − (β), then R(t) → 0 during a finite time. Therefore, if R 0 > R * , R(t) tends to a finite value for any β > β c , because R − (β) < R * < R 0 . However, if R 0 < R * , the cell shrinks even at β > β c , if β is still less than a certain valueβ (see Fig 4(c)), which is determined by the relation R 0 = R − (β) (see Fig. 4(d)). Because f (R 0 ,β) = 0, the value ofβ can be found by solving the equation, For β >β, R − (β) < R 0 , therefore the cell radius tends to a finite value (see Fig. 4(e)).
In Fig. 5 we present the numerical solution of the ODE (15a) for the values β = 1 . As we can see, the cell radius can increase monotonically until it reaches the steady state value. This is because in the framework of our model, the volume is not conserved and influences the dynamics through the parameter δ(r, t), see equation (5b).
The decrease of the cell radius (in the language of the phase-field model, the transition of the phase u = 1 into the phase u = 0) is caused by negative terms in the righthand side of (1a),(3a), among them the term −δ(1 − u)u, which is negative at large R, and by the diffusion term, which creates an effective "surface tension" of the cell surface. The positive term −αu r p hinders the decrease of the cell radius. If β is not sufficiently large, the polarization is not strong enough to stop the collapse of the cell. In that case, the cell shrinks until it disappears (see Fig. 3, β = 0.4.) Note that at large β (see Fig. 4(e)) we have Ω(β) 1, and R − (β) 1, therefore one can approximate the expression of f (R) in (12) as: .
We can see that for sufficiently large β, R − (β) can be arbitrary small. Therefore, for arbitrary small R 0 , there exists suchβ that the cell radius tends to a finite value, if β >β.
For values of parameters indicated in Fig. 2 and R 0 = 0.5,β = 0.96 . In Fig. 3 and 5, we show the numerical plots for the parameters β = 1, and R 0 = 0.5.
Let us emphasize that the described effect is not physical: for such values of parameters, the model does not reflect the true behavior of cells.

IV. DYNAMICS OF GENERAL SHAPE INTERFACE
We employ the scaling and definitions of the previous section. Now we have to consider the azimuthal dependence of variable. One can calculate, , The nonlocality in (1b) is approximated as follows, We use expansions (7),and define the auxiliary function Λ, At the leading order one can fined, therefore similarly to the previous section one can calculate the solutions The equation for u at the order O( ) have the form, where the volume variation has the form, We apply the solvability condition, which is the orthogonality of the equation's right-hand side to the solution of the homogenous equation u 0ζ , and obtain a closed form of the interface dynamics, where is the mean local curvature of the surface r = R(θ, ϕ, t).
For an explicit expression of the curvature see Appendix A. Note that in the spherically symmetric case Λ = 1, . Therefore, equation (12) is recovered from equation (19). Notice that Ω in equation (19), depends on all of the parameters that describe the nonequilibrium molecular effects of the subcell level, see (13), and (14). Equation (19) is a closed evolutionary equation for the 3D cell interface dynamics , which is an integro-differential equation, i.e., it is nonlocal, unlike that obtained in the spherical case (15a). For the details of the application of the solvability condition in order to obtain equation (19), we refer the reader to [2] and [15], where we perform similar calculations.
Notice that in (19) the expression ΛR t corresponds to the normal velocity v n of the interface, thus that equation is a generalization of the well-known curvature flow. By a proper scaling transformation, t → a 2 t and R(t) → aR(t), the equation of motion of the cell boundary can be brought to a canonical form, v n = −2D u H −Ṽ + Ω.
In addition it suggests an answer for the unrevealed force -velocity relation for the actin network that was highlighted in [10] in the context of shape dynamics of a 2D cell.

A. Stability of the radial interface
Let our base radial solution in (15a) be perturbed R =R(t) +R(θ, ϕ, t),R R .
The linearization of (19) around that base solution yields, cos θdθ.
First, let us consider solutions satisfying the condition Let us consider the normal modeR(θ, φ, t) = T (t)Θ(θ)e imϕ . For m = 0 the integral term in (21) vanishes thus we do not have a contribution of the volume variation. For m = 0 it vanishes due to condition (22) . Applying the separation of variable method one obtain, The solution Θ is bounded if 2+µ = n(n+1); in the latter case we obtain the spherical harmonics solution with the associated Legendre function, R = e imϕ P m n (− sin θ)T n (t), n = 0, 1, 2, ..., m = 0, 1, ..., n.
All solutions with n = 0 satisfy condition (22) due to the orthogonality property of the spherical harmonics, therefore T n (t) satisfy, For n = 1, which corresponds to the spatial translation of the sphere as a whole, we find that T 1 (t) = const. Disturbances with n ≥ 2, which describe the shape distortions, decay with time.
In the case n = m = 0, which corresponds to a change of the sphere radius, the integral I = π/2 −π/2 Θ(θ) cos θdθ = 0, thus from equation (21) we obtain the equation, Let us divide both sides of the equation by T Θ. We can see that both sides of the obtained equality are functions only of t: (23) Multiplying both sides of (23) by cos θ, integrating over θ from −π/2 to π/2, and dividing by I = 0, we find that The analysis of the sign of the expression in the righthand side of (24) confirms the result obtained in Section III: solution R = R − is unstable and solution R = R + is stable with respect to the radius change. T 0 can grow, but whenR(t) approaches its stationary values, the derivative d(f (R)/R)/dr at R =R(t) becomes negative (see Fig. 3), hence the spherical cell is stable with respect to spherical disturbances. Note that these results are similar to the 2D case [2], where we find that the circular cell shape is stable concerning a small disturbance.

V. CONCLUSION
We perform the analysis of a minimal phase field model that is a 3D generalization of the 2D model developed and investigated numerically in [1], [13], and [14] (a similar analysis of the latter model was done in [2]). In this model the order parameter u is coupled with 3D polarization (orientation) vector field P of the actin network. The model is supposed to describe the 3D cell motility immersed in 3D ECM via actin based protrusion mechanism [11].
We considered the rotational symmetric case i.e., spherical shape interface, where we obtained a closed ordinary differential equation describing the evolution of the radius (15a). We found the minimum value β c for the actin creation that is compatible with the existence of a stationary cell solution (16). We found that when β c < β, the circular cell can have some stationary radius, while in the case β ≤ β c the cell shrinks until it disappears, which is meaningless in the context of cell dynamics. Also, we considered the general shape 3D cell dynamics. We found the leading order solutions, (17a)-(17c), and derived a closed integro-differential equation (19) governing the 3D cell dynamics, which includes the normal velocity of the membrane, curvature, volume relaxation rate, and a parameter Ω determined by the molecular effects of the subcell level. This result is similar to the 2D case [2].
We found an equation of motion of the cell interface that can be written in the canonical form, v n = −2D u H −Ṽ + Ω.
The stability analysis shows that the non-spherical shape and the motion of the cell as a whole cannot appear due to the development of a linear instability of the spherical cell with the spherically symmetric polarization field localized near the cell boundary. In the framework of the considered model, the transition to a non-spherical shape needs a finite-amplitude disturbance significantly changing the polarization field inside the cell. The analysis of that transition, which can be carried out only numerically, is beyond the scope of the present paper.