Abstract
Cell migration plays an essential role in cancer metastasis. In cancer invasion through confined spaces, cells must undergo extensive deformation, which is a capability related to their metastatic potentials. Here, we simulate the deformation of the cell and nucleus during invasion through a dense, physiological microenvironment by developing a phenomenological computational model. In our work, cells are attracted by a generic emitting source (e.g., a chemokine or stiffness signal), which is treated by using Green’s Fundamental solutions. We use an IMEX integration method where the linear parts and the nonlinear parts are treated by using an Euler backward scheme and an Euler forward method, respectively. We develop the numerical model for an obstacleinduced deformation in 2D or/and 3D. Considering the uncertainty in cell mobility, stochastic processes are incorporated and uncertainties in the input variables are evaluated using Monte Carlo simulations. This quantitative study aims at estimating the likelihood for invasion and the length of the time interval in which the cell invades the tissue through an obstacle. Subsequently, the twodimensional cell deformation model is applied to simplified cancer metastasis processes to serve as a model for in vivo or in vitro biomedical experiments.
Introduction
Cell locomotion is closely involved in various physiological and pathological processes. For example, migration of leukocytes is important for the inflammatory response and movement of fibroblasts and also vascular endothelial cells are essential for wound healing (Lauffenburger and Horwitz 1996). On the contrary, cell migration can play a detrimental role during cancer metastasis, where the dissemination of cancer cells initializes the invasionmetastasis cascade as introduced by Chambers et al. (2002b), Fidler (2003), Lambert et al. (2017).
The diversity of cancers exceeds 200 distinct disease entities, which have differences in the normal cells of origin and similarities in subsequent cancer metastasis. Compared to primary tumors, metastatic cancers cause the overwhelming majority of cancerassociated deaths as high as 90% (Lambert et al. 2017; Seyfried and Huysentruyt 2013; Gupta and Massagué 2006). During the metastatic spreading of tumors, cancer cells can undergo transitions between two forms of movement, which are the amoeboid mode and the mesenchymal mode to optimize their invasiveness (Paňková et al. 2010; Sahai and Marshall 2003). Moreover, Pinner and Sahai (2008) observe that cancer cells are able to move quickly (up to 15 \(\upmu \mathrm{m/min}\)) like some leukocytes and rapidly change their shapes and directions of migration in an amoeboid manner with intravital confocal microscopy technology. Amoeboid movement could happen in the absence of matrix protease (Wolf et al. 2003; Wyckoff et al. 2006) where cancer cells alternatively generate large contractile force pushing fibers of matrix away and squeeze between small paths. However, if the contractile force is insufficient to deform the stiff extracellular matrix (ECM), the matrix metalloproteases (MMP’s) will be secreted by cancer cells to degrade the ECM and thereby invade further (Kalebic et al. 1983; Wolf et al. 2013). In summary, cancer cells frequently chemically and/or mechanically ‘dig’ their ways through ECM in order to reach the distinct parts of the body.
When a single cancer cell is metastasizing through a narrow cavity, it must deform its morphology by extending its membrane into an elongated protrusion; this is often driven by external signals such as chemotaxis, durotaxis or tensotaxis. Large cell deformations will also induce changes in the nucleus morphology. Extensive deformation of the nucleus can induce damage, and reduce the nuclear envelope integrity, see for instance the work by Denais et al. (2016). However, the cancer cell is also capable of repairing its ruptured nuclear envelope and damaged DNA after the penetration. Then, the cell may be able to further promote cancer development. Thus, as noted by Denais et al. (2016), the stage of nuclear envelop rupture could represent a particularly fragile point, thereby providing an opportunity to develop new antimetastatic cancer drugs to inhibit DNA repair and increase cell death. Cell deformation during cancer metastasis has been difficult to study in detail both in vivo and vitro, and further understanding of cell deformation mechanisms is crucially important. In cases where the pore sizes are much smaller than the size of the nucleus, the nucleus mostly arrests and fails to penetrate the pore due to a defective nuclear deformability. On the contrary, with pore diameters above a threshold, e.g., 7 \(\upmu \mathrm{m}\) in the work (Wolf et al. 2013), MMPindependent migration in dense ECM relies on the hourglassshaped deformation of the nucleus. Hence, in our current work, we develop a mathematical model to investigate the correlation between the deformation of a cell and of its nucleus, and show the dynamic changes in cell mechanostructure that occur during the invasion process.
Mathematical modeling has been shown to be an important tool to quantify the relations in many biomedical processes such as wound healing, cell migration and tumor progression in various scales. Cell deformation and migration models exist in the colony scale, e.g., in the works by Rey and GarciaAznar (2013), Byrne and Drasdo (2009), and Vermolen and Gefen (2012), where the cell geometry is fixed to be circular or spherical, respectively, in two and threedimensional simulations. On a smaller scale, one looks at the deformation of individual cells, and to this extent, cellular automata models have been developed and combined with finiteelement strategies by Borau et al. (2014) and Oers et al. (2014). Other cell deformation models are based on phasefield models, like in the work by Marth and Voigt (2014), or on viscoelasticity with moving meshes as in (Madzvamuse and George 2013). A phenomenological approach to cell migration and deformation is proposed in Vermolen and Gefen (2013) and Vermolen et al. (2014), wherein the latter work cell migration and deformation have been modeled in relation with the immune response system where white blood cells migrate out of the venules and transmigrate through the venule walls to chase and engulf pathogens. Moreover, Odenthal et al. (2013) introduce a deformable cell model to describe the mechanical communication among the interacting cells and between the cell and its environment. Another deformable model regarding the interactions with emphasis on the relationship between varying matrix geometries and adhesion, contractility as well as cell velocity can be found in (Tozluoğlu et al. 2013). In terms of the nucleus deformable models, MoussaviBaygi et al. (2011) establish a coarsegrained model of the nuclear pore complex to simulate the nucleocytoplasmic transport. As the increasing attention in the cell mechanics, agentbased models are booming, see (van Liedekerke et al. 2015), where three types agentbased models are described.
Cao et al. (2016) develop a chemomechanical model to investigate the impacts of transmigration through confined interstitial spaces on the geometrical and mechanical features of the cell nuclei. In their model, the shape alterations of the cell and nucleus during the transendothelial migration driven by actomyosin contraction force can perturb genomic organization, which in turn affects the behavior of the cells. More nuclear profiles regarding chromatin deformations and nuclear envelope deformations during transmigration are further investigated. This mechanical model successfully predicts the morphological evolution when one cell transmigrate an endothelial gap (Cao et al. 2016). In comparison, our model extends the process and behavior of cell transmigration driven by a chemical/stiffness signal during cancer metastasis, whereas the most inner cellular mechanical properties are neglected for sake of simplicity.
None of the aforementioned studies, however, have taken into account the Monte Carlo uncertainty quantification in the cell deformation modeling. Our work aims at modeling the interaction between cell deformation (due to migration) and the deformation of the nucleus as well as quantitative analysis of unknown parameters by Monte Carlo simulations. We quantify the correlation between nuclear deformation relaxation and the cell’s ability to penetration through narrow passages, which is important in the context of metastatic invasion. Section 2 describes the mathematical model in terms of the equations, subsequently, the numerical method is presented in Sect. 3, which is followed by the description of the results in Sect. 4. Finally, conclusions are drawn in Sect. 5.
The mathematical model
This section introduces the model in terms of the mathematical relations. We start with the deformation of cell and its nucleus in two dimensions and extend the formalism to three spatial dimensions subsequently. Moreover, the model is applied to simplified physiological transmigration of cancer cells and six parameters are studied by Monte Carlo simulations.
The model in two dimensions
The nucleus must move in coordination with the cell cytoskeletal dynamics at the front edge and rear end (Friedl et al. 2011). To mimic this cell’s cytoskeleton, a cell is treated as a collection of 30 parallel nodal points that are located on the cell membrane and on the outer boundary of the cell nucleus. We have compared the number of nodal points N (\(N = 10, 30, 50, 100\)) and we found that if the cell is freely moving that the pattern is hardly influenced by the number of springs, whereas the CPU time increases proportionally with the number of springs. If the number of springs is very large, then the time step needs to be adjusted if the cell is in contact with an obstacle. In particular, it may happen if the resolution is too high that the nodal points on the cell boundary overtake each other when they are in (partial) contact with a rigid boundary. Taking the model in Fig. 6 as an example (no perturbation of the random walk), the CPU time and penetration time \(\tau \) are compared with various N in Table 1. The table shows that CPU time increases, whereas the cell penetration time \(\tau \) is comparable with the increase of N. Each node on the cell membrane is connected to its corresponding node on the surface of the cell nucleus. On each of the nodes on the cell membrane surface, an external signal, such as a concentration gradient in the case of chemotaxis or durotaxis, is computed. This signal determines the movement of the nodal point. Next to this signal, the migration of the nodal point is determined by its position relative to its corresponding point of the nucleus boundary via the deformation relaxation of the cell’s cytoskeleton. In this way, the deformation and migration of the cell is modeled and sketched in Fig. 1.
We consider a generic signal, of which the gradient determines the migration of the nodal points on the cell boundary membrane. This signal could be the extracellular stiffness or the concentration of a chemoattractant or a light intensity for instance. In the work by Massalha and Weihs (2017), the gelstiffnessdependent differences among cells with various metastatic potentials have been observed to be correlated with cancer invasiveness, where the metastatic cells apply a wide spectrum of traction forces (100–600 nN) for their adhesion to a stiffer gel. For the sake of presentation, we denote the intensity of the signal by \(c(t,\mathbf{x})\), where t and \(\mathbf{x}\), respectively, denote time and spatial position. The signal, as well as its gradient, can be obtained from a given relationship in which the gradient is determined either analytically or numerically. A numerical evaluation in a finiteelement framework could be carried out by for instance gradient recovery techniques or by mixed finiteelement formulations. In the present paper, we consider a chemical attractant, such as a generic growth factor that attracts the cells. For the sake of illustration, we consider a point source since this allows for a simple treatment using Green’s Fundamental solutions. To this extent, let the emitting source of the chemoattractant by positioned on \(\mathbf{x}_S\), then on an unbounded domain, we solve
Here, D and \(\gamma _S\) represent the diffusion coefficient of the chemokine and the secretion rate of the source. Moreover, \(\delta \) is the Delta Dirac function, while \(\mathbf{x} (t)\) and \(\mathbf{x}_S(t)\) denote positions of the nodal points on the cell membrane and the source. The fundamental solution to this equation in an unbounded domain is used in two spatial dimensions as follows,
In the presence of multiple cells, the superposition principle is used to construct the solution. We note that the signal can be taken as generic as one wishes. The above equation just serves as an illustration. Note that the above equation predicts negative values for the concentration if the distance between the point of observation and the source is too large. In our simulations, the distances are such that the above expression predicts nonnegative values only. Let the set \(\mathbf{x}_i(t)\) and \(\mathbf{x}^n_i(t)\), respectively, denote the nodal points on the cell boundary membrane and on the surface of the nucleus of the cell. Then, the migration of the nodal points on the cell boundary membrane is determined by
Here \(\hat{\mathbf{x}}_i\) represents the vector connecting the initial position of nodal point i on the cell boundary membrane to the initial position of point i on the cell nucleus (see Fig. 1). This vector defines the equilibrium cell shape. In this text, we only consider circular and spherical cells, however, this formulation allows to consider cells of generic shapes such as dendritic shapes. Furthermore, \(\beta \) stands for the cell’s response to external signals, and \(\alpha > 0\) denotes the cell’s deformation relaxation. Over a spectrum of cell types, the mobility of the cell boundary has a locally persistent random character (Lauffenburger and Horwitz 1996), thus the last term takes care of the randomness movements of each node, where \(\eta \) is a constant and \(\mathrm{d}{} \mathbf{W}(t)\) denotes a vector Wiener process with independent samples from a normal distribution with zero mean and variance \(\mathrm{d}t\). The above equation warrants convergence to the equilibrium cell shape if there is no external stimulus for the deformation and migration of the cell.
Next, we introduce the equation of motion for the nodal points on the surface of the nucleus. We proceed similarly to the previous treatment of the nodal points on the cell boundary membrane, where we link the positions of the nodal points on the surface of the nucleus to their counterparts on the boundary membrane as well as to the position of the midpoint of the cell nucleus. To this extent, we obtain for \(i \in \{1,\ldots ,N\}\)
Here \(\alpha ^n\), \(\mathbf{x}_c\) and \(\hat{\mathbf{x}}_i^n\), respectively, stand for the deformation relaxation of the nucleus, the position of the center of the nucleus and the vector connecting the initial position of point i on the surface of the nucleus to the initial center of the nucleus (see Fig. 1). Furthermore, the random character of the mobility of the boundary of the nucleus has been taken into account. This treatment of the points on the surface of the nucleus provides the interaction between the nucleus and the cell membrane. However, this interaction such that the deformation of the nucleus is delayed and damped with respect to the deformation of the membrane.
In order to maintain the right orientation of the cell, we also introduce the rotation matrix after rotation of an angle \(\phi \) relative to the xaxis:
which transforms a vector \(\mathbf{x} \in \mathbb {R}^2\) to
The rotation matrix \(B(\phi )\) is used to determine the new equilibrium points of the cell boundary membrane and of the surface of the nucleus. Therefore, one cell is able to converge to its initial shape as well as to its rotation as a result of its migration to simulate the cell morphological polarization, see Fig. 2 for a sketch. The angle \(\phi \) is determined such that it is closest to the current position of all the nodal points on the cell boundary membrane:
where \(\tilde{\mathbf{x}}_i\) represents the initial position of the ith node on the cell membrane surface with the cell center position at time t. After the above problem has been solved, then the angle of rotation of the cell with respect to the xaxis is known. This angle, \(\phi \), is substituted into the equations of motion for all nodal points on the cell membrane surface and on the surface of the cell nucleus, which gives for \(i \in \{1,\ldots ,N\}\):
and
Next we consider the treatment of an obstacle. Imagine that the surface or contour (in 2D) of the obstacle is given by \(\partial \Omega \) and let the unit normal vector be given by \(\mathbf{n}\), then we require that the component of the migration vector, \(\mathrm{d} \mathbf{x}_i(t)\) has no component in the normal direction of the obstacle’s surface, hence we require that the inner product of \(\mathrm{d} \mathbf{x}_i(t)\) and \(\mathbf{n}\) vanishes, that is
From this we subtract the component of \(d \mathbf{x}_i(t)\) in the direction of \(\mathbf{n}\), hence this gives the following adjustment
Note that herewith the obstacle slows down the migration of the cells. This principle is also applied if cells are colliding into each other. The current model simplifies the mechanics of the cell considerably. Inertial effects would change equations (3, 4, 8, 9) into second order equations with respect to the time derivative. The first term with the firstorder time derivative is generally associated with friction or damping. Since in most studies inertia is neglected compared to friction terms (Galle et al. 2005; Drasdo et al. 1995; Odell et al. 1980), we are faced with a system of firstorder differential equations. Note that the incorporation of more complex mechanics also increases the parameter space of the model, where input parameters are often hard to get.
Extension to three spatial dimensions
Chemotaxis migration is modeled by using the Green’s function as a solution of Eq. (1). However, compared to the 2D model, the Green’s functions in 3D changes to,
where both of \(\mathbf{x (t)}\) and \(\mathbf{x}_s(t)\) have x, y and z components.
The surface of the outer membrane and the nuclear surface are divided into mesh points. For this case, superpositions of the threedimensional Green’s Fundamental solutions are used, as well as the same principles for collision with obstacles and other cells. Further, the rotation can be imposed around all the three coordinate axes, and to this extent, the rotation matrix \(B(\phi )\), entailing a rotation about the zaxis, where \(\phi \) denotes the angle with respect to the xaxis, we extent the rotation matrix to all three coordinate axes:
Here \(B_q(\phi _q)\) denotes the rotation matrix about the q–axis (\(q \in \{x,y,z\}\)), given by
All other principles remain the same and rotation is determined using a minimization with respect to the three coordinate angles.
The application to cancer metastasis
The ECM has preexisting pores (diameter varies from 1 to 20 \(\upmu \mathrm{m}\)) or fiberlike (ranging from less than 3–30 \(\upmu \mathrm{m}\) in width) and channellike (varying from 100 to 600 \(\upmu \mathrm{m}\)) tracks (Paul et al. 2017). Furthermore, cells are viscoelastic objects such that morphological deformation happens frequently during the cancer invasion process (Mak and Erickson 2013). Thence, our model for cell and nucleus deformation is applied to a simplified process occurring during cancer metastasis in a pore and a channellike microenvironment. During the invasion, a cancer cell is normally able to squeeze obstacles like cells, tissue, capillarysized vessels and deform itself as well as its nucleus to penetrate and seed in other organs. In the model, we use a constraint cavity with varying roughnesses to simulate the microenvironment during cancer spread, which is mimicked using a trigonometric function as follows.
where y depicts the rough tube bounds and \(y_0\) is a constant used to adjust the width of a tube. Through changing \(\epsilon \) and \(\omega \), different roughnesses (varying amplitudes and frequencies) can be simulated.
The numerical method
Time integration
We describe the twodimensional case and provide information for the 3Dcase if it substantially differs from the 2Dcase. Initially the cell outer membrane surface is divided into N mesh points with respect to the cell center located on \((x_c,y_c)\) as follows
where we assume the cell to be circular in 2D with radius R. The counterparts of the mesh points on the nuclear surface are given by
where \(R^n < R\) represents the radius of the cell nucleus. The 3D spherical cell is described analogously for \(i \in \{1,\ldots ,N\}\) and \(j \in \{1,\ldots ,M\}\)
and on the nuclear surface by
These initial values can be applied to the multicell configuration similarly. To determine the positions of the nodal points on the outer membrane surface, we use an IMplicitEXplicit (IMEX) time integration to update the positions at the next time step in such a way that the linear parts are treated in an Euler backward method, whereas the nonlinear parts are treated using a forward Euler method. This treatment has been chosen to avoid the need of solving a nonlinear system using an iterative procedure. This treatment results into the following equation (for the singlecell twodimensional case) for the nodes on the outer membrane
and
for the nodes on the nuclear surface, for \(i \in \{1,\dots ,N\}\). Here, \(\Delta \mathbf{W}\) is a twodimensional Wiener process with variables from a normal distribution with zero mean and \(\Delta t \) variance. For the definition and introduction of the vector Wiener process, one can refer to Steele (2012). For the gradient of the concentration (or any other signal that triggers cell migration and deformation), we use the following IMEX convention based on the Green’s Fundamental solutions in 2D (in 3D analogously)
Cell shape
In order to compute the coordinate of the cell center of mass, we need the area or volume of the cell and nucleus. The area A(t) in 2D is computed by realizing that the cell is a polygon, which follows from
For the 3D counterpart, we divide the cell into triangles (in order to allow any finiteelement surface mesh), and compute the volume V(t) of the cell by
where \(N_{\mathrm{el}}\) denotes the number of triangles that are used to approximate the cell (or nuclear) surface, and \(j_1,~j_2\) and \(j_3\) refer to the indexes of the vertices of triangle j. Further, \(\frac{1}{2} \Delta _j\) denotes the area of the jth triangle, where we compute \(\Delta _j\) by
and the unit outward normal vector by
and hence to compute \(\Delta _j n_x\), it suffices to take the x–coordinate of \((\mathbf{x}_{j2}\mathbf{x}_{j1} ) \times ( \mathbf{x}_{j3}  \mathbf{x}_{j1} )\).
For the 3Dcase, we note that the area is computed by summing the areas of all the triangles, that is
The Monte Carlo simulations
In our model, most experimental data is difficult or even impossible to collect, therefore, we refer to other literature data or estimate the input data and thereby evaluating the quantification of the propagation of uncertainty in the variables is very important. To investigate the output influence and correlation among variables, Monte Carlo simulations are carried out based on the model of cancer metastasis. There, a cell transmigrates through a narrow rough tubular path to get from one part of the surrounding tissue to another part. Passage through the tube requires deformation of the cells’ cytoplasm and nucleus and affects the corresponding penetration time \(\tau \) which is quantified under different conditions.
Suppose the variable \(X \in \{D, \beta , \alpha , \alpha ^n\}\) follows a normal distribution \(X \sim N(\mu , \sigma ^2)\), where \(\mu \) and \(\sigma \) represent the mean of the distribution and the standard deviation. Then, the stochastic variable X could be generated by
here \(N_\mathrm{s}\) denotes the number of samples. The strength of the linear association between every variable and penetration time \(\tau \) is quantified by the correlation coefficient r given by
Note the correlation coefficient is always bounded by \([\,1, 1]\), where \(\,1\) or 1, respectively, indicates a perfect negative or positive linear correlation.
Error analysis
Numerical methods yield approximate results, where the numerical error E arises from the IMEX method and the Monte Carlo simulations. The IMEX time integration error \(E_{\mathrm{ti}}\) is defined by
here C represents a positive constant, \(\hat{\tau }\) and \(\hat{\tau }^{\Delta t}\) denote the real mean penetration time and numerical mean penetration time, respectively. The numerical result becomes accurate with the limitation of a sufficiently small time step \(\Delta t\). Furthermore, the accuracy of the Monte Carlo simulations depends on the number of samples \(N_\mathrm{s}\) and this error \(E_{\mathrm{mc}}\) is achieved by
where \(S_n\) denotes the sample standard deviation and \(\hat{\tau }_{N_\mathrm{s}}^{\Delta t}\) is the sample mean as a result of \(N_\mathrm{s}\) samples, which is \(\hat{\tau }_{N_\mathrm{s}}^{\Delta t} = \frac{\sum _{j=1}^N\tau _j^{\Delta t}}{N_\mathrm{s}}\). Here \(\tau _j^{\Delta t}\) denotes the penetration time of sample j, and the sample standard deviation is given by \(S_n = \left[ \frac{\sum _{j=1}^{N_\mathrm{s}}(\tau _j^{\Delta t}\hat{\tau }_{N_\mathrm{s}}^{\Delta t})^2}{N_\mathrm{s}1}\right] ^{\frac{1}{2}}\). Note that this error decreases with increasing number of trials. Therefore, the total error E is given by
To keep the numerical approximation as accurate as possible, the time step should be small enough and number of samples should be sufficient large. Take the Monte Carlo simulations with six parameters as an example. If we fix all the six parameters and set the random walk parameter to zero, then the computation is fully deterministic. The time step is 0.0001 h and the constant C can be estimated using Richardson error estimation by
where \(C = \frac{\hat{\tau }^{\Delta t}\hat{\tau }^{2\Delta t}}{\Delta t} = 133\), which is a mean value of ten times calculations. With 10,000 Monte Carlo (with sampling in the six parameters and using random walk) samples, the analytical error analysis can be derived as,
where 0.0687 is the sample standard deviation \(S_n\). Therefore, the total error in the Monte Carlo simulations with six parameters is bounded by 0.014 h.
The numerical simulations
First, we describe the simulations in which one cell migrates toward the gradient of an increasing stimulus along obstacles in 2D and 3D. Subsequently, this deformation model of cell and its nucleus is applied to a simplified cancer metastasis phenomenon. Furthermore, six parameters are studied and analyzed by Monte Carlo simulations.
Parameter values
Most often the experimental parameter values are not available to us, therefore estimating input values based on experimental literature is essential. For example, we use 10 \(\upmu \mathrm{m}\) in 2D and 16 \(\upmu \mathrm{m}\) in 3D for diameters of the nucleus referring to the work by Friedl et al. (2011), where the diameter of the nucleus varies from 10–20 \(\upmu \mathrm{m}\) in 2D and 5–15 \(\upmu \mathrm{m}\) in 3D. Analogously, other default input values are listed in Table 2, as well as the sources from the literature whenever possible.
Cell migration along a rigid object in 2D and 3D
One cell migrating along a rigid object in 2D
In solid tumors, cell migration shows trends in its direction according to the presence of chemotactic gradients or other external cues. Since there are many parallels existing in the mechanisms underlying the movement of cancer cells and immune cells within tissues as well as in the blood circulation (Pinner and Sahai 2008), the modeled cell can be an immune cell with a chemical source of antigen or a cancer cell with a source of oxygen or substrate/ECM stiffness.
The cell moves according to the gradient of chemokine. Snapshots at different stages of the migration are shown in Fig. 3, where the red, green and gray objects visualize the cell, nucleus and a rigid obstacle, respectively. Furthermore, the signal source location is represented by an asterisk. To pass a stiff barrier or overcome an obstacle, the migrating cell has to reshape and adapt the mechanostructure of the cytoplasm and the membrane. That is done via exerting contractile forces or withstanding the stresses from neighbor cells, which are mediated by the cell cytoskeleton (Brunner et al. 2006). According to the experimental observation of Brunner et al. (2006), one migrating cell could push a small obstacle upward by exerting forces and crawl underneath this obstacle. Given a larger obstacle in our simulation, the cell and nucleus are more likely to crawl along the rigid boundary by morphological adjustments to different extents. Ultimately, the cell and nucleus are able to return to their initial shapes due to cell polarity once the source is no longer active.
One cell moving along a rigid object in 3D
In threedimensional interstitial tissues, cells typically utilize one of two mechanisms for invasion: mesenchymal or amoeboid, respectively, involving degradation of the surrounding ECM or squeezing through subcellsize pores in the ECM; these mechanisms require, respectively, proteinases that can degrade the ECM or deformations of the cell shape. (Friedl et al. 2011). To simplify the problem and make it timeefficient to solve, we only consider the mechanical deformability of the cell in this model rather than deformability of both the environment and the cell. The degradation of the ECM is hence modeled implicitly in the \(\beta \)term in Eq. 3. Note that if the ECM decayrate would be zero, then the \(\beta \)parameter would be zero as well. Hence the \(\beta \)parameter accounts for the decay of ECM and the mobility of the node. In a future study, the decay process of the ECM could be modeled more explicitly so that the migration and deformation process of the cell can be modeled to be ratedetermined by the slowest process. We model a 3D cell with a spherical equilibrium geometry that travels over an obstacle toward a source that secretes a chemokine, e.g., for immune cells the source may be a pathogen. In Fig. 4, consecutive snapshots of a 3D cell that reaches a source and engulfs are shown. It can be seen that the cell deforms mechanostructurally and that the cell shape returns to its equilibrium spherical shape once the stimulus has been removed. Note that the cell is still attached to the obstacle once the source has disappeared. Due to this mechanical attachment and cell elasticity, the cell deforms back to its equilibrium and thereby pushes itself away from the obstacle such that there is only attachment at one point of the cell boundary to the obstacle. This is a characteristic of the current model in which steadystate adherence has been neglected. The figures illustrate how the model takes into account the hard mechanical impingement between the cell and the rigid obstacle.
In general, dimensionality does not affect the expected numerical result in this case. Furthermore, the computational time of a 2D model is much shorter as a result of the need for fewer gridpoints on the boundaries of the cell and nucleus, and thereby we use 2D model for further application and analysis in this work.
Application to cancer metastasis in 2D
There are preexisting openings (pores, fiberlike or channellike tracks) in ECM that enable cancer cells to migrate with an independence of MMP’s (Paul et al. 2017). In this section, we apply the model to the transmigration of cancer cells through pores and channels to migrate from one part to another part of the tissue without degrading ECM.
Simulation on penetration of a cell through a cavity
We initially consider a single cell penetrating through a cavity, which is formed by two circular obstacles, without secreting proteolytic enzymes and remodeling the ECM; i.e., the cell migration is assumed to utilize the amoeboid mode. The initial state is shown in Fig. 5 (topleft). The cell is attracted to an imaginary source (indicated by the blue asterisk) that releases a chemokine or a ECM stiffness signal. The migration of the cell is directed up the gradient of the chemokine, and it is limited by the presence of the two physical obstacles. Further, it can be seen that the cell is mechanically compressed as a result of its shrinkage due to its migration through the cavity and the nucleus deforms whenever the size of the pore is smaller than the size of the nucleus, see Fig. 5. As soon as the cell exits the constriction and is no longer mechanically compressed, the nucleus returns to its equilibrium circular shape. Once the source has been engulfed, the cell shape returns to its equilibrium circular shape. The model only incorporates temporary adherence to the obstacle, no permanent adherence. After disappearance of the source, only restoration of the cell shape is modeled.
Simulation on penetration of a cell through a tube channel
Cell deformation is normally studied in vitro by using microfluidic devices(Mak et al. 2013; Paul et al. 2016; MyrandLapierre et al. 2015). In the latter work, discordshaped red blood cells have been shown to be able to repeatedly deform when penetrating through microcapillaries with a diameter of 2.5 \(\upmu \mathrm{m}\) or even less. As mentioned in Sect. 2.3, there are abundant preexisting fiberlike and channellike tracks formed by the alignments of the collagen architecture in interstitial tissues and organs, which guide or inhibit cell migration (Wolf et al. 2009; Paul et al. 2017). In Fig. 6, a schematic representation of an endothelial cell wall with a channel of approximate 10 \(\upmu \mathrm{m}\) in width is depicted (Paul et al. 2017). This value is considered here to guarantee that the cell is able to penetrate through it in most cases.
Mechanical boundaries could regulate some biomedical processes and Mak et al. (2013) demonstrate that if the confined dimensional modulation of a microfluidic device has a mechanical barrier smaller than the cell nucleus, then metastatic breast adenocarcinoma cells likely deform in elongated morphological states and invade distinct sites. Here, taking mechanical boundaries into account, we use the trigonometric function (from equation (15)) to simulate the different roughnesses through changing the value of parameter \(\epsilon \) and \(\omega \). A highly rough boundary of the channel is defined if the perturbation [see Eq. (15)] has a high frequency or/and a big amplitude, which is determined by the surface of the endothelial cells. Whereas a lower frequency (also a lower amplitude) as we depict in Fig. 6 could show where each cell is located. The discrepancy between the endothelial cellular surfaces and the channel through which the cancer (or immune) cell migrates, could be a consequence of the extracellular matrix around the cells. Then the boundary of the channel can have various roughnesses, which combined with other parameters, are analyzed by using Monte Carlo simulations based on this model. Moreover, this model is incorporated with Poisseuille flow to simulate the micro blood flow referring our work (Chen et al. 2018). To investigate how the cell speed changes in the current scenario, the speed evolution with the respect of time is plotted in Fig. 8a without the perturbation of vector Wiener process. As we expected, the cell speed slows down when it starts to squeeze the opening and subsequently accelerates to move toward the emitting source. When the \(\tau \) equals approximate 0.37 hour, the instantaneous speed reaches a peak and drops to zero after the engulfment of the source and cell shape recovery. During the transmigration in the tube, the cell migrates with a speed vibrating up and down at 200 \(\upmu \mathrm{m/h}\), which is in the range 1–5 \(\upmu \mathrm{m/min}\) for the typical speed of amoeboid movement observed in vivo in the work (Pinner and Sahai 2008). Moreover, the cell speed can be controlled under various conditions, like the number of emitting sources, the diffusion coefficient, cell mobility.
Simulation on cancer metastasis
Immune cells and cancer cells similarly deform in chemically or mechanically induced locomotion. The work by Springer (1994) reports that leukocytes can attach to the wall of a blood vessel by binding to adhesion molecules of the endothelial cells, subsequently the leukocytes flatten themselves, and then squeeze through openings which are much smaller than themselves among the endothelial cells. Analogously, metastatic cells utilize similar mechanisms when intravasting into or extravasating out of blood vessels. Cancer metastasis is a multistep cascade that can be divided into the following steps, (1) escape from the primary tumor site; (2) survive transit in the bloodstream or lymphatic vessels after successful intravasation; (3) disseminate and extravasate subsequently; (4) start to proliferate and colonize secondary sites at distant organs (Chambers et al. 2002a; Kopfstein and Christofori 2006). We attempt to simulate the steps of intravasation and extravasation and several consecutive snapshots showing the shape changes of cell and nucleus are provided in Fig. 7, where a schematic diagram of a capillarysized channel is depicted. In order to get around hypoxia (or lack of nutrition) as a result of competitive growth in cancer cell colonies or as a response to a stiffness gradient, metastatic cells show migratory exploratory behavior toward regions outside the cology they reside in. This migration can be inspired by gelstiffnessdependent differences in traction forces or strain energies in Massalha and Weihs (2017). Therefore cancer cells are capable of penetrating through small openings in endothelium. This process is highly inefficient, and during this dissemination, the majority of cancer cells would die and only \(<\,0.02\%\) of them are able to seed at distant sites successfully (CeliàTerrassa and Kang 2016; Luzzi et al. 1998). Analogously, the cell speed evolution of this model is shown in Fig. 8b, and the speed is around 200 \(\upmu \mathrm{m/h}\) in the channel and reaches a peak instantaneously when the cell gets close to the source. The reason for this peak is the singularity in Eq. (2) at the position of the source, which gives a very large gradient of the concentration near the source. This peak could be regularised by either adding a timedependency (through an analytic solution or through a numerical solution of the concentration) or by replacing the chemotaxis by a factor such that the velocity stays bounded. All these approaches make the model more complication and since the objective was a construct a simple model, this has been omitted. At the final stages, the speed decreases to zero due to lack of attraction signals and the vector Wiener process. We also remark that the large variations in Fig. 8b are caused by the cell having to pass through the apertures and having to migrate along the wall of the channel. This interaction between the cell boundary and obstacle causes the switch between repulsion and migration along the tangent of the obstacle and attraction as a result of a component normal to the tangent of the boundary of the obstacle. This effect of this discontinuous switch mechanism can only be inhibited by choosing a smaller time step.
Parameter study with Monte Carlo simulations
If certain input values contain uncertainties, Monte Carlo simulations could be a way to evaluate the impacts of output. This method enables us to estimate of the impact from variables ranging from various statistical distributions like Pareto, uniform, normal, lognormal, Chisquare, exponential (Mooney 1997). Furthermore, Monte Carlo simulations have been used over a spectrum of systems, which is typically concluded in following four steps, (1) generate the input random values based on their probability distribution functions; (2) calculate samples; (3) repeat the abovementioned steps with a number of trials \(N_\mathrm{s}\); (4) calculate the mean and construct a relative frequency distribution of the simulated results (Mooney 1997; Mahadevan 1997). Furthermore, one can estimate the correlation between the various input and output parameters.
The model introduced in Sect. 4.3.2 is used in Monte Carlo simulations, with the channel boundary of 60 \(\upmu \mathrm{m}\) in length and approximately 10 \(\upmu \mathrm{m}\) in width. The transit time interval that starts once one of the cell’s boundary points enters the channel and lasts until the last point exits the channel is defined as the penetration time \(\tau \). In this section, the influences of several parameters on the penetration time \(\tau \) are investigated.
As we discussed in Sect. 3.4, the accuracy of the simulation result depends on the number of samples. To achieve an accurate approximation, the number of samples is tested that is shown in Fig. 9.
Note that the axes represent the logarithm of sample count and the mean of transit time, respectively. If the sample count in the Monte Carlo simulations is too small, then the average penetration time has not yet converged (see Fig. 9 for \(N_\mathrm{s} < 200\)). We observe that using 10,000 samples only gives very small fluctuations of the average penetration time (see Fig. 9). The result has converged sufficiently to approximate 0.356 h. However, to evaluate the uncertainty of input data quantitatively, 10,000 samples are chosen in our simulation which give acceptable computation times in the order of hour. Using equation (31), the Monte Carlo error is estimated by \(\Vert E_{\mathrm{mc}}\Vert = \Vert \hat{\tau }^{\Delta t}\hat{\tau }_{N_\mathrm{s}}^{\Delta t}\Vert \simeq \frac{S_n}{\sqrt{N_\mathrm{s}}}\).
Monte Carlo simulations on parameters D, \(\beta \), \(\alpha \), \(\alpha ^n\)
We start with the Monte Carlo simulations on four input parameters which are the diffusion coefficient of the chemokine D, cell point mobility \(\beta \), cell deformation relaxation \(\alpha \) and the nucleus deformation relaxation \(\alpha ^n\). We sample them from the normal distribution, then they can be generated by Eq. (28) with the default values in Table 3.
The mean value of each is the same as the value in Table 1 and corresponding standard deviation reflects the degree of dispersion among samples. The values have been chosen mathematically based on an extensive testing. In Fig. 10, we plot a histogram of 10,000 samples as well as a cumulative distribution function (CDF) of the estimated probability of penetration time \(\tau \). Thence, the xaxis denotes the consecutive variable penetration time \(\tau \) and the yaxis represents the frequency of occurrence or the probability \(P_n (t \le \tau )\) of the corresponding variable depending on the chart considered in Fig. 10.
Taking different roughnesses of the channel boundary into consideration, various values of \(\epsilon \) (\(\epsilon = 0, 0.5, 1.0, 1.5\,\upmu \mathrm{m}\)) and \(\omega \) (\(\omega = 0, 0.25, 0.5, 0.75\, \upmu \mathrm{m}^{1}\)) are set and compared in Fig.11. The \(\epsilon \) parameter manifests the magnitude of vertical fluctuation, i.e., the amplitude; while \(\omega \) determines the frequency of the fluctuations of boundary. A smooth boundary has a small \(\omega \) value, then one cell is able to move through it much faster than through a rough channel. In Fig.11a, we observe four cumulative distribution functions with different slopes \(f(\tau )\), which represents the probability density. Thus, any probability \(P_n\) of a time interval \([\tau \frac{\Delta t}{2}, \tau +\frac{\Delta t}{2}]\) occurring can be calculated by the formula \( P_n(\tau \frac{\Delta t}{2}\le \hat{\tau } \le \tau +\frac{\Delta t}{2}) \approx f(\tau ) \cdot \Delta t\). Conversely, with the same probability, taking \(P_n = 0.5\) for an example, we can get the information about the transit time of one cell with 50% probability in various conditions, where \(\tau _1(\epsilon = 0 \, \upmu \mathrm{m})< \tau _2(\epsilon = 0.5 \, \upmu \mathrm{m})< \tau _3(\epsilon = 1.0 \, \upmu \mathrm{m}) < \tau _4(\epsilon = 1.5 \, \upmu \mathrm{m})\). Analogously, Fig.11b showing four cumulative distribution functions of the penetration time \(\tau \) are compared under different conditions with varying \(\omega \). With 50 % probability, one cell takes penetration time \(\tau _1\) with a straight boundary \(\omega = 0 \, \upmu \mathrm{m}^{1}\), the penetration time rises to \(\tau _2\), \(\tau _3\) and \(\tau _4\) with the increase in roughnesses \(\omega = 0.25 \, \upmu \mathrm{m}^{1}\), \(\omega = 0.5 \, \upmu \mathrm{m}^{1}\) and \(\omega = 0.75 \, \upmu \mathrm{m}^{1}\), respectively. In conclusion, both the standard deviation in the arrival times and the mean arrival time increase with increasing values of \(\epsilon \) and \(\omega \). Subsequently, we fix the roughness parameter values to \(\epsilon = 1.0 \, \upmu \mathrm{m}\) and \(\omega = 0.5 \, \upmu \mathrm{m}^{1}\), then the impacts of four parameters \(D, \beta , \alpha , \alpha ^n\) on the penetration time \(\tau \) are investigated and the correlation analysis are shown in Fig. 12. Based on the results, there is some positive correlation between the penetration time \(\tau \) and both D and \(\alpha \) with correlation coefficient r equal to 0.6068 and 0.49772, respectively. Moreover, \(\beta \) has a negative linear correlation with \(\tau \), whereas, the nucleus deformation relaxation has no obvious correlation with penetration time in this situation. However, the nucleus is the stiffest cellular component, which inhibits the confined cell migration if the pore diameter in the ECM is below a critical threshold (Wolf et al. 2013; Davidson et al. 2014). Therefore, the correlation between the penetration time and nucleus stiffness is expected to be highly positive if the width of channel is smaller than a critical threshold.
Monte Carlo simulations on parameters \(\epsilon \) and \(\omega \)
We next analyze the other two parameters \(\epsilon \) and \(\omega \), which reflect the amplitude and frequency of the channel boundary. Suppose \(\epsilon \) and \(\omega \) are too large, i.e., \(\epsilon ,\omega >2\) in our simulations, the trigonometric functions probably would trap migrating cells due to sharp peaks or corners. Therefore, \(\epsilon \) and \(\omega \) are generated carefully with uniform normal distribution by the following equation,
This above equation guarantees the value of \(\epsilon \) and \(\omega \) are uniformly distributed and bounded by \((0.5, \ 1.5)\) and \((0, \ 0.6)\). Based on the 10,000 samples, Fig. 13a shows the corresponding histogram, which looks like a lognormal chart and fits from a qualitative point of view with the experimental results by Abuhattum and Weihs (2016), where the migration speeds of single preadipocytes without chemoattractants follow a lognormal distribution. A cumulative percentage of the number of occurrences regarding the cell penetration time \(\tau \) is plotted in Fig. 13b.
Analogously, scatter diagrams about \(\epsilon \) and \(\omega \) with penetration time \(\tau \) indicating their correlations are shown in Fig. 14. With the increase in roughness, one cell travels a longer time to penetrate the channel in most cases. Furthermore, the increment of \(\omega \) makes a contribution to the total travel time of one cell. This is also reflected by the correlations of \(r = 0.4310\) and \(r = 0.7100\) between the penetration time and \(\epsilon \) and \(\omega \), respectively.
Monte Carlo simulations on parameters D, \(\beta \), \(\alpha \), \(\alpha ^n\), \(\epsilon \) and \(\omega \)
To test the essential variables simultaneously, all six parameters D, \(\beta \), \(\alpha \), \(\alpha ^n\), \(\epsilon \) and \(\omega \) are analyzed by Monte Carlo simulations. The histogram of the penetration time \(\tau \) is shown in Fig. 15a which can be fitted to a lognormal distribution. Furthermore, a CDF result is shown based on a sample of 10,000 times simulations in Fig.15b.
To investigate the impacts of variables on output results and analyze the correlations of each variable with penetration time \(\tau \), a couple of scatter plots are shown in Fig. 16, respectively. Adding some control variables that are statistically distributed yields more uncertainty to the system. The increase in uncertainty generally decreases the correlation. Therefore, in current simulation of six parameters, the correlation of parameters D, \(\beta \), \(\alpha \), \(\epsilon \) and \(\omega \) with time \(\tau \) decrease slightly compared with the simulations with the variation of four parameters. The correlation between \(\tau \) and \(\alpha ^n\) is still negligible. Further, Fig. 16 shows that the roughness (\(\epsilon \) and \(\omega \)) dominantly influences the cell travel time.
Discussion and conclusions
In this work, we develop a cellbased model to describe the morphological evolution of the cell and nucleus in a phenomenological way. The cell cytoskeleton spanning between the nucleus and the cell membrane is simulated by 30 springs. As we expected, an immune cell or a single cancer cell can deform according to the specific obstacles or paths when it encounters a stiff obstacle in a 2D or 3D environment. Compared with some existing models, e.g., a model investigating the role of nucleus deformation in the cell deformation under different geometrical and fluid flow conditions (SerranoAlcalde et al. 2017) and a threedimensional model describing nucleus mechanics during cell migration and deformation (Giverso et al. 2018), one of the major advantages of our modeling is its efficiency regarding CPU time, which enables to carry out Monte Carlo simulations for evaluation of parameter sensitivity. A further merit of the current model is its simplicity. If one is able measure the velocity of points on the surface of the cell under the influence of (the gradient of) a generic (being a concentration or a stiffness for instance) signal, then the \(\beta \)parameter can be determined. If one further is able to measure the retraction speed on the boundaries of the cell and the nucleus once the signal has disappeared, then it can fit the \(\alpha \) parameters.
The uncertainties in the input values necessitate us to study the impact of uncertainty by carrying out Monte Carlo simulations. With 10,000 samples, the correlations of each variable D, \(\beta \), \(\alpha \), \(\alpha ^n\), \(\epsilon \) and \(\omega \) with cell penetration time \(\tau \) are analyzed. The results show that \(\alpha ^n\) has no significant correlation with the penetration time in current situation, where the reason probably is the low range of parametric values in our simulations. A larger range, with variations over a lognormal distribution could give a higher correlation. The use of very high values of \(\alpha ^n\) in the model when the cell is penetrating through an aperture needs more investigation. Moreover, SerranoAlcalde et al. (2017) state that a small cell nucleus does not play a crucial role in cell deformabilitybased experiments under a fluid flow. Therefore, the deformability of the nucleus could be impacted by the size of the nucleus, and thereby influence the penetration time. Whereas, other variables influence the cell penetration time \(\tau \) to varying degrees, where the correlation of roughness is the most significant.
To make the problem tractable, some assumptions are made based on the simplified biomedical phenomenon, which are: (1) the equilibrium morphology of the cell is circular in 2D and spherical in 3D, respectively; (2) the cell is not allowed to die, which means the cell cannot be removed, in any extreme narrow scenarios; (3) cell mobility is simulated by a source secreting a single cytokine evenly and continuously until it is consumed , which makes the model consist of a system of ordinary differential equations; (4) the obstacles are absolutely stiff such that they cannot deform and thereby we do not need to consider the degradation of substrate/ ECM on the obstacles. Further, the introduction of elastic obstacles also needs the inclusion of mechanical balance based on Newton’s law for the objects. Although this would be an interesting extension of the model, we omit this in the current paper since this extension enlarges the parameter space for the Monte Carlo simulations. In order to improve the model, the following aspects could be considered in future work.

Compared to a 2D model, a 3D model is more physiological, however, there is no significant qualitative difference in terms of expected numerical results. Moreover, taking the Monte Carlo simulations into account, the CPU time for simulating the 2D model is much more reasonable. However, a 3D model will still be an interesting research direction in the future.

Amoeboid and mesenchymal movement, as the two basic forms of cell locomotion, mutually transform and participate in the process of cell migration. The former is also called pseudopodia movement including lamellipodia and filopodia, which normally takes place close to the cell front as a result of cell polarization (Lauffenburger and Horwitz 1996; Lämmermann and Sixt 2009; Paul et al. 2017). Since the interconversion between the amoeboid model and the mesenchymal model due to the cytoskeleton rearrangement happens during cancer cell migration (Zhao et al. 2011), the filopodia that is an extension of active membrane of cell front and rear might be considered in future work.

In the current work, we define constant values for the cell deformation relaxation \(\alpha \) and cell mobility \(\beta \) everywhere, while they in general depend on chemokines. Therefore, to introduce surfaceresident chemical species, some surface partial differential equations can be incorporated such that it describes the evolution of the chemical signals over the membrane surface. This amounts to solving
$$\begin{aligned} \begin{aligned}&{\underline{a}}_t + \nabla _{\Gamma } \cdot (\mathbf{v} {\underline{a}})  D_a \Delta _{\Gamma } {\underline{a}} = {\underline{f}} ({\underline{a}}), \\&\mathbf{v} = \frac{\mathrm{d}}{\mathrm{d}t} {\underline{x}}(t), \qquad (t,\mathbf{x}(t)) \in \mathbb {R}^+ \times \Gamma (t). \end{aligned} \end{aligned}$$(36)This is an interesting and relevant research direction, which will be taken into consideration in future work.

A tumor is typically surrounded by a dense network of collagen fibers, which are normally utilized by motile cancer cells to guide their paths (Sahai 2007). Furthermore, mutated cancer cells are capable of remodeling the normal ECM around them, abnormal ECM or the density of fibers preferably reshapes aligned direction in a parallel arrangement, which forms an anisotropic medium and thereby has a significant impact on cell migration. If we formalize this directional dependence through the socalled orientation tensor \(\Psi \). Then we get the following revision on the response to the external signal of the migration equations:
$$\begin{aligned} \begin{aligned} \mathrm{d}{} \mathbf{x}_i(t) =\,&(\beta _0 \mathbf{I} + \beta _1 \Psi ) \nabla c(t,\mathbf{x_i}(t)) \mathrm{d}t + \alpha ( \mathbf{x}_i^n(t)\\&+ \hat{\mathbf{x}}_i  \mathbf{x}_i(t)) \mathrm{d}t +\eta \mathrm{d}{} \mathbf{W}(t), \\&i \in \{1,\ldots ,N\}, \end{aligned} \end{aligned}$$(37)where \(\beta _0\) and \(\beta _1\) are two constants and \(\Psi \) can be obtained by
$$\begin{aligned} \Psi (t,\mathbf{x}) = \begin{pmatrix} \Psi _{xx} &{} \Psi _{xy} \\ \Psi _{xy} &{} \Psi _{yy} \end{pmatrix}. \end{aligned}$$(38)For the formalism, one can refer to the work by Cumming et al. (2010) and a further application in the work (Chen et al. 2018).

We note that the relaxation parameter of the nucleus has little correlation with the transmigration time. This finding seems counterintuitive. According to the studies of SerranoAlcalde et al. (2017), the stiffness of the nucleus hardly plays a role in cell deformability experiments if the nucleus is relatively small. However for larger sizes, this deformability of the nucleus may become more important.
Over the past several decades, a significant progress has been made in medical technology and attempts have been made to investigate the complexity of cancer initiation and progression. For example, cell deformability has been shown to have certain correlations with disease states of cells and metastatic potentials (Guck et al. 2005; Mak and Erickson 2013). Nonetheless, the biological mechanisms of a multistep metastatic cancer still remain poorly understood (Lambert et al. 2017). To make a contribution, our group will continue to work on biological mathematical modeling to predict the behavior of cells in the microenvironment and aid the biological experiments for the further understanding of cancer and drug development.
References
Abuhattum S, Weihs D (2016) Asymmetry in traction forces produced by migrating preadipocytes is bounded to 33%. Med Eng Phy 38(9):834–838
Borau C, Polacheck WJ, Kamm RD, GarcíaAznar JM (2014) Probabilistic voxelfe model for single cell motility in 3D. In Silico Cell Tissue Sci 1(1):2
Brunner CA, Ehrlicher A, Kohlstrunk B, Knebel D, Käs JA, Goegler M (2006) Cell migration through small gaps. Eur Biophys J 35(8):713–719
Byrne H, Drasdo D (2009) Individualbased and continuum models of growing cell populations: a comparison. J Math Biol 58(4):657–687
Cao X, Moeendarbary E, Isermann P, Davidson PM, Wang X, Chen MB, Burkart AK, Lammerding J, Kamm RD, Shenoy VB (2016) A chemomechanical model for nuclear morphology and stresses during cell transendothelial migration. Biophys J 111(7):1541–1552
CeliàTerrassa T, Kang Y (2016) Distinctive properties of metastasisinitiating cells. Gene Dev 30(8):892–908
Chambers AF, Groom AC, MacDonald IC (2002a) Dissemination and growth of cancer cells in metastatic sites. Nat Rev Cancer 2:563
Chambers AF, Groom AC, MacDonald IC (2002b) Metastasis: dissemination and growth of cancer cells in metastatic sites. Nat Rev Cancer 2(8):563–572
Champion JA, Mitragotri S (2006) Role of target geometry in phagocytosis. Proc Natl Acad Sci USA 103(13):4930–4934
Chen J, Weihs D, Vermolen FJ (2018) Monte carlo uncertainty quantification in modelling cell deformation during cancer metastasis CMBBE2018 conference proceeding
Chen J, Weihs D, Vermolen FJ (2018) A model for cell migration in nonisotropic fibrin networks with an application to pancreatic tumor islets. Biomech Model Mechanobiol 17:367–386
Cumming BD, McElwain DLS, Upton Z (2010) A mathematical model of wound healing and subsequent scarring. J R Soc Interface 7(42):19–34
Davidson PM, Denais C, Bakshi MC, Lammerding J (2014) Nuclear deformability constitutes a ratelimiting step during cell migration in 3D environments. Cell Mol Bioeng 7(3):293–306
Denais CM, Gilbert RM, Isermann P, McGregor AL, Lindert M, Weigelin B, Davidson PM, Friedl P, Wolf K, Lammerding J (2016) Nuclear envelope rupture and repair during cancer cell migration. Science 352(6283):353–358
Drasdo D, Kree R, McCaskill JS (1995) Monte carlo approach to tissuecell populations. Physic Rev E 52(6):6635
Fidler IJ (2003) The pathogenesis of cancer metastasis: the’seed and soil’hypothesis revisited. Nat Rev Cancer 3(6):453–458
Friedl P, Wolf K, Lammerding J (2011) Nuclear mechanics during cell migration. Curr Opin Cell Biol 23(1):55–64
Galle J, Loeffler M, Drasdo D (2005) Modeling the effect of deregulated proliferation and apoptosis on the growth dynamics of epithelial cell populations in vitro. Biophys J 88(1):62–75
Giverso C, Arduino A, Preziosi L (2018) How nucleus mechanics and ECM microstructure influence the invasion of single cells and multicellular aggregates. Bull Math Biol 80:1017–1045
Guck J, Schinkinger S, Lincoln B, Wottawah F, Ebert S, Romeyke M, Lenz D, Erickson HM, Ananthakrishnan R, Mitchell D et al (2005) Optical deformability as an inherent cell marker for testing malignant transformation and metastatic competence. Biophys J 88(5):3689–3698
Gupta GP, Massagué J (2006) Cancer metastasis: building a framework. Cell 127(4):679–695
Jayaraman S, Joo NS, Reitz B, Wine JJ, Verkman AS (2001) Submucosal gland secretions in airways from cystic fibrosis patients have normal [Na+] and pH but elevated viscosity. Proc Natl Acad Sci 98(14):8119–8123
Kalebic T, Garbisa S, Glaser B, Liotta LA (1983) Basement membrane collagen: degradation by migrating endothelial cells. Science 221(4607):281–283
Kopfstein L, Christofori G (2006) Metastasis: cellautonomous mechanisms versus contributions by the tumor microenvironment. Cell Mol Life Sci 63(4):449–468
Lambert AW, Pattabiraman DR, Weinberg RA (2017) Emerging biological principles of metastasis. Cell 168(4):670–691
Lämmermann T, Sixt M (2009) Mechanical modes of amoeboid cell migration. Curr Opin Cell Biol 21(5):636–644
Lauffenburger DA, Horwitz AF (1996) Cell migration: a physically integrated molecular process. Cell 84(3):359–369
Luzzi KJ, MacDonald IC, Schmidt EE, Kerkvliet N, Morris VL, Chambers AF, Groom AC (1998) Multistep nature of metastatic inefficiency: dormancy of solitary cells after successful extravasation and limited survival of early micrometastases. Am J Pathol 153(3):865–873
Madzvamuse A, George UZ (2013) The moving grid finite element method applied to cell movement and deformation. Finite Elem Anal Des 74:76–92
Mahadevan S (1997) Monte carlo simulation. Mechanical EngineeringNew York and BaselMarcel Dekker, pp 123–146
Mak M, Erickson D (2013) A serial micropipette microfluidic device with applications to cancer cell repeated deformation studies. Integr Biol 5(11):1374–1384
Mak M, ReinhartKing CA, Erickson D (2013) Elucidating mechanical transition effects of invading cancer cells with a subnucleusscaled microfluidic serial dimensional modulation device. Lab Chip 13(3):340–348
Marth W, Voigt A (2014) Signaling networks and cell motility: a computational approach using a phase field description. J Math Biol 69(1):91–112
Massalha S, Weihs D (2017) Metastatic breast cancer cells adhere strongly on varying stiffness substrates, initially without adjusting their morphology. Biomech Model Mechanobiol 16(3):961–970
Mooney CZ (1997) Monte carlo simulation, vol 116. Sage Publications, Thousand Oaks
MoussaviBaygi R, Jamali Y, Karimi R, Mofrad MRK (2011) Brownian dynamics simulation of nucleocytoplasmic transport: a coarsegrained model for the functional state of the nuclear pore complex. PLoS Comput Biol 7(6):e1002049
MyrandLapierre ME, Deng X, Ang RR, Matthews K, Santoso AT, Ma H (2015) Multiplexed fluidic plunger mechanism for the measurement of red blood cell deformability. Lab Chip 15(1):159–167
Odell G, Oster G, Burnside B, Alberch P (1980) A mechanical model for epithelial morphogenesis. J Math Biol 9(3):291–295
Odenthal T, Smeets B, Van Liedekerke P, Tijskens E, Van Oosterwyck H, Ramon H (2013) Analysis of initial cell spreading using mechanistic contact formulations for a deformable cell model. PLoS Comput Biol 9(10):e1003267
Paňková K, Rösel D, Novotnỳ M, Brábek J (2010) The molecular mechanisms of transition between mesenchymal and amoeboid invasiveness in tumor cells. Cell Mol Life Sci 67(1):63–71
Paul CD, Hung WC, Wirtz D, Konstantopoulos K (2016) Engineered models of confined cell migration. Annu Rev Biomed Eng 18:159–180
Paul CD, Mistriotis P, Konstantopoulos K (2017) Cancer cell motility: lessons from migration in confined spaces. Nat Rev Cancer 17(2):131–140
Pinner S, Sahai E (2008) Imaging amoeboid cancer cell motility in vivo. J Microsc 231(3):441–445
Rey R, GarciaAznar JM (2013) A phenomenological approach to modelling collective cell movement in 2D. Biomech Model Mechanobiol 12(6):1089–1100
Sahai E (2007) Illuminating the metastatic process. Nat Rev Cancer 7(10):737–749
Sahai E, Marshall CJ (2003) Differing modes of tumour cell invasion have distinct requirements for rho/rock signalling and extracellular proteolysis. Nat Cell Biol 5(8):711–719
Savinell JM, Lee GM, Palsson BO (1989) On the orders of magnitude of epigenic dynamics and monoclonal antibody production. Bioproc Biosyst Eng 4(5):231–234
SerranoAlcalde F, GarcíaAznar JM, José GómezBenito M (2017) The role of nuclear mechanics in cell deformation under creeping flows. J Theor Biol 432:25–32
Seyfried TN, Huysentruyt LC (2013) On the origin of cancer metastasis. Crit Rev Oncog 18(1–2):43
Springer TA (1994) Traffic signals for lymphocyte recirculation and leukocyte emigration: the multistep paradigm. Cell 76(2):301–314
Steele JM (2012) Stochastic calculus and financial applications. Springer, Berlin
Tozluoğlu M, Tournier AL, Jenkins RP, Hooper S, Bates PA (2013) Matrix geometry determines optimal cancer cell migration strategy and modulates response to interventions. Nat Cell Biol 15(7):751
Van Liedekerke P, Palm MM, Jagiella N, Drasdo D (2015) Simulating tissue mechanics with agentbased models: concepts, perspectives and some novel results. Comput Part Mech 2(4):401–444
van Oers RFM, Rens EG, LaValley DJ, ReinhartKing CA, Merks RMH (2014) Mechanical cellmatrix feedback explains pairwise and collective endothelial cell behavior in vitro. PLoS Comput Biol 10(8):e1003774
Vermolen FJ, Gefen A (2012) A semistochastic cellbased formalism to model the dynamics of migration of cells in colonies. Biomech Model Mechanobiol 11(1):183–195
Vermolen FJ, Gefen A (2013) A phenomenological model for chemicomechanically induced cell shape changes during migration and cellcell contacts. Biomech Model Mechanobiol 12:301–323
Vermolen FJ, Mul MM, Gefen A (2014) Semistochastic celllevel computational modeling of the immune system response to bacterial infections and the effects of antibiotics. Biomech Model Mechanobiol 13(4):713–734
Wolf K, Mazo I, Leung H, Engelke K, Von Andrian UH, Deryugina EI, Strongin AY, Bröcker EB, Friedl P (2003) Compensation mechanism in tumor cell migration. J Cell Biol 160(2):267–277
Wolf K, Alexander S, Schacht V, Coussens LM, von Andrian UH, van Rheenen J, Deryugina E, Friedl P (2009) Collagenbased cell migration models in vitro and in vivo. In: Seminars in cell & developmental biology. Elsevier
Wolf K, Te Lindert M, Krause M, Alexander S, Te Riet J, Willis AL, Hoffman RM, Figdor CG, Weiss SJ, Friedl P (2013) Physical limits of cell migration: control by ecm space and nuclear deformation and tuning by proteolysis and traction force. J Cell Biol 201(7):1069–1084
Wyckoff JB, Pinner SE, Gschmeissner S, Condeelis JS, Sahai E (2006) Rockand myosindependent matrix deformation enables proteaseindependent tumorcell invasion in vivo. Curr Biol 16(15):1515–1523
Zhao P, Zhang W, Wang SJ, XiaoLing Y, Tang J, Huang W, Li Y, Cui HY, Guo YS, Tavernier J et al (2011) HAb18G/CD147 promotes cell motility by regulating annexin IIactivated RhoA and Rac1 signaling pathways in hepatocellular carcinoma cells. Hepatology 54(6):2012–2024
Author information
Affiliations
Corresponding author
Ethics declarations
Funding
This study is funded by China Scholarship Council.
Conflict of interest
The authors declare that they have no conflict of interest.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Electronic supplementary material
10237_2018_1036_MOESM1_ESM.avi
Video1 One cell migrating along a rigid obstacle in 2D simulation. The corresponding consecutive snapshots are shown in Fig.3 of manuscript. (video 2.91 MB)
10237_2018_1036_MOESM2_ESM.avi
Video2 One cell migrating along a rigid obstacle in 3D simulation. The corresponding consecutive snapshots are shown in Fig.4 of manuscript. (video 3.65 MB)
10237_2018_1036_MOESM3_ESM.avi
Video3 One cell penetrating a cavity made of two rigid obstacles in 2D simulation. The corresponding consecutive snapshots are shown in Fig.5 of manuscript. (video 2.52 MB)
10237_2018_1036_MOESM4_ESM.avi
Video4 One cell penetrating through an endothelial cell wall in 2D simulation. The corresponding consecutive snapshots are shown in Fig.6 of manuscript. (video 6.18 MB)
10237_2018_1036_MOESM5_ESM.avi
Video5 One cell intravasting and extravasating of a blood or lymphatic vessel in 2D simulation. The corresponding consecutive snapshots are shown in Fig.7 of manuscript. (video 7.55 MB)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
About this article
Cite this article
Chen, J., Weihs, D., Van Dijk, M. et al. A phenomenological model for cell and nucleus deformation during cancer metastasis. Biomech Model Mechanobiol 17, 1429–1450 (2018). https://doi.org/10.1007/s1023701810365
Received:
Accepted:
Published:
Issue Date:
Keywords
 Cell deformation
 Nucleus deformation
 Monte Carlo simulations
 Cancer metastasis
 Cellbased model