Integrating Actin and Myosin II in a Viscous Model for Cell Migration

This article presents a mathematical and computational model for cell migration that couples a system of reaction-advection-diffusion equations describing the bio-molecular interactions between F-actin and myosin II to a force balance equation describing the structural mechanics of the actin-myosin network. In eukaryotic cells, cell migration is largely powered by a system of actin and myosin dynamics. We formulate the model equations on a two-dimensional cellular migrating evolving domain to take into account the convective and dilution terms for the biochemical reaction-diffusion equations, with hypothetically proposed reaction-kinetics. We employ the evolving finite element method to compute approximate numerical solutions of the coupled biomechanical model in two dimensions. Numerical experiments exhibit cell polarization through symmetry breaking which are driven by the F-actin and myosin II. This conceptual hypothetical proof-of-concept framework set premises for studying experimentally-driven actin-myosin reaction-kinetic network interactions with generalizations to multi-dimensions.


INTRODUCTION
Cell migration is fundamental in many biological processes and plays a key role in wound healing, immune response, development of embryos, inflammation, cancer invasion among others [1][2][3][4][5][6]. According to experimental observations, cell migration in eukaryotic cells is powered by actinmyosin system [7]. The actin cytoskeleton and its corresponding motor proteins play crucial role in cell movement. Actin is a polymer that can exist in two forms: F-actin and G-actin forms [1,5,6]. It converts from inactive state (G-actin) to active state (F-actin) through a process called actin polymerization and conversely from active state to inactive state through actin depolymerization [1,6,8]. The regulatory proteins are responsible for these actin polymerization-depolymerization processes [1,8]. Actin polymerization promoting proteins such as nucleator proteins help in creating new actin filaments [8]. On the other hand, actin depolymerization factor such as cofilin is capable of binding to the actin filament causing it to disintegrate and form G-actin monomers [1,8]. Actin filaments have two distinct ends: a plus end called barbed end which is a fast growing end and a minus end which is a slow growing end [7,8]. Actin filaments can assemble structures forming networks and bundles through interaction with motor proteins [9,10]. These structures produce cell protrusions called lamellipodia [6,7,11]. Actin cytoskeleton is therefore the main structure that contributes actively to force generation and therefore drives cell movement [6,12].
Actin filament is responsible for force generations that drive cell migration through two major processes, namely: (i) rapid polymerization of actin network at the cell periphery through the growth of lamellipodia [5,6,[13][14][15]. This leads to expansion of the plasma membrane and thus to the development of a contact area with the substrate and (ii) the development of stress fibers and networks that are contractile due to the action of motor proteins that tend to slide actin filaments relative to each other [6,15]. The most abundant motor proteins is myosin II [9,16]. Myosin II is known to bind to actin filaments in the cortex, crosslinking and contracting the actin filaments, causing cortical tension and mechanical resistance, which drives the cells' overall behavior. These active forces from polymerization of actin and contraction of stress fibers are eventually transmitted to substrates through adhesion sites, hence providing the necessary forces required for cell propulsion [1,12]. Similarly, new actin polymerization occurs between the cortex and the membrane, giving rise to outward pressure-driven cell membrane pseudopod formation, while in other regions of the cell, cortical tension pulls along the rest of the cell body (through contractility). All these processes when combined appropriately lead to cell movement and deformation. Cells are however not limited to crawling as a means of migrating. It has been observed that cells can swim while freely suspended in a fluid [17][18][19]. Infact, neutrophils can migrate with little to no adhesion [20].
For many centuries, experimental biology has occupied the minds of many researchers in the quest to understand the complexity of cell motility. In recent decades, mathematical and computational modeling has rapidly become an essential research technique that has greatly contributed to the understanding of the subject of cell motility [4]. The advantage of computational modeling is that it can overcome intrinsic experimental limitations and therefore allow for virtual experiments which sometimes are impossible to carry out in reality. This allows experimentalists to probe and pause new experimental hypotheses. It is a fact that cell motility involves a large number of proteins that interact together in a complex way [21]. Proposing an accurate model to account for the vast molecular interactions involved in cell motility is therefore a non-trivial activity. It is also largely known that the interaction of actin with its associated proteins is usually a major factor in the derivation of models for cell motility [4]. Many models are built on the concept that cell motility is composed of the following stages: protrusion, adhesion and contraction [1]. The cell pushes out the front, then it assembles tight adhesions to the surface at the leading edge and weakens such adhesions at the rear and finally the cell develops contractions that pull the weakly adherent rear toward the strongly adhered front [4]. Hence, these involve biomechanical interactions which describe the dynamical behavior between intra-and extra-cellular processes to enable the cell to migrate.
There have been different strategies of modeling cell migration in the last several decades. The first modeling efforts were directed at quantifying actin treadmill and using thermodynamics to understand the nature and magnitude of the polymerization force [4]. These early works introduced fundamental ideas that are still used in developing complex models [4]. In [22], "Polymerization Brownian ratchet" model was proposed to describe actin polymerization as rigid mechanism which elongates by rectifying the Brownian motion of the membrane. According to this model, when the end of an actin filament comes into contact with a membrane, the membrane would diffuse away and therefore create a gap sufficient for monomers to be added. These ratchet models used differential-difference equations [4,22]. Later in [23], an improvement was made to this model to consider the filaments as elastic springs whose behavior is a function of the bending modulus of the filament and the angle it makes with a load at its tip. The thermal fluctuations of actin filaments displaces the actin filaments from the membrane and creates a gap for the elongation of the filament [23]. This model was able to predict an optimal angle between the actin filament and the load for the effective force transmission. In [24], ratchet models underwent further development when it was suggested that some of the actin filaments are attached to the surface they push. This model, the 'tethered ratchet' model, explains the mechanism by assuming that the filament attach to the surface transiently, dissociating fast and growing freely until getting capped and losing contact with the surface altogether [24]. In [25][26][27], actin polymerization and depolymerization were treated as stochastic processes.
In [28], a mathematical model that describes key details of actin dynamics in protrusion associated with cell motility was developed. This model was based on the dendritic-nucleation hypothesis for lamellipodial protrusion in non-muscle cells such as Keratocytes. An output of the model was a relationship between the protrusion velocity and the number of filament barbed ends pushing the membrane. They observed that this relationship has a local maximum: too many barbed ends deplete the available monomer pool and too few are insufficient to generate protrusive force. Their result suggested that to achieve rapid motility, some tuning of parameters affecting actin dynamics must be operating in the cell [28].
Continuum models have also been developed to study cell motility. These models can be classified as moving boundary problems since a motile cell is a typical example of such a class of partial differential equations. In many cases, the models describing the biochemical interactions between the molecular species are coupled to a momentum balance law that accounts for the forces driving the migration of the cell. Examples of moving boundary problems for cell migration include (but are not limited to): a two-phase fluid model for cytosol and the actin network in [29], a one-dimensional viscoelastic model of the cytoplasm and active stress generation in [30], a one-dimensional model for the actin distribution and its reaction in [31], a two-dimensional elastic continuum model in [32] and a cytomechanical model that couples a force balance mechanical equation for actin network to a reaction-diffusion equation for actin [33]. This cytomechanical model was later extended in [34] where they used a cartesian coordinate system.
Another strategy is to use phase-field models as described in these works [35][36][37][38][39][40][41][42][43]. In [44], a Keratocyte cell was modeled as a two dimensional sheet with fixed area. The shape of the cell membrane was determined by the interactions of various forces including the surface tension, the bending forces and the pressure that constraints the cell area. In [35,36], a phasefield model was used to simulate vesicles' deformation and tumbling while in [37,38], the phase-field was used to study the three dimension deformation of a vesicle membrane using an energetic framework. An extension of this model to multicomponent vesicles was studied in [39]. In [41], a phasefield model for the Keratocyte cell that couples actin flow and adhesion mechanism during cell migration was presented. In the model, both myosin II contraction and actin polymerization were treated as active stresses. The adhesion sites could switch between the gripping mode and the slipping mode and their dynamics were integrated with actin flow. They also included tension and bending forces at the membrane. The various forces involved in cell migration are then considered and translated into velocity which then evolve the cell shape. Most recently, Moure and Gomez [43] gave a comprehensive review of phase-field models of individual and collective cell migration. They showed individual cell migration in confined and fibrous environments that highlight the mechanochemical interplay between the cell and the extracellular environment. Phase-fields have been adopted as one plausible numerical approach for handling moving boundary problems by embedding the sharp moving interface into a higher space-dimension which is fixed and then defining a phase-field function that describes the location of the cell (typically, 1 inside the cell) and (0 outside). The advantages of employing phase-fields include (i) the reformulation of a moving boundary problem onto a stationary domain and (ii) the method allows for topological changes in cell shape (e.g., cell division). On the other hand, phase-fields have disadvantages in that (i) by embedding the model into a higher dimension, the problem is now solved on a higher space dimension with one further partial differential equation for the phase-field function which is used to track the moving cell within the stationary and rather large space domain, and that extra modeling is required if topological changes are to be prohibited (for example by including a volume constraint term in the momentum balance law for example). Furthermore, refined meshes must be employed to resolve the partial differential equations as the small interface width order parameter is taken small. In this study, we will resort to modeling the sharp interface rather than the diffuse-interface formulation, for the moving boundary problem which is embedded naturally in the two-space dimension where a momentum balance law model is posed.
There exists many other noteworthy studies that have been carried out to model cell motility. Recently, a multiscale computational model coupling fluid mechanics, solid mechanics, and pattern formation was developed to simulate fully 3D, pseudopod-driven motility of amoeboid cells through a fluid filled porous medium [20]. The results and analysis presented showed a strong coupling between cell deformability and ECM properties. In [45], both two and three dimensions (2D and 3D) cell motility models were simulated using an activatorinhibitor system coupled with an evolving surface finite element method. In [46], a phase-field model for a 3D amoeboid cell was studied, including an activator-inhibitor system to describe the cell biochemistry, transport equations to describe cytosolic biochemistry dynamics and hydrodynamic drag to describe adhesive forces. Simulations were performed in 2D for cells navigating around obstacles on a substrate and in 3D for cells in rigid periodic cylindrical fibrous networks. In [47], a 2D model of motile cells was developed using an immersed-boundary method that resolves cell deformation, internal and external fluid flow and a reaction-diffusion system in the entire volume of the cell while in [48], a 2D model also using an immersed-boundary method in which the cytoskeleton is represented as a dynamic network of springs immersed in a fluid was presented. In [49], a deformable cell driven under a prescribed axisymmetric oscillating force in an unbounded medium was considered. A phase field approach to simulate migration of neutrophils in 3D in response to external cues was presented in [50]. In [51], a moving mesh finite element method for the approximate solution of coupled bulk-surface reaction-diffusion equations on an evolving two dimensional domain was studied. They used a moving mesh partial differential equation and generated bulk and surface meshes. They applied the method to model the twoway interaction of a migrating cell with an external chemotactic field. In [52,53], a Filament Based Lamellipodium Model (FBLM) which is a two-dimensional, anisotropic, two-phase model derived from microscopic description was studied. This model describes the actin network in terms of two transversal families of parallel filaments stabilized by cross-links and substrate adhesion. In [54], the FBLM was extended to investigate the effects of myosin polymers. The basic assumption here being that pairs of crossing actin filaments may be connected by myosin filaments.
Recently, Mackenzie et al. [55] investigated mass conservation property of a fully discrete finite element approximation of an ALE reformulation of a bulk-surface system on an evolving domain applied to a problem on cell polarization and chemotaxis on a moving domain. In this work, moving mesh methods were employed to move the internal mesh given the knowledge of the surface boundary through Dirichlet boundary conditions. In their work, the flow velocity generated by the force balance equation is assumed to be continuous across the boundary, thereby giving rise to global deformation of the mesh. Applications to cell polarization include the recent work by Cusseddu et al. [56] where the bulk-surface finite element method is used to generate numerical simulations over simple and complex geometries for the bulk-surface wave pinning model in 3-dimensions. In [57], a three-dimensional model for chemotactic motion of amoeboid cells was proposed. It accounted for the interactions between the extracellular substances, the membrane-bound proteins and the cytosolic components involved in the signaling pathway that lead to cell motility. Most recently, Moure and Gomez [58] studied the role of the cell nucleus in 2D cell migration where they used a computational model of fish keratocytes.
Our model is highly non-linear and is defined on a continuously evolving domain representing the cell. The nonlinearity of the model makes it difficult to obtain analytical solutions. Numerical methods are therefore a good choice in solving these models. Numerical methods for partial differential equations consist of two parts: a space discretization to transform the system of partial differential equations into a system of ordinary differential equations and a time discretization to transform the system of ordinary differential equations into a system of linear or non-linear algebraic equations depending on the time discretisation scheme applied. Finally, techniques from numerical linear algebra can be employed to solve the resulting system of linear equations. For the system of non-linear equations, a linearization for the non-linear terms is required in order to use numerical linear algebra. Space discretization methods include (but are not limited to) finite differences [59,60], finite elements [61][62][63], boundary elements [64] among other methods. The finite element method is well-known to handle complex and evolving cellular domains and can be easily generalized into multidimensions. It is a robust method that has been widely used to model cell motility and other problems involving reaction-diffusion equations in both stationary and continuously deforming domains [65][66][67][68][69][70][71][72][73].
Several time discretization have been used to obtain solutions for partial differential equations on both stationary and evolving domains. Fully explicit methods require very small time steps which result in computations that are expensive especially when it comes to multi-dimensions. Ruuth [74] presented different IMplicit-EXplicit (hence IMEX) schemes for solving reaction-diffusion systems on stationary domains, and their generalizations to growing domains were undertaken in [69].
The key essence of these schemes is that an implicit scheme is applied to approximate the diffusive term and an explicit scheme is used to approximate the reaction kinetics. The IMEX schemes presented in [69,74] include a first order semi-implicit backward differentiation formula (1-SBDF) which applies a backward differentiation formula to the diffusive term, Crank-Nicolson, Adams-Bashforth (CNAB) which applies Crank-Nicolson to the diffusive term and second order Adams-Bashforth to the reaction terms, Crank-Nicolson Leap Frog (CNLF) which applies something similar to Crank-Nicolson to the diffusive term and a leap frog to the reaction terms, the second order semi-implicit backward differentiation formula (2-SBDF), which applies a second order formula to the diffusive term, the third order semi-implicit backward differentiation formula (3-SBDF) which applies a third order formula to the diffusive term and the first order backward Euler finite difference scheme (1-SBEM) which treats both the diffusive term and linear part of the reaction term implicitly and non-linear part of the reaction semiimplicitly. We note that the 1-SBDF, the 2-SBDF and the 3-SBDF schemes are known to give strong decay of high frequency error components while unfortunately the CNAB and CNLF schemes are known to give a weak damping of high frequency error components [74]. From numerical experiments, the 2-SBDF is recommended as a good scheme to many two dimensional problems [69,74]. Recently, Madzvamuse and Chung [73] used a fully implicit scheme to solve a system of bulk-surface coupled reaction-diffusion equations. This scheme requires some special linearization techniques as shown in [73]. For this work, it is sufficient to apply the second order semi-implicit backward differentiation formula (2-SBDF) to discretize in time, the system of reaction-diffusion equations on the migrating cell, while a first order forward Euler's formula is used to update the moving mesh.
The focus of this work is to incorporate both actin and myosin II as well as convective terms in a reaction-diffusion model for cell motility where the models are posed on a moving domain with a sharp moving interface boundary. Our model involves a reaction-diffusion equation that considers myosin II activity in driving cell motion. We consider generalized reaction kinetics to model polymerization and depolymerization processes, where the kinetics take into account positive feedback from actin filaments. We formulate a system of reaction-diffusion equations to describe the actomyosin spatiotemporal dynamics during cell migration that in turn drives the actin filament structure here described by a viscous mechanical model in the absence of cell membrane and adhesion forces, hence our model is simpler than its phase-field counterpart [41]. Since the actomyosin system is always contained inside the cell and does not cross the cell boundary, we consider zero-flux boundary conditions. These boundary conditions ensure that biochemical concentrations do not flow across the boundary allowing for self-organization of the molecular species within the cell. The proposed framework is an alternative numerical or computational approach where phase-fields have been employed to offer numerical solutions on stationary domains which embed the migrating cell domain in a higher space-dimension [41,75]. This entails that we do not need to solve an extra equation for the phase-field function that is used to track the moving boundary. Hence, our main contribution is the development of an evolving finite element method for solving the moving boundary problem posed on its physical domain that is moving. We note that in the literature, different numerical methods have been adopted but the use of evolving finite elements for solving models that couple biomechanics with biochemistry are very few [71,72]. Furthermore, there is no need to employ refined mesh around the interface to resolve the numerical solution of the moving boundary problem within the small region encompassing the moving interface [76]. Our model is able to show expansion, contraction and deformation of the cell following actomyosin activity, with no substantial changes in shape. Our results are consistent with biological observations involving the migration of Keratocytes which are known to migrate with minimal changes in shape. Our simpler model captures the symmetry breaking process from a unit-circle shape to a polarized shape with a well-defined front and back in the absence of membrane and adhesion forces [56].
Hence, the structure of this article is as follows: in section 2, we present the model equations for cell migration and carry out non-dimensionalization of the model system. In section 3, we present a detailed moving finite element method (also now known as the evolving finite element method) and discretize in time to arrive at fully discrete equations. We also describe how we solve the discrete system of equations. We present various numerical results from our simulations in section 4. Finally, we conclude and outline future research directions in section 5.

A VISCOUS MODEL FOR CELL MIGRATION
We model the network of actin filaments in the cell as a viscous gel with active stresses generated from the action of actin and myosin II. The viscous model for cell migration considered in this work comprises of a system of reaction-advection-diffusion equations coupled to a force balance mechanical equation. Our current model does not include cellular adhesion and membrane forces that is brought about by tension and bending, which were considered in [41,75]. While adhesion plays an important role in cell migration, it has been shown that motility can occur with nearly no adhesion for example during cell swimming. Infact, neutrophils can migrate with little to no adhesion [20].

Model Formulation
Let t ⊂ R 2 be a simply connected, bounded and continuously deforming domain representing the cell shape at time t ∈ (0, T], T > 0 and ∂ t be the boundary of the cell with normal n = (n 1 , n 2 ) at a point x(t). At any point x(t) = (x(t), y(t)) ∈ t , let ρ m = ρ m (x(t), t) be the myosin II concentration, ρ a = ρ a (x(t), t) be the F-actin concentration in t and β = β(x(t), t) be the actin flow velocity which we assume is identical to the mesh/domain velocity of the cell. It must be noted that other works exist [51,55] where the flow velocity is different from the mesh velocity. We note that points x(t) move with velocity β(x(t), t), and that this velocity also describes the change of shape as it is also the velocity of the boundary points, through the continuity of boundary conditions. The model we consider is a simpler version of previous models for cell migration [41,75,77] and takes into account biochemical and mechanical structure of the cell. It is a viscous model that couples a force balance mechanical equation to reaction-advection-diffusion equations for actin and myosin II as summarized below.
To close the system, we apply zero-flux and stress-free boundary conditions given by and the terms in the reaction-diffusion system given by The terms ∇ · (ρ m β) and ∇ · (ρ a β) denote the convective and dilution processes due to cell movement and growth and these were not considered in some previous studies of this nature, for example [33,77]. Next, we introduce the quantities Dρ m Dt and Dρ a Dt , known as the material derivatives of ρ m and ρ a and these will be used throughout the sequel. The material derivatives are defined as respectively [62]. The reaction kinetics for the actin equation, f (ρ a , ρ cyt a ), depends on the F-actin ρ a and G-actin ρ cyt a concentrations [41,75] and that G-actin is assumed to be well mixed inside the cell [8]. Diffusion for myosin II D m (ρ a ) is assumed to depend on the F-actin concentration in such a way that diffusion is reduced when F-actin concentration increases [75]. We note that in the absence of F-actin, myosin II diffuses with the constant rate D 0 m . The constant D a is a positive diffusion coefficient for F-actin [41,75]. The terms σ myo (x, t), σ ν (x, t), and σ poly (x, t) are the myosin II driven contractile, viscous, and actin generated stresses, respectively, and are given by [41] where 1 2 ∇β(x, t) + (∇β(x, t)) T represents the rate of strain tensor, I is the identity tensor in two dimensions and δ(l) labels the cell periphery. The constants ν 0 , η 0 m , and η 0 a are the shear viscosity coefficient, myosin II contraction coefficient, and Factin protrusion coefficient, respectively. Table 1 shows a list of all the model parameters used. We have prescribed zeroflux boundary conditions to model spontaneous random cell migration and set initial conditions for both actin and myosin II. The polymerization stress only acts in the cell periphery. In order to describe this stress, we will assume an initial domain of a disk and specify that σ poly (x, t) only acts in the region which is some distance l from the centre of the disk domain. With this we define the function δ(l) to be of the form 1 if the point x is such its distance from the origin of the disk is more than l, 0 otherwise. This means that far from the cell periphery, only the contractile stress acts on the cell and close to the cell periphery both the contractile and the actin generated polymerization stresses act on the cell.

Remark
It must be noted that the model system (2.1a-2.1c), although posed on a two-dimensional deforming domain, generalizes immediately to three-dimensional volumes. Similarly, the application of the evolving finite element method to the model system also generalizes to three-dimensional volumes.

The Non-dimensionalized Model
We introduce non-dimensional rescaled variables as follows: Total actin density 800 [41] where R is the scaling factor for length. For notational simplicity, we drop all hats and the resulting non-dimensionalized viscous model for cell migration with non-dimensionalized parameters given in Table 2 therefore reads with initial conditions given by and boundary conditions where ρ a , ρ m , and β are the dependent variables for this model. The initial domain is now a disk with radius r = 1. Since the polymerization force is assumed to work only in the periphery of the cell [8], we let it act only in the region l ≥ 0.8. The delta function δ(l) ensures that polymerization force only act in this region.

THE EVOLVING FINITE ELEMENT METHOD FOR SIMULATING THE VISCOUS MODEL
We employ the evolving (moving grid) finite element method to derive a numerical scheme for the non-dimensionalized viscous model for cell migration. We use the evolving finite elements [61] to discretize the viscous model in space and apply a second order semi-implicit backward differentiation formula (2-SBDF) [69,74] to discretise the system of reactiondiffusion equations in time. To update the mesh after solving the equations in each timestep, we apply a first order forward Euler method. Convergence analysis for the reaction-diffusion system was undertaken on stationary domains (in the absence of the viscous model) and second-order convergence was established (results not shown). Demonstrating convergence for the full viscous model is an open problem due to the non-linear coupling between the mechanical and biochemical models. The evolving finite element method is an efficient method that is able to deal with complex and irregular geometries and has been widely used for growing and deforming domains [67][68][69]72].
From numerical experiments, the 2-SBDF is recommended as an efficient numerical time integrating scheme for reactiondiffusion equations. For more illustration on these, see for example [74]. Further information on the implementation of the evolving finite element method and time-integration can be found in [80]. The evolving finite element method therefore involves the following steps: derivation of weak formulation of the partial differential equations, the finite element spatial discretization to obtain a system of semi-discrete equations and a temporal discretization to obtain fully discrete system of equations.

Weak Formulation of the Reaction-Advection-Diffusion Equations
Reaction-advection-diffusion equations for myosin II and F-actin are, respectively, given by

1b)
where ρ m = ρ m (x(t), t) and ρ a = ρ a (x(t), t) are the myosin II and F-actin concentrations, respectively, and β = β(x(t), t) = (β 1 , β 2 ) is the flow velocity of the actin network. Here, we note that the reaction kinetics of F-actin only depends on ρ a variable and no reaction kinetics for the myosin II equation.
In order to obtain the weak formulation, we will rearrange (3.1a) and (3.1b). We apply product rule for gradient to the advection (3.2b) We recall that the quantities ∂ρ m ∂t + β · ∇ρ m and ∂ρ a ∂t + β · ∇ρ a are called material derivatives of ρ m and ρ a and are written as Dρ m Dt = ∂ρ m ∂t + β · ∇ρ m and Dρ a Dt = ∂ρ a ∂t + β · ∇ρ a . Now, using this definition for material derivatives above, we write (3.2a) and (3.2b) as respectively. In order to obtain the weak formulations, we multiply (3.3a) and (3.3b) by sufficiently smooth functions ψ 1 (x, t), ψ 2 (x, t) and integrate using Green's formula [81,82] in the domain t and use the boundary conditions ∂ρ a ∂n = ∂ρ m ∂n = 0. This yields We further use the product rule for the time derivatives in the Equations (3.4a) and (3.4b) and write Finally Reynold's transport theorem [83] gives (3.6b)

Weak Formulation of the Force Balance Equation
The force balance equation on an evolving domain t representing the cell is given by 7) where σ ν (x, t), σ myo (x, t), and σ poly (x, t) are the viscous, myosin II driven, and F-actin generated stresses, respectively. In order to write the weak formulation of the force balance above, we first decouple the stresses into x and y directions as follows. The force balance equation in x and y directions is therefore We then multiply the decoupled equations (3.9a) and (3.9b) by sufficiently smooth functions ψ 3 , ψ 4 , respectively, use Green's formula to integrate in the domain and apply the stress free boundary condition given. The boundary terms will vanish and the weak formulation thus reads: find for sufficiently smooth functions ψ 3 , ψ 4 , where f = δ(l)η 2 ρ a − η 1 ρ m . We find it convenient to rewrite the right hand sides of equations (3.10a) and (3.10b) such that we have the derivatives of the shape function ψ 3 and ψ 4 instead of function f which depends on the solution values. Applying Green's Theorem [82] to equations (3.10a) and (3.10b) gives where dS is the surface element and n = (n 1 , n 2 ) is the outward unit normal to the boundary.

Weak Formulation of the Coupled Problem
The weak formulation for the coupled problem is summarized as follows: (3.13) and similarly for the force balance equation, we have for all sufficiently smooth functions {ψ k } 4 k=1 ∈ H 1 ( t ), where f = δ(l)η 2 ρ a − η 1 ρ m .

Finite Element Discretisation
Numerical solutions for the weak formulations are defined in an infinite dimensional space H 1 ( t ). The essence of the finite element method is to seek solutions in a finite dimensional space.
We let h,t be the computational domain which is a polyhedral approximation to t . We define T h (t) to be a triangulation of h,t made up of non-degerate rectangular elements K i such that T h (t) = i K i . We call each K i an element of the mesh T h (t) where h is the diameter of the largest element. For the mesh T h (t), we require that it is made up of a finite number of elements and the elements must intersect along a complete edge, or at a vertex or not at all. The space discretisation is carried out using quadrilateral elements and we seek piece-wise linear approximation of the solution. We note that points in h,t now move with velocity β h = (β h 1 , β h 2 ) and that D Dt now stands for the material derivative with respect to the discrete flow velocity β h . We define the finite element space V h (t) by (3.16) We will seek finite element solutions of the viscous model in this space. The discretised version of (3.12-3.15) therefore reads: find 18) and similarly for the force balance equation, we have We define a basis function for the space V h (t) by φ i (x, t) ∈ V h (t) for i = 1, 2, ..., N h such that where x j (t) is the jth nodal point of the mesh and N h is the total number of degrees of freedom of the nodes. We seek finite element approximations of the form We note that the shape function φ j is a function of time t. We will make use of the following Lemma: Lemma: Transport property of the basis functions: The finite element space on the discretized domain is a space of continuous piece-wise linear functions whose nodal basis functions have the remarkable property on element K for all φ i where the derivative denotes the material derivative [84].
This gives respectively, for all i = 1, 2, ..., N h . The parameter a mon (t) represents the well mixed actin monomers concentration at time t and corresponds to the integral h,t ρ cyt,h a d h,t . Now, integrating over h,t yields the semi-discrete equations for the reaction-advection-diffusion equations as

Semi-discrete Equations for the Force Balance Equation
In (3.19) and (3.20), we substitute {ψ h l (x, t)} 4 l=3 by φ i (x(t), t), i = 1, 2, ..., N h and β h 1 , β h 2 , ρ h m , ρ h a by their corresponding finite element approximations and integrate over h,t . The semidiscrete equation in x-direction will be for all i = 1, 2, ..., N h and can be rearranged as for all i = 1, 2, ..., N h . We note that (3.27) can be written as a 11 (t)U 1 (t) + a 12 (t)U 2 (t) + · · · + a 1N h (t)U N h (t) + b 11 (t)V 1 (t) where a ij (t), b ij (t), and F 1 i (t) are integrals over h,t given by This means we can split the left hand side of the above system of equations into two parts: one with U j (t) and the other with V j (t). Similarly, the y-direction of the force balance equation is for all i = 1, 2, ..., N h , and expanded as with c ij (t), d ij (t), and F 2 i (t) integrals over h,t given by These systems of equations can be written more compactly in matrix-vector form as follows . . .

Semi-discrete Equations for the Coupled Problem
The semi-discrete equations for the reaction-advection-diffusion and the force balance equation are of the form.

Fully Discrete Model
To obtain fully discrete equations for the coupled problem, we discretise the time interval (0, T] into a finite number N of uniform sub-intervals with sub-interval size τ = t n+1 − t n and write t n = τ n. From the semi-discrete equations we employ the 2-SBDF as follows + 2 −eM(t n )ρ(t n ) + k 4 a mon (t n )H(t n ) + k 3 a mon (t n )L(ρ(t n )) − −eM(t n−1 )ρ(t n−1 ) + k 4 a mon (t n−1 )H(t n−1 ) +k 3 a mon (t n−1 )L(ρ(t n−1 )) .

One-Step Backward Euler Scheme
We will need solutions at the last two time steps in order to implement the 2-SBDF. We therefore implement a onestep backward Euler method to solve for ρ(t 1 ) and ω(t 1 ) and then proceed with 2-SBDF scheme. A one-step backward Euler scheme for (3.33) and (3.34) is given by.

Displacements of the Nodes
We now describe how the nodes of the mesh will be displaced to obtain a new nodal location. The position of any new node will be a function of its current position and the amount of displacement it has achieved. Let t n+1 = t n + τ and consider the points x(t n ) ∈ h,t n and x(t n+1 ) ∈ h,t n+1 in the respective domains. We can define a first order linear approximation of the flow velocity as follows: This means that the new domain can be approximated by where ξ (x(t n )) is the solution of the force balance equation at time t n . We note that the mesh is updated at each time step and that τ ξ (x(t n )) is the displacement from point x(t n ) to x(t n+1 ). Thus, the new nodal position will be given as the sum of the current node and its displacement. It must be observed that (3.42) is simply an update of the mesh from one time level to the other, hence a first order forward Euler scheme is sufficient.

NUMERICAL RESULTS
The resulting numerical scheme is a coupled system of linear equations. To implement the discretization, we use deal.II library which is an open source and efficient finite element library written in C++ [85]. The coupled systems of linear equations are solved using a preconditioned conjugate gradient (PCG) and generalized minimal residual (GMRES) methods [86][87][88].
Let t = nτ , where τ denotes the time-step size and n the number of time-steps. We consider the unit disk 0 as the initial domain and subdivide this domain using quadrilateral elements. We would like to make a note that our numerical method was validated for convergence and stability for the case of a stationary non-migrating cell (results not shown). We showed that the 2-SBDF method used to solve the system of reaction-diffusion equations is a second order.
We only exhibit simulations of an evolving cell. We used a finite element mesh with 5185 quadrilateral elements and applied a 2-SBDF scheme with time-step τ = 10 −3 to compute solutions to final time t = 4. The parameters used for the simulations are displayed in Table 2. The solution for this model is in the form of F-actin and myosin II concentrations and the speed of the cell |(β h 1 , β h 2 )|. F-actin and myosin II solutions are the solution of the reaction-advection-diffusion equations (2.4b) and (2.4c) while speeds of the cell come from solution of the force balance equation (2.4a). We begin simulations on a unit disk to represent the cell at initial time with zero initial speed. Initial conditions are random perturbations about ρ a = ρ m = 1.0 as shown in Figure 1.
We observe that the initial conditions chosen determine the dynamics of F-actin and cell shape. Choosing a perturbation about ρ a = 1.0 as the initial condition for ρ a variable leads to a uniform expansion of the cell. Figure 2 shows one set of solution set with initial data for the concentrations as random perturbation about ρ m = 1 and ρ a = 1 and all other parameters as in Table 2. In Figure 2C, we observe that the solution β h has its highest value around the points where x 2 + y 2 = 0.8. These correspond to the points where δ(l) is discontinuous. Figure 2D shows change in area with time of the evolving cell. We varied the initial condition for ρ a variable. Choosing a nonzero concentrations of ρ a = 1.0 only in one half of the cell leads to symmetry breaking where the cell identifies its front and rear. This leads to an irregular expansion of the cell and hence in a directed migration of the cell toward the direction with high F-actin concentrations. Figures 3A-E is a solution set when considering initial data for myosin II as random perturbation about ρ m = 1, initial condition for ρ a non-zero only in one half of the cell and all other parameters as in Table 2. Figures 3E,F show the initial mesh and mesh at the end of the simulation, respectively.
Actin changes from its active state (F-actin) to inactive state (G-actin) and vice-versa through polymerization and depolymerization processes and hence the total amount of actin is conserved at all time as shown in Figure 4A. F-actin assembles together causing expansive stress on the cell. We varied the parameter for total amount of actin, ρ tot a , in the cell by considering three cases: ρ tot a = 10, increasing the parameter  to ρ tot a = 16 and decreasing the parameter to ρ tot a = 5. We observed that the more the total amount of actin, the more the expansion of the cell. Figures 4B-D show change in the area of the evolving cell with time as a result of varying the parameter ρ tot a . In our model, myosin II only diffuses inside the cell and therefore no reaction kinetics involving myosin II variable.

CONCLUSION AND FUTURE DIRECTIONS
In this work, we considered an alternative approach to classical phase-fields approach by formulating the model equations on a sharp interface moving boundary problem for a single cell migrating on a two-dimensional moving domain. The problem involves a system of reaction-diffusion equations posed on a growing domain, where the velocity of the moving actin filaments is identical to the flow velocity of the domain and this velocity is modeled by a momentum balance equation within a viscous framework. In this approach the flow of the actin network is treated as a Newtonian fluid. We then employed the evolving finite element method to provide numerical approximate solutions of the model equations on the physical Lagrangian moving domain. This method does not allow for topological changes in cell shape unlike the phase-field model and also does not require a sufficiently refined mesh close to the sharp moving boundary interface. Furthermore, the method renders itself naturally to three-dimensional extensions.
Actin filaments and myosin II are the main sources of stresses in the cell and are responsible for driving cell migration. We varied the parameter for total actin and observed a linear relationship between the cell expansion and total amount of actin. The more the total amount of actin, the more the expansion. A decrease in the total amount of actin beyond a certain threshold leads to cell shrinking. Myosin II is responsible for cell contraction. By varying the contraction coefficient for myosin II, we observed contraction of the cell. The initial condition also played a role in the dynamics of cell migration. We considered two sets of initial conditions for the F-actin and myosin II variables. A random perturbation about ρ a = 1 led to uniform expansion of the cell where the periphery of the cell expanded or contracted uniformly. By considering a non-zero initial concentration of F-actin only in one half of the cell, we observed a directed growth of the cell where the cell expanded in the direction with more actin concentration and began to migrate in that direction. The findings and conclusion from our work is therefore as follows: in the absence of advection of actin and myosin II and domain evolution, the biochemical model for F-actin and myosin II will reach steady state and that the actomyosin system is responsible for driving cell migration. Some of the parameters and variables that are important in the dynamics of cell migration are: the initial conditions for F-actin, the total amount of actin inside the cell, the contraction and polymerization coefficients.
In this paper, no rigorous numerical convergent analysis has been undertaken due to the complex nature of the coupling between mechanics and biochemistry during domain growth and this aspect forms part of our current studies. Furthermore, extensions of this work include the introduction  Table 2, (C) area of the cell with an increased total amount of actin ρ tot a = 16 while all other parameters in Table 2 held constant and (D) decreasing area of the cell with a reduced total amount of actin ρ tot a = 5 while all other parameters in Table 2 held constant.
of a mechanism for volume constraint or conservation, incorporation of cell adhesion and membrane forces, cell-to-cell interactions, extensions to many cells, cell migration through obstacles and formulating experimentally driven actomyosin reaction dynamics, just to mention a few.