A Formalism for Modelling Traction forces and Cell Shape Evolution during Cell Migration in Various Biomedical Processes
AA Formalism for Modelling Traction forces andCell Shape Evolution during Cell Migration in
Various Biomedical Processes
Q. Peng, F.J. Vermolen, D. Weihs
Abstract
The phenomenological model for cell shape deformation and cell migration [42, 6] isextended with the incorporation of cell traction forces and the evolution of cell equilibriumshapes as a result of cell differentiation. Plastic deformations of the extracellular matrixare modelled using morphoelasticity theory. The resulting partial differential differentialequations are solved by the use of the finite element method. The paper treats variousbiological scenarios that entail cell migration and cell shape evolution. The experimentalobservations in [24], where transmigration of cancer cells through narrow apertures isstudied, are reproduced using a Monte Carlo framework.
Cells may attain various shapes and sizes, for example, stem cells can differentiate and adoptthe shape and functionality of many different cell types in our body: fan-like keratocytes, hand-shaped nerve growth cones and spindle-shaped fibroblasts [27, 32]. It has been recognized thatcell geometry influences cellular activities like cell growth and death, cell mobility and adhesionto the direct environment [2, 18, 25, 27, 36]. The shape of a mobile cell is determined by itsboundaries, which dynamically vary with a local balance between retraction and protrusion[10]. There are multiple constituent elements affecting the cell shape, for instance, the cy-toskeleton and the cell-substrate adhesions, which have been studied in depth in the past years.However, it is still a great challenge to understand the mechanisms that determine the globalcell morphology in the context of its function [18, 27].Signalling molecules play an important role in cell migration and cell shape. During woundhealing, chemotaxis is one of the most important cues for migration of immune cells and fi-broblasts in inflammatory and proliferative phases [19, 8, 11, 28]. Metastasis of cancer cellscan be induced by nutrients and oxygen, since tumour growth requires an adequate supplyof oxygen and nutrients. Under most pathological circumstances, oxygen and nutrients aresupplied though the local blood vasculature [44, 38]. Commonly, signalling molecules are ac-tivated at the plasma membrane, and de-activated in the cytoplasm. On the other hand, theconcentration of signalling molecules determines the cytoskeletal dynamics [27].In wound healing, cells migrate and change shape in both the epidermis and the dermislayers. Re-epithelialization is the most essential part for the skin to re-establish its barrierfunction [37, 39, 12]. However, the mechanisms of re-epithelialization are poorly understood.In the early stage of the epidermis closure in a wound, the basement membrane between theepidermis and dermis extends slightly over the ends of the incised dermis, creating an “extensionmembrane” (or so-called epidermal tongue) [31]. The mechanism of the occurrence of theepidermal tongue is still unclear. A possible explanation is that, the suprabasal cells (whichlie upon the layer of basal cells) form the tongue by migrating over the leading basal cells1 a r X i v : . [ q - b i o . CB ] D ec nd dedifferentiating to basal cells (which are adhered to the basement membrane between theepidermis and dermis) to form new leaders [37, 31, 43, 33]. When epidermal epithelial cellsare ”crawling” and ”climbing up” to reestablish the epidermis, they elongate and flatten [37].In the dermis, it has been widely documented that the differentiation of fibroblasts is one ofthe key events during wound healing. Differentiation changes the spindle-shaped fibroblast todendritic-shaped myofibroblasts. Subsequently, cells’ mechanobiology is modified considerablyas well. The differentiated myofibroblasts exert much larger forces on the extracellular matrix(ECM) than fibroblasts [28]. Excessive numbers of myofibroblasts will result in contractures,which are morbid and pathological macro-scale contractions. Usually, contractures concur withdisabilities and dysfunction, and have a grave impact on patients’ daily life.Cancer metastasis has been reported as the main reason of death in cancer patients [25].During the migration of a cancer cell to its destination, especially migrating through a narrowand stiff cavity, it has to deform to adapt to the obstacles. More invasive cancer cells, appearbe more pliable and dynamic both internally [13] and externally [15, 7, 40], and thus able toadjust their cytoskeleton and morphology, which might provide a possible diagnosis for cancer.In addition to that, cancer cells are observed to apply a significantly larger traction force on thesubstrate, compared to benign cells [25], yet the specific mechanisms the induce these increasedforces are still poorly understood.Mathematical modelling has been proven to be an important tool to have a deeper insightinto many biological processes that are potentially difficult to control in experiments, for ex-ample, wound healing and tumour growth. Depending on the scale of the observed domain,continuum models and agent-based models are widely used. Continuum models have the ad-vantage of modelling a larger scale, however, the model neglects the individual cellular activityand cells are not tracked [42]. Agent-based model is suitable to model cellular activities ofevery cell, for instance, cell migration and cell deformation. Hence, an agent-based model isselected in this manuscript and this work is an extension of [6] and [42]. In [6], a model of thedeformation of both the cell and the nucleus is developed. Furthermore, a parameter sensitivityanalysis is carried out on the basis of Monte Carlo simulations. However, the study in [6] doesnot consider the traction forces applied by the cell and the impact on the substrates. Com-pared to the work of Vermolen and Gefen [42], we use finite-element methods to solve all thepartial differential equations, rather than Green’s functions. Therefore, a more precise solutionis delivered. Furthermore, we implement a more intricate approach to model the traction forcesapplied by cells in various applications. In addition to circular projections of cells in [6] and[42], we model elliptic and hypocycloid-shaped cells in this manuscript.This manuscript is structured as follows: Section 2 explains the agent-based model of cellmigration, in the form of a set of partial differential equations. Possible applications of thismodel and the corresponding numerical results are exhibited in Section 3. Finally, conclusionsare shown in Section 4. In this manuscript, the phenomenological model of cell deformation is extended from the workin [6, 42], in particular, in two dimensions. With essential biological assumptions and simplifi-cations, the model mainly describes the impact of extracellular components on the cell defor-mation and displacement. Subsequently, more applications can be developed, for instance, celldifferentiation and cell repulsion. Different from the work in [45], the model in this manuscriptneglects the intracellular environment, hence, it is not capable to present the Poisson’s effectof the cell.The cell membrane is split into finite line segments by the nodal points, and the centreposition of cell is determined by the mean of all the positions of the nodal points. The equilib-2 igure 2.1: A schematic of the distribution of the nodal points on the cell membrane. In ourmodel, we assume the equilibrium shape of the cell is maintained by a collection of springs. rium shape of the cell is kept by a collection of springs, which connects each nodal point on thecell membrane to the centre of the cell, respectively; see Figure 2.1. For each nodal point, thedisplacement is determined by various mechanisms of directed motion and random turning [10],which will be discussed in details in the following contents. Regarding different applications ofthis model, there will be some model adjustments.
We assume that cell migration is mainly driven by chemotaxis (or mechanotaxis) , which iscommonly observed in wound healing and cancer cell metastasis with various types of sig-nalling molecules. In wound healing, immune cells are directed to chase a bacterium or virus;high concentration of Transforming Growth Factor-beta (TGF-beta) induces the migrationof (myo)fibroblasts towards the wound from the uninjured skin [8, 11]. Cancer metastasis istriggered by cancer cell proliferation and cell migration. On a cellular level, chemotaxis orhaptotaxis is an essentional cue for cancer cell migration and hence for the dissemination oftumours [44]. Numerous studies indicate that the availability of oxygen and nutrients is one ofmost crucial factors for the growth of tumours [34]. During tumour growth, the concentrationof oxygen and nutrients depletes in the vicinity of the tumour [38]. Therefore, cancer cells havea tendency to migrate towards regions with higher concentrations of oxygen and nutrients.A point source, which is modellled by a Dirac Delta distribution, is represented as thesource of a generic signal. Together with the reaction-diffusion equation, the concentration ofthe signal is determined by: ∂c ( x , t ) ∂t + ∇ · ( v c ( x , t )) − ∇ · ( D ∇ c ( x , t )) = kδ ( x ( t ) − x s ) , (2.1)where c ( x , t ) is the concentration of the signalling molecule, D is the diffusion rate, k is thesecretion rate of the signal source, x s is the position of the source, and v is the displacementvelocity of the substrate that results from the cellular forces exerted on their surroundings. Thevelocity is computed by solving the balance of the momentum, which will be discussed in theSection 2.2.Initially, we assume there is no signalling molecules over the computational domain, that is, c ( x ,
0) = 0 . As a boundary condition, we use the following Robin condition ∂c∂ n + κ s c = 0 , which deals with a balance between the diffusive flux across the boundary and the flux betweenthe boundary and the region far away from the domain of computation. The symbol κ s ,3hich is non-negative, represents the mass transfer coefficient. Note that as κ s → κ s → ∞ represents the case that c → In wound healing, (myo)fibroblasts exert forces on their direct environment, i.e. extracellularmatrix, which result into contraction of the tissue [8, 11, 16, 20]. For cancer cells, Massalhaand Weihs [25] indicate that the metastatic cells exert traction forces ranged from 100 − nN on the gel, of which the Young’s modulus ranged from 2 . − . kP a . Furthermore, for stiffersubstrates, the cancer cells remain rounded with changing area and they exert large tractionforces with large magnitudes to its direct environment. Hence, the model includes passiveconvection of the substrate, which can provide a more realistic model in various applications.As the cell membrane is broken into multiple line segments by nodal points, point forces areimplemented here to depict the forces exerted by the cell, which are applied on the midpointof each line segment; see Figure 2.2 as an example of a square-shape cell. Among differentapplications, the force direction may differ. For example, if the cell encounters an obstacle,a repulsive force will be exerted to resist the compression of the cell; in wound contraction,(myo)fibroblasts exert pulling forces on the extracellular matrix (ECM). Figure 2.2: An example of pulling point forces, tractions applied by a cell, in which the cellmembrane is split by four nodal points.
Morphoelasticity is widely used in the biological modelling to describe elastic growth, forinstance, the growth of tumours [14], the seashell growth [35] and large contractions in woundhealing [19, 3] etc. In wound healing, morphoelasticity describes the phenomena when the de-formation of the skin is so large that the deformations are plastic. Conservation of momentum,combined with the evolution equation for the effective Eulerian strain, results into the following4odelling equations [19]: ρ [ D v Dt + v ( ∇ · v )] − ∇ · σ = f , in Ω, D (cid:15) Dt + (cid:15) skw( L ) − skw( L ) (cid:15) + [tr( (cid:15) ) −
1] sym( L ) = − α (cid:15) , in Ω, v ( x , t ) = , on ∂ Ω, (2.2)where ρ is the density of the extracellular matrix, L = ∇ v and α is a non-negative constant.Note that if α = 0, then as soon as the force f = , then the tissue will gradually recover toits original shape and volume. Here, D y Dt = ∂ y ∂t + v ∇ · y is material derivative where y is anytensor field and v is the migration velocity of any point within the domain of computation. Inorder to have a fixed boundary, we use a homogeneous Dirichlet boundary condition for thevelocity. From a mechanical point of view, we treat the computational domain as a continuouslinear isotropic domain. Further, as a result of the presence of liquid phases in the tissue,the mechanical balance is also subject to viscous, that is friction, effects. Therefore, we useKelvin–Voigt’s viscoelastic dashpot model, of which the stress tensor reads as σ = σ elas + σ visco = E ν s { (cid:15) + tr( (cid:15) )[ ν s − ν s ] I } + µ sym( L ) + µ tr(sym( L )) I , (2.3)where ν s is the Possion’s ratio of the substrate, (cid:15) is the strain tensor, µ and µ are the shearand bulk viscosity respectively. The morphoelasticity model solves non-linear equation andboth velocity v and strain tensor (cid:15) are unknowns. The deformation of the domain is actuallydetermined by the strain tensor. The displacement of the domain can be approximated byintegrating the velocity over time: u ( t ) ≈ (cid:82) t v ( s ) ds .In the application of wound healing, (myo)fibroblasts are the cells pulling the ECM andcausing the contractions. The temparory force f of each (myo)fibroblast reads as f = N (cid:88) j =1 P ( x , t ) n ( x ( t )) δ ( x ( t ) − x j ( t ))∆Γ j , (2.4)where N is the number of nodal points on the cell membrane, P ( x , t ) is the magnitude of theforce, n ( x ( t )) is the unit inward pointing normal vector (towards the cell centre) at x ( t ), x j ( t )is the midpoint of line segment j , and ∆Γ j is the length of line segment j .Further, we consider (myo)fibroblasts colliding with each other, then the cells exert repellingforces on the other one. Here, cells are not allowed to intersect each other; see Figure 2.3 as ascheme. Suppose for (myo)fibroblasts i , there are N m = { j , . . . , j m } line segments of cell mem-brane mechanically contacting with other (myo)fibroblasts. Then, on line segments j ∈ N m ,the (myo)fibroblast exert repelling force, while on the rest of the line segments, (myo)fibroblastreleases pulling forces on the ECM. Hence, the temporary force of the (myo)fibroblast i is givenby f i = N (cid:88) j =1 ,j / ∈ N m P ( x , t ) n ( x ( t )) δ ( x ( t ) − x ij ( t ))∆Γ i,j − (cid:90) l m Q ( d ( x ) , t ) n ( x ( t )) δ ( x ( t ) − x ij ( t )) d Γ i , (2.5)where l m is the portion of the (myo)fibroblast membrane mechanically contacting with othercell, Q ( d ( x ) , t ) is the force magnitude per length and d ( x ) is the penetration depth. According5o two-dimensional Hertz theory [30, 22], for each elastic body, the total force magnitude˜ Q ( d ( x ) , t ) is linearly proportional to the penetration depth d ( x ):˜ Q ( d ( x ) , t ) = π d ( x ) E ∗ , (2.6)where E ∗ is the total equivalent Young’s modulus derived by1 E ∗ = 1 − ν E + 1 − ν E . Here, ν i and E i with i ∈ { , } represent the Poisson ratio and Young’s modulus of two elasticbodies, respectively. In particular, if two bodies have the same elastic characteristics (i.e. ν = ν = ν and E = E = E ), then E ∗ = E − ν ) . We assume that the magnitudes of the repulsive force, which is exerted on the boundarysegments of cell i that are in contact with another cell, are identical. In other words, Q ( d ( x ) , t )is given by Q ( d ( x ) , t ) = ˜ Q ( d ( x ) , t ) (cid:107) l m (cid:107) = π d ( x ) E ∗ / (cid:107) l m (cid:107) , (2.7)where (cid:107) l m (cid:107) is the total length of the portion of the membrane of (myo)fibroblast i mechanicallycontacting with other (myo)fibroblast. Subsequently, the total temporary force is f = (cid:80) N C i =1 f i ,where N C is the number of (myo)fibroblasts. Figure 2.3: When cells collide with each other, they will deform and exert repulsive forces.The dashed curves show the equilibrium shape of cells, and the black curve is the overlappingmembrane of both cells.
Another application is metastasis and invasion of cancer cell. Usually, in vitro, a microtubeexperiment is conducted [24]. To simplify the model, here we only consider one cancer cellgoing through a microtube; see Figure 2.4. Similar to the case when (myo)fibroblasts collide,we assume that l m is the portion of cell membrane mechanically contacting with the wall of themicrotube. Here, we will exclude the pulling force. In other words, the force released by thecancer cell is only the repelling force exerted on the wall of the microtube: f = − (cid:90) l m Q ( d ( x ) , t ) n ( x ( t )) δ ( x ( t ) − x j ( t )) d Γ . (2.8)The magnitude of the force here follows the same definition as in Eq (2.7), where d ( x ) is theradius subtracting the distance from the cell membrane to the cell centre.6 igure 2.4: The cell is compressed by the wall of the microtube. To inhibit any changes from itsequilibrium shape, it exerted the repelling force on the wall of the microtube, which is proportionto the compressed distance h . According to [6, 42], the cell cytoskeleton is depicted as a collection of springs between the centreposition of the cell and the nodal points on the cell membrane. Therefore, the equilibrium shapeof the cell is kept by these springs, regardless of the original cell shape. In this manuscript,we consider circular, elliptic and star-shape cells as two-dimensional projections. Combiningchemotaxis (or mechanotaxis), passive convection and random walk, the displacement of thenodal point j on the cell membrane is given by d x j = β ( CSI (Ω C ) , A (Ω C )) ∇ c ( x , t ) (cid:107)∇ c ( x , t ) (cid:107) + γ dt + E c ( x c ( t )+ ˆ x j − x j ( t )) dt + v dt + σ rw d W ( t ) , (2.9)Here, E c represents the cell elasticity; ˆ x = x c ( t ) − ˜ x j ( t ) is the vector connecting the equilibriumposition of nodal point i on the cell membrane to the cell centre, x c is the central position ofthe cell and ˜ x j represents the equilibrium position of the nodal point j corresponding to thecell centre x c (see Figure 2.1); γ is a small positive constant to prevent the denominator beingzero; v is the velocity of the substrate determined by Eq (2.2); σ rw is the portion of randomwalk, and d W ( t ) is a vector-Wiener process, which accounts for random walk. Furthermore, β = β ( CSI (Ω C ) , A (Ω C )) is the weight of chemotaxis (or mechanotaxis), where we define theCell Shape Index (CSI) of cell Ω C by CSI (Ω C ) = 4 πA (Ω C ) l ( ∂ Ω C ) , where A (Ω C ) is the projected area of the cell, l ( ∂ Ω C ) is the circumference of the cell membrane.According to [18], reduction of cell volume and deformation of cell shape reduce the mobilityof cell. For simplicity, we propose a linear relation here: β ( CSI (Ω C ) , A (Ω C ) = β × µ m ( CSI (Ω C ) /CSI (Ω C ) + A (Ω C ) /A (Ω C )) / , where β is the maximal response from the cell to the signal, µ m is the mobility reductioncoefficient, and CSI (Ω C ) and A (Ω C ) represent the CSI and volume of the equilibrium cell.In order to maintain the right orientation of the cell, we introduce a matrix after rotationof an angle φ , as in [6]: B ( φ ) = (cid:18) cos( φ ) − sin( φ )sin( φ ) cos( φ ) (cid:19) , (2.10)such that φ can be computed from˜ φ = arg min φ ∈ [0 , π ) (cid:32) N (cid:88) i =1 (cid:107) B ( φ ) ˜ x i ( t ) − x i ( t ) (cid:107) (cid:33) . (2.11)7ence, the displacement of nodal point j is adapted to d x j = β ( CSI (Ω C ) , A (Ω C )) ∇ c ( x , t ) (cid:107)∇ c ( x , t ) (cid:107) + γ dt + E c ( x c ( t ) + B ( ˜ φ ) ˆ x j − x j ( t )) dt + v dt + σ rw d W ( t ) . (2.12)If there is an obstacle encountered by the cell, adjusting the displacement is necessary.Denote ∂ Ω ob as the boundary of the obstacle, which is possibly another cell or the wall of themicrotube. For the nodal point colliding the obstacle, it cannot pass over the boundary of theobstacle. Hence, for the displacement of nodal point j , the normal direction of the boundaryof the obstacle in d x j ( t ) is vanished. To rephrase it, we adjust the displacement of the nodalpoint if it collides the obstacle by d x j ( t ) ← d x j ( t ) − ( d x j ( t ) , n ob ( x j ( t ))) n ob ( x j ( t )) , if x j ( t ) ∈ ∂ Ω ob , (2.13)where n ob ( x ) is unit pointing normal vector (outward the cell centre).Furthermore, if a cell tries to go through a microtube, which mimic the cancer cell metastasisand invasion, not only the direction of the cell migration is limited (since the cell cannot passover the microtube), but also the microtube will slow down the velocity of the cell as a resultof friction. Note that the magnitude of the friction is proportional to the repelling force. Inour model, we simply distract part of the velocity in the tangential direction of the obstacle.Hence, the displacement of the nodal point which collides the wall of the microtube is given by d x j ( t ) ← d x j ( t ) − µ f (cid:107) f ( x j ( t )) (cid:107) × ( d x j ( t ) , τ ob ( x j ( t ))) τ ob ( x j ( t )) , if x j ( t ) ∈ ∂ Ω ob , (2.14)where µ f is the cell friction coefficient, f ( x j ( t )) is the repelling force exerted by the cell and τ ob ( x ) is the tangential direction of the obstacle boundary ∂ Ω ob .This model provides a simple computational framework to describe the dynamics of thecell shape under multiple circumstances. However, the model does not describe the Poissoneffect of the cell if the cell is compressed, since the model mainly considers the extracellularenvironment. Hence, in this manuscript, the cell length will not be investigated. We exhibit several possible applications in this section, namely, cells migrating as a result ofchemotactic signals, cells differentiating to another phenotype, cells repelling each other andone cell migrating through a microtube. Some parameters are the same in every application.If there is no specification, the parameter values are shown in Table 3.1.
The basic application is that cell migrates towards the concentration gradient of a signallingmolecule, which can be oxygen, growth factors or virus. The displacement is mainly determinedby the gradient of the concentration of the signal. Subsequently, the closest part of the cellto the emitting source will develop a ”nose”; see Figure 3.1. All parameter values have beendocumented in Tables 3.1 and 3.2. The ’nose’ behaviour (or so-called a triangular tailed shape)has also been observed in biological experiments, which is due to the locations of the adhesionsites over the cell membrane [27]. In the simulation we accommodate for the engulfment of thechemical source by switching off the chemical signla once the cell physically contacts it. Oncethis signal has been switched off, the concentration gradient flattens as a result of diffusionprocesses and therewith the cell recesses back to its equilibrium (original) shape and volume.At that moment, the cell geometry is no longer determined by chemotaxis.8 able 3.1: Parameter values used in all the applications
Parameter Description Value Units Source E s Substrate elasticity 100 kg/ ( µm · min ) [21] E c Cell elasticity 5 kg/ ( µm · min ) [5] µ f Cell friction coefficient 0 . − [1] ν s Poisson’s ratio of the ECM 0 . − [19] ν c Poisson’s ratio of(myo)fibroblast and cancer cell 0 . − [41] k Secrete rate of the signal 2 . kg/ ( µm · min ) [28] κ s Parameter in Robin’s boundarycondition to solve Eq (2.1) 100 − [28] µ Shear viscosity of the ECM 33 . − [28] µ Bulk viscosity of the ECM 22 . − [28] β Maximal mobility of points oncell membrane 10 min − Estimated in thisstudy N Number of nodal points on thecell membrane 40 − Estimated in thisstudy µ m The coefficient of cell mobilityreduction 1 − Estimated in thisstudy σ rw Weight of random walk 1 − Estimated in thisstudy α Degree of permanentdeformation in Eq (2.2) 0 . − Estimated in thisstudy
From t = 0, the cells are attracted towards the centre of the computational domain,whichis the location of the source of signalling molecules. Due to the difference of the gradientof the concentration of the signal over the domain, cells are deformed into a droplet shape,where the ”nose” points in the direction of the point source. As the diffusion of the signalproceeds, the ”nose” gradually disappears and cells recover to the equilibrium shape. Toevaluate the cell geometry quantitatively, we provide the evolution of the CSI and cell volume asa function of time in Figure 3.2. These quantities are of interest from a clinical point of view [18].Resulting from the displacement of the direct environment, the volume of the cell decreases.The permanent volume changes of the cell are imposed by the permanent displacements fromthe morphoelastic model. Furthermore, the cells are tethered within a rigid deformed structure,hence, it makes cells deform as well. We note that cells have already recovered to the originalshape but not the volume. It can be implied that a stiffer cell deforms less compared to a softerone. Cell differentiation is a process of a cell changing from one phenotype to another one, forexample, a stem cell differentiates into various phenotypes, like blood cells and nerve cellsetc. In this manuscript, we mainly focus on the cellular differentiation in wound healing. Inthe proliferative phase of wound healing,, some regular fibroblasts (which are spindle-shaped[4, 29]) differentiate into myofibroblasts (which are dendritic-shaped [17, 9]), which pull theskin ever harder and cause the permanent contractions. It is commonly recognized that highconcentration of TGF-beta induces the fibroblast-to-myofibroblast differentiation [18, 8, 9]. Inthis section, since we mainly want to present a model of differentiation, only differentiation fromfibroblasts to myofibroblasts is exhibited as an example: the signal is TGF-beta and initially,there are four regular fibroblasts in unwounded region, which are simulated by ellipses.We assume that the two phenotypes of cells have the same volume when they are in the9 a) t=0 (b) t=0.2min (c) t=0.5min(d) t=1.5min (e) t=2.5min (f) t=5min
Figure 3.1: The figure shows the positions of cells at consecutive times. The red circles showthe original shape and position of cells, and the red dots in the centre is the point source of ageneric signal.Figure 3.2: The cell shape index and relative ratio of cell volume of all cells in Figure 3.1 areshown in the plot. The solid curves represent the cell volume, and the dashed curves are the cellshape index. Different colours of curves indicate different stiffness of the cells. able 3.2: Parameter values estimated in the application of cell migrating toward the signalsource Parameter Description Value Units R Cell radius 5 µm ∆ t Time step 0 . minD Diffusion rate of the signal 200 µm /minx Length of computational domain inx-coordinate 120 µmy Length of computational domain iny-coordinate 120 µmx w Length of the subdomain in the centreof computational domain inx-coordinate 40 µmy w Length of the subdomain in the centreof computational domain iny-coordinate 40 µm equilibrium status. Here, for ellipse and hypocycloid, there are two parameters to determineeach shape: long (denoted by a e ) and short axis’s (denoted by b e ) determine the ellipse; theradius of basis circle (de noted by a h ) and rotating circle (denoted by b h ) determine the hypocy-cloid. Note that the hypocycloid-shaped cell may not be realistic, and it is mainly to show thatthe model is capable to model the differentiation of cells, in which mostly the cellular skeletonand geometry are altered. To develop a smooth differentiation process, we introduce a functionsuch that each parameter changes over time: (cid:40) R a ( ω ) = a h ω + a e (1 − ω ) ,R b ( ω ) = b h ω + b e (1 − ω ) , (3.1)where R a ( ω ) and R b (Ω) represent two parameters to determine the shape, and ω = 1 − exp {− λ ω ( t − t ω ) } . Here, λ ω is the parameter of the exponential distribution and t ω is thetime point when the fibroblast starts differentiating.Figure 3.3 presents the cells positions at different time. As regular fibroblasts approachto the region with high concentration of TGF-beta, some of them start differentiating intomyofibroblasts gradually. Subsequently, they exert larger forces on the ECM and contractionsare developed in the wound, which is marked with red curves as a subdomain. The parametervalues of the simulations are from Table 3.1 and Table 3.3. Cells will deform when they encounter each other or an obstacle. On the contacting surface,cells will exert a repelling force (as it is shown in Eq 2.7) to recover to its equilibrium shape.We consider two cells colliding with each other and adjust the displacement of nodal pointson the cell membrane by Eq (2.13). Here, cells are not allowed to intersect with each other.Hence, initially, cells are located with a small distance between each other. In Figure 3.4, wepresent the cell positions at different times, and cells deform due to mechanical contact (hardimpingement). The parameter values are from Table 3.1 and Table 3.2.
Metastasis is a difficult phenomena to study due to its large variation in spatiotemporal scales.Hence, studying the mechanics of one single cell is essential since the individual cell needs11 a) t=0 (b) t=0.5min (c) t=4min(d) t=4.5min (e) t=5.5min (f) t=7min
Figure 3.3: The figure shows the positions of cells at consecutive times by the black curves. Thered circles show the original shape and position of cells, and the red dots in the centre is the pointsource of TGF-beta which triggers the differentiation from regular fibroblasts to myofibroblasts.Table 3.3: Parameter values estimated in the application of cell differentiation
Parameter Description Value Units a e Length of long axis in elliptic cell 6 . µmb e Length of short axis in elliptic cell 4 µma h Radius of the basis circle to drawhypocycloid-shape cell 20 / √ µmb h Radius of the rotating circle to drawhypocycloid-shape cell 5 / √ µmλ ω Parameter in the exponentialdistribution to compute ω − ∆ t Time step 0 . minD Diffusion rate of the signal 233 . µm /minx Length of computational domain inx-coordinate 120 µmy Length of computational domain iny-coordinate 120 µmx w Length of the subdomain in the centreof computational domain inx-coordinate 40 µmy w Length of the subdomain in the centreof computational domain iny-coordinate 40 µm a) t=0 (b) t=0.2min (c) t=1.5min(d) t=3min (e) t=5min (f) t=6min Figure 3.4: The figure shows the positions of cells in blue and green at consecutive times whentwo cells collide. The red circles show the original shape and position of cells, and the red dotsin the centre is the point source of a signal.
13o break out from the tumour and invade though the ECM. To achieve that, Mak et al. [24]developed an active microfluidic system with several features to mimic the metastasis andinvasion of the cancer cell. In this section, we will use our model to reproduce results of theexperiment in [24].
Following the settings in the experiment, we define a microtube with a varying width: a 15 µm larger channel (LC) and 3 . µm subnucleus barriers with length 10 µm (SNB10). Since the mainreason for the active migration of the cell is not evident in [24], we keep on using randomwalk, with either chemotaxis or fixed velocity of the cell. Rather than having a periodic settingof SNBs, we have one SNB in the middle connected with two LCs and run the simulationsrespectively. In this manuscript, we only run the simulations regarding SNB10 in [24].We consider the reduction of the cell mobility caused by cell shape and cell volume [18],which is explained in Eq (2.12), the repelling forces exerted by the cell on the obstacles in Eq(2.13) and the friction between the cell and the wall of the microtube in Eq (2.14).The position and shape of the cell are shown in Figure 3.5, which indicates how the cellmigrates through the microtube. Since the repelling force on the wall of the microtube isincluded, we investigate the results regarding the cell velocity, pressure and the cell shapeindex over time; see Figure 3.6. The parameter values are taken from Table 3.1 and 3.4. Table 3.4: Parameter values used in the application of cell going through a microtube
Parameter Description Value Units Source R Cell radius 9 µm [24]∆ t Time step 0 . min Estimated in thisstudy D Diffusion rate of the signal 874 . µm /min Estimated in thisstudy x Length of computational domaincontaining SNB10 in x-coordinate 400 µm Estimated in thisstudy y Length of computational domaincontaining SNB10 in y-coordinate 400 µm Estimated in thisstudy
Initially, there is a short distance before the cell enters the microtube, here the cell en-counters no distraction. Therefore, the the cell travels at maximal speed and the cell is notcompressed in the beginning. Next to this, the gradient of the signal results into the ’nose’behaviour and hence, the cell shape index changed. As the cell enters the wider part of themicrotube, it slows down due to the friction, and the cell is compressed, therefore, the cellstarts exerting pushing forces on the wall of the tube. In the LC part, the cell shape indexstays stable around 0 .
95. Further, the cell approaches the SNB, which is much more narrowthan the LC, the cell suffered more from the friction and the compression from the microtube.As a consequence, the minimal cell velocity, the cell shape index and the maximal pressure isrecorded when the cell is in the SNB. After that, cell moves further towards the signal sourcethrough the LC again. Hence, the cell velocity and cell shape index increase again, while thecell pressure reduces. According to [1] and [26], we manage to keep the cell velocity and pres-sure in a reasonable range: 6 − µm/min and the maximal pressure that a cell can handle isaround 12 kP a . 14 a) t=0 (b) t=0.35min(c) t=3.5min (d) t=21min(e) t=24.5min (f) t=28min(g) t=35min (h) t=49min Figure 3.5: The figure shows the positions and shapes of cells at consecutive times by the bluecontours when it travels through a microtube. The red circles show the original shape and positionof cells, and the red dots in the end of the microtube is the point source of a signal. a) Cell Velocity (b) Cell Pressure (c) Cell Shape Index Figure 3.6: The cell velocity, pressure and shape index over time when the cell migrates throughthe microtube, where there is a µm subnucleus barrier. The simulation mimics the experimentin [24]. In [24], the displacement of the cell is categorised as four phases: • Phase 1: the cell enters the microtube via the LC, and slows down in particular when itis approaching the SNB; • Phase 2: the cell is compressed strongly to enter the SNB; • Phase 3: the cell fails to migrate monotonically forward when it is in the SNB; • Phase 4: the cell enters the LC again and continues to migrate monotonically.Hence, in the simulations, we try to collect the data and distinguish these different phases.Different from [24] that the microtube is designed periodically (such that the sample can becollected multiple times with one individual cell), one cell is supposed to go through one setof the microtube in each simulation. To reproduce the experimental results, Monte Carlosimulations are conducted to estimate the probability of the occurrence of phase 3 and thetime cost for each phase, with different aforementioned reasons of active migration. The inputvalues for the Monte Carlo simulations are shown in Table 3.5.
Table 3.5: Parameter values used in the application of cell going through a microtube
Parameter Description Distribution Source µ f Friction coefficient for the cell goingthrough the microtube U (0 . , .
06) [1] µ m The coefficient of cell mobilityreduction U (0 . ,
1) Estimated in thisstudy
We run the simulations with four assumptions of the main mechanism provoking the ac-tive cell displacement: chemotaxis, fixed velocity with 10 µm/min , velocity generated from(6 , µm/min and (6 , µm/min in horizontal direction according to [26]. The number ofsamples and the Monte Carlo error of the occurrence of phase 3 collected from the Monte Carlosimulations of each aforementioned category are displayed in Table 3.6. Figure 3.7 illustratesthe probability of the occurrence of phase 3, which is the stage when the cell stops monotonicforward migration. The mechanism that makes the cell move forward is not clear in [24], hencethis could be reason for a mismatch between the experimental and simulated results. However,the results are significantly close regarding the probability of the occurrence of phase 3. Coin-cidently, according to our current simulations, with either chemotaxis or velocity 10 µm/min ,the probability of the occurrence of phase 3 is the same in the first 3 digits.16 able 3.6: Monte Carlo simulations of various models, in which the main mechanisms of cellactive displacement differ Mechanism of CellActive Displacement Number of Samplesfrom Monte CarloSimulations Monte Carlo Errorof the Occurrenceof Phase 3Simulation 1
Chemotaxis 1400 9 . × − Simulation 2
Fixed velocity v = 10 1390 9 . × − Simulation 3
Fixed velocity v ∈ (6 ,
15) 1360 1 . × − Simulation 4
Fixed velocity v ∈ (6 ,
20) 1378 1 . × − Figure 3.7: The probability of the occurrence phase 3 in [24] and from Monte Carlo simula-tions by implementing different mechanisms of cell active displacement (see Table 3.6 for moreinformation). The parameter values are taken from Table 3.5. a) Chemotaxis (b) Fixed velocity of 10 µm/min (c) Random velocity generated from (6 , µm/min (d) Random velocity generated from (6 , µm/min Figure 3.8: The figure shows the time cost of each phase in [24] and from Monte Carlo simula-tions by using the model. Red dots with the error bar represent the experimental data from [24]and the box plots are the data collected from the simulations. v = 10. The reason is that the microtube restricts the displacement of the cell inthe vertical direction, therefore, the cell mainly migrates in the horizontal direction. In general,the results between any simulation and the experiment differ more, compared with the resultsof the occurrence of phase 3, in particular, phase 1 and phase 3. Therefore, to investigate thepossible reasons of mismatching results in phase 1 and phase 3 in phase time, we reran thesimulation with the same settings as Simulation 4 in Table 3.6, except for the cell stiffnessmodified to E c = 1. The results are shown in Fig 3.9. With a softer cell, the simulation datain phase 3 match better with the experimental data. However, now a discrepancy betweensimulation and experiments results in phase 4 instead. Figure 3.9: The time cost of each phase from the Monte Carlo simulation, in which the cellstiffness is E c = 1 and the cell velocity is randomly generated from (6 , . Red dots with errorbar represent the experimental data from [24] and the box plots are the data collected from thesimulation. There are several possible reasons causing the discrepancy. Firstly, for phase 3, we onlyobtain valid data when the cell moves non-monotonically, which results into the reduction ofthe sample size. Secondly, the length of LC is not stated clearly in [24], therefore, we couldonly estimate that from the scale in the figures. Thirdly, it is not clear if the velocity of activemigration of the cell is constant, while in our simulation, the velocity can change over time,depending on the gradient of the chemotactic signal. Despite all these uncertainties, we stillmanaged to reproduce the results which are close to the experimental results. Fourthly, thetransaction of each phase from [24] to our simulations may cause a mismatch of the durationtime of each phase. Fifthly, it has been observed in [23] that after the first time moving throughthe narrow channel, cells deform easier to move faster through the following narrow channels,which may indicate that the cell characteristics change regarding its geometry.
A phenomenological model for cell shape evolution and migration has been developed, withprimary focus on the mechanics of the extracellular environment. Furthermore the impact ofpassive convection, due to local displacements within the extracellular matrix, on the evolution19f the cell shape has been taken into account. A morphoelastic model has been used in thecurrent study to incorporate permanent deformations of the extracellular tissue. The modelcan be applied to mimic several microscopic biological observations such as cell deformationand migration during wound contraction and cancer metastasis. To validate the model, theexperimental setup in [24] has been modeled. This experiment entailed cell migration throughmicrotubes with different widths and with a varying width over the length. The model is ableto reproduce the most important trends that were observed in the experimental data despitesome experimental uncertainties such as the determination of which phase a cell is in during thetransmigration process. Furthermore, the current model provide a basis that can be expandedto describe more experimentally observed phenomena in cell geometry.
References [1] T. E. Angelini, A. C. Dunn, J. M. Urue˜na, D. J. Dickrell, D. L. Burris, and W. G.Sawyer. Cell friction. Faraday Discussions, 156:31, 2012. doi: 10.1039/c2fd00130f. URL https://doi.org/10.1039/c2fd00130f .[2] E. L. Barnhart, K.-C. Lee, K. Keren, A. Mogilner, and J. A. Theriot. An adhesion-dependent switch between mechanisms that determine motile cell shape. PLoS Biology, 9(5):e1001059, may 2011. doi: 10.1371/journal.pbio.1001059. URL https://doi.org/10.1371%2Fjournal.pbio.1001059 .[3] M. Ben Amar, M. Wu, M. Trejo, and M. Atlan. Morpho-elasticity of inflammatory fibrosis:the case of capsular contracture. Journal of the Royal Society Interface, 12(111):20150343,2015.[4] R. Chaudhari. Myofibroblasts: Functions, evolution, origins, and the role in disease. SRMJournal of Research in Dental Sciences, 6(4):234, 2015. doi: 10.4103/0976-433x.156219.URL https://doi.org/10.4103/0976-433x.156219 .[5] J. Chen, D. Weihs, and F. J. Vermolen. A model for cell migration in non-isotropic fibrinnetworks with an application to pancreatic tumor islets. Biomechanics and Modeling inMechanobiology, pages 1–20, 2017.[6] J. Chen, D. Weihs, M. V. Dijk, and F. J. Vermolen. A phenomenological model forcell and nucleus deformation during cancer metastasis. Biomechanics and Modeling inMechanobiology, 17(5):1429–1450, May 2018. doi: 10.1007/s10237-018-1036-5. URL https://doi.org/10.1007/s10237-018-1036-5 .[7] S. E. Cross, Y.-S. Jin, J. Rao, and J. K. Gimzewski. Nanomechanical analysis of cells fromcancer patients. Nature Nanotechnology, 2(12):780–783, Dec. 2007. doi: 10.1038/nnano.2007.388. URL https://doi.org/10.1038/nnano.2007.388 .[8] B. D. Cumming, D. McElwain, and Z. Upton. A mathematical model of wound healingand subsequent scarring. Journal of The Royal Society Interface, 7(42):19–34, 2009.[9] V. D. Desai, H. C. Hsia, and J. E. Schwarzbauer. Reversible modulation of myofibroblastdifferentiation in adipose-derived mesenchymal stem cells. PLoS ONE, 9(1):e86865, Jan.2014. doi: 10.1371/journal.pone.0086865. URL https://doi.org/10.1371/journal.pone.0086865 . 2010] H. Ebata, A. Yamamoto, Y. Tsuji, S. Sasaki, K. Moriyama, T. Kuboki, and S. Kidoaki.Persistent random deformation model of cells crawling on a gel surface. Scientific Reports,8(1), Mar. 2018. doi: 10.1038/s41598-018-23540-x. URL https://doi.org/10.1038/s41598-018-23540-x .[11] S. Enoch and D. J. Leaper. Basic science of wound healing. Surgery (Oxford), 26(2):31–37,2008.[12] P. Friedl and D. Gilmour. Collective cell migration in morphogenesis, regeneration andcancer. Nature Reviews Molecular Cell Biology, 10(7):445–457, July 2009. doi: 10.1038/nrm2720. URL https://doi.org/10.1038/nrm2720 .[13] N. Gal and D. Weihs. Intracellular mechanics and activity of breast cancer cells correlatewith metastatic potential. Cell Biochemistry and Biophysics, 63(3):199–209, May 2012.doi: 10.1007/s12013-012-9356-z. URL https://doi.org/10.1007/s12013-012-9356-z .[14] A. Goriely and D. Moulton. Morphoelasticity: a theory of elastic growth. New Trendsin the Physics and Mechanics of Biological Systems: Lecture Notes of the Les HouchesSummer School: Volume 92, July 2009, 92:153, 2011.[15] J. Guck, S. Schinkinger, B. Lincoln, F. Wottawah, S. Ebert, M. Romeyke, D. Lenz,H. M. Erickson, R. Ananthakrishnan, D. Mitchell, J. K¨as, S. Ulvick, and C. Bilby.Optical deformability as an inherent cell marker for testing malignant transformationand metastatic competence. Biophysical Journal, 88(5):3689–3698, May 2005. doi:10.1529/biophysj.104.045476. URL https://doi.org/10.1529/biophysj.104.045476 .[16] E. Haertel, S. Werner, and M. Sch¨afer. Transcriptional regulation of wound inflammation.In Seminars in Immunology, volume 26, pages 321–328. Elsevier, 2014.[17] I. Hirahara, Y. Ogawa, E. Kusano, and Y. Asano. Activation of matrix metalloproteinase-2 causes peritoneal injury during peritoneal dialysis in rats. Nephrology DialysisTransplantation, 19(7):1732–1741, May 2004. doi: 10.1093/ndt/gfh262. URL https://doi.org/10.1093/ndt/gfh262 .[18] K. Keren, Z. Pincus, G. M. Allen, E. L. Barnhart, G. Marriott, A. Mogilner, and J. A.Theriot. Mechanism of shape determination in motile cells. Nature, 453(7194):475–480,May 2008. doi: 10.1038/nature06952. URL https://doi.org/10.1038/nature06952 .[19] D. Koppenol. Biomedical implications from mathematical models for the simulation ofdermal wound healing. PhD-thesis at the Delft University of Technology, the Netherlands,2017.[20] B. Li and J. H.-C. Wang. Fibroblasts and myofibroblasts in wound healing: force generationand measurement. Journal of tissue viability, 20(4):108–120, 2011.[21] X. Liang, V. Crecea, and S. A. Boppart. Dynamic optical coherence elastography: areview. Journal of innovative optical health sciences, 3(04):221–233, 2010.[22] S. Liu, A. Peyronnel, Q. Wang, and L. Keer. An extension of the hertz theory for 2d coatedcomponents. Tribology Letters, 18(4):505–511, Apr. 2005. doi: 10.1007/s11249-005-3611-z.URL https://doi.org/10.1007/s11249-005-3611-z .[23] M. Mak and D. Erickson. A serial micropipette microfluidic device with applications tocancer cell repeated deformation studies. Integrative Biology, 5(11):1374–1384, Sept. 2013.doi: 10.1039/c3ib40128f. URL https://doi.org/10.1039/c3ib40128f .2124] M. Mak, C. A. Reinhart-King, and D. Erickson. Elucidating mechanical transition effectsof invading cancer cells with a subnucleus-scaled microfluidic serial dimensional modulationdevice. Lab Chip, 13(3):340–348, 2013. doi: 10.1039/c2lc41117b. URL https://doi.org/10.1039/c2lc41117b .[25] S. Massalha and D. Weihs. Metastatic breast cancer cells adhere strongly on varying stiff-ness substrates, initially without adjusting their morphology. Biomechanics and Modelingin Mechanobiology, 16(3):961–970, Dec. 2016. doi: 10.1007/s10237-016-0864-4. URL https://doi.org/10.1007/s10237-016-0864-4 .[26] C. P. McCann, P. W. Kriebel, C. A. Parent, and W. Losert. Cell speed, persistence andinformation transmission during signal relay and collective migration. Journal of CellScience, 123(10):1724–1731, Apr. 2010. doi: 10.1242/jcs.060137. URL https://doi.org/10.1242/jcs.060137 .[27] A. Mogilner and K. Keren. The shape of motile cells. Current Biology, 19(17):R762–R771, Sept. 2009. doi: 10.1016/j.cub.2009.06.053. URL https://doi.org/10.1016/j.cub.2009.06.053 .[28] Q. Peng and F. Vermolen. Agent-based modelling and parameter sensitivity analy-sis with a finite-element method for skin contraction. Biomechanics and Modeling inMechanobiology, pages 1–27, 2020.[29] S. H. Phan. Biology of fibroblasts and myofibroblasts. Proceedings of the AmericanThoracic Society, 5(3):334–337, Apr. 2008. doi: 10.1513/pats.200708-146dr. URL https://doi.org/10.1513/pats.200708-146dr .[30] V. L. Popov. Contact Mechanics and Friction. Springer Berlin Heidelberg, 2010. doi:10.1007/978-3-642-10803-7. URL https://doi.org/10.1007/978-3-642-10803-7 .[31] L. Ritti´e. Cellular mechanisms of skin repair in humans and other mammals. Journal of CellCommunication and Signaling, 10(2):103–120, May 2016. doi: 10.1007/s12079-016-0330-1.URL https://doi.org/10.1007/s12079-016-0330-1 .[32] P. Robey. “mesenchymal stem cells”: fact or fiction, and implications in their therapeuticuse. F1000Research, 6, 2017.[33] P. Rousselle, F. Braye, and G. Dayan. Re-epithelialization of adult skin wounds: Cellularmechanisms and therapeutic strategies. Advanced Drug Delivery Reviews, 146:344–365,June 2019. doi: 10.1016/j.addr.2018.06.019. URL https://doi.org/10.1016/j.addr.2018.06.019 .[34] E. T. Roussos, J. S. Condeelis, and A. Patsialou. Chemotaxis in cancer. Nature ReviewsCancer, 11(8):573–587, July 2011. doi: 10.1038/nrc3078. URL https://doi.org/10.1038/nrc3078 .[35] S. Rudraraju, D. E. Moulton, R. Chirat, A. Goriely, and K. Garikipati. A computationalframework for the morpho-elastic development of molluskan shells by surface and volumegrowth. PLoS computational biology, 15(7):e1007213, 2019.[36] M. Saeed and D. Weihs. Finite element analysis reveals an important role for cellmorphology in response to mechanical compression. Biomechanics and Modeling inMechanobiology, 19(3):1155–1164, Dec. 2019. doi: 10.1007/s10237-019-01276-5. URL https://doi.org/10.1007/s10237-019-01276-5 .2237] K. Safferling, T. S¨utterlin, K. Westphal, C. Ernst, K. Breuhahn, M. James, D. J¨ager,N. Halama, and N. Grabe. Wound healing revised: A novel reepithelialization mechanismrevealed by in vitro and in silico models. The Journal of Cell Biology, 203(4):691–709, Nov.2013. doi: 10.1083/jcb.201212020. URL https://doi.org/10.1083/jcb.201212020 .[38] D. W. Siemann and M. R. Horsman. Modulation of the tumor vasculature and oxygenationto improve therapy. Pharmacology & Therapeutics, 153:107–124, Sept. 2015. doi: 10.1016/j.pharmthera.2015.06.006. URL https://doi.org/10.1016/j.pharmthera.2015.06.006 .[39] A. J. Singer and R. A. Clark. Cutaneous wound healing. New England Journal of Medicine,341(10):738–746, Sept. 1999. doi: 10.1056/nejm199909023411006. URL https://doi.org/10.1056/nejm199909023411006 .[40] V. Swaminathan, K. Mythreye, E. T. O ' Brien, A. Berchuck, G. C. Blobe, and R. Superfine.Mechanical stiffness grades metastatic potential in patient tumor cells and in cancer celllines. Cancer Research, 71(15):5075–5080, June 2011. doi: 10.1158/0008-5472.can-11-0247.URL https://doi.org/10.1158/0008-5472.can-11-0247 .[41] W. R. Trickey, F. P. Baaijens, T. A. Laursen, L. G. Alexopoulos, and F. Guilak. Determina-tion of the poisson ' s ratio of the cell: recovery properties of chondrocytes after release fromcomplete micropipette aspiration. Journal of Biomechanics, 39(1):78–87, Jan. 2006. doi:10.1016/j.jbiomech.2004.11.006. URL https://doi.org/10.1016/j.jbiomech.2004.11.006 .[42] F. J. Vermolen and A. Gefen. A phenomenological model for chemico-mechanically inducedcell shape changes during migration and cell–cell contacts. Biomechanics and Modelingin Mechanobiology, 12(2):301–323, May 2012. doi: 10.1007/s10237-012-0400-0. URL https://doi.org/10.1007/s10237-012-0400-0 .[43] F. J. Vermolen and E. Javierre. A finite-element model for healing of cutaneous woundscombining contraction, angiogenesis and closure. Journal of Mathematical Biology, 65(5):967–996, Nov. 2011. doi: 10.1007/s00285-011-0487-4. URL https://doi.org/10.1007/s00285-011-0487-4 .[44] R. C. Wek and K. A. Staschke. How do tumours adapt to nutrient stress? The EMBOJournal, 29(12):1946–1947, June 2010. doi: 10.1038/emboj.2010.110. URL https://doi.org/10.1038/emboj.2010.110 .[45] J. Zhao, F. Manuchehrfar, and J. Liang. Cell–substrate mechanics guide collectivecell migration through intercellular adhesion: a dynamic finite element cellular model.Biomechanics and Modeling in Mechanobiology, 19(5):1781–1796, Feb. 2020. doi: 10.1007/s10237-020-01308-5. URL https://doi.org/10.1007/s10237-020-01308-5https://doi.org/10.1007/s10237-020-01308-5