A computational study of amoeboid motility in 3D: the role of extracellular matrix geometry, cell deformability, and cell–matrix adhesion

Abstract

Amoeboid cells often migrate using pseudopods, which are membrane protrusions that grow, bifurcate, and retract dynamically, resulting in a net cell displacement. Many cells within the human body, such as immune cells, epithelial cells, and even metastatic cancer cells, can migrate using the amoeboid phenotype. Amoeboid motility is a complex and multiscale process, where cell deformation, biochemistry, and cytosolic and extracellular fluid motions are coupled. Furthermore, the extracellular matrix (ECM) provides a confined, complex, and heterogeneous environment for the cells to navigate through. Amoeboid cells can migrate without significantly remodeling the ECM using weak or no adhesion, instead utilizing their deformability and the microstructure of the ECM to gain enough traction. While a large volume of work exists on cell motility on 2D substrates, amoeboid motility is 3D in nature. Despite recent progress in modeling cellular motility in 3D, there is a lack of systematic evaluations of the role of ECM microstructure, cell deformability, and adhesion on 3D motility. To fill this knowledge gap, here we present a multiscale, multiphysics modeling study of amoeboid motility through 3D-idealized ECM. The model is a coupled fluid‒structure and coarse-grain biochemistry interaction model that accounts for large deformation of cells, pseudopod dynamics, cytoplasmic and extracellular fluid motion, stochastic dynamics of cell-ECM adhesion, and microstructural (pore-scale) geometric details of the ECM. The key finding of the study is that cell deformation and matrix porosity strongly influence amoeboid motility, while weak adhesion and microscale structural details of the ECM have secondary but subtle effects.

1 Introduction

Cellular locomotion is an essential characteristic inherent to all organisms. A variety of locomotion techniques can be found in nature. Of particular interest is the amoeboid mode of motility, which is commonly observed in cells such as Dictyostelium discoideum (Levine and Rappel 2013). Amoeboid motility is also prevalent within the human body, with examples including the extravasation of leukocytes through blood vessel walls and their subsequent migration through the surrounding tissue toward sites of inflammation (Gaylo et al. 2016), fibroblast migration during reconstruction of damaged tissue (Petrie and Yamada 2015), epithelial cell migration for wound healing (Scianna 2015), and positioning of cells during fetal development (Raz and Mahabaleshwar 2009). Metastasizing cancer cells have also been observed using the amoeboid phenotype as they migrate away from a tumor and through connective tissue (Wells et al. 2013; Paul et al. 2017).

Cellular motility, in general, relies on four key biophysical processes: pseudopod generation, cell-substrate adhesion, degradation and remodeling of the surrounding tissue, and pseudopod retraction (Lämmermann and Sixt 2009). Pseudopods are slender, outward-extending protrusions of the cell body which grow, bifurcate, and retract in a dynamic fashion (Van Haastert 2011). They are formed when signaling cascades within the cell result in the activation of regulating proteins such as WASP or ARP2/3, which polymerize cytosolic G-actin monomers into F-actin filaments (Bray 2000). These filaments are then cross-linked directly beneath the cell’s plasma membrane, thereby generating an outward protrusive force. Transmembrane adhesive proteins, such as integrins, form connections with corresponding surface receptors, resulting in discrete bonds between the cell and the substrate (Harunaga and Yamada 2011). In the rear of the cell, myosin-II motor proteins generate membrane contraction, simultaneously breaking existing bonds and pulling the cell forward. In addition, cells in vivo are subjected to a complex, 3D environment created by the surrounding tissue. This extracellular matrix (ECM) resembles a heterogeneous and porous structure composed of diverse protein fibers embedded in a gel-like polysaccharide fluid (Yamaguchi et al. 2005; Stevens and George 2005; Wolf et al. 2009). Cells often migrate by significantly degrading and remodeling the ECM through proteolytic activity (Paul et al. 2017).

The relative contributions of the aforementioned processes vary across different motility phenotypes, while the environments that cells are subjected to also play a major role. Many cells, in particular the mesenchymal variants, migrate using strong cell-ECM adhesions and significant ECM remodeling (Friedl and Wolf 2003). Strong cell-ECM adhesion is also critical for migration in 2D, as it allows the cell to pull on the substrate and generate the traction forces necessary for migration. In contrast, amoeboid cells migrate using weak or nearly no cell-ECM adhesion, and without degrading and remodeling the ECM. Instead, they often choose pre-existing passages and pores in the ECM (Sahai 2005; Lämmermann et al. 2008; Wirtz et al. 2011; Joyce and Pollard 2009). Many cells also possess an adaptive quality in their migration behavior. For instance, when blocked from degrading matrix tissue by protease inhibitors, mesenchymal cancer cells were shown to convert to the amoeboid mode, effectively skirting the therapeutic treatment and retaining their invasive properties (Wolf et al. 2003).

A large volume of work conducted over several decades exists on cellular motility over 2D substrates (see Zaman 2013). In contrast, amoeboid motility is fundamentally 3D in nature. One potential difference between 2D versus 3D motility arises in terms of the roles of matrix stiffness and confinement. Matrix stiffness is well known to affect cell motility over 2D substrates, since a stiffer substrate allows for stronger adhesion (Brábek et al. 2010). In contrast, confinement in a 3D matrix allows sufficient physical contacts between the cell and the ECM. Physical contacts in 3D, though, depend on the geometry and microstructure of the ECM. The ECM microstructure, as characterized by its scaffold architecture (amorphous versus regular lattice, matrix porosity, fiber alignment etc.), also varies from organ to organ (Stevens and George 2005; Shieh 2011). How the ECM microstructure affects cell motility in 3D is not well understood. ECM microstructure is also important from a tissue engineering perspective. The internal architecture of engineered scaffolds varies greatly depending on their bio-application. Although it has been recognized that variation in microstructures of engineered scaffolds greatly affects cell proliferation (Giannitelli et al. 2014), their role on 3D cell motility has not been systematically addressed. Our first objective in this study is to address how cellular motility depends on different architectures of the ECM scaffold.

Cell deformation is essential for amoeboid motility, as it allows for the formation of pseudopods. In the absence of strong adhesion, amoeboid cells must utilize their deformability to squeeze through narrow pores within the matrix. In order to navigate through these complex ECM pores and passages, cells must dynamically change their shapes and alter their pseudopod dynamics, thereby adapting to the ECM microstructure. Cell deformability has been proposed as a biomarker for metastatic potential of cancer cells (Xu et al. 2012), and metastatic cells were found to be significantly softer than benign cells (Cross et al. 2007). Cell softening has been shown to result in higher migration speeds for fibroblasts and breast cancer cells, also correlating with the increased proliferation of malignant cells (Fritsch et al. 2010). When different ECM microstructures are considered, it is not known whether cell deformation differently influences cellular motility under varying ECM geometries. Moreover, it is unknown how ECM microstructures affect pseudopod dynamics. Therefore, our second objective is to study the influence of cell deformation on amoeboid motility under varying ECM geometry.

Cell-ECM adhesion itself is an important biophysical component of cell motility, although amoeboid cells can migrate with weak or no adhesion present. Amoeboid-type cancer cells migrate using very little adhesive capability (Liu et al. 2015; Paluch et al. 2016), while neutrophils lacking specific integrins showed no significant differences migrating in 3D when compared to wild-type leukocytes (Mandeville et al. 1997). Integrins can impede invasion, as seen in melanoma, breast, and colon cancer cells (Friedl and Wolf 2003). There are also cases of mutant cells with a reduced ability to form adhesions (Niewöhner et al. 1997), thereby affecting cell migration. A systematic study of the role of adhesion in amoeboid motility is lacking. In mesenchymal cells, focal adhesions form at the cell’s leading edge and break at the cell rear; thus, there is a gradient in adhesion which may help pull the cell forward (Morley et al. 2014). For the amoeboid cell, which has been described as having a diffuse and low expression of adhesion (Wolf et al. 2003), the mechanism is not as clear. Thus, our third objective is to quantify the role of weak adhesion in the 3D motility of amoeboid cells.

To that end, here we undertake a computational modeling study of pseudopod-assisted amoeboid motility of deformable cells through 3D, microstructure-resolved ECM. The model considered here is a multiscale, multiphysics, coupled fluid–structure-biochemistry interaction model to simulate motility of individual amoeboid cells. It combines (1) cell deformation, (2) pseudopod dynamics, (3) cytoplasmic and extracellular fluid motion, (4) 3D ECM microstructure (geometric) details, and (5) stochastic adhesion dynamics mediated by receptor-ligand interactions. In this approach, nano-scale biomolecular reactions are coarse-grained using a dynamic pattern formation model, but the processes at larger scales, such as cell deformation, are resolved with high fidelity.

Several prior works have studied cell locomotion in 2D geometries. For example, Bottino and Fauci used the immersed-boundary method, modeling the actin cytoskeleton as a linear spring network capable of protrusion and contraction, while dynamic adhesion between the cell and substrate was considered using elastic links (Bottino and Fauci 1998). Vanderlei et al. also used the immersed-boundary method for deformable cells and coupled it with an internalized polarization-controlling reaction–diffusion system (Vanderlei et al. 2011). Shao et al. used a phase-field model with discrete adhesion sites characterized by gripping or slipping modes (Shao et al. 2012). Variation of protrusive, contractile, and adhesive forces suggested a collective effect on cell migration. Sakamoto et al. used the finite-element method for an axisymmetric cell (Sakamoto et al. 2014). Amoeboid cells were found to have weaker adhesion capabilities, and more frequent transitions between elongation and retraction stages. Cirit et al. studied the interplay between adhesion and protrusion at the leading edge of a cell (Cirit et al. 2010). The densities of nascent and stable adhesions, and myosin were modeled using coupled ODEs. It was found that optimal leading-edge protrusion occurred at intermediate ECM densities. Lim et al. considered adhesion-free, bleb-based migration in a microchannel, modeling the plasma membrane and elastic cortex as two distinct elastic networks connected through Hookean adhesive springs (Lim et al. 2013). Migration speed was shown to reach maximum magnitude at an intermediate confinement. Hecht et al. simulated a chemotactic 2D amoeboid cell lacking adhesion, with protrusion driven by reaction–diffusion equations, while subjected to unbounded medium as well as obstacle and maze-like geometries (Hecht et al. 2011). The importance of geometry on migration was noted, since cells solely dependent on a chemotactic signal for guidance became trapped. Copos et al. studied spatiotemporal adhesion patterns of a 2D cell on a flat substrate, using a force-balance approach and linear springs to model adhesion (Capos et al. 2017). Amoeboid cells were shown to have extension and contraction events from nonuniformly spaced adhesions. Schlüter et al. examined the dynamics of a rigid adhesive cell migrating on a 2D substrate composed of movable cylindrical fibers, while considering matrix stiffness and orientation (Schlüter et al. 2012). Results showed cells preferred stiffer matrices over softer ones, and cell persistence increased with fiber orientation.

In comparison, 3D models of amoeboid motility in ECM geometries are less in number. Najem and Grant used a phase-field model of a neutrophil subjected to a chemotactic signal but excluding fluid interaction (Najem and Grant 2013). Elliott et al. developed a finite-element model for a chemotactic amoeboid cell using a reaction–diffusion system for pseudopod dynamics, in addition to a 2D model of cell motility in a moveable obstacle field (Elliott et al. 2012). Moure and Gomez used a diffuse-domain approach, considering actin and myosin dynamics within the cytosol, as well as activator dynamics on the cell membrane using a reaction–diffusion equation (Moure and Gomez 2017; Moure and Gomez 2018). Despite such recent progress, there is still a lack of studies on systematic evaluations of the role of varying ECM microstructures, cell deformability and cell-ECM adhesion in 3D motility of amoeboid cells. The present study seeks to fill in this knowledge gap.

2 Model

The present model is a multiscale, multiphysics model coupling fluid mechanics, solid mechanics, coarse-grain biochemistry, receptor-mediated adhesion dynamics, and microstructure-resolved ECM. A finite-element method is used to compute cell deformation, while a finite-volume/spectral method is used for intra- and extracellular fluid motion. Protein biochemistry is modeled using a reaction–diffusion system and solved using the finite-element method. A stochastic Monte-Carlo method is used to predict adhesive bond formation between the cell and substrate. The ECM microstructures are modeled as embedded objects using an immersed-boundary method (IBM). The IBM is also used to couple the cell membrane with the surrounding fluid. A schematic of different model components and the numerical mesh is shown in Fig. 1. Descriptions about the different components are given below. To be noted, there are no external cues or gradients used to drive the cell’s migration.

Fig. 1

Numerical model for 3D amoeboid cell motility through ECM. a, b Finite-element mesh on the deforming cell surface. c Surface mesh on the ECM microstructures. d Eulerian mesh discretizing the entire computational domain that is comprised of the extracellular space (extracellular fluid and ECM scaffold) and the cytoplasm. The “immersed” cell membrane separates the cytoplasm and extracellular space. e A 2D slice showing a sample instantaneous flow field in and around the cell. Color contours represent fluid velocity magnitudes and arrows represent velocity vectors. Pink lines represent adhesive bonds between cell and ECM surface

2.1 Coarse-graining of biochemistry

Pseudopod protrusions are generated by protein reactions which are nano-scale processes (Lämmermann and Sixt 2009; Alberts et al. 2007; Lenz 2008). In order to simulate cell movement over a distance of several cell lengths, we coarse-grain the nano-scale protein dynamics using a system of reaction–diffusion (RD) equations (Meinhardt 1999; Maini et al. 1997). However, the model must retain the essential dynamics of pseudopods observed experimentally: pseudopods continually generate and grow, bifurcate into daughter pseudopods, meander over the cell surface, and finally retract (Hecht et al. 2011; Levine et al. 2006; Xiong et al. 2010). Furthermore, the growth of actin filaments is a locally positive-feedback process due to signaling molecules, but is also restricted by actin capping proteins. Hence, the process is described as a short-range self-enhancement and long-range inhibition. Such processes can be modeled as a Turing instability (Hecht et al. 2011; Levine et al. 2006; Hecht et al. 2010; Meinhardt 1982). In this model, intracellular proteins are grouped into two categories: activators, which represent proteins like Arp2/3, WASP, actin filaments, etc., and inhibitors, which represent capping and severing proteins. A set of nonlinear RD equations are then written in terms of activator (a₁) and inhibitor (a₂ and a₃) concentrations on the cell surface as given by Eqs. (1–3), where the dot represents the time derivative, ∆_s is the surface Laplacian, D₁ and D₃ are activator and inhibitor diffusivities, and r₁ and r₃ are reaction rates (Elliott et al. 2012). A random noise ε(x,t) is used to perturb the system into producing instabilities.

$$\dot{a}_{1} = D_{1} \Delta_{S} a_{1} + \frac{{(r_{1} + \varepsilon )\left( {\frac{{a_{1}^{2} }}{{a_{2} }} + k_{1} } \right)}}{{(S_{3} + a_{3} )\left( {1 + s_{1} a_{1}^{2} } \right)}} - r_{1} a_{1}$$

(1)

$$a_{2} = \frac{1}{S}\int\limits_{S} {a_{1} {\text{d}}S}$$

(2)

$$\dot{a}_{3} = D_{3} \Delta_{S} a_{3} - r_{3} a_{3} + k_{2} a_{1}$$

(3)

The above RD system is solved on the deforming cell surface using a finite-element method. Figure 1 shows the FEM mesh. Each of the activator and inhibitor concentrations is initiated to unity across the cell surface. Once the reaction–diffusion system is started, and through the influence of the random noise ε, a Turing instability soon develops. Depending on the choice of the parameters, the solution yields various patterns of activator concentration representing the Turing instability (Maini et al. 1997). Numerical experiments were performed to determine the specific range of parameters for which bifurcating patterns are observed (Fig. 2, Table 1). The Turing instabilities lead to growth, bifurcation, and decay of such localized activator concentration, making it an effective means of representing the coarse-grain dynamics of actin filament growth. The activator concentration is then used to model the pseudopod-generating protrusive force as f_P = ξa₁n where ξ is the force per actin filament, and n is the unit vector normal to the cell surface. When a deformable cell is considered, this model recreates the experimentally observed bifurcating pseudopod dynamics.

Fig. 2

Example of a Turing pattern used to generate pseudopods. Shown here is the activator concentration (scaled by initial concentration) over time. The activator concentration grows in a small area, meanders on the cell surface, and then bifurcates into two daughter patches, one of which subsequently dies, and the cycle repeats

Table 1 Physical parameters used in the model

2.2 Modeling a deformable cell

Amoeboid cells utilize their deformability to migrate, and often assume extremely deformed shapes. Cells are viscoelastic; in order to conveniently resolve large viscoelastic deformation, we model the cell as a viscous drop surrounded by a hyperelastic, zero-thickness membrane. The viscous fluid represents the cytoplasm, while the composite structure of the lipid bilayer and cortex is combined into a single membrane. The membrane cortex primarily exerts a resistance against shearing deformation while the lipid bilayer acts against bending and cell surface area dilation. For the shearing deformation and area dilation, the following strain energy function is used (Skalak et al. 1973),

$$W_{E} = \frac{{G_{S} }}{4}\left[ {\left( {I_{1}^{2} + 2I_{1} - 2I_{2} } \right) + {\text{CI}}_{2}^{2} } \right]$$

(4)

Here G_S is the membrane shear elastic modulus, CG_S is the area dilation modulus, $I_{1} = \varepsilon_{1}^{2} + \varepsilon_{2}^{2} - 2$ and $I_{2} = \varepsilon_{1}^{2} + \varepsilon_{2}^{2} - 1$ are the strain invariants, and ε₁ and ε₂ are the principal stretch ratios. The stresses generated in the membrane can be evaluated from the above constitutive law as

$$\tau_{1} = \frac{1}{{\varepsilon_{2} }}\frac{{\partial W_{E} }}{{\partial \varepsilon_{1} }}\;{\text{and}}\;\tau_{2} = \frac{1}{{\varepsilon_{1} }}\frac{{\partial W_{E} }}{{\partial \varepsilon_{2} }}.$$

(5)

The bending resistance is modeled following Helfrich’s formulation for bending energy,

$$W_{B} = \frac{{E_{B} }}{2}\int\limits_{S} {(2k - c_{0} )^{2} } {\text{d}}s,$$

(6)

where E_B is the bending modulus, k is the mean curvature, c₀ is the spontaneous curvature, and S is the surface area (Zhongcan and Helfrich 1989). FEM is used to compute the membrane stresses. The details of the method are given in our previous works (Yazdani and Bagchi 2013; Yazdani and Bagchi 2012; Cordasco et al. 2014).

2.3 Cell–matrix adhesion

Cell–matrix adhesion is assumed to occur by formation of adhesive bonds resulting from receptor-ligand-type interactions. Amoeboid adhesion has been described as diffusely distributed across the cell membrane (Friedl and Wolf 2003; Titus and Goodson 2017). To address this requirement in our model, adhesion is restricted to about 10% of the cell surface. Numerically, bonds are allowed to form on only 10% of the cell surface nodes. We assume that microvilli are located at these nodes, and bonds form at the microvilli tips. Microvilli are assumed to be unstretchable, but the bonds behave as stretched springs under force loading following a Hookean model (Hammer and Apte 1992)

$$f_{b} = k_{b} \left( {l - l_{0} } \right)$$

(7)

where f_b is the force in each bond, k_b is the spring constant, and l and l₀ are the stretched and unstretched bond lengths. Reaction kinetics dictating the rates of bond formation and breakage are given as (Hammer and Apte 1992)

$$K_{f} = K_{f}^{0}$$

(8)

$$K_{r} = k_{r}^{0} \exp \left[ {\frac{{ - [\sigma_{t} \left( {l - l_{0} } \right)]^{2} }}{{K_{B} T}}} \right]$$

(9)

where K_f and K_r are the rates of bond formation and breakage, K ⁰_f and K ⁰_r are the rates corresponding to the unstretched bond, σ_t is a parameter that relates reactivity of the molecules to the intermolecular distance, K_B is the Boltzmann constant and T is the absolute temperature. Each microvillus location is checked for potential bond formation if it is less than 1 µm from an ECM surface. Whether a bond forms depends upon a Monte Carlo model where the probability of the bond formation at an unbounded node is given by (Hammer and Apte 1992)

$$P_{f} = 1 - \exp \left( { - K_{f} \Delta t} \right)$$

(10)

where ∆t is the time interval. If $P_{f}$ is greater than a random number $N_{f}$ (0 < N_f < 1), a bond will form. Alternatively, if a bond exists, its probability of breakage $P_{r}$ is determined as

$$P_{r} = 1 - \exp \left( { - K_{r} \Delta t} \right)$$

(11)

If $P_{r}$ exceeds a random number N_r (0 < N_r < 1), the bond will break. This method imparts a stochastic nature to the cell–matrix adhesion dynamics.

2.4 Cytoplasmic and extracellular fluid flow

The cytoplasmic and extracellular fluids are assumed to be incompressible and Newtonian. Because inertia is negligible, fluid motion is governed by the unsteady Stokes equation and the incompressibility condition.

$$\rho \frac{{\partial \varvec{u}}}{\partial t} = - \nabla p + \nabla \cdot \mu \left[ {\nabla \varvec{u + }\nabla \varvec{u}^{T} } \right]$$

(12)

$$\nabla \cdot \varvec{u} = 0$$

(13)

where $\varvec{u}$ and $p$ represent the fluid velocity and pressure fields. The densities and viscosities of the cytoplasmic and extracellular fluids are assumed to be identical and equal to those of water. The model is not constrained by this fact, however, as the influence of viscosity difference was considered in our previous work (Campbell and Bagchi 2017). The extracellular fluid is further assumed to be otherwise stagnant, until cell movement causes fluid motion. The cell membrane is ‘‘immersed’’ within the surrounding fluid. Protrusive force, adhesive force, membrane deformation, and fluid motion in the cytoplasm and extracellular fluid are coupled together using the continuous forcing immersed-boundary method (IBM) (Peskin 2002). In this method, the cell interior and exterior rely on a single set of governing equations, with each respective zone identified by an indicator function. The cell membrane, which serves as the interface between the two fluids, is included through the use of a source term in the Stokes equation, taking the forces acting on the cell membrane into account as

$$\rho \frac{{\partial \varvec{u}}}{\partial t} = - \nabla p + \nabla \cdot \mu \left[ {\nabla \varvec{u} + \nabla \varvec{u}^{T} } \right] + \mathop \smallint \limits_{S}^{ } \left( {\varvec{f}_{E} + \varvec{f}_{B} + \varvec{f}_{P} + \varvec{f}_{A} } \right)\delta {\text{d}}S$$

(14)

where $\delta$ is the 3D Dirac delta function, S represents the cell surface, $\varvec{f}_{E}$ and $\varvec{f}_{B}$ are the membrane forces due to shearing deformation and area dilation, and bending, respectively, $\varvec{f}_{P}$ is the protrusive force, and $\varvec{f}_{A}$ is the adhesive force. The delta function is only nonzero at the location of the membrane. Once the fluid velocity is solved via Eq (14), the membrane velocity $\varvec{u}_{m}$ can be computed by interpolating the fluid velocity from the surrounding Eulerian nodes utilizing the Delta function as

$$\varvec{u}_{m} \left( {\varvec{x^{\prime}}} \right) = \mathop \smallint \limits_{S}^{ } \varvec{u}\left( \varvec{x} \right)\delta \left( {\varvec{x} - \varvec{x^{\prime}}} \right){\text{d}}\varvec{x}$$

(15)

The membrane points $\varvec{x}_{m}$ are then advected as ${\text{d}}\varvec{x}_{m} /{\text{d}}t = \varvec{u}_{m}$ to new positions, resulting in an updated cell configuration. Because the above method directly couples fluid motion with membrane deformation, the fluid drag on cell is directly and accurately resolved.

A projection-based method for the unsteady Stokes equations is used to obtain velocity and pressure fields. The semi-implicit Crank-Nicolson method is used for time discretization, while spatial derivatives are treated using the finite-volume method in nonperiodic directions, and Fourier series in periodic direction.

2.5 Microstructure-resolved ECM

The ECM is modeled as a 3D, fluid-filled porous medium with resolved pore-scale (microstructure) details. Four different types of ECM scaffolds are considered as shown in Fig. 3: (1) a network of randomly oriented slender rods, (2) a 3D lattice made of orthogonal slender rods, (3) an array of aligned slender rods, and (4) an array of spheres. While none of these geometries represents the exact microstructure of ECM, some physiological justifications for their selection can be given. A diverse assortment of ECM geometries can be found within the body as well as in artificial tissue scaffolds. The microstructure of the ECM is composed of fibers (collagen, elastin, etc.) and other cells (fibroblasts, leukocytes, mast cells, etc.). Although an individual collagen fiber is only a few nanometers in diameter, they often exist in bundles, thereby increasing their cross-sectional dimensions by several orders. Connective tissue can then be described as a loose, 3D heterogeneous, isotropic network characterized by cross-linked collagen bundles and fibers, which bears similarity to the randomly oriented rods geometry (Wolf et al. 2009; Charras and Sahai 2014; Wolf and Friedl 2011). Random fibers represented by slender rods were considered in 2D modeling by Schlüter et al. (Schlüter et al. 2012), as well as a highly porous, periodic, 3D variant in Moure and Gomez (2017, 2018). The 3D lattice, constructed of orthogonal rods, is used to represent engineered scaffolds (Giannitelli et al. 2014; Kim et al. 2013; Pedersen et al. 2010). Aligned filaments can be found in several parts of the body, such as in collagen networks (Wolf et al. 2009) or between muscle or nerve fibers and blood vessels (Paul et al. 2017). The dermis features dense and tightly aligned collagen bundles which can also be correlated to the aligned rod geometry considered here (Ploetz et al. 1991). Breast tumors were noted to contain bundled collagen fibers which protrude outward in a radial fashion at the tumor-stroma interface (Paul et al. 2017). Basement membrane features a dense, 2D cross-linked sheet-like geometry (Wolf and Friedl 2011). The sphere array can represent the bone scaffold which resembles sponge-like porous structures (Thavornyutikarn et al. 2014).

Fig. 3

ECM geometries considered. a An array of randomly oriented slender rods; b 3D lattice of orthogonal slender rods; c array of aligned rods; d array of spheres

The rods and spheres which compose the scaffolds are assumed to be rigid, and the space in between them is filled with the extracellular fluid. The amoeboid cell moves through this fluid-filled space while adhering to the nearby rod or sphere surfaces. The cell movement generates flow in the surrounding fluid, which is solved accurately as described in §2.4. The fluid must satisfy the no-slip condition on the surface of the ECM microstructures. For this, we use the sharp-interface ghost-node immersed-boundary method (GNIBM) (Balogh and Bagchi 2017). The basic premise of the method is to enforce the no-slip condition along a rigid surface that does not coincide with the rectangular Eulerian mesh used to solve the fluid motion by enforcing a constraint at an adjacent location outside the fluid domain. The methodology can be applied to accurately solve the fluid motion in any arbitrary complex geometry, and the details can be found in Balogh and Bagchi (2017).

2.6 Computational domain and resolution

The computational domain is a cube with each side about 19 times the cell radius. The domain is discretized using a rectangular mesh having 360³ points (Fig. 1). The cell surface is discretized by 20,480 FEM triangles. The ECM surface is discretized using triangular mesh that consists of ~ (1.6–4.2) × 10⁶ triangles, depending on the geometry considered. A dimensionless time step of 0.001 is used.

2.7 Model parameters

Physical parameters are listed in Table 1. The initial undeformed shape of the cell is taken to be spherical. The simulations are performed using dimensionless variables, which are also used to present the results. Lengths are scaled by the cell radius R, time is scaled by R²/D₁, and velocities are scaled by D₁/R, where D₁ is the activator diffusivity. An important dimensionless parameter is the ratio of the protrusive force to the cell membrane elastic force $\alpha = \xi /{\text{RG}}_{S}$.

Membrane stiffness varies in cells. For instance, immune cells are relatively softer than fibroblasts (Petrie and Yamada 2015). Stiffness also varies in malignant and drug-treated cells (Xu et al. 2012; Cross et al. 2007; Fritsch et al. 2010; Wakatsuki et al. 2000). We consider the range α = 1–7, with the lower values corresponding to a relatively stiff cell and the higher values representing a very deformable cell. The parameter that defines cell-ECM adhesion is the ratio of bond constant to membrane elastic force, σ = k_b l₀/RG_S. Since amoeboid motility occurs with weak or no adhesion, we consider a range σ = 0–3. The parameter defining the extracellular space is the porosity (void fraction) ϕ, defined as the ratio of the extracellular fluid volume to total volume. The effect of porosity is studied for a spherical array matrix by varying ϕ from 0.54 to 0.92. Otherwise, ϕ is kept constant at 0.68 for the orthogonal lattice, aligned rods, and array of spheres, while it is 0.88 for the random rods network. This range corresponds to physiological conditions in addition to those encountered within tissue engineering (Kasper et al. 2012).

Extensive validation of the cell deformation model has been done in our previous works (Yazdani and Bagchi 2013; Yazdani and Bagchi 2012; Balogh and Bagchi 2017). This includes cell aspiration in a micropipette, and cell deformation in externally applied flow. Validation of the reaction–diffusion model was presented in Campbell and Bagchi (2017).

3 Results

3.1 General motility behavior

Figure 4 presents a time sequence of an amoeboid cell migrating through the orthogonal lattice geometry as predicted by our simulations (see Supplementary Material for a movie). It should be noted that no chemotactic gradient is present in the model. A localized region of high activator concentration first develops on the initially undeformed cell. This results in a local protrusive force on the cell membrane, generating a de novo pseudopod. The de novo pseudopod extends outward until the underlying Turing instability causes the activator patch to bifurcate into two separate patches, resulting in two daughter pseudopods. Soon after, a reduction in activator concentration leads to the retraction of one pseudopod. The remaining pseudopod bifurcates again, and the cycle repeats. Using the sequence of pseudopod bifurcation, growth, and retraction, the cell squeezes through the fluid-filled spaces in the ECM as the resultant protrusive force pulls it forward. Significant deformation of the highly dynamic cell is evident. As the pseudopods meander and bifurcate over the cell membrane, a pseudo-random motion ensues.

Fig. 4

Simulation results: a sequence of images of a migrating amoeba (in color) through the 3D orthogonal lattice. The extracellular solid phase is shown in gray. Color contours on cell surface represent activator concentration (scaled by initial concentration). Arrows indicate direction of motion. b Trajectories of a few cells are shown. c Instantaneous cell speed is shown for two representative cases

Figures 4b, c show 3D trajectories of a few cells and their time-dependent speed, respectively. Random-like motility is evident from these plots as well. Interestingly, as is also evident from the sequence plot, loops and nearly 180^o turns can be seen in trajectories which are the result of a cell circumventing the scaffold structure. Deviations in trajectory can also occur due to obstruction by the ECM and adhesion. Large ‘spikes’ are observed in the velocity–time plot, suggesting that cell movement is impulsive rather than continuous. They occur when pseudopods are generated, which are usually short-lived with a lifetime of less than ~ 0.5 dimensionless time. Smaller spikes appear when a nascent pseudopod directly hits an obstacle surface and is terminated before fully extending.

Migration behavior can be quantified through several means, such as the mean-squared displacement (MSD), and velocity auto-correlation function (ACF). For a random walk (Brownian) motion, the MSD would follow a ballistic behavior (~ t²) at short time scales and a diffusive behavior (~ t) at longer time scales (Dieterich et al. 2008). In contrast, the MSDs predicted by our simulations as shown in Fig. 5 reveal a directional, super-diffusive behavior with an exponent of ~ 1.8 at intermediate time scales. This result is in agreement with recent experimental studies (Dieterich et al. 2008; Wu et al. 2014) and suggests that 3D cell migration does not follow a simple random walk behavior, even in the absence of external cues or signals. The velocity auto-correlation can additionally be utilized to further examine cell dynamics. For a random walk, the ACF would record as nearly zero, while a persistent random walk, which is characterized by a succession of uncorrelated movements of duration equal to a persistence time, would exhibit an exponential decay (Wu et al. 2015). Yet as seen in Fig. 5, the ACF predicted by our simulations shows a slowly decaying trend before decreasing significantly at larger times. This is also the case observed experimentally in Dieterich et al. (2008), where an exponential decay (i.e., the persistent random walk) fails to adequately capture the data trends, indicating that 3D cell migration cannot be characterized as such. The persistent and correlated cell trajectories are due to focused and oriented pseudopod generation as will be discussed later.

Fig. 5

a Mean-squared displacement (MSD) and b velocity auto-correlation function (ACF) predicted from several simulations

The influence of cell deformability on its trajectory is shown in Fig. 6a. Membrane deformability is seen to significantly alter cell behavior. Smaller values of α, which correspond to increased membrane stiffness, result in a cell incapable of significant migration through the ECM, with smaller distances traveled and more frequent turning events. As deformability increases, however, trajectories begin to show more directed paths taken by the cell for longer periods of time and with less turning events (Fig. 6). Such migration behavior can be directly correlated to pseudopod dynamics. For relatively larger values of α, cells assume polarized shapes, and pseudopod generation is localized primarily near the ‘leading edge’ (Fig. 6d). Pseudopod focusing allows the cell to maintain a constant direction for a long time, and the slender shape allows it to navigate through smaller ECM pores. For smaller α, protrusive forces are not strong enough to overcome membrane tension, resulting in nearly rounded cell shapes (Fig. 6c), making them incapable of significant migration through the ECM. In addition, pseudopods are generated all over the cell surface, resulting in the loss of directionality for a cell with smaller values of α.

Fig. 6

Cell trajectories under varying deformability (a α = 1, black; 3, red; 5, green; 7, blue), and varying ECM porosity (b φ = 0.54, black; 0.68, red; 0.79, green; 0.87, blue; 0.92, magenta). Pseudopod dynamics for varying deformability (c α = 1; D: α = 7; ECM is removed from the image). e Obstacle-mediated trajectory reversal; arrows indicate direction of motion

For amoeboid motility, a loss of directionality and matrix penetration in cell migration with increasing ECM confinement is also predicted by our model, as shown in Fig. 6b where cell trajectories for different values of matrix porosity φ are plotted. At higher porosity, the cell is able to migrate over a larger distance. As the porosity is decreased, however, the cell’s ability to migrate also decreases, and for sufficiently small values of φ, cell motion can be completely hindered. The presence of confinement also results in frequent directional changes in the cell trajectory. Such changes are caused by the interaction of a growing pseudopod with an ECM surface obstructing the cell motion. When a growing pseudopod collides with an ECM surface and the cell is unable to move, the activator and inhibitor concentrations in the reaction–diffusion equations are reset to unity. The absence of any force-generating activator patch causes the pseudopod to retract. In time, the reaction–diffusion system generates a new instability, and therefore a new pseudopod develops in an alternate location, allowing the cell to detour around an obstacle and continue its migration along a new direction. At times, the new pseudopod can form at a location opposite to the previous one, which causes the cell to completely reverse its direction (Fig. 6e).

A loss of directionality, although not severe, is also observed with increasing cell-ECM adhesion. Turns in trajectories appear to be more severe in the presence of adhesion than without adhesion. The mechanism by which adhesion results in a loss of directionality is shown in Fig. 7 and is explained as follows. When a pseudopod extends near an ECM surface, bonds are likely to form. If the pseudopod “collides” with the surface and is prevented from further growth, the activator concentration reverts back to the baseline level of unity. The adhesive bonds keep the local cell membrane near the obstacles, thereby reducing the chance of another pseudopod forming at the same location. Over time, a de novo pseudopod forms elsewhere on the cell, pulling it to a different direction, thereby changing the direction of cell motion. A cell with no adhesion may simply keep producing pseudopods at the same area, since the membrane could retract and grow another pseudopod. When adhesion is included, however, the membrane is kept near the obstacle, forcing pseudopods to grow in other directions.

Fig. 7

Cell-ECM adhesion-mediated directional change. a–f Sequence of such an event. Adhesive bonds keep the cell membrane close to the substrate, thereby preventing pseudopods from forming locally. Over time, de novo pseudopods form at other locations, altering the cell’s migration direction. The ECM has been removed from the image to show the mechanism

Although persistent random walk (PRW) cannot adequately represent cell movement, the persistence time can be extracted from the predicted MSDs and compared with published experimental data. Predicted persistence time was found to range from several seconds to over several hours, with many instances in the range 15–30 min. In comparison, in Ly and Lumsden (2009) which provides a survey of persistence times obtained in different experiments for 3D amoeboid migration in fiber matrices, values from several seconds to over 100 min have been noted. Specifically, for lymphocytes, the persistence time is in the range 1–10 min, while for leukocytes and macrophages it is slightly above 10 min, in agreement with our prediction. The predicted values show an increasing trend with increasing cell deformability, which can be attributed to the cell polarizing effect at higher values of α, as discussed above. For instance, a cell at α = 1 has randomly oriented pseudopods, and therefore lacks sustained direction, resulting in a small persistence time. As α increases, the cell becomes more polarized because of the focused pseudopod dynamics. Then, the cell can move in a certain direction longer, resulting a higher persistence time. However, no definitive trend of persistence time is observed with respect to bond constant and ECM microstructures.

3D migration can be further characterized by considering two quantitative results, the accumulated distance traveled and the invasion depth, which are important indicators in cancer research, since the capability of a tumor cell to travel through and invade distant tissues is a measure of its metastatic potential (Erkell and Schirrmacher 1988). The accumulated distance and maximum invasion depth are plotted in Fig. 8 against membrane deformability for several values of bond constant. Both accumulated distance traveled (L) and invasion depth (P) are observed to strongly depend on α, as increasing deformability allows the cell move farther and penetrate the matrix more effectively. Interestingly, both L and P tend to saturate at higher deformability, implying that once a cell becomes excessively deformable, it may no longer gain any advantage toward penetrating a matrix. Adhesion has an opposite effect: cells with higher values of σ show a reduced ability of matrix penetration, while cells without adhesion show greater penetration. The influence of adhesion on L and P is, however, much less compared to deformability, and is only observed at higher deformability. Considering a cell radius of 10 µm, our predicted values of L and P are 20–200 µm, and 10–80 µm, respectively. Several experimental works have reported the accumulated distance in vitro on the order of ~ 100 to 1600 µm, and invasion distance ~ 20–100 µm (Erkell and Schirrmacher 1988; Schmalstieg et al. 1986; Terry et al. 2012; Anguiano et al. 2017; Mckegney et al. 2001; Kitzing et al. 2010; Schor et al. 1983; Sapudom et al. 2015). There is, however, a strong degree of variability in experimental results, depending on cell and matrix properties, as well as the time elapsed during the experiment which is in the range of 12 h to a few days. Considering the physical time we have simulated (~ 2 h), the predicted L and P values are higher than those experimentally reported if the same time is considered. This difference is possibly due to a more confined environment than cells face in an actual matrix.

Fig. 8

a Accumulated distance traveled (L, dash lines) and invasion distance (P, continuous lines) versus cell deformability for various adhesion strengths. b, c Ratio (P/L) of invasion distance to accumulated distance for spherical lattice and orthogonal lattice, respectively, as a function of cell deformability for varying adhesion strengths. σ = 0, black; 1, red; 2, green; 3, blue

By taking the ratio of invasion distance to accumulated distance (P/L), we can define the directional persistence as done in Fritz-Laylin et al. (2017), which is displayed in Fig. 8b, c, where a magnitude of 1 indicates a highly persistent motion, while a value of 0 indicates a complete lack of persistence. Directional persistence is seen to increase with cell deformability, though a dependence on the ECM microstructure is observed. The increase is most prominent in the aligned rods and orthogonal lattice. For these matrices, the cell can penetrate more easily with increasing deformability, but can also travel greater distances through the matrix. An increase in persistence is therefore seen when the invasion distance grows faster than accumulated distance. For the random rods and spherical lattice, more subtle changes are seen: the persistence appears to peak near intermediate deformability, and then decreases as the membrane is further softened. This is a direct result of the accumulated distance increasing faster than the invasion distance, thereby leading to an almost constant directional persistence. The orthogonal lattice offers open channels for the cell to penetrate more easily, while the spherical lattice forces the cell to detour around obstacles, making penetration harder. The randomly oriented rod geometry is more porous, thereby allowing the cell to travel larger distances without a significant increase in penetration. P/L tends to slightly decrease with increasing adhesion, although a definite trend is not concluded due to the data variance.

3.2 Pseudopod dynamics

The influence of matrix geometry, cell deformability, and adhesion on pseudopod dynamics is shown in Fig. 9. Pseudopod dynamics are quantified in terms of the average number of pseudopods, and the pseudopod lifetimes. Figure 9a, b shows that the average number of active pseudopods is between 1 and 1.5. The fractional value appears because at some instances there are two active pseudopods present, while at other times only one is present. A small increase in the pseudopod number is observed with increasing cell deformability, because the reaction–diffusion system used to generate the bifurcating activator patch becomes more unstable with increasing values of α, so that at a given time there is a higher likelihood that two (daughter) pseudopods exist. The pseudopod number, however, does not show any noticeable dependence on ECM microstructure and bond strength. Figure 9c, d shows pseudopod lifetime, defined as the duration of the activator patch generating the pseudopod. The pseudopod lifetime is strongly influenced by cell deformability as it decreases with increasing α. This trend is also related to the increased instability of the activator patch with increasing deformability, as just mentioned. The porosity of the ECM affects pseudopod lifetimes; a reduced porosity results in frequent termination of the pseudopods due to cell-ECM collisions, while at higher porosity there are less cell-ECM collisions, and the pseudopod lifetime is therefore dictated mostly by the underlying Turing instability.

Fig. 9

Pseudopod dynamics. a, b Average number of pseudopods for orthogonal lattice (a), and random rod array (b). c, d pseudopod lifetime (dimensionless) for orthogonal lattice (c), and random rod array (d). Bond strength varies as σ = 0 (black), 1 (red), 2 (green), 3 (blue)

The reduced pseudopod lifetime at higher cell deformability is consistent with the increased persistence of the cell motion as discussed before, and can be understood from the observed pseudopod dynamics in the simulations (Fig. 6c, d). Once an activator patch bifurcates into daughter patches, they can meander over time away from their point of origin. A shorter lifetime also means a more frequent bifurcation, which ensures that the daughter pseudopods cannot meander farther away from their initial location. This results in pseudopods mostly localizing at the leading edge of the cell, giving it directional persistence at higher deformability.

3.3 Cell speed

The average cell speed V_C is shown in Fig. 10 for different ECM geometries, cell deformability (α), and adhesion strength (σ). Cell speed is observed to increase with increasing α. This is because a more deformable cell can more easily squeeze through the ECM pores. As noted above, higher deformation causes more focused pseudopods at the leading edge, and hence a higher resultant protrusive force, while reduced deformation results in pseudopods being generated all across the cell surface, leading to a reduced resultant force. The cell speed appears to approach a saturation at higher deformability, since the cell volume is preserved, and the cell surface area cannot dilate indefinitely. This suggests that deformability can only aid cell migration up to a certain point once confinement becomes a larger factor. As observed in each plot, the standard deviations also increase with α, since a more deformable cell can more easily navigate its movement through the matrix, frequently turning around the ECM structures.

Fig. 10

Average cell speed (dimensionless) as a function of cell deformability in different ECM and under varying adhesion strengths σ = 0 (black, □), 1 (red, Δ), 2 (green, $\nabla$), 3 (blue, Ο). Error bars represent standard deviations

In terms of different ECM microstructures, the predicted average speed is in a similar range for the 3D orthogonal lattice, aligned rods, and the spherical array. The randomly oriented rod geometry gives higher speeds. This is because the first three geometries have the same porosity (φ = 0.680), while the randomly rod geometry has higher porosity (φ = 0.88), allowing for more free space to migrate through. We have also performed additional simulations for spherical arrays by varying porosity, which showed a decreasing cell speed with decreasing porosity (Fig. 11). Together, these results suggest that amoeboid speed is strongly dependent on porosity, while microstructural details have a secondary effect for the relatively high porosity range considered in the simulations.

Fig. 11

Average cell speed for the spherical array as a function of porosity. Here α = 1 (black), 3 (red), 5 (green), 7 (blue)

This secondary influence of geometric details is rather subtle. As seen in Fig. 10, the cell speed at higher deformability is slightly higher for the orthogonal lattice compared to the aligned rods. The effect is also revealed in the manner cell speed saturates with increasing deformability. Figure 10 shows that V_C saturates at smaller values of $\alpha$ for the aligned rods, but at larger values for the spherical array, and more so for the orthogonal lattice. This is because although these three matrices have the same porosity, the effective pore diameter is different. We have calculated the effective pore diameter as the largest diameter of spheres that can be fitted within the void regions of a matrix. The ratio of the effective pore diameter to cell diameter is 1.16, 0.67 and 0.61, respectively, for the orthogonal scaffold, sphere array, and aligned rods, causing increasing confinement, which in turn results in the reduced speed and rapid saturation.

V_C as a function of adhesion strength $\sigma$ is presented in Fig. 12 for the orthogonal array, aligned rods, and sphere array. It shows that average cell speed remains nearly constant across the range of adhesion considered as compared to a non-adhesive cell. There is a slight downward trend in average speed as σ is increased in cases of high membrane deformability, which is more evident for the orthogonal lattice and sphere array compared to the aligned rods. This is because at higher deformability, the cell assumes a more elongated shape with greater surface area dilation. This leads to more cell-ECM contact and adhesive bond formation, and subsequent hindrance of the cell in the aligned rods and sphere arrays, both of which have reduced effective pore diameters smaller than the orthogonal lattice. Therefore, these results show that for amoeboid motility, weak adhesion has a secondary effect on cell speed.

Fig. 12

Average cell speed as a function of adhesion strength. a Orthogonal lattice, b aligned rods, c sphere array. α = 1 (black), 3 (red), 5 (green), 7 (blue)

3.4 Adhesion characteristics

In mesenchymal cells, focal adhesions form at the cell’s leading edge, and rupture at the cell rear; this establishes an adhesion gradient which helps pull the cell forward (Morley et al. 2014). For the amoeboid cell which has been described as having a diffuse and low expression of adhesion (Wolf et al. 2003), the mechanism is not as clear. Figure 13 shows the adhesion dynamics of amoeboid cells as observed in our simulations. It shows the distribution of bonds and contours of adhesion force on a cell surface from one simulation and at four consecutive instances. Unlike mesenchymal cells, the adhesion distribution is not continuous and is without any front-to-back gradient. Rather the distribution is highly localized in certain regions. As the cell is migrating through the matrix, these “islands” of bonds are transiently breaking and forming wherever the cell is close to the surface of the ECM microstructures. The instances chosen in the figure reveal bond distributions with a forward bias, rearward bias, central bias, and combination forward and rear bias. This specific example suggests a “tank-treading” of the adhesion sites as the cell migrates. This is further clarified by converting the 3D data to a line plot of normalized bond number density $n_{b} \left( s \right)$ and normalized adhesion force density $f_{b} \left( s \right)$ obtained by averaging over circumferential bands along the cell length $s$ from the front ($s = 1$) to the rear ($s = - 1$). The distribution moves from front to the rear of the cell over time. While in vivo cells do generate adhesive complexes which transition from the cell front to rear (Friedl and Wolf 2003), this is however more likely relevant to a mesenchymal migration type. Bonds in our simulations are transiently breaking and reforming, and the relative displacement between a bond cluster and an obstacle is nonzero. The cell is simply moving relative to the obstacle, while the bond cluster appears to slide down the cell membrane toward the rear.

Fig. 13

Contours of normalized adhesion force density and distribution of bonds (pink lines) on a cell surface shown at four instances during migration. The 3D plot is converted into a line plot in the lower panel. a–d shows bond distributions with a forward bias, rearward bias, central bias, and combination forward and rear bias, respectively. Arrows indicate direction of cell motion

Figure 14a shows the time-averaged normalized distribution of bond number density as a function of the cell length ($- 1 \le s \le 1$) for different geometries. Since the “islands” forming bonds are dynamically changing with time, this plot can be interpreted as a frequency distribution of where on the cell body the bonds are forming. The orthogonal lattice, aligned rod, and spherical array share the same distinct trend of bond density, showing a peak near the middle of the cell which diminishes toward the front and rear. This universal trend is seen across adhesion strength and deformability. This data further illustrates the distinction in the adhesion dynamics between amoeboid and mesenchymal cells, which in the latter case would exhibit a continuous front-to-back gradient (Morley et al. 2014). The centered distribution as predicted here is because of the concave nature of the cell surface that arises when in contact with the cylindrical or spherical microelements of the ECM. As the cell is moving through the matrix, it is deformed into a concave shape by the ECM scaffolds, and for the most part, the surface around the cell middle is in contact with the scaffold rather than the front and rear ends. This trend is however not observed for the randomly oriented rod array, for which the number of bonds is much lower due to higher matrix porosity.

Fig. 14

A Distribution of bond number density (time-averaged) along the cell length ($s$) from front ($s = 1$) to rear ($s = - 1$). Aligned rods (a, black), spherical array (b, blue), orthogonal lattice (c, green), and random rods array (d, red). B and C Average number of bonds in different geometry, deformability ($\alpha$) and adhesion strength ($\sigma = 1$, open square; $\sigma = 2$, filled squares)

The average number of active bonds is plotted in Fig. 14b, c for different geometries and as a function of deformability and adhesion strength. It is observed that matrix geometry plays a large role on the number of bonds formed, and on the trend in bond quantity as $\alpha$ and $\sigma$ vary. Both aligned rods and spherical array geometries show an increase in the number of bonds as $\alpha$ increases, while in the randomly oriented rods and orthogonal lattice, the bond count decreases. The change in bond count is as high as three-fold. For all geometries, bond count increases with increasing $\sigma$. The reason for differing trends in bond quantity as membrane deformability $\alpha$ is increased has to do with the cell-ECM contact area. A cell with a higher $\alpha$ is highly polarized, so depending on the geometry it can either interact with the matrix more or less frequently. The aligned rods and spherical array have smaller pore sizes, while the lattice has a larger pore size. Porosity is even larger for the randomly oriented rods. With increasing void space, the cell can navigate without much contact with the ECM. For the aligned rods and spheres, there is more cell-ECM interaction due to the smaller pore size, and therefore more bonds forming.

Figure 15a shows the net adhesive force on the cell in different geometries as a function of deformability and bond strength. It follows the similar qualitative trend of the bond counts. The predicted adhesive force is in the range 0.2–0.4 nN, which is two orders of magnitude less than that reported for adhesion-based migration, e.g., mesenchymal (Li et al. 2003). In contrast, the net protrusive force predicted here is about one order of magnitude higher than the adhesive force. Thus, the protrusive force is primarily acting against the fluid drag and friction from the ECM surface. It appears that the adhesive force in our simulations is acting against the cell motion, unlike in the mesenchymal motility. This is illustrated in Fig. 15b by plotting the time average of the normalized scalar product ($\varvec{x} \cdot \varvec{F}_{\varvec{a}}$) between the cell’s trajectory and the adhesive force vector. This quantity shows a decreasing trend with increasing cell deformability across all ECM geometry considered. Its value ranges between ~ 0 and − 0.5, implying that its influence is small but against the cell motion. The relatively small influence of the adhesive force is also seen if the instantaneous cell velocity is correlated to the instantaneous adhesion force and bond count. Figure 15C shows a time sequence of the adhesive force, bond count, and cell velocity. While bond count and adhesion force correlate well, their correlations with the cell velocity are weak.

Fig. 15

a Time-averaged adhesion force versus cell deformability for different ECM geometries and bond strengths. b Time-averaged and normalized scalar product between the instantaneous net adhesion force and cell trajectory. Line symbols in a and b are same as in Fig. 14a, b. c Time history of bond counts (blue), net adhesive force (red), and cell velocity (black) for a representative case (orthogonal lattice, $\alpha = 5, \sigma = 1$)

3.5 Influence of cell nucleus

We now consider the influence of the cell nucleus. The undeformed shape of the nucleus is spherical. The nucleus is modeled in a similar manner as the cell, that is, a viscous liquid surrounded by hyperelastic membrane. Similar models have been used by previous researchers (Jadhav et al. 2005; Marella and Udaykumar 2004; Khismatulin and Truskey 2005). The nucleus is the stiffest organelle within the cell, and consequently acts as a bottleneck to determine whether migration can be achieved (Wolf et al. 2013). The nuclear envelope provides mechanical support and protects genetic material from the cytoplasm (Denais et al. 2016). The degree of mechanical support, influenced by nuclear lamins, determines how strongly the nucleus will deform (Denais et al. 2016). Additionally, the nucleus is a highly dynamic organelle, that can also undergo remodeling. For instance, Wolf et al. (Wolf et al. 2013) experimentally found that while cell arrest occurred due to excessive nuclear deformation, the ratio of minimum deformed nuclear cross-section to undeformed cross-section was consistently near 10%, suggesting a significant nucleus deformation before subsequent rupture and cell damage occur (Denais et al. 2016). To account for the variability of nucleus stiffness, we consider as range of nucleus membrane elasticity that is 2 to 100 times greater than the cell membrane elasticity. We also vary the nucleus diameter in the range 0.5–0.75 the cell diameter. The nucleus has been observed to follow a power-law rheology: stiff at short time scales but able to soften and deform at longer time scales in a viscoelastic fashion (Dahl et al. 2005). This results in a nucleus protected from transient applications of stress, while allowing itself, and the cell, to adapt to its environment over time. Additionally, by applying cycles of deformation and relaxation to cells in a microfluidic chamber, Mak and Erickson (2013) have shown that if relaxation is not fully achieved, a cell’s initial deformation event facilitates subsequent deformation events, thereby imparting faster strain dynamics and more efficient invasion characteristics to a cell. The model implemented here does not include such complex rheology, although the overall viscoelastic behavior is retained by the combined effect of the viscous fluid and elastic membrane. However, as noted before similar models of nucleated cells (e.g., leukocytes) have been used in the past by us as well as other researchers and were validated for leukocyte dynamics (Jadhav et al. 2005; Marella and Udaykumar 2004; Khismatulin and Truskey 2005; Pappu et al. 2008).

Figure 16 shows sequences of migrating cells with a nucleus included, along with 3D trajectories. The simulations predict that cells with larger and stiffer nuclei are unable to penetrate the matrix more than one cell diameter (16B), while those with smaller and softer nucleus penetrate significantly (16A). As noted in Denais et al. (2016), amoeboid migration may become unsustainable due to the limited deformability of the cell nucleus, which can lead to cell arrest, or if deformed excessively, may rupture and lead to cell death. We have observed cases of cell arrest occurring due to the matrix topology constricting the nucleus and preventing forward cell migration. As seen in the figure, bifurcating pseudopodia have extended the cell forward, leaving the nucleus near the rear of the cell. After several iterations of pseudopod formation, however, the cell reorients by creating new pseudopods in a different direction which may or may not be successful to induce further motion depending on the nucleus properties. The observed arrested cell was unable to pull the nucleus through the constriction, instead opting to repolarize. An explanation for this failure to penetrate the obstruction is likely due to the lack of an explicit contractility model, since myosin-II-mediated contractility is known to help push the nucleus forward in less adherent immune cells and certain cancer cells in dense ECM (Petrie and Yamada 2015; Lämmermann et al. 2008). During periods of rapid migration through obstructions, the nucleus was seen to move toward the rear of the cell, which is known to occur during retrograde cortical flow (Petrie and Yamada 2015). On the other hand, streaming flow generated by protruding pseudopods was observed to bring the nucleus near the cell front, which has been observed in blebbing cells (Liu et al. 2015). Fast changes in cell polarization were also seen to set the nucleus near either the cell front or rear. On average, though, the nucleus position remains near the center half of the cell.

Fig. 16

Influence of nucleus. a, b Sequence of cell migration through a spherical matrix with smaller and softer nucleus (a D_nuc/D_cell = 0.5, G_nuc/G_cell = 2), and with larger and stiffer nucleus (b D_nuc/D_cell = 0.7, G_nuc/G_cell = 100). The nucleus is colored in dark blue. c 3D trajectories for varying nucleus size and stiffness; black, red and green represent G_nuc/G_cell = 100, 10, 2, respectively; thick and thin curves represent D_nuc/D_cell = 0.75 and 0.7, respectively. The dash circle which represents the undeformed cell is shown to indicate that the cell with stiffer and larger nucleus is unable to penetrate the matrix more than one cell diameter. d Comparison of average cell speed of the nucleated cells (D_nuc/D_cell = 0.5; thick lines) and non-nucleated cells (thin lines) with (red) and without (black) adhesion

We then analyzed the motility behavior of the nucleated cells for D_nuc/D_cell = 0.5 that are able to penetrate the matrix, and found that different measures such as MSD, ACF, penetration and time average cell speed remain nearly same as those of non-nucleated cells. Figure 16d compares the average cell speed for the two cases. This suggests that for the cells that are able to penetrate, the addition of a nucleus did not fundamentally alter cell behavior. We attribute this to pore sizes larger than the effective nucleus size.

4 Discussion and conclusion

We have simulated pseudopod-assisted migration of deformable amoeboid cells through 3D microstructure resolved, idealized ECM. The model used here is a multiscale, multiphysics, coupled fluid‒structure‒biochemistry interaction model. The key elements of the model are as follows: (1) A finite-element model for cell membrane deformation. (2) A coarse-grain pattern formation model for acto-myosin biochemistry represented in terms of nonlinear reaction–diffusion equations that is able to generate Turing instabilities. When combined with cell deformation, the model generates membrane protrusions that grow, bifurcate, and retract, mimicking experimentally observed pseudopod dynamics (Levine and Rappel 2013). (3) Cytoplasmic and extracellular fluid motion governed by the unsteady Stokes equations and solved using combined spectral and finite-volume methods. (4) Microstructure (pore-scale geometry)-resolved ECM represented by rigid scaffolds immersed in the extracellular fluid. (5) Cell-ECM adhesion represented by receptor-ligand interactions and modeled by a Monte Carlo method. (6) Coupling of different elements of the model by two variants of the immersed-boundary methods (IBM): a continuous forcing IBM for the deformable cell membrane, and a direct forcing (ghost-node) approach for the ECM scaffold. Previous computational models often neglected cell deformation (Schlüter et al. 2012; Zaman et al. 2005), microstructural details of the ECM (Sakamoto et al. 2014; Hecht et al. 2011), or intra- and extracellular fluid motion (Najem and Grant 2013; Elliott et al. 2012), and used ad hoc methods for cell-ECM adhesion. The present model is a further improvement to such efforts, and goes closer toward an exact and complete model of 3D motility.

4.1 Model deficiencies

The ECM in our work is considered rigid. The justification for this is that amoeboid cells lack the proteolytic capability needed for matrix degradation and instead use pre-existing pores and channels to migrate, which can reasonably be modeled as rigid. In absence of strong adhesion, such rigid surfaces likely provide sufficient traction on highly deformable cells. Nevertheless, a more comprehensive model of ECM is needed in future that would incorporate several important properties as discussed below.

First, the ECM is not composed of a singular material. It is a complex and heterogeneous combination of collagen and elastin fibers, adhesion proteins which bind the ECM together, viscoelastic glycosaminoglycan gels such as hyaluronan, and copious other proteins (Theocharis et al. 2016). Furthermore, the ECM transcends several length scales, from fibrils and signaling molecules on the nano-scale, to overall tissue structure on the micro and macro scale (Stevens and George 2005). A diverse set of arrangements and combinations of the aforementioned ECM components allows for a large variety of tissues to exist within the body, such as cartilage, connective tissue, bone, etc. Second, the ECM is deformable and may undergo elastic and viscoelastic deformations caused by migrating cells. Malandrino et al. observed both elastic and plastic pulling of collagen and fibrin fibers toward the body of migrating breast cancer cells (Malandrino et al. 2019). Fiber displacement during migration of neutrophilic cells in 3D collagen matrices was noted in Chen et al. (2014). MDA-MB-231 breast cancer cells were shown to be highly invasive with large forces focused on the poles due to anisotropic adhesions, resulting in high matrix contraction near the cell (Koch et al. 2012).

There are, however, some justifications for treating the ECM as rigid in studying amoeboid motility. While mesenchymal cells characterized by strong adhesion, contractility, and proteolytic activity cause permanent restructuring of the local ECM environment near the cell (Paul et al. 2017; Malandrino et al. 2019), pseudopod-driven amoeboid cells with little to no adhesion and lack of both proteolytic capability and contractility only apply minor fiber deformation as they migrate (Friedl and Wolf 2003; Chen et al. 2014). These cells often migrate through pre-existing pores and channels, without a need to further degrade the matrix. Wolf et al. found amoeboid-like cells migrating through collagen generated no structural remodeling, bundling, or destruction of the network after cell detachment, with individual fibrils returning to their original positions (Wolf et al. 2003). However, the degree to which ECM elements are deformed remains highly dependent on the cell type and properties of the ECM (Charras and Sahai 2014; Sapudom et al. 2015), and should be studied by more advanced models. While the ECM geometries considered in this work are idealized and simplified versions of in vivo tissue, they can represent scaffolds generated for the purpose of tissue engineering which often have regular lattice microstructures. Depending on the application, take bone for instance, scaffolds may be required to have appropriate mechanical strength (Thavornyutikarn et al. 2014). This implies the matrix may be mechanically superior to that of collagen, and hence, would show minor-elastic or negligible deformations.

The current framework of immersed-boundary method (IBM) can be used to include deformable ECM with realistic geometry. The extracellular fluid can be modeled as viscoelastic fluid using appropriate conservation laws (e.g., Oldroyd model). Elasticity and viscosity of the polymeric constituent can be altered in such models to study their effect on cell motility. The continuous forcing IBM that is used here to treat the deformable cell membrane can also be used to model deformable fibers, and any other cells present in the ECM as embedded (i.e., immersed) objects. Fiber deformation can be modeled following viscoelastic constitutive models and finite-element method, in the same way that cell membrane deformation is numerically treated. Therefore, fiber density, orientation, and heterogeneity in different fiber properties (e.g., length, thickness, and elasticity of elastin versus collagen fibers) can be readily considered in this approach.

Another deficiency in our model is the lack of decision-making components based on feedback control, topographical alignment, and contact guidance. The ECM can guide cells that encounter multiple pathways during migration. Mak and Erickson (Mak and Erickson 2014) ran cancer cells through a bifurcating microfluidic array while varying channel dimensionality and directionality. Cell affinity for choosing larger channels was observed to occur at high probabilities, though when sub-nucleus sized branches were made colinear with the original channel path, the probability of choosing the larger channel was reduced. These results indicate that while narrow channels are not a first preference of the cell, polarization effects can bias cells into choosing them over larger channels. Because our model lacks an explicit membrane contractility mechanism, however, attempts by the cell to squeeze through very narrow constrictions usually end in failure and an alternate path being found. Another dynamic not considered in our simulations is hydraulic pressure, which can also affect cell decision making in 3D matrix. As the cell migrates through dense tissue, a tight seal may develop between the cell and ECM. The resulting pseudopodial protrusions would then exclusively need to push the column of interstitial fluid ahead to advance, causing a hydraulic pressure (Zhao et al. 2019). Experiments have shown that cells presented with multiple avenues under these conditions will usually select the channel with the least hydraulic resistance, even if a chemical gradient directs the cell elsewhere (Zhao et al. 2019; Prentice-Mott et al. 2013). Because our ECM is not sufficiently dense to form a tight seal with the cell, we cannot consider the effect of hydraulic resistance on amoeboid migration. However, we note that the tug-of-war dynamic mentioned in Prentice-Mott et al. (2013) does occur in our simulations, though the formation and loss of competing leading edges is attributed to the stochastic nature of the cell and not by hydraulic resistance. Finally, contact guidance is frequently observed to occur, where cells orient themselves along structural elements in a persistent manner (Brábek et al. 2010). In fact, Doyle et al. (2009) demonstrated an extreme case of contact guidance using fibroblasts migrating along singular fibers. These ‘one-dimensional’ geometries were shown to prevent the formation of lateral cell protrusions, thereby increasing the efficiency of migration along a single axis. In the case of all four matrix geometries used in our simulations, in particular the aligned rods, no bias was seen during migration along fiber elements, most likely due to the weak and diffuse adhesive bonds exhibited by the amoeboid cell. Furthermore, nearby elements of the matrix allow the cell multiple options for migration. Mesenchymal cells with stronger adhesion and adhesion gradients like that used in Doyle et al. (2009) would be more likely to display contact guidance.

4.2 Comparison with experiments

The model predicts highly deformed amoeboid cell shapes as observed in experiments (Lämmermann and Sixt 2009; Van Haastert 2011), and also characteristic (pseudo) random-like migration (Friedl and Wolf 2003). Direct quantitative comparisons against experimental results can be made in terms of migration speed. Dimensionless average cell speed predicted was approximately in the range of 0.05–1.00, which scales to 0.3–6.0 μm/min using a cell radius of 10 µm and a membrane protein diffusivity of 1 µm²/s. Experimental data on amoeboid cell migration speeds in 3D matrices vary in the literature based on cell type and ECM properties as 0.07–0.6 μm/min (Wolf et al. 2003), 2–25 μm/min (Brábek et al. 2010), 4 μm/min (Petrie and Yamada 2015), 5–8 μm/min (Sarris et al. 2012), 0.1–20 μm/min (Friedl and Wolf 2003), and 1–5.5 μm/min (Lämmermann et al. 2008). Our prediction, therefore, falls within the range of these measurements. Comparisons are made in terms of persistence time and penetration distance. Predicted persistence time was found to range from several seconds to over several hours, with many instances in the range 15–30 min. In agreement, the experimental persistence time is in the range 1–10 min for lymphocytes, and slightly above 10 min for leukocytes and macrophages (Ly and Lumsden 2009). Predicted values of accumulated distance traveled and invasion depth are 20–200 µm, and 10–80 µm, respectively. Experimental works have reported the accumulated distance in vitro on the order of ~ 100 to 1600 µm, and invasion distance ~ 20–100 µm (Erkell and Schirrmacher 1988; Schmalstieg et al. 1986; Terry et al. 2012; Anguiano et al. 2017; Mckegney et al. 2001; Kitzing et al. 2010; Schor et al. 1983; Sapudom et al. 2015). This difference arises because of significantly longer time in experiments than the simulations. Further comparison is done using directional persistence, defined as the ratio of the invasion depth to the accumulated distance traveled. This quantity is predicted to be in the range ~ 0.1–0.6, which matched very well with experimental range of ~ 0.1–0.5 (Peela et al. 2016), where breast cancer cells were allowed to invade 3D hydrogels in vitro, and a mean of ~ 0.375 (Fritz-Laylin et al. 2017) where HL-60 cells followed chemoattractant in a glass chamber.

It should be noted that a quantitative comparison with experimental data is not easy because cell motility is strongly dependent on cell deformability and matrix geometric and materials properties, which are inputs to our model. Such information is often not reported. Furthermore, cells and ECM can both undergo remodeling during an experiment, which affects migration. Because of the limitations of the model, we refrain from simulating any specific cell type and specific ECM. Given such limitations, a more quantitative comparison is not possible at this point.

To establish whether cell migration can be described as a random walk or persistent random walk, we used predicted cell trajectories and extracted the mean-squared displacement (MSD) and velocity auto-correlation function (ACF). We found a super-diffusive behavior (~ $t^{1.8}$) at intermediate time scales. Thus, cell migration does not follow a random walk despite the lack of external signals or gradients. This qualitatively agrees with the experimental studies in Dieterich et al. (2008) and Wu et al. (2014). These results are further compounded by examining the ACF curves, which do not show an exponential decay, indicating that migration cannot be classified as a persistent random walk, and agreeing with the experimental study in Dieterich et al. (2008).

While a large volume of work exists on cell motility on 2D substrates, amoeboid motility is inherently 3D. Furthermore, amoeboid cells migrate without degrading and remodeling the ECM, rather using the pre-existing pores and defects therein. Additionally, these cells use none or very weak adhesion to the ECM. In the absence of matrix degradation and strong adhesion, amoeboid cells must utilize their deformability and the microstructure of the ECM to gain enough “traction” to migrate. Despite recent progress in modeling cellular motility in 3D, there is a lack of systematic evaluations of the role of ECM microstructure, cell deformability, and weak cell-ECM adhesion on 3D motility. The current study has sought to fill this knowledge gap. The key finding of this study is that cell deformability and ECM porosity strongly affect migration behavior, while microstructure geometry and cell-ECM adhesion have weaker effects.

Migration speed is shown to be a strong function of cell deformability ($\alpha$), with more deformable cells migrating at higher speeds. This trend agrees with the notion that malignant cells are softer than benign cells (Cross et al. 2007), and that cell stiffness is a potential biomarker (Xu et al. 2012). We showed that higher deformability allowed the cell to take on a polarized shape, with pseudopod generation almost exclusive to the leading edge. This allows the cell to penetrate the matrix by squeezing through narrow pores, and undergo a long time directed (or, persistent) migration. Conversely, at reduced deformability, pseudopod generation is random and distributed over the entire cell surface, which prevents any long-term directed motion. The cell also lacks sufficient compliance to squeeze through matrix pores. As such, the persistence time extracted from the predicted MSDs increased with increasing cell deformability.

The influence of cell deformability $\alpha$ on matrix penetration is quantified using three parameters: invasion distance, total accumulated distance, and their ratio, termed the directional persistence (Erkell and Schirrmacher 1988; Schmalstieg et al. 1986; Terry et al. 2012; Anguiano et al. 2017; Mckegney et al. 2001; Kitzing et al. 2010; Schor et al. 1983; Sapudom et al. 2015; Fritz-Laylin et al. 2017). The accumulated distance increases with increasing $\alpha$, due to the same reasons discussed above. The invasion distance and directional persistence, on the other hand, saturate at higher $\alpha$. This is because at high enough $\alpha$, the Turing patterns become more unstable, leading to frequent bifurcations, and hence directional change in cell motion. Thus, a very deformable cell can accumulate a larger distance traveled without significantly penetrating the matrix. Then after a certain point, additional increases in cell deformability may not help in matrix penetration. Highly deformable neutrophil cells have been noted to experience directional changes, thus supporting this notion (Henry et al. 2014). This notion is further corroborated by considering the pseudopod dynamics, which also show a strong dependence on $\alpha$. The number of active pseudopods was seen to slightly increase with increasing $\alpha$, which leads to increased cell speed. Pseudopod lifetimes, on the contrary, decreased significantly. Both these trends are due to the increasing instability of the underlying Turing patterns.

The microstructure of the ECM is defined using porosity (or conversely, confinement) and pore-scale geometry. Amoeboid cell speed is shown to strongly depend on porosity $\varphi$, with larger porosities leading to higher cell speed and directed migration, while higher confinement reduced migration capacity and instituted more frequent cell directional changes. Confinement can result in complete reversal of cell trajectory in our simulations; such behaviors have also been observed in experimental studies of leukocyte migration in 3D (Gaylo et al. 2016).

In terms of pore-scale geometry, four different scaffolds are considered: aligned rods, randomly oriented rods, orthogonal lattice, and spherical array. While none represent the exact details of the ECM, they do have some degree of resemblance to in vivo or artificial tissue. We find that the influence of such pore-scale details on motility behavior is secondary but subtle. They did not result in any significant/conclusive difference in pseudopod dynamics and persistence time. They did, however, affect the way cell speed, directional persistence, and cell-ECM adhesion depended on cell deformability. Orthogonal lattice and aligned rods offer straight passages, while the randomly oriented rod array has higher porosity, leading to increased directional persistence. The spherical lattice forces the cell to detour frequently around obstacles, resulting in reduced persistence. While the average cell speed was predicted to be in similar ranges across different microstructures for a given porosity, it was slightly lower and saturated at relatively smaller values of $\alpha$ for the aligned rods. These trends are dictated by the effective pore size. While the porosity remained constant, the effective pore size progressively decreased from orthogonal lattice to spherical lattice to aligned rods. As such, the cell speed was the lowest in the aligned rods. A major difference was seen for the average number of bond counts, which increased significantly with increasing $\alpha$ for aligned rods and spherical lattice, but showed opposite trends for the orthogonal lattice and random rod array. Since the aligned rods and spherical lattice have a much smaller pore size, cell-ECM contact, and therefore bond numbers increased as the cell became more polarized with increasing $\alpha$. In contrast, at higher pore sizes as in the case of the orthogonal lattice, and for higher porosity as for the random rods, a highly polarized cell shape caused a reduced cell-ECM contact, and hence lowered bond counts.

Cell-ECM adhesion is mediated primarily by integrin molecules. We model adhesion using stretch-dependent reaction rates, and linear spring-like behavior of bonds. The stochastic nature of bond formation and breakage is also implemented in the model. Although amoeboid cells are known to migrate with weak or no adhesion (Wolf et al. 2003), the role of adhesion, if any, is not clear. In this work, we selectively consider the A1-amoeboid mode of migration, which is characterized by low adhesion (Petrie and Yamada 2015). Experiments have shown that in considering 3D amoeboid migration, adhesive interactions between the cell and matrix are low in magnitude (Wells et al. 2013). Furthermore, while amoeboid cells have been shown to effectively migrate in 3D in the total absence of integrins, the difficult question still remains of whether cells do not need integrins to migrate, or simply that integrins are unnecessary for migration (Paluch et al. 2016). It is from this question that we consider the effect of weak adhesion. We predicted that adhesion had a secondary influence on motility. While pseudopod dynamics did not change with adhesion, a slight loss of directionality, matrix penetration, and cell speed was observed. Experimental findings support this prediction, as integrins were shown to impede invasion, as seen in types of melanoma, breast, and colon cancers (Friedl and Wolf 2003). The correlation between instantaneous cell velocity and bond count or adhesion force was observed to be weak. Also, the predicted adhesive force is two orders of magnitude less than that reported for mesenchymal phenotypes which use strong cell-ECM adhesion (Li et al. 2003). This is likely a testament to the diffuse, transient, low-affinity nature of the bonds. Unlike a continuous front-to-back adhesion gradient which is often established during strong adhesion-mediated motility (Morley et al. 2014), we predicted that for amoeboid cells, adhesion distribution is highly patchy, and focused around the gut of the cell. Perhaps the most significant observation witnessed from including adhesion in our simulations is the geometrical dependence on bond count as a function of cell deformability as described in the results section. This implies that the cell-ECM coupling can dictate the expressions of a cell as observed by various experimental studies (Huttenlocher and Horwitz 2011; Friedl and Wolf 2010), which may have implications in tissue engineering and cancer research. Along the same thought, whether an amoeboid cell requires adhesion to migrate in 3D is strongly dependent on matrix confinement (Liu et al. 2015; Friedl and Wolf 2010), and as our model can simulate virtually any matrix geometry, including adhesion serves as good generalization.

4.3 Experimentally testable predictions

This study led to some novel predictions which can be experimentally validated. One such result is the apparent geometric dependence on bond count with increasing deformability. It was found that as cell deformability increased, cells were more likely to generate bonds in the aligned rod and spherical array geometries, while the reverse was seen to happen in the orthogonal lattice and random rod geometries. This cell-ECM coupling may pose an interesting feedback problem, since as the matrix conforms the cell to its geometry, the cell likewise responds by adapting itself (Charras and Sahai 2014). Therefore, experiments may demonstrate that while one type of matrix geometry allows for more cell-ECM bonds as cell deformability is altered, another geometry may show a reduction in bonds using the same adjustments. Since cells can interpret the ECM through adhesive connections (Desgrosellier and Cheresh 2010), more or less active bonds on the cell may trigger different phenotypes to be expressed, and therefore different cell behavior.

Another novel prediction is the increasing action of adhesive bonds against the direction of amoeboid migration with increasing deformability. The effect is more pronounced in the spherical array and orthogonal lattice geometries, though adhesion is also seen to inhibit migration in the aligned rod and random rod geometries to a lesser extent. Because cells possess compensation mechanisms which allow the alteration of their migratory behavior (Wolf et al. 2003), our results suggest the sphere and lattice geometries are more hostile to migration with diffuse adhesion and may favor an amoeboid-based migration with no adhesion present. On the other hand, the aligned and random rod geometries do not penalize adhesive migration as much, and therefore may be more likely to welcome adhesion-based amoeboid migration, or even mesenchymal migration. Experiments could validate whether amoeboid cells lacking adhesive capabilities could migrate more or less efficiently in these geometries than other cell-phenotypes, and vice versa. Additionally, these predictions may aid in the design of tissue-engineered scaffolds toward the effectiveness of particular cell lines.

To conclude, we presented a multiscale, multiphysics, coupled fluid‒structure‒biochemistry interaction model to study amoeboid motility in 3D idealized ECM scaffolds by resolving pore-scale geometry of the ECM, cell deformation, intra- and extracellular fluid motion, cell-ECM interaction, and adhesion dynamics. The key finding is that cell deformation and matrix porosity strongly influence cell motility, while adhesion and microscale structural details of the scaffold have secondary but subtle effects. While several predictions can be supported by published experimental findings, a few limitations of the model are worth noting. First, while the ECM is modeled as fluid-filled porous medium, the scaffold is taken to be rigid. In reality the scaffold could be deformable. Second, both the intra- and extracellular fluids are assumed to be Newtonian, while in reality they are viscoelastic. Further, we considered only pseudopod-assisted motility, while other forms such as contractility-based motility also exist for amoeboid cells (Liu et al. 2015; Morley et al. 2014). Inclusion of these modules in our model will be considered in future.