Original Research ARTICLE
Integrating Actin and Myosin II in a Viscous Model for Cell Migration
- Department of Mathematics, School of Mathematical and Physical Sciences, University of Sussex, Brighton, United Kingdom
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.
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–6]. According to experimental observations, cell migration in eukaryotic cells is powered by actin-myosin system . 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 . 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–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–19]. Infact, neutrophils can migrate with little to no adhesion .
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 . 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 . 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 . Many models are built on the concept that cell motility is composed of the following stages: protrusion, adhesion and contraction . 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 . 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 . These early works introduced fundamental ideas that are still used in developing complex models . In , “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 , 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 . This model was able to predict an optimal angle between the actin filament and the load for the effective force transmission. In , 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 . In [25–27], actin polymerization and depolymerization were treated as stochastic processes.
In , 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 .
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 , a one-dimensional viscoelastic model of the cytoplasm and active stress generation in , a one-dimensional model for the actin distribution and its reaction in , a two-dimensional elastic continuum model in  and a cytomechanical model that couples a force balance mechanical equation for actin network to a reaction-diffusion equation for actin . This cytomechanical model was later extended in  where they used a cartesian coordinate system.
Another strategy is to use phase-field models as described in these works [35–43]. In , 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 phase-field 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 multi-component vesicles was studied in . In , a phase-field 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  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 . The results and analysis presented showed a strong coupling between cell deformability and ECM properties. In , both two and three dimensions (2D and 3D) cell motility models were simulated using an activator-inhibitor system coupled with an evolving surface finite element method. In , 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 , 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 , 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 , 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 . In , 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 two-way 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 , 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.  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.  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 , 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  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 non-linearity 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–63], boundary elements  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–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  presented different IMplicit-EXplicit (hence IMEX) schemes for solving reaction-diffusion systems on stationary domains, and their generalizations to growing domains were undertaken in . 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 semi-implicitly. 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 . From numerical experiments, the 2-SBDF is recommended as a good scheme to many two dimensional problems [69, 74]. Recently, Madzvamuse and Chung  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 . 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 . 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 . 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 .
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.
2. 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 .
2.1. Model Formulation
Let 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 = (n1, n2) 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.
with initial conditions
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 and , known as the material derivatives of ρm and ρa and these will be used throughout the sequel. The material derivatives are defined as
respectively . The reaction kinetics for the actin equation, , depends on the F-actin ρa and G-actin concentrations [41, 75] and that G-actin is assumed to be well mixed inside the cell . Diffusion for myosin II Dm(ρa) is assumed to depend on the F-actin concentration in such a way that diffusion is reduced when F-actin concentration increases . We note that in the absence of F-actin, myosin II diffuses with the constant rate . The constant Da 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 
where represents the rate of strain tensor, I is the identity tensor in two dimensions and δ(l) labels the cell periphery. The constants ν0, , and are the shear viscosity coefficient, myosin II contraction coefficient, and F-actin protrusion coefficient, respectively. Table 1 shows a list of all the model parameters used. We have prescribed zero-flux 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
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.
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.
2.2. The Non-dimensionalized Model
We introduce non-dimensional rescaled variables as follows:
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 , we let it act only in the region l ≥ 0.8. The delta function δ(l) ensures that polymerization force only act in this region.
3. 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  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 reaction-diffusion 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–69, 72]. From numerical experiments, the 2-SBDF is recommended as an efficient numerical time integrating scheme for reaction-diffusion equations. For more illustration on these, see for example . Further information on the implementation of the evolving finite element method and time-integration can be found in .
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.
3.1. Derivation of the Weak Formulation
3.1.1. Weak Formulation of the Reaction-Advection-Diffusion Equations
Reaction-advection-diffusion equations for myosin II and F-actin are, respectively, given by
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 terms and write the equations as
We recall that the quantities and are called material derivatives of ρm and ρa and are written as and . 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 . 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  gives
3.1.2. Weak Formulation of the Force Balance Equation
The force balance equation on an evolving domain Ωt representing the cell is given by
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 such that
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  to equations (3.10a) and (3.10b) gives
where dS is the surface element and n = (n1, n2) is the outward unit normal to the boundary.
3.1.3. Weak Formulation of the Coupled Problem
The weak formulation for the coupled problem is summarized as follows:
find such that
and similarly for the force balance equation, we have
for all sufficiently smooth functions , where f = δ(l)η2ρa − η1ρm.
3.2. Finite Element Discretisation
Numerical solutions for the weak formulations are defined in an infinite dimensional space . 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 Th(t) to be a triangulation of Ωh,t made up of non-degerate rectangular elements Ki such that . We call each Ki an element of the mesh Th(t) where h is the diameter of the largest element. For the mesh Th(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 and that now stands for the material derivative with respect to the discrete flow velocity βh. We define the finite element space Vh(t) by
We will seek finite element solutions of the viscous model in this space. The discretised version of (3.12–3.15) therefore reads: find
and similarly for the force balance equation, we have
for all , where
We define a basis function for the space Vh(t) by ϕi(x, t) ∈ Vh(t) for i = 1, 2, …, Nh such that
where xj(t) is the jth nodal point of the mesh and Nh 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 .
3.2.1. Semi-discrete Equations for the Reaction-Advection-Diffusion Equations
In equations (3.17) and (3.18), we substitute by ϕi(x, t), i = 1, 2, …, Nh and by their corresponding finite element approximations and use (3.23).
respectively, for all i = 1, 2, …, Nh. The parameter amon(t) represents the well mixed actin monomers concentration at time t and corresponds to the integral . Now, integrating over Ωh,t yields the semi-discrete equations for the reaction-advection-diffusion equations as
where and are the solution vectors and amon(t) is actin monomers concentration at time t. M(t) is the time-dependent global mass matrix, K(t) is the time-dependent global stiffness matrix, H(t) is the time-dependent global force vector and S(ρ(t)) and L(ρ(t)) are the time-dependent matrix and vector, respectively, which are functions of the solution vectors. These are given by.
3.2.2. Semi-discrete Equations for the Force Balance Equation
In (3.19) and (3.20), we substitute by ϕi(x(t), t), i = 1, 2, …, Nh and by their corresponding finite element approximations and integrate over Ωh,t. The semi-discrete equation in x-direction will be
for all i = 1, 2, …, Nh and can be rearranged as
for all i = 1, 2, …, Nh. We note that (3.27) can be written as
where aij(t), bij(t), and 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 Uj(t) and the other with Vj(t). Similarly, the y-direction of the force balance equation is
for all i = 1, 2, …, Nh, and expanded as
with cij(t), dij(t), and integrals over Ωh,t given by
These systems of equations can be written more compactly in matrix-vector form as follows
where aij, bij, cij and dij are functions of time t. The vectors and are the solution vectors. We let
and write (3.29) in block-vector form as
We also note from (3.27) and (3.28) that C(t) = (B(t))T and write (3.30) as
and have the following semi-discrete equation for the force balance equation.
3.2.3. 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.
3.3. 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 τ = tn+1 − tn and write tn = τn. From the semi-discrete equations
we employ the 2-SBDF as follows
For (3.35), we compute A(t) and F(ρ(t), ω(t)) at time step tn as follows
The above equations yield the following fully discrete equations
with τ as the time-step size, ρ(tn+1) and ω(tn+1) are the actin filaments and myosin II solutions at time tn+1 and ξn+1 is the velocity solutions at time tn+1.
3.3.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 one-step backward Euler method to solve for ρ(t1) and ω(t1) and then proceed with 2-SBDF scheme. A one-step backward Euler scheme for (3.33) and (3.34) is given by.
3.3.2. 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 tn+1 = tn + τ and consider the points and 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(tn)) is the solution of the force balance equation at time tn. We note that the mesh is updated at each time step and that τξ(x(tn)) is the displacement from point x(tn) to x(tn+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.
4. 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++ . The coupled systems of linear equations are solved using a preconditioned conjugate gradient (PCG) and generalized minimal residual (GMRES) methods [86–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 . 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.
Figure 1. Initial conditions for the variables. We begin simulations by considering initial concentrations of the variables, respectively, as random perturbations about (A) ρm = 1.0, (B) ρa = 1.0, and (C) stationary cell at initial time. Blue signifies lowest values while red highest values.
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 x2 + y2 = 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 non-zero 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.
Figure 2. Graphical display of numerical solution illustrating a growing cell at time t = 4. This is a case when considering initial condition as random perturbations about ρa = ρm = 1.0. Blue signifies lowest values while red signifies highest values. (A) solution for ρm variable, (B) solution for ρa variable, (C) cell speed at final time showing higher speeds as the concentration of F-actin increases, and (D) area of the growing cell as a function of time showing increase in area of the cell with time.
Figure 3. Graphical display of numerical solution illustrating a migrating cell at time t = 4 when considering a non-zero initial condition of ρa = 1 only in one half of the cell. Blue signifies lowest values while red signifies highest values. (A) Solution for ρm variable, (B) solution for ρa variable and (C) cell speed at final time. Polymerization of actin leads to cell deformation where the cell expands more in the regions with higher actin filament concentrations. (D,E) are the finite element meshes of the simulations at initial and final times, respectively, which have been enlarged for clarity.
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, , in the cell by considering three cases: , increasing the parameter to and decreasing the parameter to . 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 . In our model, myosin II only diffuses inside the cell and therefore no reaction kinetics involving myosin II variable.
Figure 4. (A) Conservation of mass for actin with total amount of actin as , (B) area of the growing cell as a function of time for the case as shown in Table 2, (C) area of the cell with an increased total amount of actin while all other parameters in Table 2 held constant and (D) decreasing area of the cell with a reduced total amount of actin while all other parameters in Table 2 held constant.
5. 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 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.
Data Availability Statement
All datasets generated for this study are included in the article.
All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.
This work (AM) was partly supported by the EPSRC grant number EP/J016780/1, the European Union Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement no. 642866, the Commission for Developing Countries, and the Simons Foundation. AM is a Royal Society Wolfson Research Merit Award Holder funded generously by the Wolfson Foundation. AM is also a Distinguished Visiting Scholar to the University of Johannesburg, Department of Mathematics, South Africa. This work (BK) was carried out when BK was a graduate researcher at University of Sussex supported by University of Sussex Chancellor's International Research Scholarship.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
10. Pollard TD, Blanchoin L, Mullins RD. Molecular mechanisms controlling actin filament dynamics in nonmuscle cells. Annu Rev Biophys Biomol Struct. (2000) 29:545–76. doi: 10.1146/annurev.biophys.29.1.545
13. Cavalcanti-Adam EA, Volberg T, Micoulet A, Kessler H, Geiger B, Spatz JP. Cell spreading and focal adhesion dynamics are regulated by spacing of integrin ligands. Biophys J. (2007) 92:2964–74. doi: 10.1529/biophysj.106.089730
15. Shemesh T, Bershadsky AD, Kozlov MM. Physical model for self-organization of actin cytoskeleton and adhesion complexes at the cell front. Biophys J. (2012) 102:1746–56. doi: 10.1016/j.bpj.2012.03.006
25. Marée AF, Jilkine A, Dawes A, Grieneisen VA, Edelstein-Keshet L. Polarization and movement of keratocytes: a multiscale modelling approach. Bull. Math. Biol. (2006) 68:1169–211. doi: 10.1007/s11538-006-9131-7
31. Mogilner A, Marland E, Bottino D. A minimal model of locomotion applied to the steady gliding movement of fish keratocyte cells. In: Maini PK, and Othmer HG, editor. Mathematical Models for Biological Pattern Formation. New York, NY: Springer (2001). p. 269–93. doi: 10.1007/978-1-4613-0133-2_12
37. Liu P, Zhang Y, Cheng Q, Lu C. Simulations of the spreading of a vesicle on a substrate surface mediated by receptor-ligand binding. J Mech Phys Solids. (2007) 55:1166–81. doi: 10.1016/j.jmps.2006.12.001
39. Lowengrub JS, Rätz A, Voigt A. Phase-field modeling of the dynamics of multicomponent vesicles: Spinodal decomposition, coarsening, budding, and fission. Phys Rev E. (2009) 79:031926. doi: 10.1103/PhysRevE.79.031926
42. Camley BA, Zhao Y, Li B, Levine H, Rappel WJ. Crawling and turning in a minimal reaction-diffusion cell motility model: coupling cell shape and biochemistry. Phys Rev E. (2017) 95:012401. doi: 10.1103/PhysRevE.95.012401
47. Vanderlei B, Feng JJ, Edelstein-Keshet L. A computational model of cell polarization and motility coupling mechanics and biochemistry. Multiscale Model Simul. (2011) 9:1420–43. doi: 10.1137/100815335
49. Farutin A, Rafaï S, Dysthe DK, Duperray A, Peyla P, Misbah C. Amoeboid swimming: a generic self-propulsion of cells in fluids by means of membrane deformations. Phys Rev Lett. (2013) 111:228102. doi: 10.1103/PhysRevLett.111.228102
51. MacDonald G, Mackenzie JA, Nolan M, Insall R. A computational method for the coupled solution of reaction-diffusion equations on evolving domains and manifolds: application to a model of cell migration and chemotaxis. J Comput Phys. (2016) 309:207–6. doi: 10.1016/j.jcp.2015.12.038
52. Ölz D, Schmeiser C. How do cells move? Mathematical Modeling of Cytoskeleton Dynamics and Cell Migration. Cell Mechanics: From Single Scale-Based Models to Multiscale Modelling. Boca Raton: Chapman and Hall/CRC (2010). doi: 10.1201/9781420094558-c5
53. Manhart A, Oelz D, Schmeiser C, Sfakianakis N. An extended Filament Based Lamellipodium Model produces various moving cell shapes in the presence of chemotactic signals. J Theor Biol. (2015) 382:244–58. doi: 10.1016/j.jtbi.2015.06.044
55. Mackenzie JA, Rowlatt CF, Insall RH. A conservative finite element ALE scheme for mass-conserving reaction-diffusion equations on evolving two-dimensional domains. (2019). arXiv[Preprint]. arXiv:191002282.
66. Masud A, Hughes TJ. A space-time Galerkin/least-squares finite element formulation of the Navier-Stokes equations for moving domain problems. Comput Methods Appl Mech Eng. (1997) 146:91–126. doi: 10.1016/S0045-7825(96)01222-4
68. Madzvamuse A, Maini PK, Wathen AJ. A moving grid finite element method for the simulation of pattern generation by Turing models on growing domains. J Sci Comput. (2005) 24:247–62. doi: 10.1007/s10915-004-4617-7
69. Madzvamuse A. Time-stepping schemes for moving grid finite elements applied to reaction-diffusion systems on fixed and growing domains. J Comput Phys. (2006) 214:239–63. doi: 10.1016/j.jcp.2005.09.012
70. Manhart A, Oelz D, Schmeiser C, Sfakianakis N. Numerical treatment of the filament-based lamellipodium model (FBLM). In: Graw F, Pahle J, Matthaus F, editors. Modeling Cellular Systems Cham: Springer (2017).
73. Madzvamuse A, Chung AH. Fully implicit time-stepping schemes and non-linear solvers for systems of reaction-diffusion equations. Appl Math Comput. (2014) 244:361–74. doi: 10.1016/j.amc.2014.07.004
76. Yang FW, Venkataraman C, Styles V, Madzvamuse A. A robust and efficient adaptive multigrid solver for the optimal control of phase field formulations of geometric evolution laws. Commun Comput Phys. (2017) 21:65–92. doi: 10.4208/cicp.240715.080716a
86. Barrett R, Berry MW, Chan TF, Demmel J, Donato J, Dongarra J, et al. Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods. Vol. 43. Philadelphia, PA: SIAM (1994) doi: 10.1137/1.9781611971538
Keywords: cell migration, viscous model, reaction-advection-diffusion equations, force balance equation, evolving finite element method
Citation: Madzvamuse A and Kiplangat BK (2020) Integrating Actin and Myosin II in a Viscous Model for Cell Migration. Front. Appl. Math. Stat. 6:26. doi: 10.3389/fams.2020.00026
Received: 02 March 2020; Accepted: 25 June 2020;
Published: 06 August 2020.
Edited by:Sebastian Aland, Hochschule für Technik und Wirtschaft Dresden, Germany
Reviewed by:Bjorn Stinner, University of Warwick, United Kingdom
Christian Engwer, University of Münster, Germany
Copyright © 2020 Madzvamuse and Kiplangat. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Anotida Madzvamuse, A.Madzvamuse@sussex.ac.uk