A Dynamic Finite Element Cellular Model and Its Application on Cell Migration

Abstract

While the tissue is formed or regenerated, cells migrate collectively and remained adherent. However, it is still unclear what are the roles of cell-substrate and intercellular interactions in regulating collective cell migration. In this chapter, we introduce our newly developed finite element cellular model to simulate the collective cell migration and explore the effects of mechanical feedback between cells and between cell and substrate. Our viscoelastic model represents one cell with many triangular elements. Intercellular adhesions between cells are represented as linear springs. Furthermore, we include a mechano-chemical feedback loop between cell-substrate mechanics and cell migration. Our results reproduce a set of experimental observation of patterns of collective cell migration during epithelial wound healing. In addition, we demonstrate that cell-substrate determined mechanics play an important role in regulating persistent and oriented collective cell migration. This chapter illustrates that our finite element cellular model can be applied to study a number of tissue related problems regarding cellular dynamic changes at subcellular level.

1. Introduction

Thanks to the accurate description of changes in material mechanics, finite element method has been widely used in the field of bioengineering to study cellular tissue related problems such as neurulation and epithelial mechanics [1, 2]. However, majority of current finite element models are only restricted on tissues undergoing changes of shapes and displacements at small scale. In addition, during the simulation, the cellular tissue is required to be remained as one integrity. These limitations restrict the traditional finite element method to be applied to study the essential physiological processes such as morphogenesis, tissue regeneration, tumor metastasis, and cancer invasion, where cells often migrate collectively as large coherent strands or tubes. Such large scale of collective cell movement is recognized as the hallmark of tissue-remodeling events. During the past decade, to overcome the limitation of traditional finite element method, dynamic finite element method such as PFEM has been developed to extend the traditional FEM to study mechanics of materials with more flexibility or undergoing larger scale of motility. The object domain (either fluid or solid) is represented as nodes tessellated by triangular mesh. The mathematical equations governing the physical rules of the mechanical property of the discretized domain defined by the mesh connecting nodes are subsequentially solved in the standard FEM. Under the analysis using dynamic finite element method, the motion of sub-domain of the object can freely move and even separate from the main domain [3]. The advancement of dynamic finite element in achieving both accurate description of material mechanics and large scale of geometric and topological changes makes it suitable to simulate the physiological processes such as wound healing and cancer invasion. During these physiological processes, cells move in collective fashion and respond with chemical and mechanical signals through cell–cell junctions and interactions between cells and their micro-environment.

In this chapter, we introduced our newly developed dynamic finite element cellular model and its application to study the influence of cell-substrate mechanics and intercellular adhesions on collective cell migration. Our model represents each cell as a mesh of triangular elements at sub-cellular level [4]. Each triangular element exhibits viscoelastic characteristic using a Maxwellian model [5]. The effects of line tension forces along the cell boundary according to the local curvature is incorporated [6]. The intercellular adhesions are modeled as elastic springs at sub-cellular scale [7]. In addition, a mechano-chemical feedback pathway including focal adhesion, proteins of Paxillin, Rac, PAK, and Merlin, which are all responsible for cell protrusion [8] is embedded in individual cell. This pathway is collaborated with another mechano-chemical pathway, which is responsible for transmitting mechanical cue through intercellular adhesions [9]. Our model is used to study collective cell migration using a simplified wound tissue. We then compare our simulation results to an in vitro study [10]. Finally, we discussed and made the conclusion that the mechanics between cell-substrate play a crucial role in guiding highly efficient collective cell migration. This guidance cue is well maintained and transmitted between cells through the intercellular adhesions.

2. Methods

2.1 Cell geometry

In our model, a cell in 2D $\mathrm{\Omega }\subset {\mathbb{R}}^2$ is represented as an oriented polygon including a number of boundary vertices $V_{\mathrm{\partial \Omega }}\equiv \left\{v_i\in \mathrm{\partial \Omega }\subset {\mathbb{R}}^2\right\}$, where the location of the vertex ${\mathit{v}}_i$ is denoted as ${\mathit{x}}_i$. The set of boundary vertices $V_{\mathrm{\partial \Omega }}$, together with a set of internal vertices $V_{\mathit{Int}}$ and a set of triangular elements $T_{\mathrm{\Omega }}\equiv \left\{{\tau }_{i,j,k}:{\mathit{v}}_i,{\mathit{v}}_j,{\mathit{v}}_k\in V_{\mathrm{\partial \Omega }}\cup V_{\mathrm{Int}}\right\}$ define the geometry of cell $\mathrm{\Omega }$ (Figure 1a). If two cells are closely in contact, a set of adhesive springs are generated between them (Figure 1a, red bars in the dashed blue box). There are several interior vertice on each cell boundary edge. They are evenly distributed along that edge. These interior vertice are the potential locations for newly generated adhesive spring to attach on. Any force applied on that interior vertex through the attached adhesive spring will be mapped onto its nearest end-node vertex of the corresponding boundary edge.

Figure 1.

The cell geometry and the chemical pathway between cell-substrate and intercellular adhesion. (a) the cell in our model is represented as following: The cell boundary is defined by an oriented polygon including a number of boundary vertices. A triangular mesh tiling up a cell is generated based on the method of farthest sampling [11]. The E-cadhesion type of intercellular adhesions between two neighboring cells are represented as elastic springs (red bars in the blue box, the dashed blue box is for a closer view). (b) each triangular element exhibits viscoelastic characteristic using a generalized Maxwell model following [5, 12, 13]. (c) the positive feedback loop between focal adhesion and cell protrusion is built up in each vertex of the triangular mesh following [8]. Such network includes the proteins of integrin, Paxillin, Rac, and PAK. The protein Merlin on the cadherin is also included to count the effects of intercellular adhesion on cell migration [9].

2.2 Viscoelasticity of the cell

Previous researches have demonstrated that the cell cytoskeleton exhibits viscoelastic characteristic [14, 15]. Following the studies of [16, 17], we assume that, during cell deformation and cell migration, linear viscoelasticity is adequate to describe the mechanical properties of the cell.

2.2.1 Strain and stress tensors

We use the strain tensor $\mathit{\epsilon }\left(\mathit{x},t\right)$ to describe the local cell deformation at $\mathit{x}$ at time $t$. $\mathit{\epsilon }\left(\mathit{x},t\right)$ takes the form of ${\epsilon }_{1,1}=\partial u_1/\partial x_1$, ${\epsilon }_{2,2}=\partial u_2/\partial x_2$, and ${\epsilon }_{1,2}={\epsilon }_{2,1}=\frac{1}{2}\left(\partial u_1/\partial x_2+\partial u_2/\partial x_1\right)$, where $\mathit{u}\left(\mathit{x},t\right)$ defined as ${\left(u_1\left(\mathit{x},t\right),u_2\left(\mathit{x},t\right)\right)}^T\subset {\mathbb{R}}^2$ is the displacement of $\mathit{x}$ at time $t$. We use the stress tensor $\mathit{\sigma }\left(\mathit{x},t\right)$ to describe the local forces at $\mathit{x}$ at time $t$. Here $\mathit{\sigma }$ is correlated with $\mathit{\epsilon }$ by a generalized Maxwell model: $\mathit{\sigma }\left(\mathit{x},t\right)={\mathit{\sigma }}_{\infty }\left(\mathit{x},t\right)+{\mathit{\sigma }}_m\left(\mathit{x},t\right)$ [5, 12, 13], where ${\mathit{\sigma }}_{\infty }\left(\mathit{x},t\right)$ is the stress of the long-term elastic element and ${\mathit{\sigma }}_m\left(\mathit{x},t\right)$ is the stress of the Maxwell elastic element. $E_{\infty }$, $E_m$, and ${\eta }_m$ denote the long-term elastic modulus, elastic modulus of the Maxwell elastic element, and viscous coefficient of the Maxwell viscous element, respectively (Figure 1b). The strain of the Maxwell elastic element ${\mathit{\epsilon }}_1\left(\mathit{x},t\right)$ and the strain of the viscous element ${\mathit{\epsilon }}_2\left(\mathit{x},t\right)$ sum up to the strain tensor $\mathit{\epsilon }\left(\mathit{x},t\right)$: ${\mathit{\epsilon }}_1\left(\mathit{x},t\right)+{\mathit{\epsilon }}_2\left(\mathit{x},t\right)=\mathit{\epsilon }\left(\mathit{x},t\right)$.

We assume that the total free energy of a cell is the summation of its elastic energy, its adhesion energy due to the contact with the substrate, its elastic energy due to the intercellular adhesions with neighboring cells, and its energy due to the forces exerting on the boundary.

2.2.2 Cell elastic energy

The elastic energy due to the deformation of the cell $\mathrm{\Omega }$ is given by

$$E_{\mathrm{\Omega }}\left(t\right)=\frac{1}{2}{\int }_{\mathrm{\Omega }}{\left({\mathit{\sigma }}_{\infty }\left(\mathit{x},t\right)+{\sigma }_a{\mathit{\delta }}_{\mathit{ij}}\left(\mathit{x}\right)\right)}^T\mathit{\epsilon }\left(\mathit{x},t\right)d\mathit{x}+\frac{1}{2}{\int }_{\mathrm{\Omega }}{\mathit{\sigma }}_m{\left(\mathit{x},t\right)}^T{\mathit{\epsilon }}_1\left(\mathit{x},t\right)d\mathit{x},$$E1

where ${\sigma }_a$ is a homogeneous contractile pressure following [6].

2.2.3 Cell adhesion energy due to the contact with the substrate

The energy due to the adhesion between the cell and the substrate is given by [6].

$$\frac{Y\left(\mathit{x},t\right)}{2}{\int }_{\mathrm{\Omega }}\mathit{u}{\left(\mathit{x},t\right)}^{2}d\mathit{x},$$E2

where $Y\left(\mathit{x},t\right)$ is the adhesion coefficient at time $t$ and is proportional to the strength of local focal adhesions [18]: $Y\left(\mathit{x},t\right)=\frac{n_{\mathit{x},t}}{n_0}{E}_{\mathit{st}}{Y}_a$, where $n_{\mathit{x},t}$ is the number of bound integrins at location $\mathit{x}$ at $t$ (more details of calculating $n_{\mathit{x},t}$ in Model of focal adhesion), $n_0$ is a normalized constant number, $E_{\mathit{st}}$ is the stiffness of the substrate and $Y_a$ is the basic adhesion constant following [18].

2.2.4 Cell adhesion energy due to intercellular adhesion

The energy due to the intercellular adhesions, which are modeled as elastic springs, is given by $\frac{1}{2}{\sum }_l{k}_l{\mathit{u}}_l{\left(t\right)}^2$, where $k_l$ is the spring constant of the spring $l$. Its orientation angle at time $t$ is denoted as ${\theta }_l\left(t\right)$. Its transformation vector $T\left(\theta \right)$ is denoted as $\left(cos\left({\theta }_l\right),sin\left({\theta }_l\right),-cos\left({\theta }_l\right),-sin\left({\theta }_l\right)\right)$. So ${\mathit{u}}_l\left(t\right)$ can be written as ${\mathit{u}}_l\left(t\right)=T\left({\theta }_l\right)\left({\mathit{u}}_{11}\left(t\right),{\mathit{u}}_{l2}\left(t\right)\right)$, where ${\mathit{u}}_{l1}\left(t\right)$ and ${\mathit{u}}_{l2}\left(t\right)$ are the displacements of the two end-node vertice ${\mathit{x}}_1$ and ${\mathit{x}}_2$ of $l$ at time $t$. The elastic force of $l$ due to displacement of $\mathrm{\Delta }l$ is applied on ${\mathit{x}}_i$ and ${\mathit{x}}_j$ as ${\mathit{f}}_l=f\left(\mathrm{\Delta }l\right){\mathit{e}}_l$ and -$f\left(\mathrm{\Delta }l\right){\mathit{e}}_l$, respectively, where $f\left(\mathrm{\Delta }l\right)$ is the magnitude of ${\mathit{f}}_l$, ${\mathit{e}}_l$ is the unit vector of the orientation of $l$.

2.2.5 Boundary and protrusion forces

Furthermore, the local forces applied on the cell boundary also contribute to the energy. Following [6], the tension force along the cell boundary is considered. In addition, we also incorporate the protrusion force on the leading edge of migrating cell. The contribution of these two forces can be written as

$${\int }_{\mathrm{\partial \Omega }}\left(\mathit{\lambda }\left(\mathit{x},t\right)+\mathit{f}\left(\mathit{x},t\right)\right)\mathit{u}\left(\mathit{x},t\right)d\mathit{x},$$E3

where $\mathit{\lambda }\left(\mathit{x},t\right)$ is the line tension force and $\mathit{f}\left(\mathit{x},t\right)$ is the protrusion force. Line tension force is written as $\mathit{\lambda }\left(\mathit{x},t\right)=-f_m\kappa \left(\mathit{x},t\right)\mathit{n}\left(\mathit{x},t\right)$, where $f_m$ is a contractile force per unit length, $\kappa \left(\mathit{x},t\right)$ is the curvature, and $\mathit{n}\left(\mathit{x},t\right)$ is the outward unit normal at $\mathit{x}$ at time $t$ [6]. Protrusion force is denoted as $f\left(\mathit{x},t\right)=-f_a\mathit{n}\left(\mathit{x},t\right)$, where $f_a$ is the protrusion force per unit length.

2.2.6 The total free energy and its dissipation

In summary, the total free energy of cell $\mathrm{\Omega }$ at time $t$ is given by

$$\begin{array}{l}{E}_{\mathrm{\Omega }}\left(t\right)=\frac{1}{2}{\int }_{\mathrm{\Omega }}\left({\left({\mathit{\sigma }}_{\infty }\left(\mathit{x},t\right)+{\mathit{\sigma }}_a{\delta }_{\mathit{ij}}\right)}^T\mathit{\epsilon }\left(\mathit{x},t\right)+{\mathit{\sigma }}_m^T\left(\mathit{x},t\right){\mathit{\epsilon }}_1\left(\mathit{x},t\right)\right)d\mathit{x}+\frac{1}{2}{\int }_{\mathrm{\Omega }}\left({\sigma }_a,{\sigma }_a,0\right)\mathit{\epsilon }\left(\mathit{x},t\right)d\mathit{x}\\ +\frac{Y\left(\mathit{x},t\right)}{2}{\int }_{\mathrm{\Omega }}\mathit{u}{\left(\mathit{x},t\right)}^{2}d\mathit{x}+\frac{1}{2}\sum _l{k}_l{\mathit{u}}_l{\left(t\right)}^2+{\int }_{\mathrm{\partial \Omega }}\left(\mathit{\lambda }\left(\mathit{x},t\right)+\mathit{f}\left(\mathit{x},t\right)+{\mathit{f}}_l\left(\mathit{x},t\right)\right)\mathit{u}\left(\mathit{x},t\right)d\mathit{x}).\end{array}$$E4

The energy dissipation of $E_{\mathrm{\Omega }}\left(t\right)$ due to cell viscosity is determined by the viscous coefficient ${\eta }_m$ and the strain of the viscous element ${\mathit{\epsilon }}_2\left(\mathit{x},t\right)$ [19]: $-{\int }_{\mathrm{\Omega }}{\eta }_m{\left(\frac{\partial {\mathit{\epsilon }}_2}{\partial t}\right)}^{2}d\mathit{x}$. The dissipation of the total free energy of the cell can be written as

$$\begin{array}{l}\frac{\partial }{\partial t}{E}_{\mathrm{\Omega }}\left(t\right)=\frac{\partial }{\partial t}(\frac{1}{2}{\int }_{\mathrm{\Omega }}\left({\left({\mathit{\sigma }}_{\infty }\left(\mathit{x},t\right)+{\mathit{\sigma }}_a{\delta }_{\mathit{ij}}\right)}^T\mathit{\epsilon }\left(\mathit{x},t\right)+{\mathit{\sigma }}_m^T\left(\mathit{x},t\right){\mathit{\epsilon }}_1\left(\mathit{x},t\right)\right)d\mathit{x}+\frac{1}{2}{\int }_{\mathrm{\Omega }}\left({\sigma }_a,{\sigma }_a,0\right)\mathit{\epsilon }\left(\mathit{x},t\right)d\mathit{x}\\ +\frac{Y\left(\mathit{x},t\right)}{2}{\int }_{\mathrm{\Omega }}\mathit{u}{\left(\mathit{x},t\right)}^{2}d\mathit{x}+\frac{1}{2}\sum _l{k}_l{\mathit{u}}_l{\left(t\right)}^2+{\int }_{\mathrm{\partial \Omega }}\left(\mathit{\lambda }\left(\mathit{x},t\right)+\mathit{f}\left(\mathit{x},t\right)+{\mathit{f}}_l\left(\mathit{x},t\right)\right)\mathit{u}\left(\mathit{x},t\right)d\mathit{x})\\ =-{\int }_{\mathrm{\Omega }}{\eta }_m{\left(\frac{\partial {\mathit{\epsilon }}_2\left(\mathit{x},t\right)}{\partial t}\right)}^{2}d\mathit{x}.\end{array}$$E5

Since ${\eta }_m\frac{\partial {\mathit{\epsilon }}_2}{\partial t}=E_m{\mathit{\epsilon }}_1$ and ${\mathit{\sigma }}_{\infty }=E_{\infty }\mathit{\epsilon }$, and ${\mathit{u}}_l\left(t\right)=T\left({\theta }_l\right)\left({\mathit{u}}_{l1},{\mathit{u}}_{l2}\right)$. Eq. (5) can be rewritten as

$$\begin{array}{l}{\int }_{\mathrm{\Omega }}{E}_{\infty }\mathit{\epsilon }\left(\mathit{x},t\right)\frac{\partial \mathit{\epsilon }\left(\mathit{x},t\right)}{\partial t}d\mathit{x}+{\int }_{\mathrm{\Omega }}{\mathit{\sigma }}_m^T\left(\mathit{x},t\right)\frac{\partial \mathit{\epsilon }\left(\mathit{x},t\right)}{\partial t}d\mathit{x}+{\int }_{\mathrm{\Omega }}\left({\sigma }_a,{\sigma }_a,0\right)\frac{\partial \mathit{\epsilon }\left(\mathit{x},t\right)}{\partial t}d\mathit{x}\\ +Y{\int }_{\mathrm{\Omega }}\mathit{u}\left(\mathit{x},t\right)\frac{\partial \mathit{u}\left(\mathit{x},t\right)}{\partial t}d\mathit{x}+\sum _{l}T{\left({\theta }_l\right)}^{T}T\left({\theta }_l\right)k_l\left({\mathit{u}}_{l1}\left(t\right),{\mathit{u}}_{l2}\left(t\right)\right)\frac{\partial \mathit{u}\left(\mathit{x},t\right)}{\partial t}\\ +{\int }_{\mathrm{\partial \Omega }}\left(\mathit{\lambda }\left(\mathit{x},t\right)+\mathit{f}\left(\mathit{x},t\right)+{\mathit{f}}_l\left(\mathit{x},t\right)\right)\frac{\partial \mathit{u}\left(\mathit{x},t\right)}{\partial t}\mathit{dx}=0.\end{array}$$E6

Denoting $\mathit{B}=\left(\begin{array}{cc}\partial /\partial x_1& 0\\ 0& \partial /\partial x_2\\ \partial /\partial x_2& \partial /\partial x_1\end{array}\right)$, then $\mathit{\epsilon }\left(\mathit{x},t\right)=\mathit{B}\mathit{u}\left(\mathit{x},t\right)$. According to Gauss’ divergence theorem, we rewrote ${\int }_{\mathrm{\Omega }}\left({\sigma }_a,{\sigma }_a,0\right)\mathit{B}\frac{\partial \mathit{u}\left(\mathit{x},t\right)}{\partial t}$ as ${\int }_{\mathrm{\Omega }}{\sigma }_a\nabla \cdot \frac{\partial \mathit{u}\left(\mathit{x},t\right)}{\partial t}d\mathit{x}$, which leads to ${\int }_{\mathrm{\partial \Omega }}{\sigma }_a\mathit{n}\left(\mathit{x},t\right)\frac{\partial \mathit{u}\left(\mathit{x},t\right)}{\partial t}d\mathit{x}$.

Denoting ${\mathit{A}}_l=\left(\begin{array}{cccc}{c}_l^2& c_l{s}_l& -c_l^2& -c_l{s}_l\\ c_l{s}_l& s_l^2& -c_l{s}_l& -s_l^2\\ -c_l^2& -c_l{s}_l& c_l^2& c_l{s}_l\\ -c_l{s}_l& -s_l^2& c_l{s}_l& s_l^2\end{array}\right)$, where $c_l=\mathit{cos}\left({\theta }_l\right)$ and $s_l=\mathit{sin}\left({\theta }_l\right)$. Eq. (6) can be rewritten as

$$\begin{array}{l}{\int }_{\mathrm{\Omega }}{\mathit{B}}^T\mathit{\sigma }{\left(\mathit{x},t\right)}^{T}d\mathit{x}+Y{\int }_{\mathrm{\Omega }}\mathit{u}\left(\mathit{x},t\right)d\mathit{x}+\sum _l{\mathit{A}}_l{k}_l\left({\mathit{u}}_{l1}\left(t\right),{\mathit{u}}_{l2}\left(t\right)\right)=\\ -{\int }_{\mathrm{\partial \Omega }}{\sigma }_a\mathit{n}\left(\mathit{x},t\right)+\mathit{\lambda }\left(\mathit{x},t\right)+\mathit{f}\left(\mathit{x},t\right)+{\mathit{f}}_l\left(\mathit{x},t\right)d\mathit{x}\end{array}$$E7

2.2.7 Stress of viscoelastic cell and its update

By using general Maxwell model, the stress $\mathit{\sigma }\left(\mathit{x},t\right)$ can be written as [12, 13, 19]:

$$E_{\infty }\mathit{\epsilon }\left(\mathit{x},t\right)+{\int }_0^t{E}_m{e}^{\frac{t-s}{{\eta }_m/E_{\infty }}}\frac{\partial \mathit{\epsilon }\left(\mathit{x},s\right)}{\partial s}\mathit{ds}={\mathit{\sigma }}_{\infty }\left(\mathit{x},t\right)+{\mathit{\sigma }}_m\left(\mathit{x},t\right).$$E8

During the time interval $\mathrm{\Delta }t=t_{n+1}-t_n$, where $t_n$ is the $n$-th time step, ${\mathit{\sigma }}_m^{n+1}\left(\mathit{x}\right)$ can be written as [19]:

$$e^{-\frac{\mathrm{\Delta }t}{{\eta }_m/E_{\infty }}}{\mathit{\sigma }}_m^n\left(\mathit{x}\right)+\frac{E_m}{E_{\infty }}{\int }_{t_n}^{t_{n+1}}{e}^{-\frac{t_{n+1}-s}{{\eta }_m/E_{\infty }}}\mathit{ds}\frac{{\mathit{\sigma }}_{\infty }^{n+1}\left(\mathit{x}\right)-{\mathit{\sigma }}_{\infty }^n\left(\mathit{x}\right)}{\mathrm{\Delta }t}.$$E9

Therefore, the stress ${\mathit{\sigma }}^n\left(\mathit{x}\right)$ at $t_n$ can be written as

$${\mathit{\sigma }}^n\left(\mathit{x}\right)={\mathit{\sigma }}_{\infty }^n\left(\mathit{x}\right)+{\mathit{\sigma }}_m^n\left(\mathit{x}\right)$$E10

2.2.8 Force balance equation for discretized time step

For each triangular element ${\mathit{\tau }}_{i,j,k}$, Eq. (7) at time step $t_{n+1}$ can be rewritten using Eq. (10) as

$$\begin{array}{l}{\int }_{{\mathit{\tau }}_{i,j,k}}{\mathit{B}}^T\left(E_{\infty }{\mathit{Bu}}^{n+1}\left(\mathit{x}\right)+e^{-\frac{\mathrm{\Delta }t}{{\eta }_m/E_{\infty }}}{\mathit{\sigma }}_m^n+{\gamma }_m{A}_m\left(E_{\infty }{\mathit{Bu}}^{n+1}\left(\mathit{x}\right)-E_{\infty }{\mathit{Bu}}^n\left(\mathit{x}\right)\right)\right)d\mathit{x}\\ +Y{\int }_{{\mathit{\tau }}_{i,j,k}}{\mathit{u}}^{n+1}d\mathit{x}+\sum _{l\in {\mathit{\tau }}_{i,j,k}}{\mathit{A}}_l{k}_l\left({\mathit{u}}_{l1}^{n+1},{\mathit{u}}_{l2}^{n+1}\right)={\mathit{F}}^{n+1}\left(\mathit{x}\right),\end{array}$$E11

where ${\gamma }_m=E_m/E_{\infty }$, $A_m=\frac{1-e^{-\frac{\mathrm{\Delta }t}{{\eta }_m/E_{\infty }}}}{\frac{\mathrm{\Delta }t}{{\eta }_m/E_{\infty }}}$, and ${\mathit{F}}^{n+1}=-{\int }_{\mathrm{\partial \Omega }}{\sigma }_a\mathit{n}\left(\mathit{x},t\right)+\mathit{\lambda }\left(\mathit{x},t\right)+\mathit{f}\left(\mathit{x},t\right)+{\mathit{f}}_l\left(\mathit{x},t\right)d\mathit{x}$. Eq. (11) leads to the following linear force-balance equation

$${\mathit{K}}_{{\mathit{\tau }}_{i,j,k}}^{n+1}{\mathit{u}}_{{\mathit{\tau }}_{i,j,k}}^{n+1}={\mathit{f}}_{{\mathit{\tau }}_{i,j,k}}^{n+1},$$E12

where ${\mathit{K}}_{{\mathit{\tau }}_{i,j,k}}^{n+1}$, ${\mathit{u}}_{{\mathit{\tau }}_{i,j,k}}^{n+1}$, and ${\mathit{f}}_{{\mathit{\tau }}_{i,j,k}}^{n+1}$ are the stiffness matrix, displacement vector, and integrated force vector of ${\mathit{\tau }}_{i,j,k}$ at time step $t_{n+1}$ (see more details of derivation of (Eq. 12) in [11]).

We can then assemble the element stiffness matrices of all triangular elements into one big global stiffness matrix ${\mathit{K}}^{n+1}$. Therefore, the linear relationship between the concatenated displacement vector ${\mathit{u}}^{n+1}$ of all cell vertice and the force vector ${\mathit{f}}^{n+1}$ on them is given by

$${\mathit{K}}^{n+1}{\mathit{u}}^{n+1}={\mathit{f}}^{n+1}.$$E13

Changes in the cell shape at time step $t_{n+1}$ can be obtained by solving Eq. (13). For vertex ${\mathit{v}}_i$ at ${\mathit{x}}_i$, its new location at next time step is then updated as ${\mathit{x}}_i^{n+1}={\mathit{x}}_i^n+{\mathit{u}}^{n+1}\left({\mathit{v}}_i\right)$.

2.3 Mechano-chemical pathway in the cell

Upon contact with the environment, cells can transfer the mechanical cues into biochemical signals, which can trigger the initiation of further cellular behaviors [20]. In our model, we considered a mechano-chemical pathway consisting of two parts, where one is to regulate the feedback loop between focal adhesion and cell protrusion and the other is to regulate the transmission of mechanical signal between adjacent cells through intercellular adhesions.

2.3.1 Model of focal adhesion

For each vertex ${\mathit{v}}_i$ in cell $\mathrm{\Omega }$, we assign a constant number of integrin ligand on it. These integrin molecules can bind or unbind with fibronectin molecules on the substrate underneath. Following [21], the numbers of bound and unbound integrin ligand molecules are determined by

$$\frac{{\mathit{dR}}_b}{\mathit{dt}}=k_f{n}_s{R}_u-k_r{R}_b,$$E14

where $R_u$ and $R_b$ are the numbers of unbound and bound integrin ligand, respectively; $k_f$ is the binding rate coefficient; $n_s$ is the concentration of fibronectin per cell vertex; $k_r$ is the unbinding rate coefficient. $k_r$ depends on the magnitude of the traction force ${\mathit{f}}_r$ applied on ${\mathit{v}}_i$. Traction force ${\mathit{f}}_r\left(\mathit{x},t\right)$ on $\mathit{x}$ at time $t$ is given by $Y\left(\mathit{x},t\right)\mathit{u}\left(\mathit{x},t\right)$ following [6]. $k_r$ is determined by $k_r=k_{r_0}\left(e^{-0.04f_r}+4e-7e^{0.2f_r}\right)$ following [22], where $k_{r_0}$ is a constant. $k_f$ is related with the substrate stiffness by $k_f=k_{f0}{E}_{\mathit{st}}^2/\left(E_{\mathit{st}}^2+E_{\mathit{st}0}^2\right)$ [16, 23], where $k_{f0}$ and $E_{\mathit{st}0}$ are constants (see Appendix for choosing $E_{\mathit{st}0}$).

2.3.2 Model of feedback loop between focal adhesion and cell protrusion

We introduced a simplified model of a positive feedback loop to control the spatial distribution of the focal adhesions, which governs the direction of cell protrusion [8, 24]. In our model, this feedback loop involves proteins of Paxillin, Rac, and PAK (Figure 1c). Upon formation of focal adhesion, Paxillin is activated by active PAK. The active Paxillin then activates Rac, which in turn triggers the activation of PAK. The activated Rac is responsible for protruding cells [25]. Since the protein Merlin on the intercellular cadherin complex also plays a role in activating Rac [9], we include Merlin in our feedback loop. The protein concentration over time is updated through a set of differential equations following [8]:

$$\frac{\mathit{dr}}{\mathit{dt}}=k_{x,r}\left(k_m\frac{C_r^2}{C_r^2+m^2}x-r\right)$$E15

$$\frac{\mathit{dp}}{\mathit{dt}}=k_{r,p}\left(r-p\right)$$E16

$$\frac{\mathit{dx}}{\mathit{dt}}=k_{p,x}\left(k_x\frac{p^2}{p^2+C_x^2}n-x\right)$$E17

where $x$, $r$, $p$, $m$ and $n$ are the concentrations of activated Paxillin, activated Rac, activated PAK, Merlin and bound integrins, respectively. $k_{x,r}$, $k_{r,p}$, $k_{p,x}$, $k_m$, $k_x$, $C_r$, $C_x$ are the parameters of corresponding rates. The level of activated Rac was used to determine the protrusion force on the leading edge of the migrating cell (see details of cell protrusion model in Appendix).

2.3.3 Model of mechanosensing through intercellular adhesion

We added Merlin in the feedback loop (Figure 1c) following a previous study reporting that Merlin on the intercellular cadherin complex regulates the Rac activity [26]. As illustrated in Figure 2a, for two adjacent cells $C_1$ and $C_2$ where $C_1$ is the leader cell and $C_2$ is the follower cell, if both cells are at static state, Merlin molecules only locate on the cadherin spring (Figure 2b top). As reported by [9], Merlin suppressed the binding of integrin. Due to such suppression, Rac turns to inactivated on the Merlin-expressed site. Once cell $C_1$ starts to migrate, tension force is generated on the cadherin spring between $C_1$ and $C_2$. Merlin is therefore delocalized from cadherin-attached site in response to the generated tension force. As a consequence, Rac is activated there to generate protrusion force to follow the leader cell (Figure 2b bottom [26]). For simplicity, we introduced the inactive Merlin phenotype along with the active Merlin phenotype on the two end vertice of one intercellular cadherin adhesion. The negative feedback loop of Merlin-Rac is modeled through a set of differential equations following [9, 26], where active Merlin and inactive Merlin can switch their phenotype, but only active Merlin can suppress the Rac activity (Figure 1c). The delocalization of Merlin was simply modeled as Merlin switching to inactive phenotype:

Figure 2.

Mechanosensing through intercellular adhesion and tissue model for collective cell migration. (a) the intercellular adhesion of the cadherin spring (red springs in the blue box) are responsible for transmitting mechanical stimulus from leader cell ($C_1$) to its follower cells ($C_2$ and $C_3$). (b) When cells are static, Merlin which inhibits the Rac activation is bound on the cadherin spring. Once leader cell migrates, stretch is generated on the cadherin spring, Merlin on the follower cell is delocated. Therefore, Rac is activated on the follower cell. (c) the size of the wound epithelial tissue is 720 $\mathit{\mu m}$$\times$ 240 $\mathit{\mu m}$. The right boundary is set as the wound edge (yellow line). Cells can migrate towards the open space on the right. Three measurements are introduced to measure the collective cell migration: (1) migration persistence $p\left(t_n\right)$, the ratio of the distance from the current position at time $t_n$ to its initial position (green line), over the length of the traversed path (red curve); (2) normalized pair separation distance $d_{i,j}\left(t_n\right)$ is the separation distance between a pair of cells at time $t_n$ (green lines) divided by the average length of the two cells’ traversed path (red curves); (3) migration direction angle $\alpha \left(t_n\right)$ is the angle between the migration direction (red arrow) and the direction towards the wound (green arrow).

$$\frac{\mathit{dm}}{\mathit{dt}}=e-\delta \left(f_t>f_{t-\mathit{thr}}\right)k_{m,e}(k_p\left(\frac{p^2}{p^2+C_e^2}+k_e\right)m$$E18

$$\frac{\mathit{de}}{\mathit{dt}}=\delta \left(f_t>f_{t-\mathit{thr}}\right)k_{m,e}\left(k_p\left(\frac{p^2}{p^2+C_e^2}+k_e\right)m\right)-e$$E19

where $m$, $e$ and $p$ are the concentrations of Merlin, inactive Merlin, and PAK, respectively. $k_{m,e}$, $k_p$, $k_e$, $C_e$ are the corresponding rate parameters. $\delta \left(x\right)$ is a kronecker function that $\delta \left(\mathit{\text{TRUE}}\right)=1$ and $\delta \left(\mathit{\text{FALSE}}\right)=0$. $f_t$ is the tension force through the cadherin spring and $f_{t-\mathit{thr}}$ is a force threshold.

2.4 Cellular tissue for collective cell migration

In our model, the collective cell migration was modeled using a wound tissue of epithelial cells. The tissue size is 720 $\mathit{\mu m}$$\times$ 240 $\mathit{\mu m}$. The epithelial cell type is set to MCF-10A, which is used in the in vitro study [10]. The corresponding epithelial-specific parameters can be found in Table 1. We arbitrarily set the right boundary of the tissue as the wound edge, and cells can migrate towards the open space to the right of the wound edge (Figure 2c). The mechano-chemical pathway was initiated first in the cells on the wound edge after they migrate. We followed a previous study [10] to divide the location of cells into four sub-regions according to their distance to the wound edge: Regions $\mathrm{I}$, $\mathrm{II}$, $\mathrm{III}$, and $\mathrm{IV}$ whose distance to the wound edge is 0–160 $\mathit{\mu m}$, 160–320 $\mathit{\mu m}$, 320–480 $\mathit{\mu m}$, and 480–640 $\mathit{\mu m}$, respectively (Figure 2c). We ran the simulation for 12 biological hours, the same experimental duration time in the in vitro study [10].

Definition	Value	Reference
Time step lapse	0.1 $\mathit{sec}$	NA^a
Cell radius	10 $\mathit{\mu m}$	[27]
Young’s modulus of cell	5 $\mathit{kPa}$	[28]
Poisson ratio of cell	0.40	[29]
Contractile pressure ${\sigma }_a$	2 $\mathit{kPa}$	[6]
Adhesion energy constant $Y_a$	0.9 $/\mathit{\mu m}$	[6]
Spring constant of cadherin spring	3.0 $\mathit{nN}/\mathit{\mu m}$	NA
Default length of cadherin spring	100 $\mathit{nm}$	[30]
Maximum length of cadherin spring	400 $\mathit{nm}$	[31]
Protrusion force constant $f_a$	2.0 $\mathit{nN}/\mathit{\mu m}$	[11]
Integrin bound rate $k_{f0}$	0.5	[22]
Integrin unbound rate $k_{r0}$	0.4	[22]
Reference substrate stiffness $E_{\mathit{st}0}$	40 kPa	NA
Rac deactivation rate $k_{x,r}$	4/min	[8]
PAK deactivation rate $k_{r,p}$	10/min	[8]
Paxillin dephosphorylation rate $k_{p,x}$	10/min	[8]
Saturation of phosphorylated Paxillin $k_x$	1	[8]
Saturation of PAK activation $k_p$	1	[8]
Saturation of Merlin $k_m$	1	NA
Merlin phosphorylation rate $k_{m,e}$	10/min	NA
Saturation of phosphorylated Merlin $k_e$	1 $e^{-3}$	NA
Force threshold of delocating Merlin $f_{t-\mathit{thr}}$	0.15 $\mathit{nN}$	NA

Table 1.

Parameters used in the model.

^aEstimated value marked as NA.

2.5 Measurements of the collective cell migration

In our model, the collective cell migration are measured using four measurements following [10]:

The migration persistence. At time step $t_n$, the length of the straight line between cell positions at $t_n$ and initial time step $t_0$ over the length of the migrating trajectory:

$$p\left(t_n\right)=\frac{\mid \mathit{x}\left(t_n\right)-\mathit{x}\left(t_0\right)\mid }{{\sum }_{k=0}^{n-1}\mid \mathit{x}\left(t_{k+1}\right)-\mathit{x}\left(t_k\right)\mid },$$E20

where $t_0$ is the initial time, $\mathit{x}\left(t_i\right)$ is the position at time step $t_i$ (Figure 2c.1).

The normalized separation distance. At time step $t_n$, the separation distance of a pair of two adjacent cells $1$-$2$, divided by the average length of their migrating trajectories:

$$d_{1,2}\left(t_i\right)=\frac{\Vert {\mathit{x}}_1\left(t_i\right)-{\mathit{x}}_2\left(t_i\right)\mid -\mid {\mathit{x}}_1\left(t_0\right)-{\mathit{x}}_2\left(t_0\right)\Vert }{\frac{1}{2}\left({\sum }_{k=0}^{i-1}|{\mathit{x}}_1\left(t_{k+1}\right)-{\mathit{x}}_1\left(t_k\right)|+{\sum }_{k=0}^{i-1}|{\mathit{x}}_2\left(t_{k+1}\right)-{\mathit{x}}_2\left(t_k\right)|\right)},$$E21

where the numerator is the separation distance between cells $i$ and $j$ at time $t_n$, and the denominator is the average path length of cells $i$ and $j$ at time $t_n$ (Figure 2c.2).

The direction angle. The angle between the cell migration direction and the direction towards the wound.

$$\alpha \left(t_n\right)=arccos\left({\mathit{u}}_c\cdot {\mathit{u}}_w\right)\cdot sgn\left(\parallel {\mathit{u}}_c\times {\mathit{u}}_w\parallel \right),$$E22

where ${\mathit{u}}_c$ is the unit vector of the cell migration direction, ${\mathit{u}}_w$ is the unit vector of direction from the cell mass center towards the wound, $sgn\left(x\right)$ is the sign of $x$ (Figure 2c.3).

3. Results

3.1 Mechanics of cell-substrate is crucial to regulate collective cell migration

3.1.1 Morphology and migration pattern under different substrate stiffness

We fist studied collective cell migration under the mechano-chemical mechanism. The trajectories of our simulation showed that cells migrate faster and more persistently on stiffer substrate (Figure 3a and b). This is compatible with the observed pattern from the in vitro study of collective cell migration (Figure 3c and d). Furthermore, the shape of cell also changes with different substrate stiffness. Cells adopted a more spherical shape on softer substrate (Figure 4a) while cells were more elongated on stiffer substrate (Figure 4c). The same pattern of cell morphology was observed in [10], where cell extended its protrusions in all directions on softer substrate (Figure 4b) while cell protruded only on the leading edge with a long tail on stiffer substrate (Figure 4d).

Figure 3.

Cell migrating trajectories. (a–b) the migrating trajectory in our simulation using two substrate stiffness: 3 and 65 kPa. (c–d) the migrating trajectory from the in vitro study using the same substrate stiffness: 3 and 65 kPa [10]. The scale bar is 100 $\mathit{\mu m}$.

Figure 4.

Cell morphology. (a, c) the cell morphology in our simulation using two substrate stiffness: 3 and 65 kPa. (b, d) the cell morphology from the in vitro study [10] using the same substrate stiffness: 3 and 65 kPa. The cell boundary is highlighted in black.

Overall, the patterns of cell trajectory and cell morphology of our simulation are consistent with that from in vitro study. This indicated that our mechano-chemical model is valid.

3.1.2 The mechanical signal has long-distance impacts on collective cell migration

We then quantified the cell migration to explore the role of mechanics of cell-substrate and cell–cell mechanics on collective cell migration using the three measurements: persistent ratio $p\left(t_n\right)$, normalized separation distance $d_{i,j}\left(t_n\right)$ and direction angle $\alpha \left(t_n\right)$.

We first examined the migrating speed of the cell. In general, cells migrate with higher speed on stiffer substrate (Figure 5a, more details of cell migration speed can be found in Appendix). In addition, cells close to the wound edge migrated with higher speed on both stiffer and softer substrate. This speed decreased gradually as the distance to the wound edge increased. On substrate with stiffness of 65 kPa, the migration speed decreased from $0.69\pm 0.01\mathit{\mu m}/\mathit{min}$ in Region $\mathrm{I}$ to $0.49\pm 0.02\mathit{\mu m}/\mathit{min}$ in Region $\mathrm{IV}$, while on substrate with stiffness of 3 kPa, the migration speed decreased from $0.38\pm 0.02\mathit{\mu m}/\mathit{min}$ in Region $\mathrm{I}$ to $0.25\pm 0.02\mathit{\mu m}/\mathit{min}$ in Region $\mathrm{IV}$ (Figure 5a). The cell migration speed of our simulation was consistent with that from the in vitro study [10]. It is easy to interpret such pattern of cell migration speed. For cells in Region $\mathrm{I}$, especially on the wound edge, there are fewer or even no cells ahead. As the distance to the wound edge increased, it was more crowded and more difficult for cells to migrate forward.

Figure 5.

Measurements of the collective cell migration. (a–c) the cell migration speed, persistence ratio and normalized separation distance of our simulation and in the in vitro study [10]. (d–g) the migration direction angle of the cells on the leading edge and more than 500 $\mathit{\mu m}$ from the wound edge in our simulation and that in the in vitro study [10] on the substrate with stiffness 65 kPa (d–e) and 3 kPa (f–g). The colors indicating simulation and experiment are shown in (a). The error bars of our simulation depict the standard deviations of four runs of simulation.

We next examined the migration persistence of the cells. As shown in Figure 5b, cells migrate more persistently on stiffer substrate. In addition, cells close to the wound edge migrated with higher migration persistence. For cells on substrate with stiffness of 65 kPa, the persistence ratio decreased from $82\pm 2\%$ in Region $\mathrm{I}$ to $58\pm 3\%$ in Region $\mathrm{IV}$, while for cells on substrate with stiffness of 3 kPa, the persistence ratio decreased from $71\pm 1\%$ in Region $\mathrm{I}$ to $55\pm 3\%$ in Region $\mathrm{IV}$ (Figure 5b). As shown in Figure 5c, collective cell migration was coordinated better on stiffer substrate.

In addition, we examined the normalized separation distance of the pairs of migrating cells. As shown in Figure 5c, the normalized separation distance increased as the distance to the wound edge increased. In our simulation, for cells on substrate at stiffness of 65 kPa, the separation distance decreased from $0.15\pm 0.02$ in Region $\mathrm{I}$ to $0.11\pm 0.02$ in Region $\mathrm{II}$ and then increased to $0.21\pm 0.03$ in Region $\mathrm{IV}$, while for cells on substrate at stiffness of 3 kPa, the separation distance decreased from $0.22\pm 0.02$ in Region $\mathrm{I}$ to $0.17\pm 0.02$ in Region $\mathrm{II}$ and then increased to $0.19\pm 0.04$ in Region $\mathrm{IV}$ (Figure 5c). This pattern of separation distance in our simulation was also observed in the in vitro study [10].

Furthermore, we examined the migration direction angle. We compared this angle for cells on the leading edge of the tissue and cells 500 $\mathit{\mu m}$ away. Since the cell migration direction is usually along the cell polarity direction [32], we also compared this direction angle to the cell polarization direction reported in [10]. As shown in Figure 5d–g, cells exhibit more accurate migration direction towards the wound on stiffer substrate (65 kPa). Only about 10 $\%$ of the cells on the leading edge had migration direction opposite to the wound (Figure 5d, 90°–270°). For cells $>$ 500 $\mathit{\mu m}$ away from the wound edge, 30 $\%$ of them had migration direction opposite to the wound (Figure 5f, 90°–270°). However, for cells on softer substrate (3 kPa), cell migration deviated more from the direction towards the wound where 35 $\%$ of the cells on leading edge had migration direction opposite to the wound direction (Figure 5e, 90°–270°), while for cells $>$ 500 $\mathit{\mu m}$ away from the wound edge, this fraction increased to 45 $\%$ (Figure 5g, 90°–270°).

These measurements implied that substrate stiffness is important to guide collective cell migration. Cells on stiffer substrate can migrate with high persistence, good coordination between cell pairs, and accurate migration direction. Our simulation suggests that the mechano-chemical feedback loop in each cell ensured it to dictate its migration direction. Furthermore, the individual cell movements were organized into a global migrative wave through intercellular adhesions.

4. Conclusions

In this chapter, we introduced our novel finite element cellular model to explore the mechanism behind collective cell migration using a simplified tissue model. This model includes a detailed mechano-chemical feedback loop, which takes into account of formation of focal adhesion and cell protrusion initiated by Rac signaling. In addition, our model incorporates the mechanical cue transmitting between the follower cell and the leader cell. We further examined the effects of cell-substrate contact and intercellular adhesions on collective cell migration.

An important result of this study is that we find the cell-substrate mechanics plays crucial role in guiding collective cell migration with higher persistence, more accurate direction, and better coordination between cell pairs (Figure 5). Previous in vitro study has shown that cells tend to have elongated shape on stiffer substrate while cells tend to have spherical shape on softer substrate [33]. This is compatible with our simulation (Figure 4a and c). We anticipate that our finite element cellular model can be applied to a broad of studies of cellular tissue problems.

In our model, cell migration is initiated and maintained by the protrusion force on the leading edge and the cell migration speed varies with the cell-substrate friction following [16].

A.1 Cell-substrate depending on substrate stiffness

The adhesion coefficient $Y\left(\mathit{x},t\right)$ of a cell vertex $\mathit{x}$ at time $t$ and set to be proportional to the strength of focal adhesions [18]: $Y\left(\mathit{x},t\right)=\frac{n_{\mathit{x},t}}{n_0}{E}_{\mathit{st}}{Y}_a$, where $n_{x,t}$ is the number of binding integrins at location $\mathit{x}$ at $t$, $n_0$ is a normalizing constant number, $E_{s}t$ is the stiffness of the substrate, $Y_a$ is the basic adhesion constant taken from [18]. In this way, the cell-substrate friction is related with the stiffness of the substrate.

A.2 Cell protrusion depends on substrate stiffness

In our model, there is a mechano-chemical pathway dictating the cell protrusion. The bound integrin initiates the activation of Rac which regulates the cell protrusion. At time $t$, the migration direction of the cell $C$ is sampled from all the boundary vertice according to their Rac concentration. One vertex ${\mathit{v}}_i$ is stochastically selected with the probability $\frac{\mathit{Rac}\left({\mathit{v}}_i\right)}{{\sum }_i\mathit{Rac}\left({\mathit{v}}_i\right)}$. The outward unit normal vector $\mathit{n}\left({\mathit{v}}_i\right)$ of ${\mathit{v}}_i$ is chosen as the cell migration direction. Any vertex ${\mathit{v}}_j$ whose outward unit normal vector $\mathit{n}\left({\mathit{v}}_j\right)$ is positively aligned with $\mathit{n}\left({\mathit{v}}_i\right)$, is treated as leading edge vertex. The protrusion force ${\mathit{f}}_a$ is then applied on each leading edge vertex ${\mathit{v}}_i$ as ${\mathit{f}}_a\left({\mathit{v}}_i\right)=f_{a}R\left({\mathit{v}}_i\right)\mathit{n}\left({\mathit{v}}_i\right)$, where $f_a$ is a constant, $R\left({\mathit{v}}_i\right)$ is the normalized Rac concentration at ${\mathit{v}}_i$.

A.3 Calibrating the cell protrusion parameter

As shown in Figure 6a–c, the cell leading edge has higher level of bound integrin, along with higher level of Rac due to the effect of positive feedback loop. The Merlin expression is also mechano-dependent. As shown of the pair of cells in the green box of Figure 6b, after the right cell migrates, the stretch force on the cadherin spring between them make the Merlin delocate from the left cell. As a result, the left cell can express Rac to protrude following the right cell. If the pair of static cell are simply in contact (Figure 6b), the Merlin is expressed on both of them. Therefore, the Rac expression is inhibited. Both of the two cells do not protrude against each other. To fit our cell protrusion model to the in vitro data, we calibrate the parameter of $E_{\mathit{st}0}$: when $E_{\mathit{st}0}=40\mathit{kPa}$, the cell migration speed of our simulation has the best match with the in vitro studies [33, 34] (Figure 6d).

Figure 6.

The cell protrusion depends on mechano-chemical process. (a–c) the spatial distribution of the normalized concentration of bound integrin, Merlin, and Rac. The black arrows indicate the migration direction. The pattern of Merlin expression depends on cell status. Green box in (b): The left cell follows the right one. Merlin is expressed only on the right cell; Orange box in (b): The two static cells are in contact. Merlin is expressed on both of them. (d) Cell migration speed of our simulation is consistent with the experimental observation [33, 34].