A moving grid finite element method applied to a mechanobiochemical model for 3D cell migration

TL;DR

This paper develops a moving grid finite element method to numerically solve a complex 3D mechanobiochemical model of cell migration, capturing cell deformation, protrusion, and contraction driven by biochemical and mechanical interactions.

Contribution

It introduces a novel moving grid finite element approach for solving highly nonlinear 3D cell migration models involving reaction-diffusion and biomechanical forces.

Findings

Numerical simulations show diverse cell deformation behaviors.

Linear stability analysis supports early migration dynamics.

The method enables studying complex biochemical and mechanical interactions.

Abstract

This work presents the development, analysis and numerical simulations of a biophysical model for 3D cell deformation and movement, which couples biochemical reactions and biomechanical forces. We propose a mechanobiochemical model which considers the actin filament network as a viscoelastic and contractile gel. The mechanical properties are modelled by a force balancing equation for the displacements, the pressure and concentration forces are driven by actin and myosin dynamics, and these are in turn modelled by a system of reaction-diffusion equations on a moving cell domain. The biophysical model consists of highly non-linear partial differential equations whose analytical solutions are intractable. To obtain approximate solutions to the model system, we employ the moving grid finite element method. The numerical results are supported by linear stability theoretical results close to bifurcation points during the early stages of cell migration. Numerical simulations exhibited show both simple and complex cell deformations in 3-dimensions that include cell expansion, cell protrusion and cell contraction. The computational framework presented here sets a strong foundation that allows to study more complex and experimentally driven reaction-kinetics involving actin, myosin and other molecular species that play an important role in cell movement and deformation.

Tables

Table 1. Table 1: Parameters for simulations in this section.

Section	4.2	4.3	4.4	4.5	4.6
Figure	2	3	5	8	10
$\psi$	200	20	150	100	100
$c$	-40	-80	-40	-80	-100
$k_a$	0.04	0.4	0.04	0.4	0.09
$k_{ma}$	0.05	0.5	0.05	0.05	0.09
$k_{am}$	0.06	0.12	0.06	0.07	0.15

Equations

$$\nabla \cdot ({\mathbf{\sigma }}_v+{\mathbf{\sigma }}_e+{\mathbf{\sigma }}_c+{\mathbf{\sigma }}_p)$$

$$\frac{\partial a}{\partial t}+\nabla \cdot (a\mathbf{\beta })-D_a\Delta a-f(a,m)$$

$$\frac{\partial m}{\partial t}+\nabla \cdot (m\mathbf{\beta })-D_m\Delta m-g(a,m)$$

$$a(\mathbf{x}(t),t)=a_0,\mathbf{u}(\mathbf{x}(t),t)$$

$$\mathbf{\beta }$$

$${\mathbf{\sigma }}_v\cdot \mathbf{n}={\mathbf{\sigma }}_e\cdot \mathbf{n}=\mathbf{n}\cdot \nabla a=\mathbf{n}\cdot \nabla m$$

$$\delta(l)=\left\{\begin{array}[]{l l}1,&\;\text{if }\sqrt{\xi_{x}^{2}+\xi_{y}^{2}+\xi_{z}^{2}}>0.8,\\ 0,&\;\text{otherwise.}\end{array}\right.$$

$$f(a,m)$$

$$g(a,m)$$

$$ϵ(\mathbf{u}):=\frac{1}{2}(\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^T)=\left(\begin{array}{ccc}\frac{\partial u}{\partial x}& \frac{1}{2}(\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y})& \frac{1}{2}(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x})\\ \frac{1}{2}(\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y})& \frac{\partial v}{\partial y}& \frac{1}{2}(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y})\\ \frac{1}{2}(\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y})& \frac{1}{2}(\frac{\partial w}{\partial y}+\frac{\partial v}{\partial z})& \frac{\partial w}{\partial z}\end{array}\right)\text{ and }\varphi (\mathbf{u}):=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z},$$

$${\mathbf{\sigma }}_v=\left(\begin{array}{ccc}{\mu }_{1+2}\frac{\partial \dot{u}}{\partial x}+{\mu }_2(\frac{\partial \dot{v}}{\partial y}+\frac{\partial \dot{w}}{\partial z})& \frac{{\mu }_1}{2}(\frac{\partial \dot{v}}{\partial x}+\frac{\partial \dot{u}}{\partial y})& \frac{{\mu }_1}{2}(\frac{\partial \dot{w}}{\partial x}+\frac{\partial \dot{u}}{\partial z})\\ \frac{{\mu }_1}{2}(\frac{\partial \dot{v}}{\partial x}+\frac{\partial \dot{u}}{\partial y})& {\mu }_2(\frac{\partial \dot{u}}{\partial x}+\frac{\partial \dot{w}}{\partial z})+{\mu }_{1+2}\frac{\partial \dot{v}}{\partial y}& \frac{{\mu }_1}{2}(\frac{\partial \dot{v}}{\partial z}+\frac{\partial \dot{w}}{\partial y})\\ \frac{{\mu }_1}{2}(\frac{\partial \dot{w}}{\partial x}+\frac{\partial \dot{u}}{\partial z})& \frac{{\mu }_1}{2}(\frac{\partial \dot{v}}{\partial z}+\frac{\partial \dot{w}}{\partial y})& {\mu }_{1+2}\frac{\partial \dot{w}}{\partial z}+{\mu }_2(\frac{\partial \dot{u}}{\partial x}+\frac{\partial \dot{v}}{\partial y})\end{array}\right),$$

$${\mathbf{\sigma }}_e=\frac{E}{1+\nu }\left(\begin{array}{ccc}\frac{\partial u}{\partial x}+{\nu }^{\prime }\varphi (\mathbf{u})& \frac{1}{2}(\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y})& \frac{1}{2}(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x})\\ \frac{1}{2}(\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y})& \frac{\partial v}{\partial y}+{\nu }^{\prime }\varphi (\mathbf{u})& \frac{1}{2}(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y})\\ \frac{1}{2}(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z})& \frac{1}{2}(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y})& \frac{\partial w}{\partial z}+{\nu }^{\prime }\varphi (\mathbf{u})\end{array}\right),$$

$${\mathbf{\sigma }}_c=\left(\begin{array}{ccc}\psi a^2{e}^{-a/a_{sat}}+cm& 0& 0\\ 0& \psi a^2{e}^{-a/a_{sat}}+cm& 0\\ 0& 0& \psi a^2{e}^{-a/a_{sat}}+cm\end{array}\right),$$

$${\mathbf{\sigma }}_p=\left(\begin{array}{ccc}\frac{p}{1+\varphi }\left(1+\frac{2}{\pi }\delta (l){\tan }^{-1}a\right)& 0& 0\\ 0& \frac{p}{1+\varphi }\left(1+\frac{2}{\pi }\delta (l){\tan }^{-1}a\right)& 0\\ 0& 0& \frac{p}{1+\varphi }\left(1+\frac{2}{\pi }\delta (l){\tan }^{-1}a\right)\end{array}\right),$$

$$\displaystyle\frac{\partial}{\partial x}\Bigg{(}D_{11}\frac{\partial\dot{u}}{\partial x}+D_{12}\bigg{(}\frac{\partial\dot{v}}{\partial y}+\frac{\partial\dot{w}}{\partial z}\bigg{)}+C_{11}\frac{\partial u}{\partial x}+C_{12}\bigg{(}\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}\bigg{)}\Bigg{)}\hskip 142.26378pt$$

$$\displaystyle+\frac{\partial}{\partial y}\Bigg{(}D_{33}\bigg{(}\frac{\partial\dot{v}}{\partial x}+\frac{\partial\dot{u}}{\partial y}\bigg{)}+C_{33}\bigg{(}\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\bigg{)}\Bigg{)}+\frac{\partial}{\partial z}\Bigg{(}D_{33}\bigg{(}\frac{\partial\dot{w}}{\partial x}+\frac{\partial\dot{u}}{\partial z}\bigg{)}+C_{33}\bigg{(}\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\bigg{)}\Bigg{)}$$

$$\displaystyle\frac{\partial}{\partial x}\Bigg{(}D_{33}\bigg{(}\frac{\partial\dot{v}}{\partial x}+\frac{\partial\dot{u}}{\partial y}\bigg{)}+C_{33}\bigg{(}\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\bigg{)}\Bigg{)}+\frac{\partial}{\partial y}\Bigg{(}D_{11}\frac{\partial\dot{v}}{\partial y}+D_{12}\bigg{(}\frac{\partial\dot{u}}{\partial x}+\frac{\partial\dot{w}}{\partial z}\bigg{)}\hskip 113.81102pt$$

$$\displaystyle+C_{11}\frac{\partial v}{\partial y}+C_{12}\bigg{(}\frac{\partial u}{\partial x}+\frac{\partial w}{\partial z}\bigg{)}\Bigg{)}+\frac{\partial}{\partial z}\Bigg{(}D_{33}\bigg{(}\frac{\partial\dot{w}}{\partial y}+\frac{\partial\dot{w}}{\partial z}\bigg{)}+C_{33}\bigg{(}\frac{\partial w}{\partial z}+\frac{\partial v}{\partial z}\bigg{)}\Bigg{)}$$

$$\displaystyle\frac{\partial}{\partial x}\Bigg{(}D_{33}\bigg{(}\frac{\partial\dot{w}}{\partial x}+\frac{\partial\dot{u}}{\partial z}\bigg{)}+C_{33}\bigg{(}\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\bigg{)}\Bigg{)}+\frac{\partial}{\partial y}\Bigg{(}D_{33}\bigg{(}\frac{\partial\dot{v}}{\partial z}+\frac{\partial\dot{w}}{\partial y}\bigg{)}+C_{33}\bigg{(}\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\bigg{)}\Bigg{)}\hskip 85.35826pt$$

$$\displaystyle+\frac{\partial}{\partial z}\Bigg{(}D_{11}\frac{\partial\dot{w}}{\partial z}+D_{12}\bigg{(}\frac{\partial\dot{v}}{\partial y}+\frac{\partial\dot{u}}{\partial x}\bigg{)}+C_{11}\frac{\partial w}{\partial z}+C_{12}\bigg{(}\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\bigg{)}\Bigg{)}$$

$$f_1=f_2=f_3=\left[\frac{p}{1+\varphi }\left(1+\frac{2}{\pi }\delta (l)\arctan a\right)+\psi a^2{e}^{-a/a_{sat}}+cm\right],$$

$$D_{11}={\mu }_1+{\mu }_2,D_{12}={\mu }_2,D_{33}=\frac{{\mu }_1}{2},$$

$$C_{11}=\frac{E(1-\nu )}{(1+\nu )(1-2\nu )},C_{12}=\frac{E\nu }{(1+\nu )(1-2\nu )}\text{ and }{C}_{33}=\frac{E}{2(1+\nu )}.$$

$${\int }_{{\mathbf{\Omega }}_t}$$

$$\displaystyle+C_{33}\bigg{(}\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\bigg{)}\Bigg{)}+\frac{\partial\hat{\phi}}{\partial z}\left(D_{33}\bigg{(}\frac{\partial\dot{w}}{\partial x}+\frac{\partial\dot{u}}{\partial z}\bigg{)}+C_{33}\bigg{(}\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\bigg{)}\right)d{{\bf\Omega}_{t}}=-\int_{{\bf\Omega}_{t}}\frac{\partial\hat{\phi}}{\partial x}f_{1}d{{\bf\Omega}_{t}}+\int_{{\bf\partial\Omega}_{t}}\hat{\phi}f_{1}n_{1}ds,$$

$${\int }_{{\mathbf{\Omega }}_t}$$

$$\displaystyle+C_{11}\frac{\partial v}{\partial y}+C_{12}\bigg{(}\frac{\partial u}{\partial x}+\frac{\partial w}{\partial z}\bigg{)}\Bigg{)}+\frac{\partial\hat{\phi}}{\partial z}\left(D_{33}\bigg{(}\frac{\partial\dot{w}}{\partial y}+\frac{\partial\dot{w}}{\partial z}\bigg{)}+C_{33}\bigg{(}\frac{\partial w}{\partial z}+\frac{\partial v}{\partial z}\bigg{)}\right)d{{\bf\Omega}_{t}}=-\int_{{\bf\Omega}_{t}}\frac{\partial\hat{\phi}}{\partial y}f_{2}d{{\bf\Omega}_{t}}+\int_{{\bf\partial\Omega}_{t}}\hat{\phi}f_{2}n_{2}ds,$$

$${\int }_{{\mathbf{\Omega }}_t}$$

$$\displaystyle+\frac{\partial\hat{\phi}}{\partial z}\left(D_{11}\frac{\partial\dot{w}}{\partial z}+D_{12}\bigg{(}\frac{\partial\dot{v}}{\partial y}+\frac{\partial\dot{u}}{\partial x}\bigg{)}+C_{11}\frac{\partial w}{\partial z}+C_{12}\bigg{(}\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\bigg{)}\right)d{{\bf\Omega}_{t}}=-\int_{{\bf\Omega}_{t}}\frac{\partial\hat{\phi}}{\partial z}f_{3}d{{\bf\Omega}_{t}}+\int_{{\bf\partial\Omega}_{t}}\hat{\phi}f_{3}n_{3}ds.$$

$${\int }_{{\mathbf{\Omega }}_t}\frac{\partial f_j}{\partial x}\hat{\varphi }d{\mathbf{\Omega }}_t=-{\int }_{{\mathbf{\Omega }}_t}\frac{\partial \hat{\varphi }}{\partial x}{f}_{j}d{\mathbf{\Omega }}_t+{\int }_{{\partial \Omega }_t}\hat{\varphi }{f}_j{n}_{j}ds,$$

$$\frac{\partial a}{\partial t}+\nabla \cdot (a\mathbf{\beta })-D_a\Delta a=f(a,m),\frac{\partial m}{\partial t}+\nabla \cdot (m\mathbf{\beta })-D_m\Delta m=g(a,m).$$

$$\frac{Da}{Dt}-D_a\Delta a+a(\nabla \cdot \mathbf{\beta })=f(a,m),$$

$$\frac{Dm}{Dt}-D_m\Delta m+m(\nabla \cdot \mathbf{\beta })=g(a,m).$$

$$\displaystyle\frac{\partial}{\partial t}\int_{{\bf\Omega}_{t}}a\hat{\psi}d{{\bf\Omega}_{t}}+\int_{{\bf\Omega}_{t}}\Big{(}D_{a}\nabla a\cdot\nabla\hat{\psi}\Big{)}d{{\bf\Omega}_{t}}=\int_{{\bf\Omega}_{t}}\bigg{(}f(a,m)\hat{\psi}+a\frac{D\hat{\psi}}{Dt}\bigg{)}\;d{{\bf\Omega}_{t}},$$

$$\displaystyle\frac{\partial}{\partial t}\int_{{\bf\Omega}_{t}}m\hat{\psi}d{{\bf\Omega}_{t}}+\int_{{\bf\Omega}_{t}}\Big{(}D_{m}\nabla m\cdot\nabla\hat{\psi}\Big{)}d{{\bf\Omega}_{t}}=\int_{{\bf\Omega}_{t}}\bigg{(}g(a,m)\hat{\psi}+m\frac{D\hat{\psi}}{Dt}\bigg{)}\;d{{\bf\Omega}_{t}},$$

$${\int }_{{\mathbf{\Omega }}_{h,t}}$$

$$\displaystyle\hskip 85.35826pt+\frac{\partial\hat{\phi}}{\partial z}\Bigg{(}D_{33}\bigg{(}\frac{\partial\dot{w}_{h}}{\partial x}+\frac{\partial\dot{u}_{h}}{\partial z}\bigg{)}+C_{33}\bigg{(}\frac{\partial w_{h}}{\partial x}+\frac{\partial u_{h}}{\partial z}\bigg{)}\Bigg{)}d{{\bf\Omega}_{t}}=-\int_{{\bf\Omega}_{h,t}}\frac{\partial\hat{\phi}}{\partial x}f_{1}d{{\bf\Omega}_{h,t}}+\int_{{\bf\partial\Omega}_{h,t}}\hat{\phi}f_{1}n_{1}ds,$$

$${\int }_{{\mathbf{\Omega }}_{h,t}}$$

$$\displaystyle\hskip 85.35826pt+\frac{\partial\hat{\phi}}{\partial z}\Bigg{(}D_{33}\bigg{(}\frac{\partial\dot{w}_{h}}{\partial y}+\frac{\partial\dot{w}_{h}}{\partial z}\bigg{)}+C_{33}\bigg{(}\frac{\partial w_{h}}{\partial z}+\frac{\partial v_{h}}{\partial z}\bigg{)}\Bigg{)}d{{\bf\Omega}_{h,t}}=-\int_{{\bf\Omega}_{h,t}}\frac{\partial\hat{\phi}}{\partial y}f_{2}d{{\bf\Omega}_{h,t}}+\int_{{\bf\partial\Omega}_{h,t}}\hat{\phi}f_{2}n_{2}ds,$$

$${\int }_{{\mathbf{\Omega }}_{h,t}}$$

$$\displaystyle+\frac{\partial\hat{\phi}}{\partial z}\Bigg{(}D_{11}\frac{\partial\dot{w}_{h}}{\partial z}+D_{12}\bigg{(}\frac{\partial\dot{v}_{h}}{\partial y}+\frac{\partial\dot{u}_{h}}{\partial x}\bigg{)}+C_{11}\frac{\partial w_{h}}{\partial z}+C_{12}\bigg{(}\frac{\partial u_{h}}{\partial x}+\frac{\partial v_{h}}{\partial y}\bigg{)}\Bigg{)}d{{\bf\Omega}_{h,t}}=-\int{{\bf\Omega}_{h,t}}\frac{\partial\hat{\phi}}{\partial z}f_{3}d{{\bf\Omega}_{h,t}}+\int_{{\bf\partial\Omega}_{h,t}}\hat{\phi}f_{3}n_{3}ds,$$

$$\frac{\partial }{\partial t}{\int }_{{\Omega }_{h,t}}{a}_h\hat{\psi }d{\Omega }_{h,t}+{\int }_{{\Omega }_{h,t}}{D}_a\nabla a_h\cdot$$

$$\frac{\partial }{\partial t}{\int }_{{\Omega }_{h,t}}{m}_h\hat{\psi }d{\Omega }_{h,t}+{\int }_{{\Omega }_{h,t}}{D}_m\nabla m_h\cdot$$

$$u_h=\sum _{j=1}^{nde}{U}_j{\varphi }_j,v_h=\sum _{j=1}^{nde}{V}_j{\varphi }_j,w_h=\sum _{j=1}^{nde}{W}_j{\varphi }_j,a_h=\sum _{j=1}^{nde}{\alpha }_j{\varphi }_j\text{ and }{m}_h=\sum _{j=1}^{nde}{\mu }_j{\varphi }_j.$$

$$\left[\begin{array}{ccc}{\mathbf{A}}^{11}& {\mathbf{A}}^{12}& {\mathbf{A}}^{13}\\ [{\mathbf{A}}^{12}{]}^T& {\mathbf{A}}^{22}& {\mathbf{A}}^{23}\\ [{\mathbf{A}}^{13}{]}^T& [{\mathbf{A}}^{23}{]}^T& {\mathbf{A}}^{33}\end{array}\right]\left\{\begin{array}{c}\frac{d\mathbf{U}}{dt}\\ \frac{d\mathbf{V}}{dt}\\ \frac{d\mathbf{W}}{dt}\end{array}\right\}+\left[\begin{array}{ccc}{\mathbf{B}}^{11}& {\mathbf{B}}^{12}& {\mathbf{B}}^{13}\\ [{\mathbf{B}}^{12}{]}^T& {\mathbf{B}}^{22}& {\mathbf{B}}^{23}\\ [{\mathbf{B}}^{13}{]}^T& [{\mathbf{B}}^{23}{]}^T& {\mathbf{B}}^{33}\end{array}\right]\left\{\begin{array}{c}\mathbf{U}\\ \mathbf{V}\\ \mathbf{W}\end{array}\right\}=\left\{\begin{array}{c}{\mathbf{F}}^1\\ {\mathbf{F}}^2\\ {\mathbf{F}}^3\end{array}\right\},$$

Full text

1 Introduction

Cell movement is critical in multicelluar organisms due to roles in embryogenesis, wound healing, immune response, cancer metastasis, tumour invasion, and other processes, therefore, understanding cell movement is of great importance to medicine and to understanding our origins [8, 9, 11, 14]. In this study we consider a mechanobiochemical model previously studied by George et al. [15, 16, 27] which we will extend to 3-dimensions as well as introducing for the first time, the role of myosin in the modelling and computational framework. The model comprises of a system of reaction-diffusion equations for cellular proteins and a viscoelastic mechanical model for cell movement and deformation. 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 for solving complex partial differential equation (PDEs), methods that have been employed to model cell motility include finite differences, phase field methods, boundary element methods (BEM), immersed boundary methods or level set methods (LSM), [1, 7, 32, 35, 44, 45, 47]. Choosing a suitable method for a particular model is a balance between the ease of application within the model’s framework and the reliability of solutions produced.

Finite differences were used in previous incarnations of the model [1, 44]. This method is very useful and easy to implement on fixed and simple domains but it is significantly more complicated to incorporate for the evolving domains and surfaces we wish to use and there are often problems with a moving boundary. Level set methods are used extensively in cell simulations and are useful when cells split and reconnect, therefore, it may be advantageous to use this method in the future when considering cell proliferation (cell division) and apoptosis (cell death) [48]. In this 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 be generalised to multidimensions with little complication, hence it is the ideal method to numerically solve our model system. 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 diffusion equations. Additionally we develop 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. Theoretical predictions of the spatiotemporal behaviour of the solutions of the model close to bifurcation points are presented in A, these identify important bifurcation parameters. In Section 3, theoretical predictions are used to validate finite element simulations since there are no analytical solutions to compare with. Further finite element simulations illustrating 3D cell movement and deformation are exhibited in Section 4. We discuss findings and conclude in Section 5.

2 A mechanobiochemical model

The model we consider and extend is inspired by contractile models of the actin cytogel [23, 34]. These models comprise of a force balance equation modelling the displacements of the cell when deformed and a reaction-diffusion equation for the concentration of the gel that in turn drives cell movement. The idea of pressure driven protrusion and the use of concentration of actin originates from Alt and Tranquillo [1]. In their model they assume movement is produced by a balance between contractile force of the actin network pulling on the membrane and pressure pushing on the membrane. This was extended by Stephanou et al. so that large deformations could be modelled which is more realistic for most cells [44]. George further extended this model by observing (and hence modelling) that higher actin concentration in a region leads to more pressure [16]. In the previous models a polar coordinate system was used and radial extension of the cell was calculated [1, 44]. Unlike this approach, we follow the work of George and study the mechanobiochemical 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 by adding the consideration of the concentration of myosin. We model the concentrations of, and interactions between, actin and myosin using two reaction-diffusion equations. The reaction-diffusion equations are coupled to a force balance 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: ${{\bf\Omega}_{t}}\subset\mathbb{R}^{3}$ with boundary ${{\bf\partial\Omega}_{t}}$, where $t\in I=[0,T_{f}],$ $T_{f}>0$. Any point ${\bf x}\in{{\bf\Omega}_{t}}$ is defined by ${\bf x}=(x(t),y(t),z(t))$. We define the displacement of ${\bf x}$ at time $t$ by ${\bf u}=(u({\bf x}(t),t),v({\bf x}(t),t),w({\bf x}(t),t))^{T}$. Let the concentration of F-actin, and bound myosin, at point ${\bf x}(t)$ be given by $a=a({\bf x}(t),t)$, and $m=m({\bf x}(t),t)$, respectively. We describe the dynamics of the actin network by the following system of the equations,

[TABLE]

where the viscoelastic and contractile properties are described by stress tensors:

• •

viscous $\boldsymbol{\sigma}_{v}({\bf u})=\mu_{1}\frac{\partial\boldsymbol{\varepsilon}}{\partial t}+\mu_{2}\frac{\partial\phi}{\partial t}{\bf I}$ where $\mu_{1}$ and $\mu_{2}$ are shear and bulk viscosities respectively, $\boldsymbol{\varepsilon}$ is the the strain tensor ($\frac{1}{2}(\nabla{\bf u}+\nabla{\bf u}^{T})$) and $\phi$ is the dilation ($\nabla\cdot{\bf u}$).

• •

elastic $\boldsymbol{\sigma}_{e}({\bf u})=\frac{E}{1+v}(\boldsymbol{\varepsilon}+\frac{\nu}{1-2\nu}\phi{\bf I})$ where $E$ is the Youngs modulus and $\nu$ is the Poisson ratio.

• •

contractile $\boldsymbol{\sigma}_{c}(a,m)=(\psi a^{2}e^{-a/a_{sat}}+cm){\bf I}$, where $\psi$ and $c$ are the contractility coefficients for $a$ and $m$, respectively, and $a_{sat}$ is the saturation coefficient of actin.

• •

pressure \boldsymbol{\sigma}_{p}({\bf u},a)=\dfrac{p}{1+\phi}\big{(}1+\frac{2}{\pi}\delta(l)\arctan a\big{)}{\bf I}. This describes two types of pressure. First is the hydrostatic pressure which is present everywhere and corresponds to the osmotic pressure in the cell which depends on the dilation $\phi$ and pressure coefficient $p$. Close to the membrane there is also polymerisation pressure caused by the polymerising actin filaments pushing on the cell membrane. This increases with increasing concentration of filaments $a$. We choose close to the membrane to mean less than 20% of the cell radius from the edge in the initial state. To define this we use $\delta(l)$ and the points $\boldsymbol{\xi}=(\xi_{x},\xi_{y},\xi_{z})\in{\bf\Omega}_{0}$. There exists a family of bijective mappings between the initial and current domains we can let $l:{{\bf\Omega}_{t}}\times I\to\mathbb{R}$ and corresponding $\hat{l}:{\bf\Omega}_{0}\times I\to[0,1]$ then $\hat{l}(\boldsymbol{\xi},t)=l({\bf x}(\boldsymbol{\xi},t),t)$. So we calculate the distance from the centroid by

[TABLE]

In the reaction-diffusion equations we have the diffusion coefficients for actin and myosin, $D_{a}$ and $D_{m}$ respectively. Because the cell is moving we introduce the flow velocity $\boldsymbol{\beta}=\frac{\partial{\bf u}}{\partial t}$. The interactions between actin and myosin are described by the reaction terms $f(a,m)$ and $g(a,m)$ respectively. We have formulated different plausible reaction kinetics in the absence of experimental data and for the sake of brevity, we only present here one such model. We refer the interested reader to consult [30]. Although we will discuss the implementation of one specific plausible model, other models can be easily incorporated and studied in a similar fashion. For illustrative purposes we consider the following hypothetical reaction kinetics

[TABLE]

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 polymerisation/depolymerisation and $a_{c}$ is the equilibrium concentration and if the concentration is above this critical value then F-actin will depolymerise at the same rate. Next, since myosin binds to actin the amount of myosin will increase due 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 will depolymerise with higher concentrations of myosin and is subject to a saturation coefficient $K$, for $a$. The negation is true for myosin since myosin is seen to accumulate.

Thus we have three connected equations: the solutions to (1b) (actin concentration) and (1c) (myosin concentration) affect the contractile and pressure parts of the force balance equation and the solution to (1a) (displacement) affects the reaction-diffusion equations through the convection terms and the changing shape of the domain.

2.1 Summary of results from linear stability theory

The linear stability analysis detailed in A reveals that parameters, in particular $\psi$ and $c$, can be varied so that particular patterns become unstable and grow. These patterns correspond to eigenfunctions of the Laplacian on the chosen volume.

3 Finite element formulation

We have formulated a very complex non-linear system of partial differential equations. It is not possible to analytically solve this system, therefore we must turn to numerical methods to produce approximations. Previous versions of this model were represented in a polar coordinate system and solved using finite differences however George [15] used a finite element formulation. We proceed similarly and describe our methods below. In particular, we employ the moving grid finite element method [3, 25, 26, 27] to compute approximate numerical solutions of the coupled viscoelastic reaction-diffusion system defined in 3D Cartesian coordinate system.

To begin, the force balance is separated into a system of three partial differential equations representing the three space dimensions. This clarifies the derivation of the weak formulation. Since $\boldsymbol{\sigma}_{v},\boldsymbol{\sigma}_{e},\boldsymbol{\sigma}_{c}$ and $\boldsymbol{\sigma}_{p}$ (as described in Section 2) are all stress tensors we can write them in matrix form. In three dimensions, strain and dilation are given by

[TABLE]

respectively. It follows then that we can write the stress tensors in three-dimensional tensor-matrix form:

[TABLE]

where $\mu_{1+2}=\mu_{1}+\mu_{2}$ and $\nu^{\prime}=\nu/(1-2\nu)$. Substituting these expressions into $\nabla\cdot(\boldsymbol{\sigma}_{v}+\boldsymbol{\sigma}_{e}+\boldsymbol{\sigma}_{c}+\boldsymbol{\sigma}_{p})={\bf 0}$ gives us three equations

[TABLE]

where

[TABLE]

3.1 Weak formulation

To find the weak formulation, we take the usual route and multiply by a test function $\hat{\phi}({\bf x},t)\in H^{1}({{\bf\Omega}_{t}})$, where $H^{1}({{\bf\Omega}_{t}})$ is a Hilbert space, and integrate over the domain. This takes into account Green’s formula and the boundary conditions. The boundary condition $\boldsymbol{\sigma}_{v}\cdot{\bf n}=\boldsymbol{\sigma}_{e}\cdot{\bf n}=0$ means that boundary term disappears during integration. The weak formulation is to find $u({\bf x},t)$, $v({\bf x},t)$ and $w({\bf x},t)\in H^{1}({{\bf\Omega}_{t}})$, $t\in I$ such that

[TABLE]

Since $\frac{\partial f_{1}}{\partial x},\frac{\partial f_{2}}{\partial y}$ and $\frac{\partial f_{3}}{\partial z}$ are difficult to evaluate, we have used identities derived from the gradient theorem to write the weak form as above. In other words we have used the identity

[TABLE]

for $j=1,2,3$, where $x$ can also be substituted by $y$ and $z$. $n_{1},\;n_{2},\;n_{3}$ are the direction cosines of the outward unit vector ${\bf n}$ normal to ${{\bf\partial\Omega}_{t}}$.

Next we want to find the weak formulation of the reaction-diffusion equations which are given as

[TABLE]

We apply the product rule and convert to the material derivative (defined as $\frac{Da}{Dt}=\frac{\partial a}{\partial t}+a(\nabla\cdot\boldsymbol{\beta})$, in [37]). This gives

[TABLE]

Now continuing as with the force balance equation, we multiply by a test function $\hat{\psi}({\bf x},t)\in H^{1}({{\bf\Omega}_{t}})$ and integrate over the domain. The terms $(D_{a}\Delta a)\hat{\psi}$ and $(D_{m}\Delta m)\hat{\psi}$ can be simplified using the divergence theorem and for the remaining part of the left hand side we can use the Reynolds transport theorem. This means the weak formulation can be written as: Find $a({\bf x},t),\;m({\bf x},t)\in H^{1}({{\bf\Omega}_{t}}),t\in I$ such that

[TABLE]

for all $\hat{\psi}({\bf x},t)\in H^{1}({{\bf\Omega}_{t}})$.

3.2 Space discretisation

We now wish to define the problem at discrete points in space. To do this, we define the computational domain $\Omega_{h,t}$ as a polyhedral approximation to ${{\bf\Omega}_{t}}$, $T_{h,t}$ the discretisation of $\Omega_{h,t}$ made up of non-degenerate elements $\kappa_{i}$ and the finite element space $V_{h}(t):=\{v_{h}\in C^{0}({{\bf\Omega}_{t}}):v_{h}|_{\kappa}\text{ is linear}\}$. Thus the space-discrete problem is to find $u_{h}({\bf x},t)$, $v_{h}({\bf x},t)$, $w_{h}({\bf x},t)$, $a_{h}({\bf x},t)$, $m_{h}({\bf x},t)\in V_{h}(t),\;t\in I$, such that

[TABLE]

and

[TABLE]

for all $\hat{\phi},\;\hat{\psi}\in 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:

[TABLE]

This means that we are left with equations which contain only simple functions and their derivatives and point values for the variables. Substituting $u_{h}$, $v_{h}$, and $w_{h}$ into the above equations and using the Galerkin formulation, take the test functions to belong the same spaces as the nodal basis functions. Hence, the force balance equations can be written in block matrix-vector form

[TABLE]

where $\{{\bf U}(t)\}=(U_{1},...,U_{nde})$, $\{{\bf V}(t)\}=(V_{1},...,V_{nde})$ , $\{{\bf W}(t)\}=(W_{1},...,W_{nde})$ and:

[TABLE]

It must be noted that in the above we have used the fact that $\hat{\phi}\in V_{h}$ and $\hat{\psi}\in V_{h}$ where $V_{h}={\phi_{1},....,\phi_{n}}$. For the sake of ease of notation and computation we we define the following block matrices and vectors

[TABLE]

Therefore the force balance equation’s semi-discrete finite element formation can be written compactly as

[TABLE]

Now considering the reaction kinetics to be as in Eqs (3), (in Section 2), we can similarly write the reaction-diffusion equations in semi-discrete form

[TABLE]

where vector operations in $\boldsymbol{\alpha}(1-\boldsymbol{\mu})/(1+K\boldsymbol{\alpha}^{2})$ are pointwise and

[TABLE]

To compute these integrals we use Gauss numerical quadrature [36]. This is done as follows. First we can consider the integrals elementwise,

[TABLE]

Then choose a numerical quadrature comprising a set of points and weights depending on the functions to be integrated. This can be written as a formula for the integral of a function $\xi$

[TABLE]

where $\bar{\bf x}_{q}$ and $w_{q}$ are the $q$th quadrature points and weights respectively. Therefore the integrals can be approximated by

[TABLE]

These are all computed in the deal.II implementation [4].

3.3 Time discretisation

Next we carry out the temporal discretisation of the system of ordinary differential equations arising from the finite element discretisation. To proceed, we split the interval into a finite number of sub-intervals $[t^{n},t^{n+1}]$ use a uniform time step $\Delta t:=t^{n+1}-t^{n}$. We can then use a modified implicit-explicit (IMEX) finite differentiation formula [22, 25, 39]. Thus the fully discrete problem is now

[TABLE]

where the superscripts n and n+1 are the computed values on the mesh at times $t^{n}$ and $t^{n+1}$ respectively. Note that we have treated some parts implicitly (e.g. diffusion) 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 a system of linear algebraic equations. When solving (25a) we see that the block matrix on the left hand side is not symmetric therefore we use the most effective solver for this which is GMRES [40]. The equations (25b) and (25c) are solved using the conjugate gradient method [20].

3.4 Nodal displacements

The displacement of the nodes of the mesh is chosen to be equal to the flow velocity therefore $\boldsymbol{\beta}:=\frac{\partial{\bf U}}{\partial t}$. Since $t^{n+1}=t^{n}+\Delta t$ and ${\bf x}(t^{n})\in{\bf\Omega}_{t^{n}},{\bf x}(t^{n+1})\in{\bf\Omega}_{t^{n+1}}$ be points in the respective domains. We can define a first order linear approximation as:

[TABLE]

This means we can define a new approximation to the domain ${\bf\Omega}_{t^{n+1}}$ such that

[TABLE]

At each step we have a new mesh with new shape functions so we must assemble new matrices ${\bf M}^{n},{\bf H}^{n},{\bf A}^{n},{\bf B}^{n},{\bf F}^{n}$ to iteratively solve the discrete coupled problem as outlined in the following algorithm

3.5 Numerical algorithm

The fully discrete problem is solved iteratively with the following algorithm:

• •

Initialise ${\bf U}^{0}$, $\boldsymbol{\alpha}^{0}$, $\boldsymbol{\mu}^{0}$ and fixed parameters

• •

WHILE $t<endtime$

• –

Assemble ${\bf M}^{n},{\bf H}^{n},{\bf A}^{n},{\bf B}^{n},{\bf F}^{n}$

• –

Solve for ${\bf U}^{n+1}$ using (25a)

• –

Compute the new domain using ${\bf U}^{n+1}$

• –

Solve for $\boldsymbol{\alpha}^{n+1}$ and $\boldsymbol{\mu}^{n+1}$ using (25b) and (25c)

• –

$t=t+\Delta t$

• •

END

We create a mesh using Gmsh [17] and implement this algorithm using deal.II [4], a C++ software library which provides tools to solve partial differential equations which are discretised with finite element methods. Unlike the majority 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 regions which will cause the cell to move. This organisation may be caused by diffusion-driven instability, or due to the movement of the cell combined with the reaction-diffusion equations. The linear stability analysis of A holds true close to critical bifurcation points, these include parameters as well as the geometrical deformation of the cell. The conditions for stability are numerous and complex, however it is still possible to choose parameters so that particular modes can be selected. When we consider longer time, and therefore far away from equilibrium, linear stability theory no longer holds but we see significant protrusions and contractions which deform the mesh into many different shapes. Parameters used are in Table 1.

4.1 Excitation of mode $w_{1,1}^{1}$

In this example, the actin and myosin concentration solutions will be the negation of each other with actin concentration highest on one side and myosin concentration highest on the opposite side. This mode is the first eigenfunction that one might hope to see for the organisation of actin and myosin in a cell because it is similar to what is often observed in a moving cell [12, 31]. $k^{2}_{1,1}=2.0816$ is also the lowest eigenvalue. Choosing parameters $\psi=20,\;c=-80,\;k_{a}=0.04,\;k_{ma}=0.05$ and $k_{am}=0.06$ and initial conditions

[TABLE]

we observe that the mode $w_{1,1}^{1}$ is selected for actin and myosin. In Figure 1 we plot the concentrations of actin and myosin at time $t=1$. Blue indicates where the concentration is low, while red indicates that concentration is high. In this case very little deformation is seen.

4.2 Cell deformation when $w_{2,1}^{0}$ is excited initially

The simple first mode is not the only organisation which makes sense or shows similarities to organisation seen in cells. The cell can protrude in more than one direction because of actin accumulation at both ends, or deform in many other ways. Also, myosin could accumulate and ”squeeze” on both sides. Therefore we continue by isolating other modes. We see that both the parameters, and the initial conditions, have an effect on which mode will grow. The first large deformation is seen when choosing initial conditions

[TABLE]

In Figure 3 we plot the concentrations of actin and myosin and the displacement $\left(|{\bf U}|=\sqrt{U^{2}+V^{2}+W^{2}}\right)$ at each point. The cell expands at the two ends where actin concentration is high and contracts in the middle where myosin concentration is high. So far the results are visually similar to the the results obtained by [16] in the absence of myosin, one difference however is that there is only a very small volume increase because the cell is contracting in the middle as well as protruding. Other results (not shown) when the excited mode for myosin is the same as the mode for actin are very similar to the previous model [30]. Therefore, we investigate whether more interesting dynamics may occur if we try to excite differing modes for the two concentrations.

4.3 Cell deformation when $w_{1,1}^{1}$ and $w_{2,1}^{0}$ are excited for actin and myosin, respectively

While the idea that actin and myosin accumulate in opposite sides is quite well founded, their concentration gradients are rarely exactly opposite. Therefore here we investigate if differing modes can be excited for actin and myosin. Choosing appropriate initial conditions, to encourage different modes to grow, we observe more irregular deformations. In Figure 3 we plot the concentrations of actin and myosin and the displacement when the initial conditions are

[TABLE]

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.

4.4 Cell deformation when $w_{2,1}^{0}$ and $w_{1,1}^{0}$ are excited for actin and myosin, respectively

In Figure 5, the initial conditions

[TABLE]

contain the same two eigenfunctions as the last example but with different orientations, ($w_{1,1}^{0}$ is a rotation of $w_{1,1}^{1}$), we observe a quite different deformation. There is high actin concentration at the top and bottom of the sphere. Without the effect 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=5$ the protrusion slows and there is a contraction at the bottom where myosin concentration is high. There is another subsequent expansion and contraction with the actin and myosin concentrations reorganising to be is higher nearer the surface except when the cell is contracting, when the opposite is true, this is displayed in Figure 6. Figure 7b shows the translation of the cell and Figure 7a shows the change in length in the $z$-direction.

4.5 Cell deformation when $w_{1,1}^{0}$ and $w_{3,1}^{0}$ are excited for actin and myosin, respectively

Next, we begin with initial conditions

[TABLE]

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 translated in the negative $z$-direction and this is plotted in Figure 8, and the change in volume is illustrated in Figure 9.

4.6 Cell deformation when $w_{1,1}^{1}$ and $w_{4,1}^{0}$ are excited for actin and myosin, respectively

In another example of mixed modes, we start with

[TABLE]

This leads to the expansion shown in Figure 10. The cell contracts inwards at areas of high myosin concentration and protrudes in the remaining areas, there are large protrusions in two opposing directions, the largest being the direction where actin was initially highest, subsequently, actin concentrates in areas of high curvature and protrudes further.

4.7 $L_{2}$ norms

In Figure 11 we plot the norm of differences between successive solutions in the case of the full system example in Section 4.4. We see an increase, or decrease, in the $L_{2}$ norm when the rate of deformation is accelerating, or decelerating, respectively. The qualitative changes in the $L_{2}$ norms are similar, but the changes in myosin and displacement appear slightly later than actin. This may suggest, in this example, that the change in actin triggers the change in the other variables.

In most of the numerical solutions in Section 4 the cell becomes deformed in such a way that means the mesh becomes highly non-conforming and therefore the numerical method breaks down. In this case, re-meshing might be 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.

5 Conclusion

Our model revolves around an equation which balances elasticity, viscosity, contractility and pressure. Connected to this are two equations for the concentrations of F-actin and bound myosin. Unlike the previous study which our model is inspired by, (that of George [16]), we used the software library deal.II [4] to implement the finite element method. The key difference between this software and the previously used ALBERTA [42] is that the elements are hexahedra rather than tetrahedra. The implementation is therefore different but we were able to appropriately replicate the previous results by George [16] of cytomechanical model on a unit two dimensional disk (results not shown). 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, for the first time, we considered a second reaction-diffusion model to describe how myosin interacts with actin and how it contributes to cell contraction during cell migration. In the absence of experimental observations, we postulated hypothetical reaction kinetics describing the interaction between actin and myosin. As a first step in understanding solution behaviour in three dimensions of the full model, linear stability analysis close to bifurcation points was carried out and appropriate key parameter values were identified. An evolving finite element method was implemented in multi-dimensions.

Our numerical simulations showed that the model extends into three dimensions and the addition of myosin allows some 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, higher actin concentration is observed.

• •

This outward movement due to actin concentration is halted in areas with high myosin concentration. Additionally, if there is low actin and high myosin, we can see negative curvature.

• •

Identifying bifurcation parameters is more complicated than in previous models but effects of parameters can still be seen. The contractility due to myosin, $c$, strongly effects the speed of the deformation while the reaction constants $k_{a},\;k_{ma}$ and $k_{am}$, and the diffusion coefficients $D_{a}$ and $D_{m}$ play a part in which mode the actin and myosin concentrations will arrange into. It is not as possible to isolate single modes just by picking parameter values, however, choosing initial conditions as a random number multiplied by the eigenfunction means it is possible 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 conservation or constraint could be a useful extension.

The modelling and computational framework presented in this article can readily be adapted to consider new experimentally driven reaction kinetics between actin and myosin or interactions between three or more molecular species, 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-diffusion equation alone could not cause patterning [16]. Without the flow term, the prescribed reactions meant the concentration of actin would always return to the homogeneous steady state of $a=a_{c}$. In our case we have two coupled reaction-diffusion equations which are well known to induce patterning in certain cases.

We perform non-dimensionalisation to reduce parameters and simplify calculations. It also allows the reaction-diffusion equations to take the form necessary to use the standard conditions for diffusion driven instability. This is investigated in detail in [30]. We choose the nondimensionalised parameters:

[TABLE]

In the above, $L$ is the typical radius of a cell. Substituting appropriately and carrying out algebraic manipulations leads to the following non-dimensionalised system (for general kinetics)

[TABLE]

where $\sigma(\tilde{a})=\tilde{\psi}\tilde{a}^{2}e^{-\tilde{a}/\tilde{a}_{sat}}$, $\nu^{\prime}=\frac{\nu}{1-2\nu}$ and \tilde{p}(\tilde{a})=\tilde{p}\big{(}1+\frac{2}{\pi}\delta(l)\arctan\tilde{a}\big{)}. $f$ and $g$ have been nondimensionalised and we choose only functions such that system has a steady state at $(a_{s},m_{s},{\bf u}_{s})=(1,1,{\bf 0})$. Given small variations $\hat{a},\hat{m}$ and $\hat{{\bf u}}$, consider the perturbation from the steady state $\tilde{a}=a_{s}+\hat{a}=1+\hat{a},\;\tilde{m}=m_{s}+\hat{m}=1+\hat{m},\;\tilde{{\bf u}}={\bf u}_{s}+\hat{{\bf u}}=\hat{{\bf u}}$. This results in the linear system

[TABLE]

We now look for solutions of the form

[TABLE]

where $\lambda$ is the growth rate, ${\bf k}$ is the wave vector, and $a^{*},m^{*}$ and ${\bf u}^{*}$ are small amplitudes. We require solutions to be non-trivial and so we obtain the stability matrix

[TABLE]

Thus $h(\lambda)=0$ (32b) is our dispersal relation and we are concerned with the solution $\lambda$. There will be instability when $\operatorname{Re}(\lambda)>0$. We exploit this relation to isolate particular patterns/modes. The unstable modes will correspond to the eigenfunctions of the Laplacian on the sphere and $k^{2}$ the associated eigenvalues.

A.1 Eigenvalues and eigenfunctions of the Laplacian in the bulk of the unit sphere

The eigenvalues on the unit sphere ${\bf\Omega}_{0}=\{(x,y,z):x^{2}+y^{2}+z^{2}\leq 1\}$ (with homogeneous Neumann boundary condition) are well known and are obtained using separation of variables [2, 29]. There are an infinite number of discrete solutions of the form

[TABLE]

where $l,m,n$ are all integers such that $|m|\leq l\leq n,\;A_{l,n}^{m}$ are constants, $J_{\alpha}(x)$ is a Bessel function of the first kind, i.e. $J_{\alpha}(x)=\sum_{j=0}^{\infty}\frac{(-1)^{j}}{j!\Gamma(1+j+\alpha)}\left(\frac{x}{2}\right)^{2j+\alpha}$ with $\Gamma(n)=(n-1)!$, $P_{l}^{m}(x)$ are associated Legendre polynomials and $j^{\prime}_{l+\frac{1}{2},n}$ are zeros of the differential of the spherical Bessel function. We can find the eigenvalues $k_{l,n}^{2}=(j^{\prime}_{l+\frac{1}{2},n})^{2}$ numerically. It follows that for each eigenvalue $\lambda_{l,n}=k_{l,n}^{2}$ there are $2l+1$ possible eigenfunctions. Figure 12 shows the eigenfunctions for some selected values of $l$, $m$ and $n$. The wave numbers are discrete.

A.2 Parameter selection

The conditions on the positivity of the roots of (32b) are numerous and the coefficients of the polynomial are burdensome. 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 ($\psi$) are particularly significant for finding unstable wavenumbers. In Figure 13a we plot the real and imaginary parts of the solution against $k^{2}$ for three different values of c. We can see that when wavenumbers $k^{2}$ are less than $\sim 8.5$, then $\operatorname{Re}(\lambda)>0$ for the three values of $c$, therefore the wavenumbers will be unstable. Additionally these wavenumbers will be oscillatory for $c=10$. There are also regions of Hopf instability and oscillatory instability for all the three values of $c$. In Figure 13b we fix other parameters and vary $\psi$ to see that, just like in [16] higher values of $\psi$ mean higher wavenumbers can be excited.

Acknowledgements

This work (LM) was supported by an EPSRC Doctoral Training Centre Studentship through the University of Sussex. 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 the Marie Sklodowska-Curie grant agreement (No 642866). AM is a Royal Society Wolfson Research Merit Award Holder, generously supported by the Wolfson Foundation. LM acknowledges the support from the University of Sussex ITS for computational purposes.