A formalism for modelling traction forces and cell shape evolution during cell migration in various biomedical processes

The phenomenological model for cell shape deformation and cell migration Chen (BMM 17:1429–1450, 2018), Vermolen and Gefen (BMM 12:301–323, 2012), is extended with the incorporation of cell traction forces and the evolution of cell equilibrium shapes as a result of cell differentiation. Plastic deformations of the extracellular matrix are modelled using morphoelasticity theory. The resulting partial differential differential equations are solved by the use of the finite element method. The paper treats various biological scenarios that entail cell migration and cell shape evolution. The experimental observations in Mak et al. (LC 13:340–348, 2013), where transmigration of cancer cells through narrow apertures is studied, are reproduced using a Monte Carlo framework.


Introduction
Cells may attain various shapes and sizes, for example, stem cells can differentiate and adopt the shape and functionality of many different cell types in our body: fan-like keratocytes, handshaped nerve growth cones and spindle-shaped fibroblasts [27,32].It has been recognized that cell geometry influences cellular activities like cell growth and death, cell mobility and adhesion to the direct environment [2,18,25,27,36].The shape of a mobile cell is determined by its boundaries, which dynamically vary with a local balance between retraction and protrusion [10].There are multiple constituent elements affecting the cell shape, for instance, the cytoskeleton 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 global cell morphology in the context of its function [18,27].
Signalling molecules play an important role in cell migration and cell shape.During wound healing, chemotaxis is one of the most important cues for migration of immune cells and fibroblasts in inflammatory and proliferative phases [19,8,11,28].Metastasis of cancer cells can be induced by nutrients and oxygen, since tumour growth requires an adequate supply of oxygen and nutrients.Under most pathological circumstances, oxygen and nutrients are supplied though the local blood vasculature [44,38].Commonly, signalling molecules are activated at the plasma membrane, and de-activated in the cytoplasm.On the other hand, the concentration of signalling molecules determines the cytoskeletal dynamics [27].
In wound healing, cells migrate and change shape in both the epidermis and the dermis layers.Re-epithelialization is the most essential part for the skin to re-establish its barrier function [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 the epidermis and dermis extends slightly over the ends of the incised dermis, creating an "extension membrane" (or so-called epidermal tongue) [31].The mechanism of the occurrence of the epidermal tongue is still unclear.A possible explanation is that, the suprabasal cells (which lie upon the layer of basal cells) form the tongue by migrating over the leading basal cells and dedifferentiating to basal cells (which are adhered to the basement membrane between the epidermis and dermis) to form new leaders [37,31,43,33].When epidermal epithelial cells are "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 of the key events during wound healing.Differentiation changes the spindle-shaped fibroblast to dendritic-shaped myofibroblasts.Subsequently, cells' mechanobiology is modified considerably as 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 with disabilities 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 narrow and stiff cavity, it has to deform to adapt to the obstacles.More invasive cancer cells, appear be more pliable and dynamic both internally [13] and externally [15,7,40], and thus able to adjust 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 the substrate, compared to benign cells [25], yet the specific mechanisms the induce these increased forces are still poorly understood.
Mathematical modelling has been proven to be an important tool to have a deeper insight into many biological processes that are potentially difficult to control in experiments, for example, 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 advantage of modelling a larger scale, however, the model neglects the individual cellular activity and cells are not tracked [42].Agent-based model is suitable to model cellular activities of every cell, for instance, cell migration and cell deformation.Hence, an agent-based model is selected in this manuscript and this work is an extension of [6] and [42].In [6], a model of the deformation of both the cell and the nucleus is developed.Furthermore, a parameter sensitivity analysis is carried out on the basis of Monte Carlo simulations.However, the study in [6] does not consider the traction forces applied by the cell and the impact on the substrates.Compared to the work of Vermolen and Gefen [42], we use finite-element methods to solve all the partial differential equations, rather than Green's functions.Therefore, a more precise solution is delivered.Furthermore, we implement a more intricate approach to model the traction forces applied 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 cell migration, in the form of a set of partial differential equations.Possible applications of this model and the corresponding numerical results are exhibited in Section 3. Finally, conclusions are shown in Section 4.

Mathematical Modelling
In this manuscript, the phenomenological model of cell deformation is extended from the work in [6,42], in particular, in two dimensions.With essential biological assumptions and simplifications, the model mainly describes the impact of extracellular components on the cell deformation and displacement.Subsequently, more applications can be developed, for instance, cell differentiation and cell repulsion.Different from the work in [45], the model in this manuscript neglects the intracellular environment, hence, it is not capable to present the Poisson's effect of the cell.
The cell membrane is split into finite line segments by the nodal points, and the centre position of cell is determined by the mean of all the positions of the nodal points.The equilib- rium shape of the cell is kept by a collection of springs, which connects each nodal point on the cell membrane to the centre of the cell, respectively; see Figure 2.1.For each nodal point, the displacement 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 of this model, there will be some model adjustments.

Concentration of Generic Signal
We assume that cell migration is mainly driven by chemotaxis (or mechanotaxis) , which is commonly observed in wound healing and cancer cell metastasis with various types of signalling 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 migration of (myo)fibroblasts towards the wound from the uninjured skin [8,11].Cancer metastasis is triggered by cancer cell proliferation and cell migration.On a cellular level, chemotaxis or haptotaxis is an essentional cue for cancer cell migration and hence for the dissemination of tumours [44].Numerous studies indicate that the availability of oxygen and nutrients is one of most crucial factors for the growth of tumours [34].During tumour growth, the concentration of oxygen and nutrients depletes in the vicinity of the tumour [38].Therefore, cancer cells have a 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 the source of a generic signal.Together with the reaction-diffusion equation, the concentration of the signal is determined by: ∂c(x, t) ∂t where c(x, t) is the concentration of the signalling molecule, D is the diffusion rate, k is the secretion rate of the signal source, x s is the position of the source, and v is the displacement velocity of the substrate that results from the cellular forces exerted on their surroundings.The velocity is computed by solving the balance of the momentum, which will be discussed in the Section 2.2.Initially, we assume there is no signalling molecules over the computational domain, that is, 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 between the boundary and the region far away from the domain of computation.The symbol κ s , which is non-negative, represents the mass transfer coefficient.Note that as κ s → 0 then the Robin condition tends to a homogeneous Neumann condition, which represents no flux (hence isolation).Whereas as κ s → ∞ represents the case that c → 0 on the boundary, which, physically, is reminiscent to having an infinite mass flow rate at the boundary into the surroundings.The Robin condition, also referred to as a mixing boundary condition, is able to deal with both these two limits and all cases between these limits.

Passive Convection of Substrate
In wound healing, (myo)fibroblasts exert forces on their direct environment, i.e. extracellular matrix, which result into contraction of the tissue [8,11,16,20].For cancer cells, Massalha and Weihs [25] indicate that the metastatic cells exert traction forces ranged from 100 − 600nN on the gel, of which the Young's modulus ranged from 2.2 − 10.9kP a.Furthermore, for stiffer substrates, the cancer cells remain rounded with changing area and they exert large traction forces with large magnitudes to its direct environment.Hence, the model includes passive convection 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 are implemented here to depict the forces exerted by the cell, which are applied on the midpoint of each line segment; see Figure 2.2 as an example of a square-shape cell.Among different applications, 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).Morphoelasticity is widely used in the biological modelling to describe elastic growth, for instance, the growth of tumours [14], the seashell growth [35] and large contractions in wound healing [19,3] etc.In wound healing, morphoelasticity describes the phenomena when the deformation 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 following modelling equations [19]: 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 = 0, then the tissue will gradually recover to its original shape and volume.Here, Dy Dt = ∂y ∂t + v∇ • y is material derivative where y is any tensor field and v is the migration velocity of any point within the domain of computation.In order to have a fixed boundary, we use a homogeneous Dirichlet boundary condition for the velocity.From a mechanical point of view, we treat the computational domain as a continuous linear 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 use Kelvin-Voigt's viscoelastic dashpot model, of which the stress tensor reads as where ν s is the Possion's ratio of the substrate, is the strain tensor, µ 1 and µ 2 are the shear and bulk viscosity respectively.The morphoelasticity model solves non-linear equation and both velocity v and strain tensor are unknowns.The deformation of the domain is actually determined by the strain tensor.The displacement of the domain can be approximated by integrating the velocity over time: u(t) ≈ t 0 v(s)ds.In the application of wound healing, (myo)fibroblasts are the cells pulling the ECM and causing the contractions.The temparory force f of each (myo)fibroblast reads as where N is the number of nodal points on the cell membrane, P (x, t) is the magnitude of the force, 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 repelling forces on the other one.Here, cells are not allowed to intersect each other; see Figure 2.3 as a scheme.Suppose for (myo)fibroblasts i, there are N m = {j 1 , . . ., j m } line segments of cell membrane 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)fibroblast releases pulling forces on the ECM.Hence, the temporary force of the (myo)fibroblast i is given by where l m is the portion of the (myo)fibroblast membrane mechanically contacting with other cell, Q(d(x), t) is the force magnitude per length and d(x) is the penetration depth.According to 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): where E * is the total equivalent Young's modulus derived by Here, ν i and E i with i ∈ {1, 2} represent the Poisson ratio and Young's modulus of two elastic bodies, respectively.In particular, if two bodies have the same elastic characteristics (i.e. . We assume that the magnitudes of the repulsive force, which is exerted on the boundary segments of cell i that are in contact with another cell, are identical.In other words, Q(d(x), t) is given by where l m is the total length of the portion of the membrane of (myo)fibroblast i mechanically contacting with other (myo)fibroblast.Subsequently, the total temporary force is f = N C i=1 f i , where N C is the number of (myo)fibroblasts.Another application is metastasis and invasion of cancer cell.Usually, in vitro, a microtube experiment is conducted [24].To simplify the model, here we only consider one cancer cell going 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 the microtube.Here, we will exclude the pulling force.In other words, the force released by the cancer cell is only the repelling force exerted on the wall of the microtube: (2.8) The magnitude of the force here follows the same definition as in Eq (2.7), where d(x) is the radius subtracting the distance from the cell membrane to the cell centre.The cell is compressed by the wall of the microtube.To inhibit any changes from its equilibrium shape, it exerted the repelling force on the wall of the microtube, which is proportion to the compressed distance h.

Cell Deformation
According to [6,42], the cell cytoskeleton is depicted as a collection of springs between the centre position of the cell and the nodal points on the cell membrane.Therefore, the equilibrium shape of 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.Combining chemotaxis (or mechanotaxis), passive convection and random walk, the displacement of the nodal point j on the cell membrane is given by Here, E c represents the cell elasticity; x = x c (t)− xj (t) is the vector connecting the equilibrium position of nodal point i on the cell membrane to the cell centre, x c is the central position of the cell and xj represents the equilibrium position of the nodal point j corresponding to the cell centre x c (see Figure 2.1); γ is a small positive constant to prevent the denominator being zero; v is the velocity of the substrate determined by Eq (2.2); σ rw is the portion of random walk, and dW (t) is a vector-Wiener process, which accounts for random walk.Furthermore, is the weight of chemotaxis (or mechanotaxis), where we define the Cell Shape Index (CSI) of cell Ω C by 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 mobility of cell.For simplicity, we propose a linear relation here: where β 0 is the maximal response from the cell to the signal, µ m is the mobility reduction coefficient, and CSI 0 (Ω C ) and A 0 (Ω 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 rotation of an angle φ, as in [6]: such that φ can be computed from φ = arg min Hence, the displacement of nodal point j is adapted to (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 the microtube.For the nodal point colliding the obstacle, it cannot pass over the boundary of the obstacle.Hence, for the displacement of nodal point j, the normal direction of the boundary of the obstacle in dx j (t) is vanished.To rephrase it, we adjust the displacement of the nodal point if it collides the obstacle by 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 metastasis and invasion, not only the direction of the cell migration is limited (since the cell cannot pass over the microtube), but also the microtube will slow down the velocity of the cell as a result of friction.Note that the magnitude of the friction is proportional to the repelling force.In our 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 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 the cell shape under multiple circumstances.However, the model does not describe the Poisson effect of the cell if the cell is compressed, since the model mainly considers the extracellular environment.Hence, in this manuscript, the cell length will not be investigated.

Applications and Numerical Results
We exhibit several possible applications in this section, namely, cells migrating as a result of chemotactic signals, cells differentiating to another phenotype, cells repelling each other and one 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.

Cells Moving towards the Point Source
The basic application is that cell migrates towards the concentration gradient of a signalling molecule, which can be oxygen, growth factors or virus.The displacement is mainly determined by the gradient of the concentration of the signal.Subsequently, the closest part of the cell to the emitting source will develop a "nose"; see Figure 3.1.All parameter values have been documented 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 adhesion sites over the cell membrane [27].In the simulation we accommodate for the engulfment of the chemical source by switching off the chemical signla once the cell physically contacts it.Once this signal has been switched off, the concentration gradient flattens as a result of diffusion processes 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.From t = 0, the cells are attracted towards the centre of the computational domain,which is the location of the source of signalling molecules.Due to the difference of the gradient of 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 signal proceeds, the "nose" gradually disappears and cells recover to the equilibrium shape.To evaluate the cell geometry quantitatively, we provide the evolution of the CSI and cell volume as a 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 from the 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 original shape but not the volume.It can be implied that a stiffer cell deforms less compared to a softer one.

Differentiation of Cells
Cell differentiation is a process of a cell changing from one phenotype to another one, for example, a stem cell differentiates into various phenotypes, like blood cells and nerve cells etc.In this manuscript, we mainly focus on the cellular differentiation in wound healing.In the 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 the skin ever harder and cause the permanent contractions.It is commonly recognized that high concentration of TGF-beta induces the fibroblast-to-myofibroblast differentiation [18,8,9].In this section, since we mainly want to present a model of differentiation, only differentiation from fibroblasts 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 the    equilibrium status.Here, for ellipse and hypocycloid, there are two parameters to determine each shape: long (denoted by a e ) and short axis's (denoted by b e ) determine the ellipse; the radius of basis circle (de noted by a h ) and rotating circle (denoted by b h ) determine the hypocycloid.Note that the hypocycloid-shaped cell may not be realistic, and it is mainly to show that the model is capable to model the differentiation of cells, in which mostly the cellular skeleton and geometry are altered.To develop a smooth differentiation process, we introduce a function such that each parameter changes over time: 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 the time point when the fibroblast starts differentiating.Figure 3.3 presents the cells positions at different time.As regular fibroblasts approach to the region with high concentration of TGF-beta, some of them start differentiating into myofibroblasts gradually.Subsequently, they exert larger forces on the ECM and contractions are developed in the wound, which is marked with red curves as a subdomain.The parameter values of the simulations are from Table 3.1 and Table 3.3.

Repulsion between Two Colliding Cells
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 points on 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, we present the cell positions at different times, and cells deform due to mechanical contact (hard impingement).The parameter values are from Table 3.1 and Table 3.2.

Cell Moving through a Microtube
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 needs to 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 and invasion of the cancer cell.In this section, we will use our model to reproduce results of the experiment in [24].

Simulation Settings
Following the settings in the experiment, we define a microtube with a varying width: a 15µm larger channel (LC) and 3.3µm subnucleus barriers with length 10µm (SNB10).Since the main reason for the active migration of the cell is not evident in [24], we keep on using random walk, with either chemotaxis or fixed velocity of the cell.Rather than having a periodic setting of SNBs, we have one SNB in the middle connected with two LCs and run the simulations respectively.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 cell migrates through the microtube.Since the repelling force on the wall of the microtube is included, we investigate the results regarding the cell velocity, pressure and the cell shape index over time; see Figure 3.6.The parameter values are taken from Table 3.1 and 3.4.Initially, there is a short distance before the cell enters the microtube, here the cell encounters no distraction.Therefore, the the cell travels at maximal speed and the cell is not compressed 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 the microtube, it slows down due to the friction, and the cell is compressed, therefore, the cell starts exerting pushing forces on the wall of the tube.In the LC part, the cell shape index stays stable around 0.95.Further, the cell approaches the SNB, which is much more narrow than 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 is recorded when the cell is in the SNB.After that, cell moves further towards the signal source through the LC again.Hence, the cell velocity and cell shape index increase again, while the cell pressure reduces.According to [1] and [26], we manage to keep the cell velocity and pressure in a reasonable range: 6 − 20µm/min and the maximal pressure that a cell can handle is around 12kP a.

Monte Carlo Simulations
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 it is 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 be collected multiple times with one individual cell), one cell is supposed to go through one set of the microtube in each simulation.To reproduce the experimental results, Monte Carlo simulations are conducted to estimate the probability of the occurrence of phase 3 and the time cost for each phase, with different aforementioned reasons of active migration.The input values for the Monte Carlo simulations are shown in Table 3.5.We run the simulations with four assumptions of the main mechanism provoking the active cell displacement: chemotaxis, fixed velocity with 10µm/min, velocity generated from (6, 15)µm/min and (6, 20)µm/min in horizontal direction according to [26].The number of samples and the Monte Carlo error of the occurrence of phase 3 collected from the Monte Carlo simulations of each aforementioned category are displayed in Table 3.6.Figure 3.7 illustrates the probability of the occurrence of phase 3, which is the stage when the cell stops monotonic forward migration.The mechanism that makes the cell move forward is not clear in [24], hence this 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. Coincidently, 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.3.6 for more information).The parameter values are taken from Table 3.5.Furthermore, the time cost of each phase is recorded.The results with different modelling settings or simulations do not show many differences, in particular between chemotaxis and fixed velocity v = 10.The reason is that the microtube restricts the displacement of the cell in the 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 results of the occurrence of phase 3, in particular, phase 1 and phase 3. Therefore, to investigate the possible reasons of mismatching results in phase 1 and phase 3 in phase time, we reran the simulation with the same settings as Simulation 4 in Table 3.6, except for the cell stiffness modified to E c = 1.The results are shown in Fig 3 .9.With a softer cell, the simulation data in phase 3 match better with the experimental data.However, now a discrepancy between simulation and experiments results in phase 4 instead.(6,20).Red dots with error bar represent the experimental data from [24] and the box plots are the data collected from the simulation.
There are several possible reasons causing the discrepancy.Firstly, for phase 3, we only obtain valid data when the cell moves non-monotonically, which results into the reduction of the sample size.Secondly, the length of LC is not stated clearly in [24], therefore, we could only estimate that from the scale in the figures.Thirdly, it is not clear if the velocity of active migration 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 still managed to reproduce the results which are close to the experimental results.Fourthly, the transaction of each phase from [24] to our simulations may cause a mismatch of the duration time of each phase.Fifthly, it has been observed in [23] that after the first time moving through the narrow channel, cells deform easier to move faster through the following narrow channels, which may indicate that the cell characteristics change regarding its geometry.

Conclusions
A phenomenological model for cell shape evolution and migration has been developed, with primary focus on the mechanics of the extracellular environment.Furthermore the impact of passive convection, due to local displacements within the extracellular matrix, on the evolution of the cell shape has been taken into account.A morphoelastic model has been used in the current study to incorporate permanent deformations of the extracellular tissue.The model can be applied to mimic several microscopic biological observations such as cell deformation and migration during wound contraction and cancer metastasis.To validate the model, the experimental setup in [24] has been modeled.This experiment entailed cell migration through microtubes with different widths and with a varying width over the length.The model is able to reproduce the most important trends that were observed in the experimental data despite some experimental uncertainties such as the determination of which phase a cell is in during the transmigration process.Furthermore, the current model provide a basis that can be expanded to describe more experimentally observed phenomena in cell geometry.

Figure 2 . 1 :
Figure 2.1: A schematic of the distribution of the nodal points on the cell membrane.In our model, we assume the equilibrium shape of the cell is maintained by a collection of springs.

Figure 2 . 2 :
Figure 2.2: An example of pulling point forces, tractions applied by a cell, in which the cell membrane is split by four nodal points.

Figure 2 . 3 :
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 overlapping membrane of both cells.

Figure 2 . 4 :
Figure 2.4:The cell is compressed by the wall of the microtube.To inhibit any changes from its equilibrium shape, it exerted the repelling force on the wall of the microtube, which is proportion to the compressed distance h.

Figure 3 . 1 :
Figure 3.1: The figure shows the positions of cells at consecutive times.The red circles show the original shape and position of cells, and the red dots in the centre is the point source of a generic signal.

Figure 3 . 2 :
Figure 3.2: The cell shape index and relative ratio of cell volume of all cells in Figure 3.1 are shown in the plot.The solid curves represent the cell volume, and the dashed curves are the cell shape index.Different colours of curves indicate different stiffness of the cells.

Figure 3 . 3 :Figure 3 . 4 :
Figure 3.3: The figure shows the positions of cells at consecutive times by the black curves.The red circles show the original shape and position of cells, and the red dots in the centre is the point source of TGF-beta which triggers the differentiation from regular fibroblasts to myofibroblasts.

Figure 3 . 5 :
Figure 3.5: The figure shows the positions and shapes of cells at consecutive times by the blue contours when it travels through a microtube.The red circles show the original shape and position of cells, and the red dots in the end of the microtube is the point source of a signal.

Figure 3 . 6 :
Figure 3.6:The cell velocity, pressure and shape index over time when the cell migrates through the microtube, where there is a 10µm subnucleus barrier.The simulation mimics the experiment in[24].

Figure 3 . 7 :
Figure 3.7: The probability of the occurrence phase 3 in [24] and from Monte Carlo simulations by implementing different mechanisms of cell active displacement (see Table3.6 for more information).The parameter values are taken from Table3.5.

Figure 3 . 8 :
Figure 3.8: The figure shows the time cost of each phase in[24] and from Monte Carlo simulations 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.

Figure 3 . 9 :
Figure 3.9: The time cost of each phase from the Monte Carlo simulation, in which the cell stiffness is E c = 1 and the cell velocity is randomly generated from(6,20).Red dots with error bar represent the experimental data from[24] and the box plots are the data collected from the simulation.

Table 3 .
1: Parameter values used in all the applications

Table 3 .
2: Parameter values estimated in the application of cell migrating toward the signal source

Table 3 .
4: Parameter values used in the application of cell going through a microtube

Table 3 .
5: Parameter values used in the application of cell going through a microtube

Table 3 .
6: Monte Carlo simulations of various models, in which the main mechanisms of cell active displacement differ