A moving grid finite element method applied to a mechanobiochemical model for 3D cell migration
AA moving grid finite element method applied to a mechanobiochemical modelfor 3D cell migration.
Laura Murphy, Anotida Madzvamuse
University of Sussex, School of Mathematical and Physical Sciences,Department of Mathematics,BN1 9QH, UK
Abstract
This work presents the development, analysis and numerical simulations of a biophysical model for 3D cell deforma-tion and movement, which couples biochemical reactions and biomechanical forces. We propose a mechanobiochem-ical model which considers the actin filament network as a viscoelastic and contractile gel. The mechanical propertiesare modelled by a force balancing equation for the displacements, the pressure and concentration forces are driven byactin and myosin dynamics, and these are in turn modelled by a system of reaction-di ff usion equations on a movingcell domain. The biophysical model consists of highly non-linear partial di ff erential equations whose analytical solu-tions are intractable. To obtain approximate solutions to the model system, we employ the moving grid finite elementmethod. The numerical results are supported by linear stability theoretical results close to bifurcation points duringthe early stages of cell migration. Numerical simulations exhibited show both simple and complex cell deformationsin 3-dimensions that include cell expansion, cell protrusion and cell contraction. The computational framework pre-sented here sets a strong foundation that allows to study more complex and experimentally driven reaction-kineticsinvolving actin, myosin and other molecular species that play an important role in cell movement and deformation. Keywords:
Mechanobiochemical model, viscoelastic, force balance equation, cell motility, moving grid finiteelements, reaction-di ff usion equations, partial di ff erential equations, moving boundary problem
1. Introduction
Cell movement is critical in multicelluar organisms due to roles in embryogenesis, wound healing, immune re-sponse, cancer metastasis, tumour invasion, and other processes, therefore, understanding cell movement is of greatimportance to medicine and to understanding our origins [8, 9, 11, 14]. In this study we consider a mechanobiochemi-cal model previously studied by George et al. [15, 16, 27] which we will extend to 3-dimensions as well as introducingfor the first time, the role of myosin in the modelling and computational framework. The model comprises of a systemof reaction-di ff usion equations for cellular proteins and a viscoelastic mechanical model for cell movement and defor-mation. Given that the model is highly nonlinear, exact analytical solutions are not possible to obtain in closed form,instead, we will seek to compute numerical approximations to these exact solutions. Numerical methods abound forsolving complex partial di ff erential equation (PDEs), methods that have been employed to model cell motility includefinite di ff erences, phase field methods, boundary element methods (BEM), immersed boundary methods or level setmethods (LSM), [1, 7, 32, 35, 44, 45, 47]. Choosing a suitable method for a particular model is a balance between theease of application within the model’s framework and the reliability of solutions produced.Finite di ff erences were used in previous incarnations of the model [1, 44]. This method is very useful and easyto implement on fixed and simple domains but it is significantly more complicated to incorporate for the evolvingdomains and surfaces we wish to use and there are often problems with a moving boundary. Level set methods areused extensively in cell simulations and are useful when cells split and reconnect, therefore, it may be advantageous Email address:
[email protected] (Laura Murphy) a r X i v : . [ q - b i o . CB ] M a r o use this method in the future when considering cell proliferation (cell division) and apoptosis (cell death) [48]. Inthis work we are not concerned with cells splitting.The finite element method is well known to easily handle complex and evolving cellular domains and can begeneralised to multidimensions with little complication, hence it is the ideal method to numerically solve our modelsystem. Finite element methods have been widely used to model cell motility [5, 6, 10, 13, 18, 24, 28, 38, 41, 46],and can be implemented in diverse ways depending on the model.Given these considerations we develop a finite element based formulation which follows the work of George [16]with the extension into multidimensions and involving solving two reaction di ff usion equations. Additionally wedevelop our numerical solver based on deal.ii rather than ALBERTA as previously done by George [16].This article is therefore structured as follows: In Section 2, we introduce the mechanobiochemical model. Theoret-ical predictions of the spatiotemporal behaviour of the solutions of the model close to bifurcation points are presentedin Appendix A, these identify important bifurcation parameters. In Section 3, theoretical predictions are used tovalidate finite element simulations since there are no analytical solutions to compare with. Further finite element sim-ulations illustrating 3D cell movement and deformation are exhibited in Section 4. We discuss findings and concludein Section 5.
2. A mechanobiochemical model
The model we consider and extend is inspired by contractile models of the actin cytogel [23, 34]. These modelscomprise of a force balance equation modelling the displacements of the cell when deformed and a reaction-di ff usionequation for the concentration of the gel that in turn drives cell movement. The idea of pressure driven protrusionand the use of concentration of actin originates from Alt and Tranquillo [1]. In their model they assume movement isproduced by a balance between contractile force of the actin network pulling on the membrane and pressure pushingon the membrane. This was extended by Stephanou et al. so that large deformations could be modelled which is morerealistic for most cells [44]. George further extended this model by observing (and hence modelling) that higher actinconcentration in a region leads to more pressure [16]. In the previous models a polar coordinate system was used andradial extension of the cell was calculated [1, 44]. Unlike this approach, we follow the work of George and study themechanobiochemical model in its physical Cartesian coordinates without any need for coordinate transformation [16].We extend of the work by George in two key ways, firstly by extending from two, to three dimensions and secondly byadding the consideration of the concentration of myosin. We model the concentrations of, and interactions between,actin and myosin using two reaction-di ff usion equations. The reaction-di ff usion equations are coupled to a forcebalance equation which describes the movement of the cell. The model equations are outlined as follows.We assume that the cell shape is a simply connected and continuously deforming domain: Ω t ⊂ R with boundary ∂ Ω t , where t ∈ I = [0 , T f ] , T f >
0. Any point x ∈ Ω t is defined by x = ( x ( t ) , y ( t ) , z ( t )). We define the displacementof x at time t by u = ( u ( x ( t ) , t ) , v ( x ( t ) , t ) , w ( x ( t ) , t )) T . Let the concentration of F-actin, and bound myosin, at point x ( t ) be given by a = a ( x ( t ) , t ), and m = m ( x ( t ) , t ), respectively. We describe the dynamics of the actin network by thefollowing system of the equations, ∇ · ( σ v + σ e + σ c + σ p ) = in Ω t , t ∈ I , (1a) ∂ a ∂ t + ∇ · ( a β ) − D a ∆ a − f ( a , m ) = Ω t , t ∈ I , (1b) ∂ m ∂ t + ∇ · ( m β ) − D m ∆ m − g ( a , m ) = Ω t , t ∈ I , (1c) a ( x ( t ) , t ) = a , u ( x ( t ) , t ) = for x ∈ Ω , (1d) β = ω n for x ∈ ∂ Ω t , t ∈ I , (1e) σ v · n = σ e · n = n · ∇ a = n · ∇ m = x ∈ ∂ Ω t , t ∈ I , (1f)where the viscoelastic and contractile properties are described by stress tensors: • viscous σ v ( u ) = µ ∂ ε ∂ t + µ ∂φ∂ t I where µ and µ are shear and bulk viscosities respectively, ε is the the straintensor ( ( ∇ u + ∇ u T )) and φ is the dilation ( ∇ · u ). 2 elastic σ e ( u ) = E + v ( ε + ν − ν φ I ) where E is the Youngs modulus and ν is the Poisson ratio. • contractile σ c ( a , m ) = ( ψ a e − a / a sat + cm ) I , where ψ and c are the contractility coe ffi cients for a and m , respec-tively, and a sat is the saturation coe ffi cient of actin. • pressure σ p ( u , a ) = p + φ (cid:0) + π δ ( l ) arctan a (cid:1) I . This describes two types of pressure. First is the hydrostaticpressure which is present everywhere and corresponds to the osmotic pressure in the cell which depends on thedilation φ and pressure coe ffi cient p . Close to the membrane there is also polymerisation pressure caused bythe polymerising actin filaments pushing on the cell membrane. This increases with increasing concentrationof filaments a . We choose close to the membrane to mean less than 20% of the cell radius from the edge in theinitial state. To define this we use δ ( l ) and the points ξ = ( ξ x , ξ y , ξ z ) ∈ Ω . There exists a family of bijectivemappings between the initial and current domains we can let l : Ω t × I → R and corresponding ˆ l : Ω × I → [0 , l ( ξ , t ) = l ( x ( ξ , t ) , t ). So we calculate the distance from the centroid by δ ( l ) = , if (cid:113) ξ x + ξ y + ξ z > . , , otherwise. (2)In the reaction-di ff usion equations we have the di ff usion coe ffi cients for actin and myosin, D a and D m respectively.Because the cell is moving we introduce the flow velocity β = ∂ u ∂ t . The interactions between actin and myosin aredescribed by the reaction terms f ( a , m ) and g ( a , m ) respectively. We have formulated di ff erent plausible reactionkinetics in the absence of experimental data and for the sake of brevity, we only present here one such model. Werefer the interested reader to consult [30]. Although we will discuss the implementation of one specific plausiblemodel, other models can be easily incorporated and studied in a similar fashion. For illustrative purposes we considerthe following hypothetical reaction kinetics f ( a , m ) = k a ( a c − a ) + k am a ( m c − m )1 + Ka , (3a) g ( a , m ) = − k ma ( a c − a ) − k am a ( m c − m )1 + Ka , (3b)where we begin with the same reaction term, k a ( a c − a ), as used in [15, 16, 27]. k a is the rate of polymerisa-tion / depolymerisation and a c is the equilibrium concentration and if the concentration is above this critical valuethen F-actin will depolymerise at the same rate. Next, since myosin binds to actin the amount of myosin will increasedue to higher concentration of actin, hence the term − k ma ( a c − a ) where k ma is the rate of binding / unbinding of myosin.Defining m c as the unstable equilibrium concentration of m , the last term in the actin equation represents that actin willdepolymerise with higher concentrations of myosin and is subject to a saturation coe ffi cient K , for a . The negation istrue for myosin since myosin is seen to accumulate.Thus we have three connected equations: the solutions to (1b) (actin concentration) and (1c) (myosin concentra-tion) a ff ect the contractile and pressure parts of the force balance equation and the solution to (1a) (displacement)a ff ects the reaction-di ff usion equations through the convection terms and the changing shape of the domain. The linear stability analysis detailed in Appendix A reveals that parameters, in particular ψ and c , can be variedso that particular patterns become unstable and grow. These patterns correspond to eigenfunctions of the Laplacianon the chosen volume.
3. Finite element formulation
We have formulated a very complex non-linear system of partial di ff erential equations. It is not possible toanalytically solve this system, therefore we must turn to numerical methods to produce approximations. Previousversions of this model were represented in a polar coordinate system and solved using finite di ff erences howeverGeorge [15] used a finite element formulation. We proceed similarly and describe our methods below. In particular,3e employ the moving grid finite element method [3, 25, 26, 27] to compute approximate numerical solutions of thecoupled viscoelastic reaction-di ff usion system defined in 3D Cartesian coordinate system.To begin, the force balance is separated into a system of three partial di ff erential equations representing the threespace dimensions. This clarifies the derivation of the weak formulation. Since σ v , σ e , σ c and σ p (as described inSection 2) are all stress tensors we can write them in matrix form. In three dimensions, strain and dilation are givenby (cid:15) ( u ) : =
12 ( ∇ u + ( ∇ u ) T ) = ∂ u ∂ x ( ∂ v ∂ x + ∂ u ∂ y ) ( ∂ u ∂ z + ∂ w ∂ x ) ( ∂ v ∂ x + ∂ u ∂ y ) ∂ v ∂ y ( ∂ v ∂ z + ∂ w ∂ y ) ( ∂ v ∂ x + ∂ u ∂ y ) ( ∂ w ∂ y + ∂ v ∂ z ) ∂ w ∂ z and φ ( u ) : = ∂ u ∂ x + ∂ v ∂ y + ∂ w ∂ z , (4a)respectively. It follows then that we can write the stress tensors in three-dimensional tensor-matrix form: σ v = µ + ∂ ˙ u ∂ x + µ ( ∂ ˙ v ∂ y + ∂ ˙ w ∂ z ) µ ( ∂ ˙ v ∂ x + ∂ ˙ u ∂ y ) µ ( ∂ ˙ w ∂ x + ∂ ˙ u ∂ z ) µ ( ∂ ˙ v ∂ x + ∂ ˙ u ∂ y ) µ ( ∂ ˙ u ∂ x + ∂ ˙ w ∂ z ) + µ + ∂ ˙ v ∂ y µ ( ∂ ˙ v ∂ z + ∂ ˙ w ∂ y ) µ ( ∂ ˙ w ∂ x + ∂ ˙ u ∂ z ) µ ( ∂ ˙ v ∂ z + ∂ ˙ w ∂ y ) µ + ∂ ˙ w ∂ z + µ ( ∂ ˙ u ∂ x + ∂ ˙ v ∂ y ) , σ e = E + ν ∂ u ∂ x + ν (cid:48) φ ( u ) ( ∂ v ∂ x + ∂ u ∂ y ) ( ∂ u ∂ z + ∂ w ∂ x ) ( ∂ v ∂ x + ∂ u ∂ y ) ∂ v ∂ y + ν (cid:48) φ ( u ) ( ∂ v ∂ z + ∂ w ∂ y ) ( ∂ w ∂ x + ∂ u ∂ z ) ( ∂ v ∂ z + ∂ w ∂ y ) ∂ w ∂ z + ν (cid:48) φ ( u ) , σ c = ψ a e − a / a sat + cm ψ a e − a / a sat + cm
00 0 ψ a e − a / a sat + cm , σ p = p + φ (cid:16) + π δ ( l ) tan − a (cid:17) p + φ (cid:16) + π δ ( l ) tan − a (cid:17)
00 0 p + φ (cid:16) + π δ ( l ) tan − a (cid:17) , where µ + = µ + µ and ν (cid:48) = ν/ (1 − ν ). Substituting these expressions into ∇ · ( σ v + σ e + σ c + σ p ) = gives us4hree equations ∂∂ x (cid:32) D ∂ ˙ u ∂ x + D (cid:18) ∂ ˙ v ∂ y + ∂ ˙ w ∂ z (cid:19) + C ∂ u ∂ x + C (cid:18) ∂ v ∂ y + ∂ w ∂ z (cid:19)(cid:33) + ∂∂ y (cid:32) D (cid:18) ∂ ˙ v ∂ x + ∂ ˙ u ∂ y (cid:19) + C (cid:18) ∂ v ∂ x + ∂ u ∂ y (cid:19)(cid:33) + ∂∂ z (cid:32) D (cid:18) ∂ ˙ w ∂ x + ∂ ˙ u ∂ z (cid:19) + C (cid:18) ∂ w ∂ x + ∂ u ∂ z (cid:19)(cid:33) = − ∂ f ∂ x ,∂∂ x (cid:32) D (cid:18) ∂ ˙ v ∂ x + ∂ ˙ u ∂ y (cid:19) + C (cid:18) ∂ v ∂ x + ∂ u ∂ y (cid:19)(cid:33) + ∂∂ y (cid:32) D ∂ ˙ v ∂ y + D (cid:18) ∂ ˙ u ∂ x + ∂ ˙ w ∂ z (cid:19) + C ∂ v ∂ y + C (cid:18) ∂ u ∂ x + ∂ w ∂ z (cid:19)(cid:33) + ∂∂ z (cid:32) D (cid:18) ∂ ˙ w ∂ y + ∂ ˙ w ∂ z (cid:19) + C (cid:18) ∂ w ∂ z + ∂ v ∂ z (cid:19)(cid:33) = − ∂ f ∂ y ,∂∂ x (cid:32) D (cid:18) ∂ ˙ w ∂ x + ∂ ˙ u ∂ z (cid:19) + C (cid:18) ∂ w ∂ x + ∂ u ∂ z (cid:19)(cid:33) + ∂∂ y (cid:32) D (cid:18) ∂ ˙ v ∂ z + ∂ ˙ w ∂ y (cid:19) + C (cid:18) ∂ v ∂ z + ∂ w ∂ y (cid:19)(cid:33) + ∂∂ z (cid:32) D ∂ ˙ w ∂ z + D (cid:18) ∂ ˙ v ∂ y + ∂ ˙ u ∂ x (cid:19) + C ∂ w ∂ z + C (cid:18) ∂ u ∂ x + ∂ v ∂ y (cid:19)(cid:33) = − ∂ f ∂ z , where f = f = f = (cid:34) p + φ (cid:32) + π δ ( l ) arctan a (cid:33) + ψ a e − a / a sat + cm (cid:35) , (7a) D = µ + µ , D = µ , D = µ , (7b) C = E (1 − ν )(1 + ν )(1 − ν ) , C = E ν (1 + ν )(1 − ν ) and C = E + ν ) . (7c) To find the weak formulation, we take the usual route and multiply by a test function ˆ φ ( x , t ) ∈ H ( Ω t ), where H ( Ω t ) is a Hilbert space, and integrate over the domain. This takes into account Green’s formula and the boundaryconditions. The boundary condition σ v · n = σ e · n = u ( x , t ), v ( x , t ) and w ( x , t ) ∈ H ( Ω t ), t ∈ I such that (cid:90) Ω t ∂ ˆ φ∂ x (cid:32) D ∂ ˙ u ∂ x + D (cid:18) ∂ ˙ v ∂ y + ∂ ˙ w ∂ z (cid:19) + C ∂ u ∂ x + C (cid:18) ∂ v ∂ y + ∂ w ∂ z (cid:19)(cid:33) + ∂ ˆ φ∂ y (cid:32) D (cid:18) ∂ ˙ v ∂ x + ∂ ˙ u ∂ y (cid:19) + C (cid:18) ∂ v ∂ x + ∂ u ∂ y (cid:19)(cid:33) + ∂ ˆ φ∂ z (cid:32) D (cid:18) ∂ ˙ w ∂ x + ∂ ˙ u ∂ z (cid:19) + C (cid:18) ∂ w ∂ x + ∂ u ∂ z (cid:19)(cid:33) d Ω t = − (cid:90) Ω t ∂ ˆ φ∂ x f d Ω t + (cid:90) ∂ Ω t ˆ φ f n ds , (cid:90) Ω t ∂ ˆ φ∂ x (cid:32) D (cid:18) ∂ ˙ v ∂ x + ∂ ˙ u ∂ y (cid:19) + C (cid:18) ∂ v ∂ x + ∂ u ∂ y (cid:19)(cid:33) + ∂ ˆ φ∂ y (cid:32) D ∂ ˙ v ∂ y + D (cid:18) ∂ ˙ u ∂ x + ∂ ˙ w ∂ z (cid:19) + C ∂ v ∂ y + C (cid:18) ∂ u ∂ x + ∂ w ∂ z (cid:19)(cid:33) + ∂ ˆ φ∂ z (cid:32) D (cid:18) ∂ ˙ w ∂ y + ∂ ˙ w ∂ z (cid:19) + C (cid:18) ∂ w ∂ z + ∂ v ∂ z (cid:19)(cid:33) d Ω t = − (cid:90) Ω t ∂ ˆ φ∂ y f d Ω t + (cid:90) ∂ Ω t ˆ φ f n ds , (cid:90) Ω t ∂ ˆ φ∂ x (cid:32) D (cid:18) ∂ ˙ w ∂ x + ∂ ˙ u ∂ z (cid:19) + C (cid:18) ∂ w ∂ x + ∂ u ∂ z (cid:19)(cid:33) + ∂ ˆ φ∂ y (cid:32) D (cid:18) ∂ ˙ v ∂ z + ∂ ˙ w ∂ y (cid:19) + C (cid:18) ∂ v ∂ z + ∂ w ∂ y (cid:19)(cid:33) + ∂ ˆ φ∂ z (cid:32) D ∂ ˙ w ∂ z + D (cid:18) ∂ ˙ v ∂ y + ∂ ˙ u ∂ x (cid:19) + C ∂ w ∂ z + C (cid:18) ∂ u ∂ x + ∂ v ∂ y (cid:19)(cid:33) d Ω t = − (cid:90) Ω t ∂ ˆ φ∂ z f d Ω t + (cid:90) ∂ Ω t ˆ φ f n ds . ∂ f ∂ x , ∂ f ∂ y and ∂ f ∂ z are di ffi cult to evaluate, we have used identities derived from the gradient theorem to write theweak form as above. In other words we have used the identity (cid:90) Ω t ∂ f j ∂ x ˆ φ d Ω t = − (cid:90) Ω t ∂ ˆ φ∂ x f j d Ω t + (cid:90) ∂ Ω t ˆ φ f j n j ds , (9)for j = , ,
3, where x can also be substituted by y and z . n , n , n are the direction cosines of the outward unitvector n normal to ∂ Ω t .Next we want to find the weak formulation of the reaction-di ff usion equations which are given as ∂ a ∂ t + ∇ · ( a β ) − D a ∆ a = f ( a , m ) , ∂ m ∂ t + ∇ · ( m β ) − D m ∆ m = g ( a , m ) . (10)We apply the product rule and convert to the material derivative (defined as DaDt = ∂ a ∂ t + a ( ∇ · β ), in [37]). This gives DaDt − D a ∆ a + a ( ∇ · β ) = f ( a , m ) , DmDt − D m ∆ m + m ( ∇ · β ) = g ( a , m ) . Now continuing as with the force balance equation, we multiply by a test function ˆ ψ ( x , t ) ∈ H ( Ω t ) and integrate overthe domain. The terms ( D a ∆ a ) ˆ ψ and ( D m ∆ m ) ˆ ψ can be simplified using the divergence theorem and for the remainingpart of the left hand side we can use the Reynolds transport theorem. This means the weak formulation can be writtenas: Find a ( x , t ) , m ( x , t ) ∈ H ( Ω t ) , t ∈ I such that ∂∂ t (cid:90) Ω t a ˆ ψ d Ω t + (cid:90) Ω t (cid:16) D a ∇ a · ∇ ˆ ψ (cid:17) d Ω t = (cid:90) Ω t (cid:18) f ( a , m ) ˆ ψ + a D ˆ ψ Dt (cid:19) d Ω t , (12a) ∂∂ t (cid:90) Ω t m ˆ ψ d Ω t + (cid:90) Ω t (cid:16) D m ∇ m · ∇ ˆ ψ (cid:17) d Ω t = (cid:90) Ω t (cid:18) g ( a , m ) ˆ ψ + m D ˆ ψ Dt (cid:19) d Ω t , (12b)for all ˆ ψ ( x , t ) ∈ H ( Ω t ). We now wish to define the problem at discrete points in space. To do this, we define the computational domain Ω h , t as a polyhedral approximation to Ω t , T h , t the discretisation of Ω h , t made up of non-degenerate elements κ i and thefinite element space V h ( t ) : = { v h ∈ C ( Ω t ) : v h | κ is linear } . Thus the space-discrete problem is to find u h ( x , t ), v h ( x , t ), w h ( x , t ), a h ( x , t ), m h ( x , t ) ∈ V h ( t ) , t ∈ I , such that (cid:90) Ω h , t ∂ ˆ φ∂ x (cid:32) D ∂ ˙ u h ∂ x + D (cid:18) ∂ ˙ v h ∂ y + ∂ ˙ w h ∂ z (cid:19) + C ∂ u h ∂ x + C (cid:18) ∂ v h ∂ y + ∂ w h ∂ z (cid:19)(cid:33) + ∂ ˆ φ∂ y (cid:32) D (cid:18) ∂ ˙ v h ∂ x + ∂ ˙ u h ∂ y (cid:19) + C (cid:18) ∂ v h ∂ x + ∂ u h ∂ y (cid:19)(cid:33) + ∂ ˆ φ∂ z (cid:32) D (cid:18) ∂ ˙ w h ∂ x + ∂ ˙ u h ∂ z (cid:19) + C (cid:18) ∂ w h ∂ x + ∂ u h ∂ z (cid:19)(cid:33) d Ω t = − (cid:90) Ω h , t ∂ ˆ φ∂ x f d Ω h , t + (cid:90) ∂ Ω h , t ˆ φ f n ds , (cid:90) Ω h , t ∂ ˆ φ∂ x (cid:32) D (cid:18) ∂ ˙ v h ∂ x + ∂ ˙ u h ∂ y (cid:19) + C (cid:18) ∂ v h ∂ x + ∂ u h ∂ y (cid:19)(cid:33) + ∂ ˆ φ∂ y (cid:32) D ∂ ˙ v h ∂ y + D (cid:18) ∂ ˙ u h ∂ x + ∂ ˙ w h ∂ z (cid:19) + C ∂ v h ∂ y + C (cid:18) ∂ u h ∂ x + ∂ w h ∂ z (cid:19)(cid:33) + ∂ ˆ φ∂ z (cid:32) D (cid:18) ∂ ˙ w h ∂ y + ∂ ˙ w h ∂ z (cid:19) + C (cid:18) ∂ w h ∂ z + ∂ v h ∂ z (cid:19)(cid:33) d Ω h , t = − (cid:90) Ω h , t ∂ ˆ φ∂ y f d Ω h , t + (cid:90) ∂ Ω h , t ˆ φ f n ds , (cid:90) Ω h , t ∂ ˆ φ∂ x (cid:32) D (cid:18) ∂ ˙ w h ∂ x + ∂ ˙ u h ∂ z (cid:19) + C (cid:18) ∂ w h ∂ x + ∂ u h ∂ z (cid:19)(cid:33) + ∂ ˆ φ∂ y (cid:32) D (cid:18) ∂ ˙ v h ∂ z + ∂ ˙ w h ∂ y (cid:19) + C (cid:18) ∂ v h ∂ z + ∂ w h ∂ y (cid:19)(cid:33) + ∂ ˆ φ∂ z (cid:32) D ∂ ˙ w h ∂ z + D (cid:18) ∂ ˙ v h ∂ y + ∂ ˙ u h ∂ x (cid:19) + C ∂ w h ∂ z + C (cid:18) ∂ u h ∂ x + ∂ v h ∂ y (cid:19)(cid:33) d Ω h , t = − (cid:90) Ω h , t ∂ ˆ φ∂ z f d Ω h , t + (cid:90) ∂ Ω h , t ˆ φ f n ds , ∂∂ t (cid:90) Ω h , t a h ˆ ψ d Ω h , t + (cid:90) Ω h , t D a ∇ a h ·∇ ˆ ψ d Ω h , t = (cid:90) Ω h , t (cid:32) I h f ( a h , m h ) ˆ ψ + a h D ˆ ψ Dt (cid:33) d Ω h , t ,∂∂ t (cid:90) Ω h , t m h ˆ ψ d Ω h , t + (cid:90) Ω h , t D m ∇ m h ·∇ ˆ ψ d Ω h , t = (cid:90) Ω h , t (cid:32) I h g ( a h , m h ) ˆ ψ + m h D ˆ ψ Dt (cid:33) d Ω h , t , for all ˆ φ, ˆ ψ ∈ V h ( t ). We can then express u h , v h , w h , a h and m h in terms of the linear basis functions: u h = nde (cid:88) j = U j φ j , v h = nde (cid:88) j = V j φ j , w h = nde (cid:88) j = W j φ j , a h = nde (cid:88) j = α j φ j and m h = nde (cid:88) j = µ j φ j . (15)This means that we are left with equations which contain only simple functions and their derivatives and point valuesfor the variables. Substituting u h , v h , and w h into the above equations and using the Galerkin formulation, take the testfunctions to belong the same spaces as the nodal basis functions. Hence, the force balance equations can be writtenin block matrix-vector form A A A [ A ] T A A [ A ] T [ A ] T A d U dtd V dtd W dt + B B B [ B ] T B B [ B ] T [ B ] T B UVW = F F F , (16)where { U ( t ) } = ( U , ..., U nde ), { V ( t ) } = ( V , ..., V nde ) , { W ( t ) } = ( W , ..., W nde ) and: A i j ( t ) : = (cid:90) Ω h , t D ∂φ i ∂ x ∂φ j ∂ x + D (cid:32) ∂φ i ∂ y ∂φ j ∂ y + ∂φ i ∂ z ∂φ j ∂ z (cid:33) d Ω h , t , A i j ( t ) : = (cid:90) Ω h , t D (cid:32) ∂φ i ∂ x ∂φ j ∂ x + ∂φ i ∂ z ∂φ j ∂ z (cid:33) + D ∂φ i ∂ y ∂φ j ∂ y d Ω h , t , A i j ( t ) : = (cid:90) Ω h , t D (cid:32) ∂φ i ∂ x ∂φ j ∂ x + ∂φ i ∂ y ∂φ j ∂ y (cid:33) + D ∂φ i ∂ z ∂φ j ∂ z d Ω h , t , B i j ( t ) : = (cid:90) Ω h , t C ∂φ i ∂ x ∂φ j ∂ x + C (cid:32) ∂φ i ∂ y ∂φ j ∂ y + ∂φ i ∂ z ∂φ j ∂ z (cid:33) d Ω h , t , B i j ( t ) : = (cid:90) Ω h , t C (cid:32) ∂φ i ∂ x ∂φ j ∂ x + ∂φ i ∂ z ∂φ j ∂ z (cid:33) + C ∂φ i ∂ y ∂φ j ∂ y d Ω h , t , B i j ( t ) : = (cid:90) Ω h , t C (cid:32) ∂φ i ∂ x ∂φ j ∂ x + ∂φ i ∂ y ∂φ j ∂ y (cid:33) + C ∂φ i ∂ z ∂φ j ∂ z d Ω h , t , A i j ( t ) : = (cid:90) Ω h , t D ∂φ i ∂ x ∂φ j ∂ y + D ∂φ i ∂ y ∂φ j ∂ x d Ω h , t , A i j ( t ) : = (cid:90) Ω h , t D ∂φ i ∂ x ∂φ j ∂ z + D ∂φ i ∂ z ∂φ j ∂ x d Ω h , t , A i j ( t ) : = (cid:90) Ω h , t D ∂φ i ∂ y ∂φ j ∂ z + D ∂φ i ∂ z ∂φ j ∂ y d Ω h , t , B i j ( t ) : = (cid:90) Ω h , t C ∂φ i ∂ x ∂φ j ∂ y + C ∂φ i ∂ y ∂φ j ∂ x d Ω h , t , B i j ( t ) : = (cid:90) Ω h , t C ∂φ i ∂ x ∂φ j ∂ z + C ∂φ i ∂ z ∂φ j ∂ x d Ω h , t , B i j ( t ) : = (cid:90) Ω h , t C ∂φ i ∂ y ∂φ j ∂ z + C ∂φ i ∂ z ∂φ j ∂ y d Ω h , t , F j ( t ) : = − (cid:90) Ω h , t f ∂φ j ∂ x d Ω h , t + (cid:90) ∂ Ω h , t n f φ j ds , F j ( t ) : = − (cid:90) Ω h , t f ∂φ j ∂ y d Ω h , t + (cid:90) ∂ Ω h , t n f φ j ds , F j ( t ) : = − (cid:90) Ω h , t f ∂φ j ∂ z d Ω h , t + (cid:90) ∂ Ω h , t n f φ j ds . It must be noted that in the above we have used the fact that ˆ φ ∈ V h and ˆ ψ ∈ V h where V h = φ , ...., φ n . For the7ake of ease of notation and computation we we define the following block matrices and vectors[ A ] : = A A A [ A ] T A A [ A ] T [ A ] T A , [ B ] : = B B B [ B ] T B B [ B ] T [ B ] T B , { U } : = UVW and { F } : = F F F . (18a)Therefore the force balance equation’s semi-discrete finite element formation can be written compactly as[ A ] { d U } dt + [ B ] { U } = { F } . (19)Now considering the reaction kinetics to be as in Eqs (3), (in Section 2), we can similarly write the reaction-di ff usionequations in semi-discrete form ∂∂ t ( M α ) + D a K α = k a a c H − k a M α + k am M α (1 − µ )1 + K α , (20a) ∂∂ t ( M µ ) + D m K µ = − k ma a c H + k ma M α − k am M α (1 − µ )1 + K α , (20b)where vector operations in α (1 − µ ) / (1 + K α ) are pointwise and M i , j = (cid:90) Ω h , t φ i φ j , K i , j = (cid:90) Ω h , t ∇ φ i · ∇ φ j and H j = (cid:90) Ω h , t φ j . (21)To compute these integrals we use Gauss numerical quadrature [36]. This is done as follows. First we can considerthe integrals elementwise, M i , j = (cid:88) ∆ k (cid:90) ∆ k φ i φ j , A i , j = (cid:88) ∆ k (cid:90) ∆ k ∇ φ i · ∇ φ j and H j = (cid:88) ∆ k (cid:90) ∆ k φ j . (22)Then choose a numerical quadrature comprising a set of points and weights depending on the functions to be inte-grated. This can be written as a formula for the integral of a function ξ (cid:90) ξ ( x ) ≈ (cid:88) q ξ (¯ x q ) w q , (23)where ¯ x q and w q are the q th quadrature points and weights respectively. Therefore the integrals can be approximatedby M i , j ≈ (cid:88) ∆ k (cid:88) q φ i (¯ x q ) φ j (¯ x q ) w q , A i , j ≈ (cid:88) ∆ k (cid:88) q ∇ φ i (¯ x q ) · ∇ φ j (¯ x q ) w q and H j ≈ (cid:88) ∆ k (cid:88) q φ j (¯ x q ) w q . (24)These are all computed in the deal.II implementation [4]. Next we carry out the temporal discretisation of the system of ordinary di ff erential equations arising from thefinite element discretisation. To proceed, we split the interval into a finite number of sub-intervals [ t n , t n + ] use auniform time step ∆ t : = t n + − t n . We can then use a modified implicit-explicit (IMEX) finite di ff erentiation formula[22, 25, 39]. Thus the fully discrete problem is now([ A ] n + ∆ t [ B ] n ) U n + = [ A ] n { U } n + ∆ t { F } n , (25a) (cid:104) M n + + ∆ tD a K n + (cid:105) α n + = M n α n +∆ t ( k a ( a c H n − M n α n ) + k am M n α n (1 − µ n )1 + K ( α n ) ) , (25b) (cid:2) M n + + ∆ tD m K n + (cid:3) µ n + = M n µ n + ∆ t ( − k ma ( a c H n − M n α n ) − k am M n α n (1 − µ n )1 + K ( α n ) ) , (25c)8here the superscripts n and n + are the computed values on the mesh at times t n and t n + respectively. Note that wehave treated some parts implicitly (e.g. di ff usion) and other parts fully explicit (e.g. reactions).Hence we have three equations all with the same form. At each time-step we assemble the matrices to obtain asystem of linear algebraic equations. When solving (25a) we see that the block matrix on the left hand side is notsymmetric therefore we use the most e ff ective solver for this which is GMRES [40]. The equations (25b) and (25c)are solved using the conjugate gradient method [20]. The displacement of the nodes of the mesh is chosen to be equal to the flow velocity therefore β : = ∂ U ∂ t . Since t n + = t n + ∆ t and x ( t n ) ∈ Ω t n , x ( t n + ) ∈ Ω t n + be points in the respective domains. We can define a first order linearapproximation as: β ( x , t n ) = x ( t n + ) − x ( t n ) ∆ t . (26)This means we can define a new approximation to the domain Ω t n + such that x ( t n + ) = x ( t n ) + ∆ t ∂ U ∂ t = x ( t n ) + ( U n + − U n ) . (27)At each step we have a new mesh with new shape functions so we must assemble new matrices M n , H n , A n , B n , F n toiteratively solve the discrete coupled problem as outlined in the following algorithm The fully discrete problem is solved iteratively with the following algorithm: • Initialise U , α , µ and fixed parameters • WHILE t < endtime – Assemble M n , H n , A n , B n , F n – Solve for U n + using (25a) – Compute the new domain using U n + – Solve for α n + and µ n + using (25b) and (25c) – t = t + ∆ t • ENDWe create a mesh using Gmsh [17] and implement this algorithm using deal.II [4], a C ++ software library whichprovides tools to solve partial di ff erential equations which are discretised with finite element methods. Unlike themajority of other finite element software, deal.II uses hexahedral and quadrilateral elements rather than triangles.
4. Numerical simulations
We now present simulations of our model. We want to see the organisation of the molecular species into regionswhich will cause the cell to move. This organisation may be caused by di ff usion-driven instability, or due to themovement of the cell combined with the reaction-di ff usion equations. The linear stability analysis of Appendix Aholds true close to critical bifurcation points, these include parameters as well as the geometrical deformation of thecell. The conditions for stability are numerous and complex, however it is still possible to choose parameters so thatparticular modes can be selected. When we consider longer time, and therefore far away from equilibrium, linearstability theory no longer holds but we see significant protrusions and contractions which deform the mesh into manydi ff erent shapes. Parameters used are in Table 1. 9 .1. Excitation of mode w , In this example, the actin and myosin concentration solutions will be the negation of each other with actin con-centration highest on one side and myosin concentration highest on the opposite side. This mode is the first eigen-function that one might hope to see for the organisation of actin and myosin in a cell because it is similar to whatis often observed in a moving cell [12, 31]. k , = . ψ = , c = − , k a = . , k ma = .
05 and k am = .
06 and initial conditions a ( x , = + w , ( x ) × ran , m ( x , = − w , ( x ) × ran , we observe that the mode w , is selected for actin and myosin. In Figure 1 we plot the concentrations of actin andmyosin at time t =
1. Blue indicates where the concentration is low, while red indicates that concentration is high. Inthis case very little deformation is seen. , is excited initially The simple first mode is not the only organisation which makes sense or shows similarities to organisation seen incells. The cell can protrude in more than one direction because of actin accumulation at both ends, or deform in manyother ways. Also, myosin could accumulate and ”squeeze” on both sides. Therefore we continue by isolating othermodes. We see that both the parameters, and the initial conditions, have an e ff ect on which mode will grow. The firstlarge deformation is seen when choosing initial conditions a ( x , = + w , ( x ) × ran , m ( x , = − w , ( x ) × ran . In Figure 3 we plot the concentrations of actin and myosin and the displacement (cid:16) | U | = √ U + V + W (cid:17) at eachpoint. The cell expands at the two ends where actin concentration is high and contracts in the middle where myosinconcentration is high. So far the results are visually similar to the the results obtained by [16] in the absence of myosin,one di ff erence however is that there is only a very small volume increase because the cell is contracting in the middleas well as protruding. Other results (not shown) when the excited mode for myosin is the same as the mode for actinare very similar to the previous model [30]. Therefore, we investigate whether more interesting dynamics may occurif we try to excite di ff ering modes for the two concentrations. , and w , are excited for actin and myosin, respectively While the idea that actin and myosin accumulate in opposite sides is quite well founded, their concentrationgradients are rarely exactly opposite. Therefore here we investigate if di ff ering modes can be excited for actin andmyosin. Choosing appropriate initial conditions, to encourage di ff erent modes to grow, we observe more irregulardeformations. In Figure 3 we plot the concentrations of actin and myosin and the displacement when the initialconditions are a ( x , = + w , ( x ) × ran , m ( x , = + w , ( x ) × ran . The cell squeezes where there is high myosin concentration and there is a protrusion in the direction of higher actin,this is also illustrated in Figure 4 where the minimum and maximum in each spatial direction are plotted. , and w , are excited for actin and myosin, respectively In Figure 5, the initial conditions a ( x , = + w , × ran , m ( x , = + w , ( x ) × ran , contain the same two eigenfunctions as the last example but with di ff erent orientations, ( w , is a rotation of w , ),we observe a quite di ff erent deformation. There is high actin concentration at the top and bottom of the sphere.Without the e ff ect of myosin one would expect it to extend in both directions in the same way as in Section 4.2,however there is high myosin at the bottom so the cell only protrudes upwards. Then at t = z -direction.10 .5. Cell deformation when w , and w , are excited for actin and myosin, respectively Next, we begin with initial conditions a ( x , = + w , ( x ) × ran , m ( x , = + w , ( x ) × ran . This leads to a protrusion in the area with highest actin concentration which is pulling the cell in the negative z -direction. At the same time there is inward movement in areas of high myosin concentration. The cell has translatedin the negative z -direction and this is plotted in Figure 8, and the change in volume is illustrated in Figure 9. , and w , are excited for actin and myosin, respectively In another example of mixed modes, we start with a ( x , = + w , ( x ) × ran , m ( x , = + w , ( x ) × ran . This leads to the expansion shown in Figure 10. The cell contracts inwards at areas of high myosin concentration andprotrudes in the remaining areas, there are large protrusions in two opposing directions, the largest being the directionwhere actin was initially highest, subsequently, actin concentrates in areas of high curvature and protrudes further.Section 4.2 4.3 4.4 4.5 4.6Figure 2 3 5 8 10 ψ
200 20 150 100 100 c -40 -80 -40 -80 -100 k a k ma k am Table 1: Parameters for simulations in this section.Figure 1: Graphical displays of the actin and myosin concentrations at time t =
1. These are numerical solutions to the full system (1) using thefinite element formulation, as described in Section 4.1. a) t = = = t , for the conditions described inSection 4.2. There is high actin at two ends, and high myosin in the middle. We then see in (c) that the cell squeezes in the middle stretches in thetwo directions of higher actin. a) t = = = = t , for the conditions described inSection 4.3. The sphere is squeezed where there is high myosin and then there is a protrusion in the area of high actin. Displacements in the x , y and z directions are shown in Figure 4. igure 4: Plot to show the minimum and maximum of x (red), y (blue) and z (green) for the example in Section 4.3 and Figure 3. The cell iscontracting in the y direction, expanding slightly in the x direction but significantly in the positive z -direction. a) t = = = = = a) t = =
20 (c) t =
29 (d) t = z -direction (b) Translation the centre of the cell in z -directionFigure 7: Plots to illustrate how the cell expands, contracts and translates in Figure 3. (Example 4.4). a) t = = = z -direction. The cell is also being squeezed inthe x - and y -direction so we do not observe a significant volume increase. (Example 4.5 and Figure 8). a) t = = = norms In Figure 11 we plot the norm of di ff erences between successive solutions in the case of the full system examplein Section 4.4. We see an increase, or decrease, in the L norm when the rate of deformation is accelerating, or decel-erating, respectively. The qualitative changes in the L norms are similar, but the changes in myosin and displacementappear slightly later than actin. This may suggest, in this example, that the change in actin triggers the change in theother variables.In most of the numerical solutions in Section 4 the cell becomes deformed in such a way that means the meshbecomes highly non-conforming and therefore the numerical method breaks down. In this case, re-meshing mightbe useful in order to allow the mesh to adapt to large deformations giving rise to a more stable numerical method.Equally important, adding a volume constraint term to the model might also prevent large expansions or contractions.These observations will form the subject of future studies.18 igure 11: Plot of the L norm of di ff erence between successive solutions for the example shown in Figure 5. There is an initial decrease due todi ff usion, increases when the deformation is accelerating and decreases as deformation decelerates. The large error at the end is due to the meshbecoming very deformed and breaking.
5. Conclusion
Our model revolves around an equation which balances elasticity, viscosity, contractility and pressure. Connectedto this are two equations for the concentrations of F-actin and bound myosin. Unlike the previous study which ourmodel is inspired by, (that of George [16]), we used the software library deal.II [4] to implement the finite elementmethod. The key di ff erence between this software and the previously used ALBERTA [42] is that the elementsare hexahedra rather than tetrahedra. The implementation is therefore di ff erent but we were able to appropriatelyreplicate the previous results by George [16] of cytomechanical model on a unit two dimensional disk (results notshown). Once this new implementation was verified we extended the model substantially in two ways: The two-dimensional formulation was extended to three dimensions. Unlike previous studies of this modelling framework, forthe first time, we considered a second reaction-di ff usion model to describe how myosin interacts with actin and howit contributes to cell contraction during cell migration. In the absence of experimental observations, we postulatedhypothetical reaction kinetics describing the interaction between actin and myosin. As a first step in understandingsolution behaviour in three dimensions of the full model, linear stability analysis close to bifurcation points was carriedout and appropriate key parameter values were identified. An evolving finite element method was implemented inmulti-dimensions.Our numerical simulations showed that the model extends into three dimensions and the addition of myosin allowssome symmetries to be broken and more striking deformations to emerge. In summary the main observations are: • There is outward movement in areas with high actin concentration, also where there is higher curvature, higheractin concentration is observed. • This outward movement due to actin concentration is halted in areas with high myosin concentration. Addition-ally, if there is low actin and high myosin, we can see negative curvature. • Identifying bifurcation parameters is more complicated than in previous models but e ff ects of parameters canstill be seen. The contractility due to myosin, c , strongly e ff ects the speed of the deformation while the reactionconstants k a , k ma and k am , and the di ff usion coe ffi cients D a and D m play a part in which mode the actin andmyosin concentrations will arrange into. It is not as possible to isolate single modes just by picking parametervalues, however, choosing initial conditions as a random number multiplied by the eigenfunction means it ispossible to choose modes. • Several examples show small translations. • Initial conditions are a highly significant factor for the progression of the solutions. • In most examples the volume is increasing slightly but not significantly. Thus a mechanism for volume conser-vation or constraint could be a useful extension. 19he modelling and computational framework presented in this article can readily be adapted to consider newexperimentally driven reaction kinetics between actin and myosin or interactions between three or more molecularspecies, for example, studies using actin, myosin, GEF, Rho, Rac and CDC42 could be productive [19, 21, 33, 43].Other useful extensions of the model could use or formulate re-meshing strategies.
Appendix A. Linear stability analysis
In the previous model, the reaction-di ff usion equation alone could not cause patterning [16]. Without the flowterm, the prescribed reactions meant the concentration of actin would always return to the homogeneous steady stateof a = a c . In our case we have two coupled reaction-di ff usion equations which are well known to induce patterning incertain cases.We perform non-dimensionalisation to reduce parameters and simplify calculations. It also allows the reaction-di ff usion equations to take the form necessary to use the standard conditions for di ff usion driven instability. This isinvestigated in detail in [30]. We choose the nondimensionalised parameters:˜ t = L D a t , ˜ a = aa c = a , ˜ m = mm c = m , d = D m D a , ˜ K = Ka c , γ = L k a D a , (A.1a)˜ k ma = k ma k a , ˜ k am = k am k a , ˜ k m = k m k a , ˜ ∆ = L ∆ , ˜ ∇ = L ∇ , ˜ u = u L , ˜ φ = φ, (A.1b)˜ ε = ε , ˜ p = p + ν E , ˜ β = β LD a , ˜ a sat = a sat a c , ˜ ψ = ψ a c + ν E ˜ µ i = µ i D a (1 + ν ) EL . (A.1c)In the above, L is the typical radius of a cell. Substituting appropriately and carrying out algebraic manipulationsleads to the following non-dimensionalised system (for general kinetics)˜ ∇ · (cid:34) ( ˜ µ ˜ ε t + µ ˜ φ t I ) + (˜ ε + ν (cid:48) ˜ φ I ) + σ (˜ a ) I + ˜ c ˜ m I + ˜ p (˜ a )1 + ˜ φ I (cid:35) = ∂ ˜ a ∂ ˜ t + ˜ ∇ · (˜ a ˜ β ) − ˜ ∆ ˜ a − γ f (˜ a , ˜ m ) = ∂ ˜ m ∂ ˜ t + ˜ ∇ · ( ˜ m ˜ β ) − d ˜ ∆ ˜ m − γ g (˜ a , ˜ m ) = . (A.2c)where σ (˜ a ) = ˜ ψ ˜ a e − ˜ a / ˜ a sat , ν (cid:48) = ν − ν and ˜ p (˜ a ) = ˜ p (cid:0) + π δ ( l ) arctan ˜ a (cid:1) . f and g have been nondimensionalised and wechoose only functions such that system has a steady state at ( a s , m s , u s ) = (1 , , ). Given small variations ˆ a , ˆ m and ˆ u ,consider the perturbation from the steady state ˜ a = a s + ˆ a = + ˆ a , ˜ m = m s + ˆ m = + ˆ m , ˜ u = u s + ˆ u = ˆ u . This resultsin the linear system ˜ ∇ · (cid:34) ( ˜ µ ˆ ε t + µ ˆ φ t I ) + (ˆ ε + ν (cid:48) ˆ φ I ) + ˆ a σ (cid:48) (1) I + c ˆ m I + ˜ p (1 − ˆ φ ) I + ˜ p π δ ( l )ˆ a I (cid:35) = , (A.3a) ∂ ˆ a ∂ ˜ t + ˜ ∇ · ( ˆ β ) − ˜ ∆ ˆ a − f a ˆ a − f m ˆ m = , (A.3b) ∂ ˆ m ∂ ˜ t + ˜ ∇ · ( ˆ β ) − d ˜ ∆ ˆ m − g a ˆ a − g m ˆ m = . (A.3c)We now look for solutions of the formˆ a ( x , t ) = a ∗ e λ t + i k · x , ˆ m ( x , t ) = m ∗ e λ t + i k · x and ˆ u ( x , t ) = u ∗ e λ t + i k · x , (A.4)where λ is the growth rate, k is the wave vector, and a ∗ , m ∗ and u ∗ are small amplitudes. We require solutions to be20 igure A.12: Analytical solutions to the eigenvalue problem on the unit sphere i.e. (A.6) for selected values of l , m , n . For l ≥ non-trivial and so we obtain the stability matrix (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) λ + k − γ f a − γ f m λ ik − γ g a λ + dk − γ g m λ ik − ik σ (cid:48) (1) − ik ˜ p π δ ( l ) − cik ˜ µ k λ + k (1 + ν (cid:48) ) − ˜ pk (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = , where k = | k | (A.5a) = ⇒ ( h ( λ ) : = ) µλ + a ( k ) λ + b ( k ) λ + c ( k ) = , (A.5b)where a ( k ) = k (1 + d ) − γ ( f a + g m ) + + ν (cid:48) − p − c − ( σ (cid:48) (1) + ˜ p π δ ( l )) , (A.5c) b ( k ) = ˜ µ ( k − γ f a )( dk − γ g m ) + (1 + ν (cid:48) + p )( k (1 + d ) − γ ( f a + g m )) (A.5d) − c ( k + γ ( − f a + g a )) + ( σ (cid:48) (1) + ˜ p π δ ( l ))( γ ( f m + g m ) − dk ) − γ ˜ µ f m g a , (A.5e)and c ( k ) = (1 + ν (cid:48) + p ) (cid:16) ( k − γ f a )( dk − γ g m ) − γ f m g a (cid:17) . (A.5f)Thus h ( λ ) = λ . There will be instabilitywhen Re( λ ) >
0. We exploit this relation to isolate particular patterns / modes. The unstable modes will correspond tothe eigenfunctions of the Laplacian on the sphere and k the associated eigenvalues. Appendix A.1. Eigenvalues and eigenfunctions of the Laplacian in the bulk of the unit sphere
The eigenvalues on the unit sphere Ω = { ( x , y , z ) : x + y + z ≤ } (with homogeneous Neumann boundarycondition) are well known and are obtained using separation of variables [2, 29]. There are an infinite number ofdiscrete solutions of the form w ml , n ( r , θ, φ ) = A ml , n J l + ( j (cid:48) l + , n r ) e im φ P ml (cos θ ) , (A.6)where l , m , n are all integers such that | m | ≤ l ≤ n , A ml , n are constants, J α ( x ) is a Bessel function of the first kind, i.e. J α ( x ) = (cid:80) ∞ j = − j j ! Γ (1 + j + α ) (cid:16) x (cid:17) j + α with Γ ( n ) = ( n − P ml ( x ) are associated Legendre polynomials and j (cid:48) l + , n are zerosof the di ff erential of the spherical Bessel function. We can find the eigenvalues k l , n = ( j (cid:48) l + , n ) numerically. It follows21 a) Plot to show maximum real (solid line) and imaginary (dotted line)parts of the solution to the dispersal relation. (b) Plot to show maximum real part of λ as ψ is varied.Figure A.13: The e ff ects of varying parameters on behaviour of solutions. that for each eigenvalue λ l , n = k l , n there are 2 l + l , m and n . The wave numbers are discrete. Appendix A.2. Parameter selection
The conditions on the positivity of the roots of (A.5b) are numerous and the coe ffi cients of the polynomial areburdensome. Therefore, we numerically find these roots and observe the real and imaginary parts.We found that contractility due to myosin ( c ), and due to actin ( ψ ) are particularly significant for finding unstablewavenumbers. In Figure A.13a we plot the real and imaginary parts of the solution against k for three di ff erent valuesof c. We can see that when wavenumbers k are less than ∼ .
5, then Re( λ ) > c , therefore thewavenumbers will be unstable. Additionally these wavenumbers will be oscillatory for c =
10. There are also regionsof Hopf instability and oscillatory instability for all the three values of c . In Figure A.13b we fix other parameters andvary ψ to see that, just like in [16] higher values of ψ mean higher wavenumbers can be excited. Acknowledgements
This work (LM) was supported by an EPSRC Doctoral Training Centre Studentship through the University ofSussex. AM acknowledges support from the Leverhulme Trust Research Project Grant (RPG-2014-149). This work(AM) has received funding from the European Union Horizon 2020 research and innovation programme under theMarie Sklodowska-Curie grant agreement (No 642866). AM is a Royal Society Wolfson Research Merit AwardHolder, generously supported by the Wolfson Foundation. LM acknowledges the support from the University ofSussex ITS for computational purposes.
References [1] Wolfgang Alt and Robert T Tranquillo. Basic morphogenetic system modeling shape changes of migrating cells: How to explain fluctuatinglamellipodial dynamics.
Journal of Biological Systems , 3(04):905–916, 1995.[2] G Arfken, H Weber, and F Harris.
Mathematical methods for physicists . Elsevier, Amsterdam, 2013.[3] M J Baines.
Moving finite elements, Monographs on Numerical Analysis . Clarendon Press, Oxford, 1994.[4] W Bangerth, T Heister, L Heltai, G Kanschat, M Kronbichler, M Maier, and B Turcksin. The deal.ii library, version 8.3.
Archive of NumericalSoftware , 4(100):1–11, 2016. ISSN 2197-8263. doi: 10.11588 / ans.2016.100.23122.[5] Konstantinos N Blazakis. Computational methods for investigating cell motility with applications to neutrophil cell migration . PhD thesis,University of Sussex, 2015.[6] Dean Bottino, Alexander Mogilner, Tom Roberts, Murray Stewart, and George Oster. How nematode sperm crawl.
Journal of cell science ,115(2):367–384, 2002.
7] Dean C. Bottino and Lisa J. Fauci. A computational model of ameboid deformation and locomotion.
European Biophysics Journal , 27(5):532–539, Aug 1998. doi: 10.1007 / s002490050163. URL https://doi.org/10.1007/s002490050163 .[8] Dennis Bray. Cell movements: from molecules to motility . Garland Science, 2001.[9] Volker Brinkmann, Ulrike Reichard, Christian Goosmann, Beatrix Fauler, Yvonne Uhlemann, David S Weiss, Yvette Weinrauch, and ArturoZychlinsky. Neutrophil extracellular traps kill bacteria. science , 303(5663):1532–1535, 2004.[10] J Chen, J Irianto, S Inamdar, P Pravincumar, DA Lee, Dan L Bader, and MM Knight. Cell mechanics, structure, and function are regulatedby the sti ff ness of the three-dimensional microenvironment. Biophysical journal , 103(6):1188–1197, 2012.[11] John Condeelis and Je ff rey E Segall. Intravital imaging of cell movement in tumours. Nature Reviews Cancer , 3(12):921, 2003.[12] GM Cooper.
The Cell: A Molecular Approach. 2nd edition.
Sinauer Associates, Sunderland (MA), 2000.[13] C Elliott, B Stinner, and C Venkataraman. Modelling cell motility and chemotaxis with evolving surface finite elements.
J. Roy. Soc. Interface ,9(76):3027–3044, 2012.[14] Peter Friedl and Darren Gilmour. Collective cell migration in morphogenesis, regeneration and cancer.
Nature reviews Molecular cell biology ,10(7):445–457, 2009.[15] Uduak Z George, Ang´elique St´ephanou, and Anotida Madzvamuse. Mathematical modelling and numerical simulations of actin dynamics inthe eukaryotic cell.
Journal of mathematical biology , pages 1–47, 2013.[16] Uduak Zenas George.
A numerical approach to studying cell dynamics . PhD thesis, University of Sussex, 2012.[17] Christophe Geuzaine and Jean-Franc¸ois Remacle. Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities.
International journal for numerical methods in engineering , 79(11):1309–1331, 2009.[18] E Gladilin, A Micoulet, B Hosseini, K Rohr, J Spatz, and R Eils. 3d finite element analysis of uniaxial cell stretching: from image to insight.
Physical biology , 4(2):104, 2007.[19] Alan Hall. Rho gtpases and the actin cytoskeleton.
Science , 279(5350):509–514, 1998.[20] Magnus Rudolph Hestenes and Eduard Stiefel. Methods of conjugate gradients for solving linear systems.
J Res Nat Bur Stand , 49(6):409–436, 1952.[21] William R Holmes, Adriana E Golding, William M Bement, and Leah Edelstein-Keshet. A mathematical model of gtpase pattern formationduring single-cell wound repair.
Interface focus , 6(5):20160032, 2016.[22] O Lakkis, A Madzvamuse, and C Venkataraman. Implicit–explicit timestepping with finite element approximation of reaction–di ff usionsystems on evolving domains. SIAM Journal on Numerical Analysis , 51(4):2309–2330, 2013.[23] MA Lewis and JD Murray. Analysis of stable two-dimensional patterns in contractile cytogel.
Journal of Nonlinear Science , 1(3):289–311,1991.[24] G MacDonald, John A Mackenzie, M Nolan, and RH Insall. A computational method for the coupled solution of reaction–di ff usion equationson evolving domains and manifolds: Application to a model of cell migration and chemotaxis. Journal of computational physics , 309:207–226, 2016.[25] A Madzvamuse. Time-stepping schemes for moving grid finite elements applied to reaction–di ff usion systems on fixed and growing domains. J Comput Phys , 214(1):239–263, 2006.[26] A Madzvamuse, A Wathen, and P Maini. A moving grid finite element method applied to a model biological pattern generator.
J ComputPhys , 190(2):478–500, 2003.[27] Anotida Madzvamuse and Uduak Zenas George. The moving grid finite element method applied to cell movement and deformation.
FiniteElements in Analysis and Design , 74:76–92, 2013.[28] Angelika Manhart, Dietmar Oelz, Christian Schmeiser, and Nikolaos Sfakianakis. Numerical treatment of the filament-based lamellipodiummodel (fblm). In
Modeling cellular systems , pages 141–159. Springer, 2017.[29] M Morimoto.
Analytic functionals on the sphere . American Mathematical Society, Providence, R.I., 1998.[30] Laura R Murphy.
Mathematical studies of a mechanobiochemical model for 3D cell migration . PhD thesis, University of Sussex, 2018.[31] Michael Murrell, Patrick W Oakes, Martin Lenz, and Margaret L Gardel. Forcing cells into shape: the mechanics of actomyosin contractility.
Nature Reviews Molecular Cell Biology , 16(8):486–498, 2015.[32] Matthew P Neilson, John A Mackenzie, Steven D Webb, and Robert H Insall. Modeling cell movement and chemotaxis using pseudopod-based feedback.
SIAM Journal on Scientific Computing , 33(3):1035–1057, 2011.[33] Catherine D Nobes and Alan Hall. Rho, rac, and cdc42 gtpases regulate the assembly of multimolecular focal complexes associated withactin stress fibers, lamellipodia, and filopodia.
Cell , 81(1):53–62, 1995.[34] George F Oster, James D Murray, and GM Odell. Formation of microvilli. Technical report, Los Alamos National Lab., NM (USA), 1985.[35] C Pozrikidis. Numerical simulation of cell motion in tube flow.
Annals of Biomedical Engineering , 33(2):165–178, 2005.[36] William H Press, Saul A Teukolsky, William T Vetterling, and Brian P Flannery.
Numerical recipes 3rd edition: The art of scientificcomputing . Cambridge university press, 2007.[37] J.N. Reddy.
An introduction to the finite element method.
McGraw-Hill, New York, 1993.[38] Boris Rubinstein, Ken Jacobson, and Alex Mogilner. Multiscale two-dimensional modeling of a motile simple-shaped cell.
MultiscaleModeling & Simulation , 3(2):413–439, 2005.[39] S Ruuth. Implicit-explicit methods for reaction-di ff usion problems in pattern formation. J Math Biol , 34(2):148–176, 1995.[40] Youcef Saad and Martin H Schultz. Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems.
SIAMJournal on scientific and statistical computing , 7(3):856–869, 1986.[41] Y Sakamoto, S Prudhomme, and MH Zaman. Modeling of adhesion, protrusion, and contraction coordination for cell migration simulations.
Journal of mathematical biology , 68(1-2):267–302, 2014.[42] A Schmidt, KG Siebert, D K¨oster, O Kriessl, and CJ Heine. Alberta-an adaptive hierarchical finite element toolbox.
URL http: // , 2007.[43] Cory M Simon, Emily M Vaughan, William M Bement, and Leah Edelstein-Keshet. Pattern formation of rho gtpases in single cell woundhealing. Molecular biology of the cell , 24(3):421–432, 2013.[44] A Stephanou, MAJ Chaplain, and Philippe Tracqui. A mathematical model for the dynamics of large membrane deformations of isolated broblasts. Bulletin of mathematical biology , 66(5):1119, 2004.[45] Wanda Strychalski, David Adalsteinsson, and Timothy C Elston. Simulating biochemical signaling networks in complex moving geometries.
SIAM Journal on Scientific Computing , 32(5):3039–3070, 2010.[46] Melda Tozluo˘glu, Alexander L Tournier, Robert P Jenkins, Steven Hooper, Paul A Bates, and Erik Sahai. Matrix geometry determines optimalcancer cell migration strategy and modulates response to interventions.
Nature cell biology , 15(7):751–762, 2013.[47] Charles W Wolgemuth and Mark Zajac. The moving boundary node method: A level set-based, finite volume algorithm with applications tocell motility.
Journal of computational physics , 229(19):7287–7308, 2010.[48] Feng Wei Yang, Chandrasekhar Venkataraman, Vanessa Styles, and Anotida Madzvamuse. A robust and e ffi cient adaptive multigrid solverfor the optimal control of phase field formulations of geometric evolution laws. Communications in Computational Physics , 21(1):65–92,2017., 21(1):65–92,2017.