The moving grid finite element method applied to cell movement and deformation

Finite Elements in Analysis and Design 74 (2013) 76–92

Contents lists available at SciVerse ScienceDirect

Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel

The moving grid finite element method applied to cell movement and deformation Anotida Madzvamuse n, Uduak Zenas George School of Mathematical and Physical Sciences, Department of Mathematics, University of Sussex, Brighton BN1 9QH, United Kingdom

art ic l e i nf o

a b s t r a c t

Article history: Received 23 September 2012 Received in revised form 27 February 2013 Accepted 5 June 2013 Available online 7 July 2013

In this article we present a novel application of the moving grid finite element method [1] for solving a cytomechanical model that describes actin dynamics in order to generate cell movement and deformation. The cytomechanical model describes both the mechanical and biochemical properties of the cortical network of actin filaments and its concentration. Actin is a polymer that can exist either in filamentous form (F-actin) or in monometric form (G-actin) [2] and the filamentous form is arranged in a paired helix of two protofilaments [3]. By assuming slow deformations of the cell, we validate numerical results by comparing qualitatively with predictions from linear stability theory close to bifurcation points. Far from bifurcation points, the mathematical model and computational algorithm are able to describe and generate the complex cell deformations typically observed in experiments. Our numerical results illustrate cell expansion, cell contraction, cell translation and cell relocation as well as cell protrusions. A key model bifurcation parameter identified is the contractile tonicity formed by the association of actin filaments to the myosin II motor proteins. The robustness, generality and applicability of the numerical methodology allows us to tackle similar problems in developmental biology, biomedicine and plant biology where similar mechanisms are routinely used. & 2013 Elsevier B.V. All rights reserved.

Keywords: Moving grid finite element method Moving boundary problem Lagrangian kinematics Cell movement Cell deformation Cell motility Reaction–diffusion equations

1. Introduction The migratory behavior of cells plays a crucial role in many biological events (e.g. physiological and pathological process [4]) such as immune response [5–7], wound healing, development of tissues [8], embryogenesis [9,10], inflammation and the formation of tumor metastasis [7,11]. According to experimental observations, actin plays an important role in cell movement and deformation. Actin is a polymer that can exist in two forms; in monomeric form as G-actin or in polymeric form (i.e. in filamentous form) as F-actin [12]. The mechanical properties of actin network and the biochemical dynamics of actin (i.e. the concentration of F-actin, its random movement via diffusion, convective and dilution effects due to shape movement and its kinetics of polymerization and de-polymerization from F-actin to G-actin) inside the cell play an active role in cell movement and deformation [12–15]. In order to describe the dynamics of movements and changes in cell geometry of a eukaryotic cell we will consider a

n

Corresponding author. Tel.: +44 1273873259. E-mail addresses: [email protected] (A. Madzvamuse), [email protected] (U. Zenas George). 0168-874X/$ - see front matter & 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.finel.2013.06.002

cytomechanical model that is derived based on the actin dynamics. The model is highly complex and consist of a system of nonlinear equations defined on a complex evolving geometry representing the cell. We propose to apply a novel numerical method: the moving grid finite element method (MGFEM) [1,17] to compute numerical approximate solutions of the cytomechanical model. A determining factor in choosing a numerical method is usually the ease in which a method can be applied to a problem at hand and the ability of the method to produce highly accurate solutions in comparison with other available methods. The key factors being accuracy, efficiency and robustness of the numerical method Bindu [16,17,20,25]. To the best of our knowledge, this article is the first to present the use of finite elements to solve the cytomechanical model system. To-date, numerical approximation of the cytomechanical model has been done by Alt and Tranquillo [15] and Stéphanou et al. [12] based on finite differences. Rather than working in the original cartesian coordinates, the authors transform the model system into a polar coordinate system which makes it difficulty to describe the evolution of a cell as it moves outside the origin of the polar coordinate system. The cell is not able to move or deform away from the origin [12]. The origin remains a fixed point, the cell is only able to deform along the boundary but is not able to move freely on the substrate as is observed in experiments [12]. In order

A. Madzvamuse, U. Zenas George / Finite Elements in Analysis and Design 74 (2013) 76–92

to remedy this drawback, we propose to employ the MGFEM thereby enabling us to simulate accurately cell movement and deformation in complete agreement with experimental observations. Other alternative methods that could be used to overcome difficulties in solving the cytomechanical model using [18,19] finite differences include (but are not limited to) boundary element methods (BEM) [22,23] and level set methods (LSM) [24]. Since the cytomechanical model is a second order inhomogeneous and nonlinear system of partial differential equations, the application of the BEM is much more difficult compared to FEM [26]. The LSM will not be implemented here but could be implemented in the future by coupling the method to the MGFEM in order to facilitate modeling splitting and reconnecting cells. In this article, we choose to implement a finite element formulation by employing the MGFEM [1,17,29] for two advantageous reasons: (i) the MGFEM solves the model equations in their physical cartesian coordinates without the need for transformations or mappings and (ii) the MGFEM easily handles complex and evolving cellular domains. As a result, the cell is able to move and deform freely on the substrate in complete agreement with experimental observations. In view of these key advantages, we develop a MGFEM of the cytomechanical model that is defined on a two-dimensional (2D) cartesian coordinate system. The MGFEM is a novel extension of the finite element method to moving boundary problems introduced by Madzvamuse [27] to study partial differential equations posed on complex evolving domains. The essential component underlying the MGFEM is that the object function is represented by piecewise basis functions on a finite element mesh in which the nodal positions vary with time [1,17,29]. It is a highly accurate, efficient and robust numerical method which has been used successfully to compute numerical solutions of reaction–diffusion systems on evolving and continuously deforming domains with complex geometries (see [1,17,21,28] for details). The novelty of the MGFEM is its ability to allow the prescription of the nodal displacement of the computational grid points of the finite element mesh during the evolution of the domain. It is able to accommodate a priori description of the grid displacement [27], a computed grid displacement (which are solutions of a partial differential equation of an unknown variable) or a grid displacement generated by minimizing the residual of the differential equation [29,30] (as developed by [30] in the moving finite element). The difference between the moving finite element and MGFEM is the derivation of the grid displacement [21]. Therefore, this article is structured as follows: In Section 2 we present the cytomechanical model. In Section 3 we derive the weak formulation for the cytomechanical model and then present the finite element discretization of the model system. We also describe the method of computation of the domain evolution and the techniques of implementation of the MGFEM. Since the cytomechanical model is a highly complex nonlinear system of partial differential equations, analytical solutions in closed form can not be derived. To this end, it is difficult to compare numerical solutions to exact analytical solutions. Therefore, in order to validate numerical results, we call upon linear stability theory and show that close to bifurcation points, the numerical method computes solutions in close agreement with those predicted by linear theory. Details of linear stability theory are presented in Appendix B. Linear theory also provides us with plausible ranges for the parameter values. Once validated, numerical results obtained from the implementation of the MGFEM are shown in Section 4. Finally, in Section 5 we conclude and discuss future research strands.

2. A cytochemical model for cell movement and deformation Following George et al. [31], let Ωt ⊂R2 be a simply connected bounded continuously deforming domain representing the Eukaryotic

77

cell shape at time t∈I ¼ ½0; T f , T f 4 0 and ∂Ωt be the boundary describing the cell. At any given point, x ¼ ðxðtÞ; yðtÞÞ∈Ωt ⊂R2 , let a ¼ aðxðtÞ; tÞ be the F-actin concentration and u ¼ ðuðxðtÞ; tÞ; vðxðtÞ; tÞÞT be a vector of displacement of the elements of the actin network at position x∈Ωt at time t∈I. As a result of cell movement we define β to represent the flow velocity of the elements of the actin network and ωn the normal velocity of the boundary. Furthermore denote n ¼ ðn1 ; n2 Þ the outward unit vector normal to the boundary. We note that by defining u to be a vector of displacement of the elements of the actin network at position x∈Ωt at time t∈I, the following system of equations describing the actin dynamics as given by Stéphanou et al. [12] models both plasticity and viscoelasticity of the actin network: for all t∈I, ∇  ðrv þ re þ rc þ rp Þ ¼ 0

in Ωt ;

∂a −Da Δa þ ∇  ðaβÞ−ka ðac −aÞ ¼ 0 ∂t uðxðtÞ; tÞ ¼ 0 aðxðtÞ; tÞ ¼ a0 β ¼ ωn

∀x∈Ω0 ;

in Ωt ;

ð2:1bÞ ð2:1cÞ

∀x∈Ω0 ;

ð2:1dÞ

∀x∈∂Ωt ;

ðrv þ re Þ  n ¼ n  ∇a ¼ 0

ð2:1aÞ

ð2:1eÞ ∀x∈∂Ωt ;

ð2:1fÞ

where Ω0 is the initial cell domain at time t¼0. In (2.1a), rv , re , rc and rp are the viscous, elastic, contractile and pressure component stress tensors. These are defined as 8 rc ¼ sðaÞI ¼ ψa2 e−a=asat I; > > > > ∂ϵ ∂ϕ > > > > rv ¼ μ1 ∂t þ μ2 ∂t I; > <  E
ν ð2:2Þ re ¼ ϵþ ϕI ; > > 1 þ ν 1−2ν > >   > > > p 2 > > 1 þ δðlÞ arctan a I; : rp ¼ 1þϕ π where ϵ ¼ 12 ð∇u þ ∇uT Þ represents the strain tensor, I is the identity tensor, ϕ ¼ ∇  u is the dilation, and μ1 and μ2 are the shear and bulk viscosities of the actin network respectively. Finally E and ν are the Young's modulus and Poisson ratio respectively. The function sðaÞ represents the contractile activity of the actomyosin network [31]. In the above formulation, 2asat represents the saturation concentration of F-actin. In order to describe the pressure forces rp acting within the cell we make the following assumptions. We assume that the initial domain is a unit disk which we denote as Ω0 and that there exists a family of bijective mappings that maps the points ξ ¼ ðξx ; ξy Þ of the initial domain to point x on the current domain Ωt . Let l : Ωt  I-R and it's corresponding function on the initial domain be ^l : Ω0  I-½0; 1 where ^lðξ; tÞ is the distance between the centroid and the point ξ in the initial domain with ^lðξ; tÞ ¼ lðxðξ; tÞ; tÞ. The function pðaÞ≔p=1 þ ϕð1 þ ð2=πÞδðlÞ arctan aÞ describes the pressure force in different regions governed by the heavy-side function: 8 > 1 if the point ðx; tÞ with lðxðξ; tÞ; tÞ ¼ ^lðξ; tÞ is such that > < qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi δðlÞ ¼ the distance ξ2x þ ξ2y 4 0:8 in the initial domain; > > : 0 elsewhere: ð2:3Þ This signifies that far from the membrane, only the osmotic component p=ð1 þ ϕÞ of the pressure force exists within the network and depends q onffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi the dilation ϕ. In the vicinity of the membrane (i.e. where ξ2x þ ξ2y 4 0:8), a polymerization-induced pressure which depends on the local actin density reinforces the osmotic stress and pushes the membrane out at regions where the filaments are not firmly linked to the membrane [31]. δðlÞ is used here to specify the differences in pressure in the cell at the vicinity

78

A. Madzvamuse, U. Zenas George / Finite Elements in Analysis and Design 74 (2013) 76–92

Table 1 Dimensional parameters and their values as used in the mathematical model except where it is specified otherwise. Parameter

Meaning

Value

Ref.

E ν Da ka asat μ1 μ2 l0 ac

Young's modulus of the actin network Poisson's ratio of the actin network Diffusion coefficient of actin Polymerisation rate of the actin network Saturation concentration of F-actin Shear viscosity of the actin network Bulk viscosity of the actin network Specifies the vicinity of the membrane F-actin concentration at the chemical equilibrium

1.5 dyn/cm2 0.3 0.012 cm2/s 0.03 s−1 1.4 normalized density 96:15 dyn s=cm2 250 dyn s=cm2 80% of the cell radius

Estimated Estimated Stéphanou et al. [12] Estimated Stéphanou et al. [12] Bausch et al. [32] Bausch et al. [32] Estimated

1.0 normalized density

Derived from model

of the membrane and further away from it. We have assumed that the polymerization of actin, which happens predominately at the vicinity of the membrane, induces a pressure that reinforces the osmotic pressure at the vicinity of the membrane. Eq. (2.1b) describes the conservation of F-actin within the cell. It describes the concentration of F-actin, its random movement via diffusion, convective and dilution effects due to the movement of the cell and its kinetics of polymerization and depolymerization. In (2.1b), ac is a constant parameter representing the F-actin concentration at the chemical equilibrium. It differentiates the states of polymerization and depolymerization. The polymerization rate is given by ka and Da is a positive diffusion coefficient for F-actin. The reaction kinetics ka ðac −aÞ, models the actin dynamics of polymerization and depolymerization. The actin polymerization and depolymerization are considered to occur at the same rate ka. Thus the state of polymerization or depolymerization depends on the local value of F-actin concentration relative to the concentration at the chemical equilibrium ac [55]. On the boundary, we take the flow velocity of the actin network β ¼ ωn ≔∂u=∂t. Here u is the displacement of the actin network at any given point in space and is computed as the displacement solution obtained from the force balance equation (2.1a). In the interior of the domain, we assume that β ¼ ∂x=∂t≔∂u=∂t where ∂x=∂t is the mesh velocity. Thus we have a Lagrangian description of the domain evolution. For initial conditions, we define the initial domain Ω0 , t¼0, to be a unit disk. We prescribe the initial conditions for F-actin density a0 to correspond to random small perturbations around the homogeneous steady state of the cell. This state biologically corresponds to the cell condition right after mitosis (i.e. cell division) with the cell having a perfectly circular shape. The steady state of F-actin is assumed to occur as a result of an equilibrium in the kinetics of polymerization/depolymerization of actin within the cell and there is no flux of actin from the exterior. We assume that at the initial time, the cell is unstrained from its original position. The boundary conditions in (2.1f) specify a zero flux boundary for the reaction–diffusion equation (2.1b) and a stress free boundary for the force balance equation (2.1a). Remark 2.1. We note that in a previous work by Stéphanou et al. [12], the displacement uðxðtÞ; tÞ was taken to be the displacement of the elements of the actin network from the original unstrained position. In this work, we compute the displacement uðxðtÞ; tÞ continuously given the previous positions of x∈Ωt , at time t∈I. An advantage of this modification is that the cell is no longer rigid but can translate, expand and contract. Model parameter values are chosen close to those available in the literature. For example, in Stéphanou et al. [12], a diffusion coefficient Da in the range 0:00962–0:134 cm2 =s and an actin saturation concentration asat of 1.1 normalized density were used. Viscosity of Fibroblast cells is given in Bausch et al. [32] to be

200 dyn/cm2. Here we choose the value of the shear viscosity to be 96:15 dyn s cm2 and the value of the bulk viscosity of the actin network μ2 ¼ 250 dyn s=cm2 to be larger than that of the shear modulus. This is because F-actin are more resistant to compression than shear [33]. The rest of the dimensional model parameter values are shown in Table 1. It must be noted that we have not assigned parameter values for pressure coefficient p and the contractile coefficient ψ, these will be determined mathematically by use of linear stability as illustrated in Appendix B.

3. The moving grid finite element method (MGFEM) The essential component underlying the moving grid finite element method (MGFEM) is that the object function is represented by piecewise basis functions on a finite element mesh in which the nodal positions vary with time [1,17,29]. The way in which these nodes are moved and how the solution evolves is derived in the weak form of the model system. Unlike previous studies [1,17] where the nodal movement was provided by a separate definition, here nodal movement is an unknown quantity to be computed simultaneously with the actin concentrations. The MGFEM is an extension of the finite element method to moving boundary problems. Thus from Sections 3.1–3.80, we describe the finite element discretization of the coupled system (2.1) defined on a moving cell and in Section 3.9, we present the method for the cell domain evolution. We note that the MGFEM applied to reaction–diffusion equations has already been established [21] but this is not the case for the force balance equation. Our aim is to present in detail the MGFEM applied to the cytomechanical model by deriving a MGFEM of the force balance equation and then couple it to that of the reaction–diffusion equation. Reynolds transport theorem. Here we state the Reynolds transport theorem. This theorem shall be useful in the derivation of the weak formulation for the reaction–diffusion equation. Theorem 3.1 (Reynolds transport theorem). Let gðx; tÞ be a scalar function defined on Ωt and β be a flow velocity field then d dt

Z 

Z Ωt

g dΩt ¼

Ωt

 Dg þ g∇  β dΩt Dt

ð3:1Þ

A proof of this theorem can be found in [27]. To proceed, we decouple the force balance equation into a system of two partial differential equations for ease of implementation of the weak formulation. We substitute the values of rv , re , rc and rp as defined in (2.2) into the force balance equation (2.1a). This necessitates the expression of the stress tensors in tensormatrix forms. It follows that the dilation ϕ by definition is

A. Madzvamuse, U. Zenas George / Finite Elements in Analysis and Design 74 (2013) 76–92

ϕ≔∇  u ¼ ∂u=∂x þ ∂v=∂y and the stain tensor ϵ is 0 1 ∂u 1 ∂v ∂u ∂x 2 ð ∂x þ ∂y Þ 1 T A ϵðuÞ ¼ ð∇u þ ð∇uÞ Þ ¼ @ 1 ∂v ∂u ∂v 2 2 ð ∂x þ ∂y Þ ∂y

∈ R22 :

By using the above representation of the strain tensor ϵ and the dilation ϕ, we represent the stress tensors rc , rp , re and rv (defined in (2.2)) in two-dimensional tensor-matrix forms as members of R22 as follows: ! ψa2 e−a=asat 0 rc ¼ ; 0 ψa2 e−a=asat 0 1 μ1 ∂v_ ∂u_ ðμ1 þ μ2 Þ ∂∂xu_ þ μ2 ∂∂yv_ 2 ð ∂x þ ∂y Þ A; rv ¼ @ μ1 ∂v_ ∂u_ μ2 ∂∂xu_ þ ðμ1 þ μ2 Þ ∂∂yv_ 2 ð ∂x þ ∂y Þ 0 ∂u ∂v E′ @ ðν′ þ 2Þ ∂x þ ν′ ∂y re ¼ ∂v ∂u 2 ∂x þ ∂y and

0

rp ¼ @

p 1þϕ ð1

∂v ∂x

∂v ν′ ∂u ∂x þ ðν′ þ 2Þ ∂y

þ 2π δðlÞ arctan aÞ

A;

0 p 1þϕ ð1

0

1

þ ∂u ∂y

þ 2π δðlÞ arctan aÞ

1 A:

Substituting the matrix tensor representation of the stress tensors into (2.1a) we obtain the force balance equation in matrix tensor form which is then decoupled into a system of two equations. The decoupled system of differential equations of the force balance equation is thus        ∂ ∂u_ ∂v_ ∂ ∂u_ ∂v_ ∂ ∂u ∂v D11 þ D12 þ D33 þ þ C11 þ C12 ∂x ∂x ∂y ∂y ∂y ∂x ∂x ∂x ∂y þ

   ∂ ∂u ∂v ∂f C33 þ ¼− 1; ∂y ∂y ∂x ∂x

ð3:2aÞ

and        ∂ ∂u_ ∂v_ ∂ ∂u_ ∂v_ ∂ ∂u ∂v D33 þ þ D12 þ D22 þ C12 þ C22 ∂x ∂y ∂x ∂y ∂x ∂y ∂y ∂x ∂y    ∂ ∂u ∂v ∂f þ C33 þ ¼− 2; ð3:2bÞ ∂x ∂y ∂x ∂y where f 1 ≡f 2 ¼

   p 2 1 þ δðlÞ arctan a þ ψa2 e−a=asat ; 1þϕ π



D11 ¼ D22 ¼ μ1 þ μ2 ; C11 ≡C22 ¼

D12 ≡D21 ¼ μ2 ;

Eð1−νÞ ; ð1 þ νÞð1−2νÞ

C12 ≡C21 ¼

D33 ¼ μ1 =2; Eν ; ð1 þ νÞð1−2νÞ

and

79

t∈I such that      Z ∂w1 ∂u_ ∂v_ ∂u ∂v ∂w1 ∂u_ ∂v_ þ D12 þ C11 þ C12 þ þ D11 D33 ∂x ∂y ∂x ∂y ∂y ∂x ∂y Ωt ∂x    Z ∂w1 ∂u ∂v ∂f þ dΩt ¼ þ C33 w1 1 dΩt ; ð3:3aÞ ∂y ∂x ∂y ∂x Ωt and      Z ∂w2 ∂u_ ∂v_ ∂u ∂v ∂w2 ∂u_ ∂v_ þ D22 þ C12 þ C22 þ þ D12 D33 ∂x ∂y ∂x ∂y ∂y ∂x ∂x Ωt ∂y    Z ∂w2 ∂u ∂v ∂f þ þ dΩt ¼ C33 w2 2 dΩt ; ð3:3bÞ ∂y ∂x ∂x ∂y Ωt for all w1 ðx; tÞ, w2 ðx; tÞ∈H 1 ðΩt Þ, t∈I. The integrals on the right hand side of the weak formulation (3.3) are difficult to evaluate. To show this we write them out explicitly as 8    R R ∂f ∂ p 2 > > 1 þ δðlÞ arctan a dΩt > Ωt w1 1 dΩt ¼ Ωt w1 > > ∂x 1 þ ϕ π ∂x > > > > R ∂ > 2 −a=asat > > Þ dΩt ; < þ Ωt w1 ∂xðψa e    ð3:4Þ R R ∂f 2 ∂ p 2 > > δðlÞ arctan a dΩt dΩ w ¼ w 1 þ > t Ωt 2 ∂y 1þϕ > Ωt 2 ∂y > π > > > > R ∂ > 2 −a=asat > > þ Ωt w2 ðψ a e Þ dΩt ; : ∂y since f 1 ≡f 2 ¼ p=1 þ ϕð1 þ ð2=πÞδðlÞ arctan aÞ þ ψa2 e−a=asat . In view of this, we state the following identities which can be derived using the gradient and divergence theorems and are useful in the sequel in rewriting the weak form in a computationally efficient form. Let rðx; tÞ and gðx; tÞ be scalar functions of class C 0 ðΩt Þ defined in Ωt ⊂R2 and also let n ¼ ðn1 ; n2 Þ denote the outward unit vector normal to ∂Ωt for time t∈I. Then the following identities (derived from the gradient theorem) hold [16]: 8R R ∂r R ∂g > > < Ωt r ∂x dΩt ¼ − Ωt g ∂x dΩt þ ∂Ωt n1 rg ds; ð3:5Þ R ∂g R ∂r R > > dΩt ¼ − Ωt g dΩt þ ∂Ωt n2 rg ds; : Ωt r ∂y ∂y where ds is the element of arclength. Using the identities (3.5) we can rewrite the weak form as follows: Find uðx; tÞ, vðx; tÞ∈H 1 ðΩt Þ, t∈I such that      Z ∂w1 ∂u_ ∂v_ ∂u ∂v ∂w1 ∂u_ ∂v_ þ D12 þ C11 þ C12 þ þ D11 D33 ∂x ∂y ∂x ∂y ∂y ∂x ∂y Ωt ∂x    Z ∂w1 ∂u ∂v ∂w1 þ dΩt ¼ − f 1 C33 dΩt þ ∂y ∂x ∂y ∂x Ωt Z þ n1 f 1 w1 ds; ð3:6aÞ ∂Ωt

3.1. Weak formulation of the force balance equation

and      Z ∂w2 ∂u_ ∂v_ ∂u ∂v ∂w2 ∂u_ ∂v_ þ D22 þ C12 þ C22 þ þ D12 D33 ∂x ∂y ∂x ∂y ∂y ∂x ∂x Ωt ∂y    Z ∂w2 ∂u ∂v ∂w2 þ dΩt ¼ − f 2 þ C33 dΩt ∂y ∂x ∂x ∂y Ωt Z þ n2 f 2 w2 ds; ð3:6bÞ

To derive the weak formulation of the force balance equation, we multiply the system of partial differential Eqs. (3.2a) and (3.2b) by the test functions w1 ðx; tÞ and w2 ðx; tÞ∈H 1 ðΩt Þ, t ∈I respectively, and using Green's formula [34] we integrate the partial differential equations over the domain Ωt and apply the boundary condition (2.1f). We note that the boundary condition (2.1f) implies that the boundary term vanishes after carrying out the integration over the domain. The weak formulation is thus: Find uðx; tÞ, vðx; tÞ∈H 1 ðΩt Þ,

for all w1 ðx; tÞ, w2 ðx; tÞ∈H 1 ðΩt Þ, t∈I. Here n1 and n2 are the direction cosines of the outward unit vector, n normal to ∂Ωt for time t∈I. The test functions we use are piecewise linear basis functions and their spatial derivatives are easy and straight forward to compute, the weak formulation as expressed in (3.6) is a lot easier to compute compared to (3.3). We say that uðx; tÞ, vðx; tÞ are weak solutions of the force balance equation (2.1a) if uðx; tÞ, vðx; tÞ∈H 1 ðΩt Þ and (3.11) holds. Here we have assumed that

C33 ¼

E : 2ð1 þ νÞ

∂Ωt

80

A. Madzvamuse, U. Zenas George / Finite Elements in Analysis and Design 74 (2013) 76–92

Z

aðx; tÞ is known and is the solution of the reaction–diffusion equation (2.1b). Remark 3.1. We would like to note that we do not actually compute the derivatives of the delta function but we transfer these derivatives to the test function by using identities (3.5). 3.2. Weak formulation of the reaction–diffusion equation We first re-write the reaction–diffusion equation for the actin biochemical dynamics in the following amenable form Da −Da Δa þ að∇  βÞ−ka ðac −aÞ ¼ 0 Dt

ð3:7Þ

where

denotes the material derivative of the actin concentration a. To obtain the weak formulation we multiply (3.7) by a test function w3 ðx; tÞ∈H 1 ðΩt Þ; t∈I, then we integrate by parts and apply the boundary conditions (see (2.1f)):  Z  Da þ a w3 ð∇  βÞ þ Da ∇a  ∇w3 þ ka a w3 dΩt w3 ∂t Ωt Z ¼ ka ac w3 dΩt : ð3:8Þ Ωt

Using the product rule, the differential equation (3.8) can be rewritten as  Z  Dðaw3 Þ Dw3 −a þ aw3 ð∇  βÞ þ Da ∇a  ∇w3 dΩt ∂t ∂t Ωt Z Z þ ka aw3 dΩt ¼ ka ac w3 dΩt : ð3:9Þ Ωt

Ωt

Using the Reynolds transport theorem [35], we rewrite (3.9) such that the weak formulation reads: Find aðx; tÞ∈H 1 ðΩt Þ; t∈I such that Z Z Z d aw3 dΩt þ ðDa ∇a  ∇w3 þ ka aw3 Þ dΩt ¼ ka ac w3 dΩt dt Ωt Ωt Ωt Z Dw3 dΩt ; ∀w3 ∈H 1 ðΩt Þ: þ a ð3:10Þ Dt Ωt

3.3. Weak formulation of the coupled problem The weak formulation of the coupled problem (2.1) is thus: find aðx; tÞ, uðx; tÞ, vðx; tÞ∈H 1 ðΩt Þ, t∈I such that      Z ∂w1 ∂u_ ∂v_ ∂u ∂v ∂w1 ∂u_ ∂v_ þ D12 þ C11 þ C12 þ þ D11 D33 ∂x ∂y ∂x ∂y ∂y ∂x ∂y Ωt ∂x    Z ∂w1 ∂u ∂v ∂w1 þ þ dΩt ¼ − f 1 C33 dΩt ∂y ∂x ∂y ∂x Ωt Z þ n1 f 1 w ds; ð3:11aÞ ∂Ωt



Ωt

a Ωt

Dw3 dΩt ; Dt

ð3:11cÞ

for all w1 ðx; tÞ, w2 ðx; tÞ, w3 ðx; tÞ∈H 1 ðΩt Þ, t∈I. 3.4. Finite-dimensional subspaces Let Ωh;t , t∈I be a bounded domain triangulated by T h;t . Each triangular partition is known as an element S. We use barycentric coordinates as the local coordinate system on the elements of the triangulation. Let P 1 be a finite dimensional function space defined on S, where S is a reference element: ( ) 3

S≔ ðλ0 ; …; λ2 Þ∈R3 ; 0 ≤λk ≤1; ∑ λk ¼ 1 k¼1

Da ∂a ¼ þ ð∇aÞ  β Dt ∂t

Z

þ

    ∂w2 ∂u_ ∂v_ ∂u ∂v ∂w2 ∂u_ ∂v_ þ D22 þ C12 þ C22 þ þ D12 D33 ∂x ∂y ∂x ∂y ∂y ∂x ∂y ∂x    Z ∂w2 ∂u ∂v ∂w2 þ þ dΩt ¼ − f 2 C33 dΩt ∂y ∂x ∂x ∂y Ωt Z þ n2 f 2 w2 ds; ð3:11bÞ ∂Ωt

and Z Z Z d aw3 dΩt þ ðDa ∇a  ∇w3 þ ka aw3 Þ dΩt ¼ ka ac w3 dΩt dt Ωt Ωt Ωt

then there exist a one-to-one mapping from S to S. Also let x1 ðtÞ; x2 ðtÞ; x3 ðtÞ denote the vertices of the element S. Then the following parameterization using barycentric coordinates over the element S can be defined [36], 3

xðλ1 ; λ2 ; λ3 ; tÞ ¼ ∑ λk xk ðtÞ:

ð3:12Þ

k¼1

Let wðxðtÞÞ : Ω-R, t∈I, be a finite element function on an element S defined by a finite dimensional function space P 1 defined on a reference element S and the mapping λS from the reference element S to S such that wðxðtÞÞ ¼ wðλS ðxðtÞÞÞ, where w is defined on S. We define the space X h ðtÞ⊂H 1 ðΩt Þ: n o X h ðtÞ ¼ w∈C 0 ðΩt Þ∩H 1 ðΩt Þ; w∈P 1 for all S∈T h;t ; t∈I :

3.5. Finite element spatial discretization of the cytomechanical model First let us approximate the continuously evolving cell domain with a finite polygonal domain as follows: at each time t, t∈I, we discretize Ωt into a finite unstructured triangular partition Ωh;t of non-overlapping triangles, where h is the maximum size of the largest triangle. Each triangular partition is known as an element S and the set of these finite triangular elements is called a mesh and we denote this mesh by T h;t . For the finite element mesh we require that no vertex of any element lie on the interior of a side of another element and Si ∩Sj ¼ ∅ if i≠j. Thus Ω h;t ¼

⋃ SðtÞ: SðtÞ∈T h;t

We assume that the triangulation T h;t for all t∈I is shape regular. Definition 3.1 (Shape regularity (Johnson [37])). Let hS ¼ be the longest side of S (otherwise known as the diameter of S), ρS the diameter of the circle inscribed in S and h ¼ max hS : S∈T h;t

A triangulation T h;t is shape regular if there exist a positive constant ζ independent of h such that ρS ≥ζ hS

∀S∈T h;t ; t∈I:

ð3:13Þ

The constant ζ is a measure of the minimum angle in any S∈T h;t , t∈I. The regularity condition (3.13) specifies that the elements S are not allowed to be arbitrary thin. We note that all angles of the triangular elements stay away from 1801 and are thus shape-consistent. Next, we discretize the problem (3.11) using the classical Galerkin method. The discretized problem of the coupled system(PS) of (3.11) reads: find ah ðx; tÞ,

A. Madzvamuse, U. Zenas George / Finite Elements in Analysis and Design 74 (2013) 76–92

uh ðx; tÞ, vh ðx; tÞ∈X h ðtÞ, t∈I such that   Z ∂w1h ∂u_ ∂v_ ∂u ∂v D11 h þ D12 h þ C11 h þ C12 h ∂x ∂y ∂x ∂y Ωh;t ∂x   1 ∂wh ∂u_ h ∂v_ h D33 þ þ ∂y ∂y ∂x   Z 1 ∂wh ∂w1h ∂uh ∂vh þ C33 þ dΩh;t ¼ − dΩh;t f1 ∂y ∂y ∂x ∂x Ωh;t Z þ n1 f 1 w1h ds;

have the following property; DφSi   ¼ 0: Dt S

ð3:14aÞ



Ωh;t

∂y

D D φ ðxðtÞÞ ¼ φ ðλS ðxðtÞÞ ¼ 0 Dt i Dt i because λS ¼ ðλ1 ; λ2 ; λ3 Þ satisfies 0 ≤ λk ≤ 1; k ¼ 1; 2; 3 for a point ðx; y; tÞ in S. □



∂u_ h ∂v_ ∂u ∂v þ D22 h þ C12 h þ C22 h ∂x ∂y ∂x ∂y   2 ∂wh ∂u_ h ∂v_ h D33 þ þ ∂x ∂y ∂x    Z ∂w2h ∂w2h ∂uh ∂vh C33 þ dΩh;t ¼ − dΩh;t f2 þ ∂x ∂y ∂x ∂y Ωh;t Z þ n2 f 2 w2h ds; ∂w2h

ð3:17Þ

Proof. The basis functions φSi ðx; tÞj ¼ 1; 2; 3 on S satisfy φSi ðxðtÞÞ ¼ φ i ðλS ðxðtÞÞÞ, where λS is the barycentric coordinates for a point ðx; y; tÞ in S. And the material derivative of φSi ðxðtÞÞ on an element S yields

∂Ωh;t

Z

81

D12

3.6. Semi-discrete model of the force balance equation

ð3:14bÞ

∂Ωh;t

and Z Z d ah w3h dΩh;t þ ðDa ∇ah  ∇w3h þ ka ah w3h Þ dΩh;t dt Ωh;t Ωh;t Z Z Dw3h dΩh;t ; ¼ ka ac w3h dΩh;t þ a Dt Ωh;t Ωh;t

"

A ðtÞ

We seek to find the finite element numerical approximation ah ðx; tÞ, uh ðx; tÞ, vh ðx; tÞ∈X h ðtÞ⊂X ðΩt Þ expressed as linear combinations of the linear nodal basis functions φi of the form

i¼1 nde

vh ðx; tÞ ¼ ∑ Vi ðtÞφi ðx; tÞ:

nde

uh ðx; tÞ ¼ ∑ Ui ðtÞφi ðx; tÞ;

A12 ðtÞ 22

T

A ðtÞ

#( dU

ðtÞ

dt dV ðtÞ dt

"

) þ

B11 ðtÞ 12

B ðtÞ

T

B12 ðtÞ 22

B ðtÞ

#(

UðtÞ VðtÞ

(

) ¼

F1 ðtÞ

)

2

F ðtÞ

ð3:18Þ ð3:14cÞ

represent piecewise linear finite element nodal basis functions satisfying that ( 1 if i ¼ j; φi ðxj ; tÞ ¼ 0 if i≠j:

ah ðx; tÞ ¼ ∑ ai ðtÞφi ðx; tÞ;

A11 ðtÞ 12

for all w1h ðx; tÞ, w2h ðx; tÞ, w3h ðx; tÞ∈X h ðtÞ. Let nde represent the total number of degrees of freedom of the nodes for the finite element discretization. Also let the set n o φi ; φi ∈X h ðtÞ⊂H 1 ðΩt Þ; i ¼ 1; …; nde ð3:15Þ

nde

In order to write the semi-discrete finite element approximation of the force balance equation in block matrices form, we substitute w1h , w2h and w3h by the basis function φj , ðj ¼ 1; …; ndeÞ, uh and vh by their corresponding values as given in (3.16) into (3.14a) and (3.14b), then the resulting semi-discrete system of algebraic equations can be written in compact-matrix-vector form as

and

i¼1

ð3:16Þ

i¼1

We take the number of nodes nde in Ωt to be large enough, typically in the orders of tens of thousands such that ah, uh and vh are a good approximation of a, u and v respectively. Transformation of local basis functions to a reference simplex. We denote the local basis functions on a reference element S by φ i with i ¼ 1; 2; 3. The local basis function φi ; i ¼ 1; 2; 3 on an element S∈T h;t can be obtained by φi ðxðtÞÞ ¼ φ i ðλS ðxðtÞÞÞ; where λS is the barycentric coordinates for a point ðx; y; tÞ in an element S. Lemma 3.1 (Transport property of basis functions (Dziuk and Elliott [36])). The finite element space on the discretized domain is a space of continuous piecewise linear functions whose nodal basis functions

where

fUðtÞ VðtÞg

are solutions of the finite element semi-discrete (3.18). In

(3.18) we denote by UðtÞ ¼ ðU1 ðtÞ; …; Unde ðtÞÞ, VðtÞ ¼ ðV1 ðtÞ; …; Vnde ðtÞÞ, dU=dtðtÞ ¼ ðdU1 =dtðtÞ; …; dUnde =dtðtÞÞ and dV=dtðtÞ ¼ ðdV1 =dtðtÞ; …; dVnde =dtðtÞÞ vectors of the solutions and their derivatives. Akl ðtÞ, Bkl ðtÞ, ðk; l ¼ 1; 2Þ and FðtÞ≔ðF1 ðtÞ; F2 ðtÞÞT are, respectively, the time dependent stiffness matrices and the generalized force vector. Their respective entries are defined by  Z  ∂φ ðxðtÞÞ ∂φj ðxðtÞÞ ∂φ ðxðtÞÞ ∂φj ðxðtÞÞ A11 þ D33 i dΩh;t ; D11 i ij ðtÞ≔ ∂x ∂y ∂x ∂y Ωh;t  Z  ∂φ ðxðtÞÞ ∂φj ðxðtÞÞ ∂φ ðxðtÞÞ ∂φj ðxðtÞÞ þ D22 i dΩh;t ; D33 i A22 ij ðtÞ≔ ∂x ∂y ∂x ∂y Ωh;t  Z  ∂φ ðxðtÞÞ ∂φj ðxðtÞÞ ∂φ ðxðtÞÞ ∂φj ðxðtÞÞ þ D33 i dΩh;t ; D12 i A12 ij ðtÞ≔ ∂x ∂y ∂y ∂x Ωh;t with:

  ∂φ ðxðtÞÞ ∂φj ðxðtÞÞ ∂φ ðxðtÞÞ ∂φj ðxðtÞÞ þ C66 i dΩh;t ; C11 i ∂x ∂y ∂x ∂y Ωh;t  Z  ∂φ ðxðtÞÞ ∂φj ðxðtÞÞ ∂φ ðxðtÞÞ ∂φj ðxðtÞÞ þ C22 i dΩh;t ; C33 i B22 ij ðtÞ≔ ∂x ∂y ∂x ∂y Ωh;t  Z  ∂φ ðxðtÞÞ ∂φj ðxðtÞÞ ∂φ ðxðtÞÞ ∂φj ðxðtÞÞ þ C33 i dΩh;t ; C12 i B12 ij ðtÞ≔ ∂x ∂y ∂y ∂x Ωh;t

Z B11 ij ðtÞ≔

and

Z ∂φj ðxðtÞÞ dΩh;t þ n1 f 1 ðxðtÞÞφj ðxðtÞÞ ds; ∂x Ωh;t ∂Ωh;t Z Z ∂φj ðxðtÞÞ dΩh;t þ f 2 ðxðtÞÞ n2 f 2 ðxðtÞÞφj ðxðtÞÞ ds: F2j ðtÞ ¼ − ∂y Ωh;t ∂Ωh;t F1j ðtÞ≔−

Z

f 1 ðxðtÞÞ

Appendix A illustrates how the normals, n1 and n2, to the cell boundary are computed and validated. For convenience's sake, we denote " 11 # " 11 # A ðtÞ A12 ðtÞ B ðtÞ B12 ðtÞ ½AðtÞ≔ 12 T ; ; ½BðtÞ≔ B12 ðtÞT B22 ðtÞ A22 ðtÞ A ðtÞ ( 1 ) ( )

( dU ðtÞ ) UðtÞ



F ðtÞ dW dt ðtÞ ≔ dV and FðtÞ ≔ WðtÞ ≔ : ; VðtÞ dt ðtÞ F2 ðtÞ dt

82

A. Madzvamuse, U. Zenas George / Finite Elements in Analysis and Design 74 (2013) 76–92

Then we can write the semi-discrete finite element model for the force balance equation (3.18) in the following compact form:



dW ðtÞ þ ½BðtÞ WðtÞ ¼ FðtÞ: ð3:19Þ ½AðtÞ dt Remark 3.2. Some entries of the block matrices ½A and ½B are non-symmetric matrices, hence ½A and ½B are non-symmetric block matrices. 3.7. Semi-discrete model of the reaction–diffusion equation To proceed to write the semi-discrete finite element approximation of the reaction–diffusion equation in matrix-vector form, we substitute wh by the basis function φj , ðj ¼ 1; …; ndeÞ and ah by its corresponding value as given in (3.16) into (3.14c). Upon applying the transport property of basis functions (see Lemma 3.1 for details) we obtain the following semi-discrete system of algebraic equations written in compact-matrix-vector form: d ðMðtÞfaðtÞgÞ þ ðDa KðtÞ þ ka MðtÞÞfaðtÞg ¼ ka ac HðtÞ; dt

ð3:20Þ

where fag ¼ ða1 ; …; ande Þ are actin concentration solutions of the semi-discrete (3.20). M, K and H are the time dependent mass and stiffness matrices, and force vector, whose entries are defined respectively as Z Mij ðtÞ≔ φi ðxðtÞÞφj ðxðtÞÞ dΩh;t ; Z Kij ðtÞ≔ Z Hj ðtÞ≔

Ωh;t

Ωh;t

∇φi ðxðtÞÞ  ∇φj ðxðtÞÞ dΩh;t

3.9. Computation of the evolution of the domain In order to track the evolution of the domain accurately, we specify a Lagrangian kinematic description of the domain. We assume that the displacement of the boundary of the domain occurs only in the direction that is normal to its non-deformed position. The displacement of the domain boundary at any time t corresponds to the solutions of the force balance equation Wðx; tÞ ¼ ðUðx; tÞ; Vðx; tÞÞ at the boundary which in turn are dependent on the local concentration of the actin filament aðx; tÞ at any particular space and time. The displacement of the interior nodes _ of the mesh is chosen to be equal to the flow velocity β (i.e. β ¼ x). On the boundary we assume that β ¼ ωn ≔∂W=∂t, where W is the displacement solution of the force balance equation. A limitation of the Lagrangian kinematic description is that the minimum angle condition for the mesh elements can be violated for non-constant mesh velocity and large mesh velocity. A remedy is to introduce a monitor function to check the quality of the mesh (say after a couple of mesh movements) and re-mesh when the elements are no longer favourable. It must be noted that remeshing will introduce interpolation errors. In our case, the displacement solutions of the force balance equation are small in value hence the mesh velocity is small. For the time interval we considered, no mesh refinement was required. Developing adaptive time- and mesh-refinement strategies form part of our current studies for these types of models. 3.10. A schematic numerical algorithm The fully discrete coupled problem is solved iteratively. We present in Algorithm 3.1 the numerical algorithm for the method.

and

Algorithm 3.1. Fully discrete scheme. Ωh;t

φj ðxðtÞÞ dΩh;t :

Initialize W0 and a0 . Set the time-step size to τ. FOR n ¼ 1; …; T F with time-step size τ DO Assemble Mðt n Þ, Hðt n Þ, ½Aðt n Þ, ½Bðt n Þ and Fðt n Þ.

3.8. Fully discrete scheme of the coupled problem

Solve for fWgnþ1 in (3.21a).

We discretize the time interval ð0; T f  into a finite number of uniform subintervals ½t n ; t nþ1  such that 0 ¼ t 0 o t 1 …o t n o t nþ1 ≤T f . The size of each time interval τ≔t nþ1 −t n . In order to obtain a fully discrete finite element model of the coupled problem, we consider a modified implicit finite differentiation formula for the time integration of the semi-discrete finite element model of the force balance equation and the reaction– diffusion equation as given in (3.19) and (3.20) respectively. Thus the fully discrete finite element model of the coupled problem is given by ð½Aðt n Þ þ τ½Bðt n ÞÞfWgnþ1 ¼ ½Aðt n ÞfWgn þ τFðtn Þ;

ð3:21aÞ

and ½ð1 þ τ ka ÞMðt nþ1 Þ þ τðDa Kðt nþ1 ÞÞfagnþ1 ¼ Mðt n Þfagn þ τ ka ac Hðt n Þ;

ð3:21bÞ n

nþ1

n

nþ1

and fag , fag are the displacement solutions where fWg , fWg and the actin concentrations at time tn and t nþ1 respectively. For the numerical implementation, the initial data fWg0 and fag0 are interpolated on the initial mesh Ωh;0 . All integral are evaluated using Gauss numerical quadrature [16]. At each time t, t∈I, we assemble the finite elements to obtain the system of linear algebraic (3.21). This system (3.21a) is solved using a generalized minimal residual method (GMRES) [38] while that of (3.21b) is solved using a conjugate gradient method (CG) with diagonal preconditioner [39,40].

Compute the new domain given fWgnþ1 Assemble Mðt nþ1 Þ and Kðt nþ1 Þ on the new domain. Solve for fagnþ1 in (3.21b). END FOR

4. Numerical results Let t ¼ nτ, where τ and n denotes the time-step size and number of time-steps respectively. The initial domain we consider is a unit disk Ω0 . We assume that at time 0 ot s ≪1 the domain evolution is negligible and based on this assumption we are able to compare the finite element solutions to predictions from linear stability theory. The perturbations used in the initial conditions will be dependent on the eigenmode we seek to excite. We subsequently show that the finite element scheme gives numerical results that are consistent with those predicted by linear stability theory. In all our simulations rnd denotes a randomly generated number between 10−3 and 10−5 and a tilde will be used to denote non-dimensional parameter values. 4.1. Validating numerical results using linear stability theory close to bifurcation points 4.1.1. Excitation of the eigenmode w1;1 We present the numerical results for actin concentration ah and the displacement solutions uh ¼ ðuh ; vh Þ in Fig. 1(b) and (c)

A. Madzvamuse, U. Zenas George / Finite Elements in Analysis and Design 74 (2013) 76–92

83

w1.1 0.5 0 –0.5 1

1

0 –1 –1

0

Fig. 1. Surface plots of the numerical results for actin concentration ah and displacement solution uh . Parameter values used in the numerical simulations are selected such 2 that the lowest non-zero wavenumber k1;1 is excited. (a) Predicted solution for w1;1 from linear theory. (b) Numerical result for actin concentration ah, and (c) Numerical result for displacement solution uh .

respectively. These results are consistent with those predicted by linear stability theory as given in Fig. B1(a). The initial conditions used were 1:0009 þ rndn sin ðxÞ, a finite element mesh with 2113 nodes was used and τ ¼ 1:0228  10−2 .

4.1.2. Excitation of the eigenmode: −w0;2 Let us fix p~ ¼ 1:646 and ψ~ ¼ 38:24 respectively. For these parameter choices, the dispersion relation again isolates the low2 est non-zero wavenumber k1;1 . The initial conditions of the actin concentration are chosen to be equal to ac þ rndn cos ðxÞ. We found that by choosing the perturbations in the initial conditions to be rndn cos ðxÞ, we encouraged the excitation of the eigenmode −w0;2 , where w0;2 is the lowest eigenmode. The numerical simulation results are shown in Fig. 2 for the actin concentrations ah and the displacement solutions uh . We note that all numerical results for uh are plotted using its magnitude. For the sake of the comparison of the eigenmode w0;2 to the numerical solutions we present in Fig. 2(a) a plot of the eigenmode w0;2 . In Fig. 2(c), we observe that the displacement solutions uh has its highest value around the points where x2 þ y2 ¼ 0:8. These points correspond to the points where δðlÞ is discontinuous. It is clear that the method is able to provide accurate solutions as predicted by linear theory close to the bifurcation points.

Fig. 2. Surface plots of the numerical results for the actin concentration ah and the displacement solution uh showing the replication of the eigenmode −w0;2 . (a) Surface plot of the eigenmode w0;2 . (b) and (c) Display the simulation results of actin concentration ah and displacement solutions uh respectively at time t ¼ 3:61  10−3 . (d) Displays the displacement solutions uh at time t¼ 0.024 having no oscillations at the points of discontinuity. The initial conditions used for actin concentrations are ac þ rndn cos ðxÞ.

4.2. Validating numerical results using linear stability theory far from bifurcation points 2

Far from the bifurcation points of k1;1 , we have the possibility of exciting mixed or higher modes. Below we present the numerical results for actin concentration ah and displacement solutions uh showing the excitation of mixed and higher modes.

4.2.1. Excitation of mixed and higher modes The excitation of mixed modes and higher modes is possible if the value of ψ~ is allowed to vary while all other parameters value are fixed (see Table B4 for possible values of ψ~ and the corresponding number of modes the dispersion relation isolates). We present in Figs. 3 and 4 the simulation results for actin concentration ah and the displacement solutions showing the excitation of mixed and higher modes respectively. We present in Table 2 the value of ψ used in the numerical simulations. A comparison of the predicted solutions from linear theory (Fig. 3 (a) and (d)) with the numerical results (Fig. 3(b), (c) and (e), (f)) show that the numerical scheme reproduces results in close agreement to those predicted by linear stability theory for the mixed and higher modes for the actin concentration.

4.3. Numerical results illustrating cell expansion, contraction and non-uniform deformation Now that we have validated our numerical results by use of linear stability theory close to bifurcation points, we are in position to illustrate the robustness and applicability of the numerical method. By choosing appropriate model parameter values far away from bifurcation points, cell expansion and contraction (uniform) can be observed as well as non-uniform cell deformation as illustrated in the next. The simulation results showing the actin concentration solutions with uniform expansions are shown in Fig. 5 for ~ ¼ ð70:366; 0:433Þ. In Fig. 5(e) we present a plot of the area ðψ~ ; pÞ of the cell against the number of time-steps taken to show the rate at which the area is increasing. For the sake of numerical experiment, we allow p to take a negative value. By choosing p~ ¼ −0:433 and ψ~ ¼ 70:366 we observe the cell contracting uniformly at positions that are equidistant from the centroid of the cell domain. And the concentration of actin increases from the boundary to the centroid. We show in Fig. 6 the actin concentration solutions at a time (a)

84

A. Madzvamuse, U. Zenas George / Finite Elements in Analysis and Design 74 (2013) 76–92

Fig. 3. Surface plots of the numerical results of the actin concentration ah and displacement solutions uh having mixed modes excited together with the reproduced solutions predicted from linear theory. (a) and (d) Predicted solutions from linear theory for w0;2 þ 2w2;1 and w3;1 −w0;2 −3w2;1 respectively. The numerical results reproduce the following modes: (b) and (c) w0;2 þ 2w2;1 , (e) and (f) w3;1 −w0;2 −3w2;1 . ψ is chosen as given in Table 2 such that mixed modes are excited (refer to the list given in Table B4 for more detail on the eigenmodes isolated by a dispersion relation depending on the value of ψ~ and ψ).

Fig. 4. Surface plots of the numerical results of the actin concentration ah and displacement solutions uh having higher modes excited together with the reproduced solutions predicted from linear theory. (a) and (d) Predicted solutions from linear theory for w4;1 and w2;1 respectively. The numerical results reproduce the following modes: (b) and (c) w4;1 , and (e) and (f) w2;1 . ψ is chosen as given in Table 2 such that higher modes are excited (refer to the list given in Table B4 for more detail on the eigenmodes isolated by a dispersion relation depending on the value of ψ~ and ψ).

t ¼ 1:2046  10−3 , (b) t ¼ 3:6137  10−3 , (c) t¼ 0.0361 and (d) t¼0.5180. The initial conditions of actin concentration are still equal to ac þ rndn cos ðxÞ. Non-uniform cell deformations occur if all parameter values in Table B1 stay unchanged but ψ~ is increased or by choosing

parameter values ψ~ and p~ far away from bifurcation points (refer ~ Below we give the to Fig. B2 for a parameter space plot of ðψ~ ; pÞ). parameter values and the simulation results for both cases were non-uniform cell deformations occur. Here the initial conditions of actin concentration are given by ac þ rndn cos ðxÞ.

A. Madzvamuse, U. Zenas George / Finite Elements in Analysis and Design 74 (2013) 76–92

85

Table 2 Values of ψ and the initial conditions used for the excitation of mixed and higher modes presented in Fig. 3. Figure

Value of ψ (dyn/cm2)

Initial conditions for the actin concentration

No. of nodes

Time-step size τ

Fig. 3(b) and (c) Fig. 3(e) and (f)

3.1  102 2.0  102

1:0 þ j0:1 cos xj

8321 8321

1.0228  10−3 1.0228  10−3

Fig. 4(b) and (c) Fig. 4(e) and (f)

3.67  102 1.1  102

8321 2113

1.0228  10−3 1.0228  10−2

1:0 þ j0:1 sin 2 yj 1:0 þ j0:1 sin x sin yj 1:0 þ j∑20 i ¼ 1 0:1 sin x sin yj

Fig. 5. Graphical display of the simulation results of the actin concentration. Blue denotes the lowest values and red the highest. These results were obtained at the following times: (a) t ¼ 1:2046  10−3 , (b) t ¼ 3:6137  10−3 , (c) t ¼0.0361, and (d) t¼ 0.5180. The numerical value of the contractile tonicity ψ~ ¼ 70:366 and that of the pressure coefficient p~ ¼ 0:433. (e) Plot of the area of the cell against the number of time-steps taken. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

 If all model parameter values are kept fixed as given in Table





B1, taking increasing values of ψ~ results in non-uniform cell deformations. For example taking ψ~ ¼ 1:04  105 , we observed from the numerical results that the cell deforms non-uniformly and rapidly making it impossible to capture the solutions and the numerical algorithm fails. If we choose ψ~ and p~ arbitrarily such that the parameter space ~ is very far from bifurcation points (for example, ψ~ ¼ 7:8  ðψ~ ; pÞ 103 and p~ ¼ 0:026), we observe non-uniform cell contraction with actin localized around the centroid of the cell (see Fig. 7). In Fig. 7(d) we present a plot of the cell area. A plot of the index of polarity against the number of time-step taken is given in Fig. 7(e). We recall that the simulation results presented in Fig. 7(e) was obtained for p~ ¼ 0:026 and ψ~ ¼ 7:8  103 . Now if we decrease only ψ~ such that ψ~ ¼ 1:3  103 , we observe cell deformations which agree qualitatively to those observed experimentally [41]. It can be observed that actin concentration is high at the periphery of the cell boundary and is highest at regions were protrusions occur. In Fig. 8, we present the graphical display of

the simulation results of the actin concentration. A plot of the index of polarity against the number of time-steps taken is given in Fig. 8(g). We can infer from these results that the dynamics of the cell domain, the distribution of actin filaments are related to the pressure coefficient p and the contractile tonicity ψ. 4.3.1. Comparison of Fig. 8(f) with experimental observations The distribution of F-actin in cells has been well studied and epifluorescence images are available in the literature. Fig. 9(a) shows an epifluorescence image of a Swiss 3T3 fibroblast cell obtained from [41]. The fibroblast was stained by the authors for F-actin after being allowed to spread for 1 h. The epifluorescence image shows that F-actin is predominant at the cell periphery and is highest at the regions were protrusions occur. The simulation results presented in Fig. 9(b) are thus consistent with the experimental observations. Remark 4.1 (Key model parameters). Based on the linear stability analysis and the numerical simulation results we are able to identify key parameters that control cell deformations with respect to the model problem. The contractile tonicity ψ is the bifurcation parameter

86

A. Madzvamuse, U. Zenas George / Finite Elements in Analysis and Design 74 (2013) 76–92

Fig. 6. Graphical display of the simulation results of the actin concentration. Blue denotes the lowest values and red the highest. These results were obtained at the following times: (a) t ¼ 1:2046  10−3 , (b) t ¼ 3:6137  10−3 , (c) t ¼0.0361, and (d) t¼ 0.5180. The numerical value of the contractile tonicity ψ~ ¼ 70:366 and that of the pressure coefficient p~ ¼ −0:433. (e) Plot of the area of the cell against the number of time-steps taken. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)

Fig. 7. Graphical display of the simulation results of the actin concentration. Blue denotes the lowest values and red the highest. These results were obtained at the following times: (a) t ¼ 1:2046  10−3 , (b) t ¼ 3:6137  10−3 , and (c) t¼ 0.0722. The numerical value of the contractile tonicity ψ~ ¼ 7:8  103 and that of the pressure coefficient p~ ¼ 0:026. (d) Plot of the area of the cell against the number of time-steps taken. (e) Plot of the index of polarity against the number of time-steps for simulation results with contractile tonicity ψ~ ¼ 7:8  103 and the pressure coefficient p~ ¼ 0:026. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)

A. Madzvamuse, U. Zenas George / Finite Elements in Analysis and Design 74 (2013) 76–92

87

1

Index of polarity

0.95 0.9 0.85 0.8 0.75 0.7 0.65 0

50 100 150 Number of time steps

200

Fig. 8. Graphical display of the simulation results of the actin concentration. Blue denotes the lowest values and red the highest. The numerical value of the contractile tonicity ψ~ ¼ 1:3  103 and the pressure coefficient p~ ¼ 0:026. (g) Plot of the index of polarity against the number of time-steps for simulation results with contractile tonicity ψ~ ¼ 1:3  103 and the pressure coefficient p~ ¼ 0:026. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)

Fig. 9. (a) Epifluorescence image of a Swiss 3T3 fibroblast cell [41]. F-actin is predominant at the cell periphery and is highest at the region were protrusions occur. (b) A simulated cell deformation for contractile tonicity ψ~ ¼ 1:3  103 and the pressure coefficient p~ ¼ 0:026. The distribution of F-actin in the simulated cell is in qualitative agreement with experimental observations. Reproduced with permission from Senju and Miyata [41], Oxford University Press.

and determines the transition from stable to unstable state. The actin saturation concentration asat and the pressure coefficient were also found to play a key role in cell deformations.

4.4. Applications to cell dynamics Here we choose parameters for the cell such that they are consistent with those available in the literature such that Young's

modulus E of the actin filament is 4000 dyn/cm2 [42], diffusion coefficient Da of the actin filament is 1:6  10−10 cm2 =s [43,44], the polymerisation rate ka is 66=s [45], the values of the poisson ratio ν, shear and bulk viscosities μ1 and μ2 respectively are as given in Table 1. A cell of radius 10 μm (0.001 cm) is assumed. We also assume that the pressure coefficient is 800 dyn/cm2 and the contractile tonicity is 6:9  104 dyn=cm2 . We note that for these choice of dimensional parameter values the equivalent nondimensional parameter values would have a dispersion relation

88

A. Madzvamuse, U. Zenas George / Finite Elements in Analysis and Design 74 (2013) 76–92

Fig. 10. (a)–(c) are graphical displays of the numerical results of the actin concentration ah with τ ¼ 1:0228  10−3 . Blue signifies the lowest values and red the highest values. (d) Plot of the area of the cell against the number of time-steps taken. It shows the area of the cell increasing with the number of time-step taken. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.) 2

that isolates more than 7 wavenumbers inclusive of k4;1 . We used an initial perturbation for actin equal to jrnd sin x sin yj in order to bias the excitation of the mode w4;1 . The numerical simulation result gives protrusion on four fronts and is similar to the result given in Fig. 10.

5. Conclusion and discussion Understanding cell movement and deformation in cell motility, developmental biology and biomedical engineering promises to unearth new research strands currently unchartered. In many cases, the models derived are highly complex and nonlinear with no closed form analytical solutions. Furthermore, growth evolution (cell membrane movement and deformation) is an added parameter that makes it even more difficult to treat the models analytically. It remains therefore, to propose and develop novel numerical methods capable of providing approximate solutions. The numerical method should be robust, efficient, and should be able to deal easily with domain and surface evolution. In this article, we proposed the moving grid finite element method [1] to numerically compute solutions to a system of nonlinear partial differential equations describing cell movement and deformation. Our main contributions are:

 a novel application of the MGFEM to solve a cytomechanical model that allows cells to move and deform arbitrarily.

 by using linear stability theory close to bifurcation points, we 

validated numerical results where linear solutions are analytically derived. far away from bifurcation points, we generate and observe complex cell deformations typically observed in experiments. Our numerical results show cell expansion, cell contraction, cell translation and relocation and cell protrusions.

 the contractile tonicity formed by the association of actin filaments to the myosin II motor proteins and the pressure coefficient are identified as key model bifurcation parameters of the model.

The robustness, generality and applicability of our numerical methodology allows us to tackle similar models in other areas such as plant biology, biomedical engineering and cancer biology.

Acknowledgments This work (A.M.) is partly supported by the following grants: EPSRC (EP/H020349/1 and EP/J016780/1), the LMS Grant (R4P2), the British Council through its UK-US New Partnership Fund (PMI2) and the Royal Society Grant: (R4N9). U.G. would like to acknowledge financial support from Akwa Ibom State University Technology Nigeria in support of her DPhil studies.

Appendix A. Computing normals Here we present the technique for the computation of n on the polygonal unit circle corresponding to the boundary of Ωh;0 given in (3.18). The exact normals of a unit circle can be calculated analytically and are known [46]. As the number of boundary nodes increases (via global refinement of the mesh) the mesh height h decreases and we show in Table A1 that the numerically computed unit normals converges to the exact normal of a unit circle. Given a line L with end points P 1 ¼ ðx1 ; y1 Þ and P 2 ¼ ðx2 ; y2 Þ. Then a vector N ¼ ðN 1 ; N 2 Þ normal to the line L is [47] N1 ¼ y1 −y2 ;

and

N2 ¼ x2 −x1 :

A. Madzvamuse, U. Zenas George / Finite Elements in Analysis and Design 74 (2013) 76–92

Table A1 L2 errors of computed normal vectors to a circle. Mesh

No. of global refinement

Mesh size

∥nexct −ncomp ∥L2 ð∂ωh;0 Þ

1 2 3 4

0 1 2 3

1.414 0.707 0.353 0.177

0.8284 0.3045 0.1087 0.0385

89

deformation. To proceed, we first non-dimensionalize the system of equations in order to reduce the number of parameters making the mathematical analysis afterwards more amendable.

B.1. Nondimensionalization Let the length scale L be a typical radius of a cell. Substituting the following dimensionless quantities into (2.1a) and (2.1b), 8 2 asat a u ~ ~ ~ ~ ~ ~ > > > t ¼ tka ; a ¼ ac ; u ¼ L ; ∇ ¼ L∇; Δ ¼ L Δ; a sat ¼ ac ; < p~ ¼ p 1þν ϕ~ ¼ ϕ; ϵ~ ¼ ϵ; μ~ i ¼ μi ka 1þν ψ~ ¼ ψ a2c 1þν E ; E ; E ; > > β Da ~ > ~ ; d ¼ Da ¼ 2 ; :β ¼ ka L

ka L

ðB:1Þ results in 8  

> ∇  μ1 ϵt þ μ2 ϕt I þ ðϵ þ ν′ϕIÞ þ ψa2 e−a=asat I þ > > n o <   p 1 þ 2π δðlÞ arctanðaÞ I ¼ 0; ∇  1þϕ > > > : ∂a −dΔa þ ∇  ðaβÞ ¼ 1−a; ∂t

Fig. A1. (a) is a cross-section of a polygonal circle showing the outward pointing unit vector normal to a line segment. (b) shows the two outward pointing normal vectors that occur as a result of the intersection of two line segments at a nodal position. The resultant normal n is computed by applying the parallelogram law of forces.

ðB:2Þ

where ν′ ¼ ν=1−2ν, ν≠0:5. In the above, we have dropped the tildes without any loss of generality. Non-dimensionalized parameters are shown in Table B1 and these are determined from the dimensional values shown in Table 1.

B.2. Linear stability analysis on a unit disc The unit normal vector n ¼ ðn1 ; n2 Þ to L is then obtain by dividing the normal vector N by its magnitude such that y1 −y2 n1 ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðy1 −y2 Þ2 þ ðx2 −x1 Þ2

and

x2 −x1 n2 ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : ðy1 −y2 Þ2 þ ðx2 −x1 Þ2

Given a boundary of a closed surface with nodal points Pi ði ¼ 1; …; ndeÞ having coordinates ðxi ; yi Þ, where nde denotes the total number of nodes, then the outward pointing unit vector normal to any line segment ½i; i þ 1 (see for example Fig. A1(b)) is yi −yiþ1 ; n1 ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðyi −yiþ1 Þ2 þ ðxiþ1 −xi Þ2

and

xiþ1 −xi n2 ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : ðyi −yiþ1 Þ2 þ ðxiþ1 −xi Þ2

At the intersection of two line segments (i.e. at a nodal position). The resultant normal n is computed by applying the parallelogram law of forces (see Fig. A1). In Table A1 we show that the computed normals converge to the exact normal in the L2 norm as the number of mesh points increases. The mesh with no global refinement corresponds to a macro triangulation having 4 line segments that form a diamond shape. The mesh with 1, 2, 3 global refinements of the macro triangulation have polygonal boundaries that comprises of 8, 16 and 32 line segments. We note that we have incorporated nodal projection of boundary nodes during mesh refinements so that the resulting mesh approximates the unit circle better as the number of global refinement increases.

We restrict linear stability analysis to be carried out on a static unit disk. The assumption is that for time t ¼ t 0 þ nτ, t 4 t 0 with n small, the deformed cell domain is still very close to a unit disk. Hence linear analysis is valid on Ωt provided t≪1 and Ωt ≈Ω0 . Here τ is a small time step. The steady state solution of (B.2) is as ¼ 1, us ¼ 0. Consider the ^ and stability of the steady state to small perturbations a ¼ as þ a, ^ where a^ and u^ are small variations from the steady u ¼ us þ u, state. Substituting these into the nonlinear system (B.2) and neglecting all but the linear terms results in the following linear system of partial differential equations: (   ∇  ðμ1 ϵt þ μ2 ϕt IÞ þ ðϵ þ ν′ϕIÞ þ s′ð1ÞaI þ pð1−ϕÞI þ p 2π δðlÞaIÞ ¼ 0; ∂a ∂t −dΔa

þ ∇  ðβÞ þ a ¼ 0;

ðB:3Þ where s′ð1Þ ¼ ∂sðaÞ=∂aja ¼ as . Let λ and k represent the growth rate (also known as an eigenvalue) and the wave vector respectively and assume that an and un are constants. It follows that the linearized system (B.3) Table B1 Non-dimensional parameters and their values as used in the nondimensionalized mathematical model except where it is specified otherwise. These parameter values were obtained from their dimensional counterpart by using (B.1). Corresponding dimensional parameter values are given in Table 1. Parameter Meaning

Non-dimensional value

d asat μ1 μ2 ψ

0.4 1.4 2.5 6.5 62.4

Appendix B. Linear stability analysis of the model equations Here we present linear stability analysis of the cytomechanical model. Specific details of linear stability theory can be found in [31]. The primary purpose of linear stability theory in this article are two fold: (i) to present linear solutions that validate numerical solutions around bifurcation points and (ii) to determine mathematically plausible parameter values for the pressure coefficient p and the contractile coefficient ψ that will give rise to cell

p l0

Diffusion coefficient of actin Saturation concentration of F-actin Shear viscosity of the actin network Bulk viscosity of the actin network Contractility coefficient of the actin network Pressure coefficient of the actin network Specifies the vicinity of the membrane

0.26 0.8

90

A. Madzvamuse, U. Zenas George / Finite Elements in Analysis and Design 74 (2013) 76–92

provided λ satisfies

admits non-trivial solutions of the form aðx; tÞ ¼ an expðλ t þ ik  xÞ

and

ðB:4Þ

uðx; tÞ ¼ un expðλ t þ ik  xÞ;

2

2

2

k ½μλ2 þ bðk Þλ þ cðk Þ ¼ 0;

ðB:5Þ

where   2 2 2 bðk Þ ¼ μdk þ 1 þ ν′ þ μ−s′ð1Þ−p− pδðlÞ ; π

Table B2 Possibilities of stable and unstable modes to exist. 2

Possible conditions ðk ≠0Þ

Types of mode

and

Negative

cðk Þ ¼ d½1 þ ν′−pk þ ½1 þ ν′−p;

2

2

2

2

Stable modes for all k 4 0. Unstable modes will exist

Positive

2

2

Unstable modes will exist

Positive

2

2

Unstable modes will exist

Positive

bðk Þ 40 and cðk Þ 40 bðk Þ 40 and cðk Þ o 0 bðk Þ o 0 and cðk Þ4 0 bðk Þ o 0 and cðk Þo 0

2

Sign of ReðλÞ

Table B3 Zeros of the derivative of Bessel functions of the first kind: j′m;n m ¼0,…,5; n ¼1, …,4. n

j′0;n

j′1;n

j′2;n

j′3;n

j′4;n

j′5;n

1 2 3 4

0.00000 3.83170 7.01558 10.17346

1.84118 5.33144 8.53632 11.70600

3.05424 6.70613 9.96947 13.17037

4.20119 8.01524 11.34592 14.58585

5.31755 9.28240 12.68191 15.96411

6.41562 10.51986 13.98719 17.31284

Table B4 A display of the values of ψ~ (and their corresponding dimensional values of ψ ) required by the dispersion relation in order to isolate at least two unstable wavenumbers. Value of ψ~

Value of ψ (dyn/cm2)

Wavenumbers isolated (i.e. unstable)

62.439

72.045

95.327

109.993

k0;1 ; k1;1 ; k2;1

131.869

152.157

174.767

201.654

270.094

311.647

317.758

366.643

2 2 2 2 k0;1 ; k1;1 ; k2;1 ; k0;2 2 2 2 2 2 2 k0;1 ; k1;1 ; k2;1 ; k0;2 ; k3;1 ; k4;1 2 2 2 2 2 2 2 2 k0;1 ; k1;1 ; k2;1 ; k0;2 ; k3;1 ; k4;1 ; k1;2 ; k5;1 2 2 2 2 2 2 2 2 2 k0;1 ; k1;1 ; k2;1 ; k0;2 ; k3;1 ; k4;1 ; k1;2 ; k5;1 k2;2

2

2

2

2

k0;1 ; k1;1 2

2

2

2

where k2 is the solution of the eigenvalue problem Δw ¼ −k w, w≠0 with homogeneous Neumann boundary condition i.e. n^  ∇w ¼ 0, where w are time-independent eigenfunctions of the linear system of differential equation (B.3). Solving (B.5) results in the dispersion relation qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 2 2 −bðk Þ 7 b ðk Þ−4μk cðk Þ 2 λðk Þ ¼ : ðB:6Þ 2 2μk From (B.5), it is obvious that the dispersion relation is indetermi2 nate when k ¼ 0. Thus in the linear stability analysis we shall only 2 2 consider k 4 0. Instability will occur for some wavenumber k 4 0 2 if the corresponding λðk Þ has ReðλÞ 4 0. Below we use Routh– Hurwitz stability criterion [48] to determine a sufficient condition on the coefficients of the polynomial (B.5) that would result in 2 2 2 Re λðk Þ 4 0 for some k 4 0. Instability can occur for some k 4 0 if one of the following conditions is satisfied (see Table B2 for details): 2

2

1. bðk Þ 4 0 and cðk Þ o 0 or 2 2 2. bðk Þ o 0 and cðk Þ 4 0 or 2 2 3. bðk Þ o 0 and cðk Þ o 0. In order to determine the value of k2, we use the separation of variables method to solve for an eigenfunction w, w≠0 and a wavenumber k2 such that the eigenvalue problem is well defined. We present the result in Remark B.1 below.

2

Fig. B1. (a) A surface plot of the vibration mode w1;1 ðr; θÞ ¼ J 1 ðj′1;1 rÞ cos θ which corresponds to the lowest non-zero wavenumber k1;1 on a unit disk. (b)–(f) Surface plots of selected eigenmodes for the bands of wavenumbers displayed in Table B4. The xy-axis give the x- and y-coordinates of the unit disk and the z-axis gives the magnitude of the eigenmodes.

A. Madzvamuse, U. Zenas George / Finite Elements in Analysis and Design 74 (2013) 76–92

91

B.2.1. Parameter space ðψ; pÞ Here we keep all parameters given in Table B1 fixed except the pressure coefficient p~ and the contractile tonicity ψ~ . We make a ~ showing spaces that lie in a plot o f the parameter space ðψ~ ; pÞ region of Hopf instability ðReðλÞ o0 and ImagðλÞ 4 0Þ, oscillatory instability ðReðλÞ 4 0 and ImagðλÞ 4 0Þ or Turing instability ðReðλÞ 4 0 and ImagðλÞ ¼ 0Þ [52–54]. This plot is presented in 2 Fig. B2 for the case where bðk Þ has δðlÞ ¼ 1. Note that for δðlÞ ¼ 0 the parameter space plot is identical to that shown in Fig. B2. In Fig. B2 we observe that there exists Hopf, oscillatory and Turing instability regions where all have a dispersion relation that isolates 2 k1;1 as the only non-zero wavenumber. The ash area (H) signifies the region of Hopf instability with a dispersion relation that 2 isolates k1;1 as the onlyunstable non-zero wavenumber. The region (dark blue) preceding that of Hopf instability represents the region where the uniform steady state is always stable for all wavenumbers k2. Immediately after the Hopf instability region we have a 2 purple region (OS1) where an oscillatory instability exists for k1;1 . The green region is a Turing instability region and has a dispersion 2 relation that isolates k1;1 as the only unstable non-zero wavenumber. In the purple region (OS2) an oscillatory instability exists for 2 k2;1 . The light blue region is also a Turing instability region and has a dispersion relation that isolates the first two non-zero wave2 2 numbers k1;1 and k2;1 . Fig. B2. We present the parameter space plot showing the regions where instability exist for δðlÞ ¼ 1. The ash area (H) signifies the region of Hopf instability 2 with a dispersion relation that isolates k1;1 as the only unstable non-zero wavenumber. The region (dark blue) preceding that of Hopf instability represents the region where the uniform steady state is always stable for all wavenumbers k2. Immediately after the Hopf instability region we have a purple region (OS1) where 2 an oscillatory instability exists for k1;1 . The green region is a Turing instability 2 region and has a dispersion relation that isolates k1;1 as the only unstable non-zero 2 wavenumber. In the purple region (OS2) an oscillatory instability exists for k2;1 . The light blue region is also a Turing instability region and has a dispersion relation that isolates the first two non-zero wavenumbers. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)

Remark B.1. Consider a scalar function w such that Δw ¼ 2 k w; w≠0, on a disk Ω0 ¼ fðx; yÞ : x2 þ y2 ≤R2 g with homogeneous Neumann boundary condition. Then the discrete eigenvalues are given by [49,50] 2

km;n ¼ ðj′m;n =RÞ2 ;

m ¼ 0; 1; 2; ⋯; n ¼ 1; 2; ⋯;

ðB:7Þ

with corresponding eigenfunctions, wm;n given by ( wm;n ðr; θÞ ¼

J 0 ðj′0;n r=RÞ

if m ¼ 0;

J m ðj′m;n r=RÞc cos mðθ−θ1 Þ

if m 4 0:

n ¼ 1; 2; ⋯;

for constants c and θ1 . Here Jm is the mth Bessel function of the first kind and j′m;n is the nth positive zero of J′m (except for j′0;1 ¼ 0), where J′m denotes the derivative of the Bessel function Jm with respect to r. Values of j′m;n are obtained from the ‘Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables’ by Abramowitz and Stegun [51]. For the sake of completeness, in Table B3, we present some values of j′m;n for m ¼ 0; …; 5 and 2 n ¼ 1; …; 4. We consider the lowest non-zero wavenumber k1;1 ≔ 2 ðj′1;1 =RÞ and its corresponding vibration mode w1;1 ðr; θÞ ¼ J 1 ðj′1;1 r=RÞ cos θ. The wavenumber k2 are discrete. In Table B4, we give the values of ψ required by the dispersion relation to isolate a couple of wavenumbers. In Fig. B1 we show possible single and mixed modes that could evolve for bands of wavenumbers given in Table B4.

References [1] A. Madzvamuse, A.J. Wathern, P.K. Maini, A moving grid finite element method applied to a model biological pattern generator, J. Comput. Phys. 190 (2003) 478–500. [2] H. Chen, B.W. Bernstein, J.R. Bamburg, Regulating actin-filament dynamics in vivo, Trends Biochem. Sci. 25 (2000) 19–23. [3] R. Ananthakrishnan, J. Gusk, J. Käs, Cell mechanics: recent advances with a theoretical perspective, Recent Res. Dev. Biophys. 5 (2006) 39–69. [4] C. Le Clainche, M.-F. Carlier, Regulation of actin assembly associated with protrusion and adhesion in cell migration, Physiol. Rev. 88 (2008) 489–513. [5] A.J. Ridley, M.A. Schwartz, K. Burridge, R.A. Firtel, M.H. Ginberg, G. Borisy, J. T. Parsons, A.R. Horwitz, Cell migration: integrating signal from front to back, Science 302 (2003) 1704–1709. [6] F. Fleischer, R. Ananthakrishnan, S. Eckel, H. Schmidt, J. Käs, V. Svitkina, B. Michael, Actin network architecture and elasticity in lamellipodia of melanoma cells, New J. Phys. 9 (2007) 420. [7] D.A. Lauffenburgery, A.F. Horwitz, Cell migration: a physically integrated molecular process, Cell 84 (1996) 359–369. [8] F. Xue, D.M. Janzen, D.A. Knecht, Contribution of filopodia to cell migration: a mechanical link between protrusion and contraction, Int. J. Cell Biol. (2010). [9] A. Stéphanou, Spatio-Temporal Dynamics of the Cell: Characterization from Images and Computer Simulations, Lambert Academic Publishing, 2010. [10] R.A.F. Clark, The Molecular and Cellular Biology of Wound Repair, Plenum Press New York, 1996. [11] A. Stéphanou, E. Mylona, M. Chaplain, P. Tracqui, A computational model of cell migration coupling the growth of focal adhesions with oscillatory cell protrusions, J. Theor. Biol. 253 (2008) 701–716. [12] A. Stéphanou, M. Chaplain, P. Tracqui, A mathematical model for the dynamics of large membrane deformations of isolated fibroblasts, Bull. Math. Biol. 66 (2004) 1119–1154. [13] F. Binamé, G. Pawiak, P. Roux, U. Hibner, What makes cells move: requirements and obstacles for spontaneous cell motility, Mol. Biosyst. 6 (2010) 648–661. [14] G. Charras, E. Paluch, Blebs lead the way: how to migrate without lamellipodia, Nat. Rev. Mol. Cell Biol. 9 (2008) 730–736. [15] W. Alt, R.T. Tranquillo, Basic morphogenetic system modeling shape changes of migrating cells: how to explain fluctuating lamellipodial dynamics, J. Biol. Syst. 3 (1995) 905–916. [16] J.N. Reddy, An Introduction to the Finite Element method, McGraw-Hill, 1993. [17] A. Madzvamuse, Time-stepping schemes for moving grid finite elements applied to reaction–diffusion systems on fixed and growing domains, J. Comput. Phys. 214 (2006) 239–263. [18] E.J. Crampin, E.A. Gaffney, P.K. Maini, Reaction and diffusion on growing domains: scenarios for robust pattern formation, Bull. Math. Biol. 61 (1999) 1093–1120. [19] T. Liszka, J. Orkisz, The finite difference method at arbitrary irregular grids and its application in applied mechanics, Comput. Struct. 11 (1980) 83–95. [20] K.W. Morton, D.F. Mayers, Numerical Solution of Partial Differential Equations, Cambridge University Press, 1994.

92

A. Madzvamuse, U. Zenas George / Finite Elements in Analysis and Design 74 (2013) 76–92

[21] A. Madzvamuse, P.K. Maini, A.J. Wathen, A moving grid finite element method for the simulation of pattern generation by turing models on growing domains, J. Sci. Comput. 24 (2005) 247–262. [22] C.A. Brebbia, Boundary Element Methods, Springer, Berlin, Heiddberg, New York, 1981. [23] S.L. Crouch, A.H. Starfield, Boundary Element Methods in Solid Mechanics, George Allen & Unwin, London, 1983. [24] J.A. Sethian, Level Set Methods: Evolving Interfaces in Geometry, Fluid Mechanics, Computer Vision and Materials Sciences, Cambridge University Press, 1996. [25] O.C. Zienkiewicz, R.L. Taylor, J.Z. Zhu, The Finite Element Method: Its Basis and Fundamentals, Elsevier Butterworth-Heinemann, London, 2005. [26] T. LaForce, Pe281 Boundary Element Method Course Notes, 2006. [27] A. Madzvamuse, A Numerical Approach to the Study of Spatial Pattern Formation, DPhil Thesis, Exeter college, Trinity term, University of Oxford (2000). [28] A. Madzvamuse, P.K. Maini, Velocity-induced numerical solutions of reaction– diffusion systems on fixed and growing domains, J. Comput. Phys. 224 (2007) 100–119. [29] M.J. Baines, Moving Finite Elements, Monographs on Numerical Analysis, Clarendon Press, Oxford, 1994. [30] K. Miller, R.N. Miller, Moving finite elements, SIAM J. Numer. Anal. 18 (1981) 1019–1032. [31] U. George, A. Stéphanou, A. Madzvamuse, Mathematical modelling and numerical simulations of actin dynamics in the eukaryotic cell, J. Math. Biol. (2012) 1–47, http://dx.doi.org/10.1007/s00285-012-0521-1. [32] A.R. Bausch, F. Ziemann, A.A. Boulbitch, K. Jacobson, E. Sackmann, Local measurements of viscoelastic parameters of adherent cell surfaces by magnetic bead microrheometry, Biophys. J. 75 (1998) 2038–2049. [33] E.L. Barnhart, K.-C. Lee, K. Keren, A. Mogilner, J. Theriot, An adhesiondependent switch between mechanisms that determine motile cell shape, PLoS Biol. 9 (2011). [34] S. Larsson, V. Thomée, Partial Differential Equations with Numerical Methods, Springer-Verlag, Berlin, Heidelberg, New York, 2003. [35] D.J. Acheson, Elementary Fluid Dynamics. Oxford Applied Mathematics and Computing Science, Clarendon Press, Oxford, 1990. [36] G. Dziuk, C.M. Elliott, Finite elements on evolving surfaces, IMA J. Numer. Anal. 27 (2007) 262–292. [37] C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element Method, University of Cambridge, 1987. [38] Y. Sadd, M.H. Schultz, Gmres: a generalized minimal residual algorithm for solving non-symmetric linear systems, SIAM J. Sci. Stat. Comput. 7 (1986) 856–869.

[39] Y. Saad, Iterative Methods for Sparse Linear Systems, Society for Industrial and Applied Mathematics, 2003. [40] M.R. Hestenes, E. Stiefel, Methods of conjugate gradients for solving linear systems, J. Res. Natl. Bur. Stand. 49 (1952). [41] Y. Senju, H. Miyata, The role of actomyosin contractility in the formation and dynamics of actin bundles during fibroblast spreading, J. Biochem. 145 (2) (2009) 137–150. [42] M.J. Dayel, O. Akin, M. Landeryou, V. Risca, e.a. Mogilner, A, In silico reconstitution of actin-based symmetry breaking and motility, PLoS Biol. 7 (2009) e1000201. [43] F. Lanni, B.R. Ware, Detection and characterization of actin monomers, oligomers, and filaments in solution by measurement of fluorescence photobleaching recovery, Biophys. J. 46 (1984) 97–110. [44] J.R. Simon, A. Gough, E. Urbank, F. Wang, F. Lanni, Analysis of rhodamine and fluorescein-labeled f-actin diffusion in vitro by flourescence photobleaching recovery, Biophys. J. 54 (1988) 801–815. [45] N. Watanabe, Inside view of cell locomotion through single-molecule: fast f-/ g-actin cycle and g-actin regulation of polymer restoration, Proc. Jpn. Acad. 86 (2010). [46] D.G. Zill, M.R. Cullen, Advanced Engineering Mathematics, Jones and Bartlett Publishers International, Inc., 2000. [47] J. Burkardt, Computational Geometry Lab: Barycentric Coordinates in Triangles, 2009, /http://people.sc.fsu.edu/  jburkardt/presentations/cg_lab_bary centric_trianglesS (accessed October 2009). [48] L. Edelstein-Keshet, Mathematical Models in Biology, Society for Industrial and Applied Mathematics, 2005. [49] U. George, A Numerical Approach to Studying Cell Dynamics, Ph.D. Thesis, University of Sussex, 2011. [50] E.C. Zachmanoglou, D.W. Thoe, Introduction to Partial Differential Equations with Applications, Dover Publications, Inc., New York, 1986. [51] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, US Department of Commerce, 1968. [52] P. Moreo, E.A. Gaffney, J.M. Garca-Aznar, M. Doblaré, On the modelling of biological patterns with mechanochemical models: insights from analysis and computation, Bull. Math. Biol. 72 (2010) 400–431. [53] L. Yang, M. Dolnik, A.M. Zhabotinsky, I.R. Epstein, Pattern formation arising from interactions between turing and wave instabilities, J. Chem. Phys. 117 (2002) 7259–7265. [54] M.C. Cross, P.C. Hohenberg, Pattern formation outside of equilibrium, Rev. Mod. Phys. 65 (1993) 851–1112. [55] A. Stéphanou, P. Tracqui, Cytomechanics of cell deformations and migration: from model to experiment, C.R. Biol. 325 (2002) 1–14.