Amoeboid cells propel by generating pseudopods that are finger-like protrusions of the cell body that continually grow, bifurcate, and retract. Pseudopod-driven motility of amoeboid cells represents a complex and multiscale process that involves bio-molecular reactions, cell deformation, and cytoplasmic and extracellular fluid motion. Here we present a 3D model of pseudopod-driven swimming of an amoeba suspended in a fluid without any adhesion and in the absence of any chemoattractant. Our model is based on front-tracking/immersed-boundary methods, and it combines large deformation of the cell, a coarse-grain model for molecular reactions, and cytoplasmic and extracellular fluid flow. The predicted shapes of the swimming cell from our model show similarity with experimental observations. We predict that the swimming behavior changes from random-like to persistent unidirectional motion, and that the swimming speed increases, with increasing cell deformability and protein diffusivity. The unidirectionality in cell swimming is observed without any external cues and as a direct result of a change in pseudopod dynamics. We find that pseudopods become preferentially focused near the front of the cell and appear in greater numbers with increasing cell deformability and protein diffusivity, thereby increasing the swimming speed and making the cell shape more elongated. We find that the swimming speed is minimum when the cytoplasm viscosity is close to the extracellular fluid viscosity. We further find that the speed increases significantly as the cytoplasm becomes less viscous compared with the extracellular fluid, resembling the viscous fingering phenomenon observed in interfacial flows. While these results support the notion that softer cells migrate more aggressively, they also suggest a strong coupling between membrane elasticity, membrane protein diffusivity, and fluid viscosity.

## I. INTRODUCTION

Cell motility is ubiquitous in nature, existing in several different forms. One such form is pseudopod-driven motility, which is present in many single and multicellular organisms. Pseudopodia are finger-like protrusions of the cell body that continually grow, bifurcate, and retract.^{1,2} For many amoeboid cells such as *Dictyostelium discoideum* (*Dicty*), pseudopod-driven motility is the primary form of motility. This form of motility is also responsible for many physiological processes in the body. Examples include extravasation of neutrophils through blood vessel walls towards regions of infected tissues, migration of epithelial cells during wound healing, and embryonic cell migration during fetal development.^{3,4} A similar pseudopod-driven motility is also observed in metastatic cancer cells.^{5,6}

A pseudopod is generated when a cell’s external receptors are activated by some stimulus. With this signal, the cell utilizes specific proteins like Arp2/3 to activate nucleation sites, converting G-actin monomers into F-actin filaments.^{7} The growing actin filaments then exert a protrusive force, pushing the membrane outward. In order to crawl on a substrate, the cell must also adhere and exert a force while the rear of the cell retracts. This retraction is accomplished using Myosin ii proteins that generate a tensile force, thus breaking the adhesion sites near the cell posterior. Although crawling is the most common mode of amoeboid motility, it was recently observed that amoeboid cells and neutrophils can swim while suspended in a fluid using a similar pseudopod-assisted mechanism.^{8–10} Unlike in crawling, cell-substrate adhesion is not present in swimming. Pseudopod protrusions generated in a swimming *Dicty* were shown to bifurcate and retract as they moved towards the back of the cell. This capability to swim resonates with Purcell’s scallop theorem,^{11} which implies that moving cellular projections could create a non-reciprocal motion in time, thus generating a drag force sufficient for motion on the microscale.

Deformability of the cell is critical for its migration^{1,8–10,13} and has been proposed as a biomarker for the metastatic potential of cancer cells.^{14} Metastatic cells were found to be significantly softer than benign cells.^{15} Cell softening has been shown to result in higher migration speeds for fibroblasts and breast cancer cells and to correlate with increased proliferation of malignant cells.^{16} Two important parameters, among others, that control cell deformation are membrane elasticity and cytoplasmic rheology. Membrane tension must be overcome by the protrusive force to push the membrane outward. The membrane curvature and tension may in turn dictate localization of membrane proteins such as Ras and affect the actin polymerization rate and cell motility.^{17,18} Rheology of the cytoplasm and extracellular fluid also plays a critical role.^{19,20} A streaming flow occurs in a growing pseudopod as it extends against the drag exerted by the extracellular fluid.

While a significant advancement has been achieved in modeling amoeboid motility, a full 3D modeling remains a challenge.^{21} Many models consider only the biochemical processes but not cell deformation and fluid flow.^{22,23} Several models that consider cell deformation are two-dimensional.^{24–28} Many studies are limited to rather simple cell shapes.^{12,27,29} In the following, we discuss a few noteworthy modeling studies that are relevant for the current work. Vanderlei *et al.*^{30} developed a 2D model of a motile cell using an immersed-boundary method that resolves cell deformation, internal and external fluid flows, and a reaction-diffusion system in the entire volume of the cell. This model, thus, allows advection-diffusion of biochemical substances within the cell as it moves and deforms. Bottino and Fauci^{31} developed a 2D model of amoeboid cell deformation and migration using an immersed-boundary method in which the cytoskeleton is modeled as a dynamic network of springs immersed in a fluid. Their model was able to generate protrusive and contractile forces and the attachment-detachment cycle in a cell crawling over a substrate. Elliott *et al.* developed a 3D model for a deformable cell with bifurcating pseudopods.^{32} Their model and that of several others,^{12,13,25,26,28,29} however, do not include the cytoplasm and extracellular fluid, which as noted above play a critical role in cell motility. Farutin *et al.*^{3} and Wu *et al.*^{34} considered a deformable cell driven under a prescribed axisymmetric oscillating force. While the presence of the internal and external fluids is considered, the model does not represent the complex pseudopod dynamics and random cell movement. Najem and Grant^{35} used a phase field approach to simulate migration of neutrophils in 3D in response to external cues but neglected the presence of fluids. Moure and Gomez^{36} developed a phase-field approach to simulate cell movement in 2D with and without obstacles and in 3D through a fibrous network. Their model includes consideration for cytosolic fluids, cytoplasmic rheology, continuum models for acto-myosin dynamics, and membrane-signaling dynamics using an activator reaction-diffusion.

In this article, we present a 3D model of pseudopod-driven swimming of amoeboid cells using front-tracking/immersed-boundary methods. Our model combines coarse-grain biochemistry, cell deformation, and cytoplasmic and extracellular fluid motion. The cell is freely suspended in a fluid; thus, the cell-substrate adhesion is not present. Note that for amoeboid migration, a strong cell-substrate adhesion is not always needed. Under a low expression of adhesion molecule integrins, tumor cells are known to migrate using the amoeboid mode, while integrin adhesion is nearly dispensable for neutrophils.^{37,38} A swimming amoeba, therefore, also serves as a model for cell migration under a weak adhesion.

## II. MODEL

Direct numerical simulation (DNS) of pseudopod-driven swimming of amoeboid cells is considered in 3D. In a swimming amoeba, the protrusive force from the actin filaments deforms the membrane working against the membrane tension and the drag from the surrounding fluids. The cell membrane evolves in time following complex physical laws governing its mechanics. Reaction-diffusion of the membrane proteins also occurs on this deforming surface. The extension of the membrane and movement of the cell generate fluid motion both inside and outside the cell, which in turn determines the local fluid drag. The problem then becomes a complex moving interface (fluid/structure interaction) problem where the location of the evolving interface is not known *a priori* and is obtained as a solution along with the flow variables and protein concentrations.

### A. Cell membrane model

The cell membrane is composed of a lipid bilayer and an underlying cortex. The thickness of the membrane is on the order of nanometers and hence much smaller compared with the size of the whole cell which is on the order of tens of micrometers. Then, a continuum description of the cell can be used by neglecting the molecular details of the membrane. By this assumption, the cell is represented as a viscous drop made of cytoplasm surrounded by a zero-thickness (i.e., 2D) hyperelastic membrane which is suspended in the extracellular fluid. The membrane is assumed to possess resistance against shearing, area dilation, and bending. The first property is due to the cortex, while the latter two are primarily due to the lipid bilayer. Constitutive laws for the surface energy functions are introduced to obtain the stresses generated in the membrane due to deformation. For the shearing deformation and area dilation, the following strain energy function that was developed for the composite lipid bilayer/cortex of an erythrocyte membrane is used:

Here *G*_{S} is the membrane shear elastic modulus, *CG*_{S} is the area dilation modulus, $I1=\epsilon 12+\epsilon 22\u22122$ and $I2=\epsilon 12\epsilon 22\u22121$ are the strain invariants of the Green strain tensor, and *ε*_{1} and *ε*_{2} are the principal stretch ratios.^{39} The Green strain tensor is defined as **E** = (**F**^{T}·**F** − **I**), where $F=\u2202x/\u2202X$ is the deformation gradient of the current configuration *x* relative to the original configuration **X** of the membrane. The parameter *C* controls the amount of area dilation and is fixed to be 1; thus membrane dilation is allowed. We note that the plasma membrane is nearly inextensible allowing only about 2%-3% change in the surface area. However, experimental studies have shown 20%-30% change in the surface area of Dicty.^{40} Such an increase in the surface area can happen by endo- and exocytosis of the cell membrane. Extensive fusion of vesicles to the membrane of Dicty cells was experimentally observed.^{40} Additionally, many motile cells, e.g., leukocytes, exhibit a large number of small-scale wrinkles and microvilli which may act as a reservoir for additional area.^{41} These structures are too small to be directly resolved in our model. Therefore, we assume an extensible area in order to allow for the utilization of such fine structures as the membrane reservoir. In our simulations, once the initial transient period that lasts for a dimensionless time of ∼1 is passed, the surface area fluctuates about a mean with a variation of 7%-20% as reported in the aforementioned experiment.^{40} For the 2D membrane, the principal elastic stresses (or, tensions) arising from shearing and area dilation are given by

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

where *E*_{B} is the bending modulus, *κ* is the mean curvature, *c*_{0} is the spontaneous curvature, and *S* is the surface area.^{42} The spontaneous curvature serves to describe the effect of molecular difference in the two leaflets of the membrane bilayer. The two leaflets typically differ in their chemical composition and are exposed to chemically different surroundings.^{43} Due to this difference, the membrane prefers to curve in a certain manner which is quantified by the spontaneous curvature, which, in turn, determines the resting shape of the membrane.

### B. Reaction-diffusion of membrane proteins

There are several different modes that exist in the literature for the migration of deformable cells. The choice of a particular model must be dictated by the experimental observation and the fundamental mechanism behind the specific form of cellular motility considered here which is the pseudopod-driven motility of Dicty. Excellent reviews are given in Refs. 44 and 45 discussing the unique characteristics observed in the pseudopod-driven motility of Dicty, and the requirements for a suitable model are as follows:

Chemotactic cells are highly sensitive to small changes in the concentration gradient, suggesting an internal mechanism for signal amplification. In other words, pseudopod formation is an excitable process. The growth of the actin filaments creates free ends that promote further actin polymerization establishing the self-enhancing mechanism. Such an extreme sensitivity can be modeled through an appropriate pattern formation mechanism.

Actin filaments are formed in spatially isolated patches which suggest that the self-enhancing reaction occurs in a few localized “hot spots” or “active patches.” Experiments have shown “ordered patterns”

^{13,46,67,68,72}of membrane protrusions in motile cells, further supporting the idea of an intracellular pattern formation mechanism.The growth of the filaments is checked by filament-severing and capping proteins and is also limited due to the finite amount of monomeric actin.

^{7}Thus, the self-enhancing localized reaction must be accompanied by a long-range inhibition reaction.

Such short-range self-enhancement and long-range inhibition reactions can be modeled using Turing-like patterns.^{25,44,45} Further, Turing patterns naturally establish a directional polarity which is another requirement for motility. The polarized pattern remains even after the signal is removed which enable the cells to migrate even in the absence of an external cue as observed in experiments.^{47,72} We therefore use the pattern generation capability of Turing instabilities to generate pseudopods in our model cell, as in many previous studies.^{1,12,13,26,28,32,44,46}

Not all studies have used Turing-like patterns for cell motility. For example, Farutin *et al.*^{3} and Wu *et al.*^{34} considered a deformable cell driven under a prescribed force. As such, the model does not represent the complex pseudopod dynamics. In Ref. 31, the cell was propelled by a specified extension-retraction cycle. It is unclear whether this model could represent the specific pseudopod dynamics as described above. There are also many mechanochemical models that consider major biochemical reactions and couple protein concentrations to membrane deformation. For example, in the work of Yang *et al.*,^{48} a steady-state profile of membrane lipids was used to model a protrusive pressure that propelled the cell towards a chemoattractant. In the work of Liu *et al.*,^{49} protein concentrations are coupled to a Helfrich-like free energy, and the evolution of the membrane is driven by energy minimization. However, these models predicted stationary cell shape, rather than the dynamically changing pseudopods that are the hallmark of Dicty motility. As such, as noted above, Turing models seem to be very suitable to predict the dynamic nature of pseudopods and cell shape.

There are a number of Turing models that can generate biological patterns:^{45,50} the FitzHugh-Nagumo model,^{51} the Schnackenberg model,^{52} the Thomas model,^{53} the Lengyel-Epstein model,^{54} the Degn-Harrison model,^{55} and the Gray-Scott model,^{56} to name a few. In general, these models produce diverse types of patterns such as, stripes, spots of different shapes, and spirals. Again the choice of a specific model is dictated by the characteristics of the pseudopods. In addition to satisfying the three requirements (i-iii) as stated above, the model must also satisfy the following three requirements:

- (iv)
The resulting patterns must generate a few spots (1–4) instead of filling the entire domain, as most of the above models do.

^{1,2,8–10} - (v)
The patterns must be of nearly circular shape rather than stripes, spirals, or spots of other shapes, as the pseudopods are nearly cylindrical protrusions.

- (vi)
The inherent characteristic of many reaction-diffusion systems is to create steady-state patterns.

^{44,45}In contrast, pseudopods are highly dynamic. They continually grow, move along the cell body, bifurcate, and then decay.^{1,2,8–10,67,68}The specific model must create spots that exhibit this dynamic behavior, instead of a steady-state pattern.^{46}

These criteria eliminate many reaction-diffusion models from consideration.^{44,50} Meinhardt^{44} introduced a Turing model that apparently satisfied the aforementioned requirements. This model was used in the past for pseudopod-driven motility.^{28,32} We therefore select this model in our work. A detailed justification for choosing this model can be found in Ref. 44. This model utilizes the behavior of several competing species which we call activators and inhibitors. Because each species diffuses at different rates, dynamic patterns emerge as the activator and inhibitors continuously dominate one another in the system. In this model, the following system of nonlinear reaction-diffusion equations is solved on the evolving cell surface:

where the variables *a*_{1}, *a*_{2}, and *a*_{3} are the concentration of activators, *global* inhibitors, and *local* inhibitors, respectively, the · (dot) represents a time derivative, *S* represents the deforming cell surface, and Δ_{S} is the surface Laplacian. Physically, the activators can be thought as actin filaments, and the inhibitors can be thought as capping and filament severing proteins. An explanation of different parameters can be found in Refs. 44 and 50: *r*_{1} is the production and decay rates of the activators, *r*_{3} is the decay rate of the local inhibitors, *k*_{1} represents a baseline concentration of the activators, *k*_{2} is the production rate of the local inhibitors, *s*_{1} controls the saturation level of the activators so that multiple maxima (and, hence, multiple pseudopods) can coexist, and *s*_{3} is the baseline concentration of the local inhibitors. The important parameters in these reaction-diffusion equations are *D*_{1} and *D*_{3} which represent the surface diffusivity of the activators and local inhibitors, respectively. A random noise *ε*(**x**,*t*) generated from a stochastic differential equation representing the Ornstein-Uhlenbeck process is used to perturb the system into producing instabilities. No explicit external cue, such as a chemoattractant gradient, is present in our model. The cell is only exposed to and reacts solely from random noise.

In pattern formation, the inhibitors must also increase (decrease) as the activator increases (decreases) so that the latter does not grow unbounded. This is also biologically realistic: As more filaments are generated (*a*_{1} increasing), more capping and severing proteins are attached to them to check the unbounded growth of filaments.^{4,7} With more actin (high *a*_{1}) being assembled, the chances of disassembly grow higher and therefore the concentration of the inhibitors *a*_{2} and *a*_{3} also increases. Since *a*_{2} is the global inhibitor, it must change following the net (global) change in *a*_{1}; hence, the integration in Eq. (5) is over the entire cell surface. In contrast, the local inhibitor *a*_{3} changes following the local concentration of *a*_{1}. The term *r*_{1} + *ε* in Eq. (4) is the activator growth rate with some added noise to introduce temporal and spatial variations in the growth rate, and $a12/a2$ ensures that the activator concentration grows fast and nonlinearly but is also checked as the global inhibitor concentration *a*_{2} increases. Additional check on the activator growth is achieved by the terms in the denominator: $1+s1a12$ acts to reduce the growth when the activator concentration becomes too high, while *s*_{3} + *a*_{3} acts to reduce the activator growth when there is a high concentration of local inhibitor *a*_{3}.

The Turing instabilities of the surface reaction-diffusion equations lead to the growth and bifurcation of concentrated regions of activators. Following several prior studies,^{28–34,36} the protrusive force is directly related to the activator concentration as

where the parameter *ξ* represents the force per actin filament. The specific form of the protrusive force chosen here is biologically relevant for pseudopod-driven motility. Pseudopods are extended perpendicular to the cell surface; therefore, we set the protrusive force to be in the outward normal direction. In Ref. 13, the authors used a linear coupling between the activator and actin polymerization force as a generality because the complexities between signaling and polymerization made it impossible to model their dependence exactly. In addition, the authors state that other forms of the protrusive force, including nonlinear relationships, produced similar results. The form of the protrusive force as used here is, therefore, a common form utilized in many studies as noted above.

Indeed, there could be other choices to propel the cell. For example, the protein concentration gradient on the cell surface could be coupled to the membrane force in the same way that the gradient of surfactant results in viscous fingering phenomenon in liquid drops, which is known as the Marangoni effect.^{57} Following this path, protein concentrations are coupled to the membrane Helfrich energy function.^{48,49,58} However, minimization of this energy yields steady-state shapes, rather than the dynamic patterns of pseudopods as discussed before. Further, this formulation requires explicit relations between protein concentrations, bending rigidity, and membrane elasticity, which are often unknown and they introduce additional unknown parameters.^{48,49} Equation (7), on the other hand, introduces only one additional parameter ξ, which as noted above, has a physical meaning. In an alternative approach,^{17,59} the membrane spontaneous curvature was coupled to the surface protein concentration by letting the latter to concentrate at or avoid high curvature areas. While cell deformation can be recreated by this model, it is unknown whether the dynamic pseudopod bifurcation can be recreated. Our model uses the Helfrich energy function for membrane bending, albeit in an explicit form as a bending force density as discussed later in Sec. II D. This force density is directed normal to the membrane. Thus, in principle, one could envision that the protrusive force given by Eq. (7) is a modification to the Helfrich force caused by the surface proteins. Note further that locally altering the membrane shear elasticity cannot generate membrane protrusion because the shear elastic force is tangential to the membrane.

### C. Cytoplasmic and extra-cellular fluids

The cytoplasm and extracellular fluid are assumed to be incompressible and Newtonian. Since inertia is negligible, the fluid motion is governed by the Stokes equations and the incompressibility condition,

where *u* and *p* represent the fluid velocity and pressure fields both in the cytoplasm and in the extracellular fluid. We assume that the densities are the same for both fluids, but the viscosities are different. The viscosity *μ* = *μ*_{0} for the extracellular fluid and *μ* = *λμ*_{0} for the cytoplasm. The parameter *λ* is therefore the ratio of cytoplasmic to extracellular fluid viscosity.

The coupling between the protrusive force, membrane deformation, and fluid motion in the cytoplasm and extracellular fluid is done using the immersed-boundary method.^{60} In this approach, a single set of governing equations is written for both the interior and exterior of the cell, while an indicator function is used to differentiate between the two fluids. The presence of the membrane, which is the interface between the two fluids, is accounted for by introducing a source term in the Stokes equation which includes the forces acting on the cell membrane as

where *S* represents the cell surface. The 3D Dirac delta function *δ* appearing in Eq. (10) is zero everywhere except at the location of the membrane and is used to couple the protrusive and membrane forces to the fluid motion. Note that for a freely suspended swimming cell, the *net* force is zero. Since we model the protrusive force using the explicit form as Eq. (7), this force by itself does not satisfy the force-free condition. However, in the present formulation, it is the *net* force that must be zero. The zero net force condition is satisfied by the force balance implied in Eq. (10) which includes the membrane bending and elastic forces, the protrusive force, as well as the fluid pressure difference and the viscous drag.

### D. Computation of membrane tensions

A finite-element method is used to compute the membrane tension **f**_{E} arising from shearing deformation. The cell surface is discretized using 20 480 Delaunay triangles. Each triangle is assumed to remain flat upon deformation (Fig. 1). The displacement field **v** is assumed to vary linearly within each triangular element and is expressed in terms of linear shape functions *N*_{i} as **v** = *N*_{i}**v**_{i} where the index *i* = 1, 2, 3 denotes the vertices of the triangle. The finite-element shape functions can be evaluated by knowing the coordinates of the vertices and by letting, for example, *N*_{1} = 1 at the vertex 1 and 0 at vertices 2 and 3. Once the shape functions are known, the deformation gradient **F** is obtained. Subsequently, $\epsilon 12$ and $\epsilon 22$ which are the eigenvalues of **F**·**F**^{T} and the stress tensor *τ* are obtained. Having the stress tensor in each element evaluated, we obtain the resultant elastic force at each node as

where **N** is the vector of the shape functions, **P** = *ε*_{1}*ε*_{2}*τ*·**F**^{−T} is the first Piola-Kirchhoff stress tensor, and *S*_{m} is the area of each of the *m* triangles surrounding the node. Additional details can be found in our prior studies.^{61,62}

For numerical computation of the bending resistance on the triangulated cell surface, it is more convenient to use an expression of bending force density that can be derived from the bending energy expression as follows:

where *κ*_{g} is the Gaussian curvature, Δ_{LB} is the Laplace-Beltrami operator, and **n** is the normal vector.^{42} As noted before, the spontaneous curvature is an intrinsic property of the membrane, but it is difficult to measure experimentally. Often, its value is obtained by iterative calculations in mathematical models aimed to generate the experimentally observed resting shapes of cells, such as the biconcave discocyte shape of erythrocytes. Even in this case, different values can be used to predict the same resting shape of erythrocytes depending on the stresses in the cell membrane under the resting condition.^{63} For the specific cell of interest Dicty, there is no information available on the value of *c*_{0}. Therefore, for simplicity, we set it to be zero so that a bilayer patch is planar under the resting condition. The mean and Gaussian curvatures in Eq. (12) are evaluated at each vertex using a quadratic surface fitted to the vertex and its nearest neighboring vertices. Then, using the Gauss theorem, Δ_{LB}*κ* is approximated on a small surface patch *dS* as $1/dS\u222bl\u2207S\kappa \u22c5nldl$, where *l* denotes the patch boundary, $\u25bf$_{S} is the surface gradient, and **n**_{l} is the unit normal to the boundary *l*. The gradient $\u25bf$_{S} on a surface triangle can be obtained either by a linear interpolation of the surface and *κ* or using Loop’s subdivision method. The bending force density **f**_{B} is then coupled to the bulk flow by Eq. (10). Additional details of this computation are given in Ref. 61. In general, the bending forces are small compared with the protrusive forces, as long as pseudopods are not significantly extended. However, as seen in experiments^{1,8,13,47,72} and also in our simulations, long slender pseudopods with high curvature at the tips are generated. As per Eq. (12), the bending force density is proportional to *κ*^{3} and would be significant when pseudopods are extended. Further, the constitutive law used for the membrane shear deformation, Eq. (1), is known to cause buckling, based on our own study and others.^{62} Thus, from a numerical standpoint also the bending force is necessary.

### E. Reaction-diffusion system

For the nonlinear reaction-diffusion equations (4)–(6), we utilize the finite element method given in Ref. 32. The same surface mesh used to solve for membrane deformation is utilized. We describe the method using a generic reaction-diffusion equation written in terms of a species *a* on an evolving surface as

where *D* is the diffusion coefficient and *f*(*a*) is a reaction term. As customary in finite-element theory, Eq. (13) is multiplied by an arbitrary function *ϕ* and integrated over the cell surface to yield

where we utilized the fact that the net flux vanishes on the closed cell boundary. Using the shape function *N*, we express *a* and *ϕ* for each surface element as *a*^{e} = *N*_{i}*a*_{i} and *ϕ*^{e} = *N*_{i}*ϕ*_{i}. Equation (14) then becomes

where the integral is taken over each triangular element, the summation is over all elements, and *S*^{e} is the element area. Since *ϕ* is arbitrary, Eq. (15) yields

for every element, where $Me=\u222bSeNeTNedS$ is the mass matrix, $Ke=\u222bSeBeTDBedS$ is the stiffness matrix, $Fe=\u222bSeNeTf(ae)dS$ is the forcing matrix, and *B*^{e} is the approximation for the surface gradient matrix. Because linear triangular elements are used, the mass and stiffness matrices can be found directly without integration. Equation (16) is calculated over each element and globally assembled leading to the final form as

A forward difference approximation for the time derivative and a semi-implicit method for the remaining terms are used. The above formulation is used to advance the activator *a*_{1} and the local inhibitor *a*_{3} on the deforming cell surface.

### F. Flow solver and immersed-boundary method

The flow domain is a triply periodic cubic domain with lengths nineteen times the cell radius for most simulations. The domain is discretized using a fixed (Eulerian) rectangular mesh of 360^{3} nodes (Fig. 1). The governing equations for the fluid motion are solved on this mesh using a projection method for the time integration and a collocated arrangement for the variables in space. A combination of spectral and finite-difference schemes is used to compute the spatial derivatives. The delta function used to spread the singular membrane and protrusive forces to the surrounding fluid is numerically approximated with a cosine function spanning over four Eulerian points around the cell boundary as

where Δ is the Eulerian grid spacing, *x*_{i} is an Eulerian point, and $xi\u2032$ is a Lagrangian node on the cell surface. Once the fluid velocity is obtained, the membrane velocity **u**_{m} is computed by interpolating the fluid velocity from the surrounding Eulerian nodes as

Then, the membrane nodes **x**_{m} are advected as d**x**_{m}/d*t* = **u**_{m}, resulting in a new location and deformed shape of the cell.

We now describe the means by which the viscosity contrast *λ* between the cytoplasm and extracellular fluid is updated. The immersed-boundary/front-tracking approach permits multiple fluids to be modeled by a single governing equation, while the distinction is made by the use of an indicator function *I*(*x*,*t*) that is unity inside the cell and zero outside. Using this function, the viscosity field *μ*(*x*,*t*) is evolved as the deformable cell changes its shape and acquires a new position as

It can be shown that the indicator function follows a Poisson equation as

### G. Validation

Extensive validation of the cell membrane model has been done in several of our previous studies.^{61–63} It includes numerical experiments, such as cell aspiration in a micropipette, and deformation in shear flow. For brevity, we avoid repeating these studies here and point the readers to the aforementioned references.

Validation of the reaction-diffusion equations is now presented. First we consider pure diffusion of a species *a* on a rigid sphere following $\u2202a/\u2202t$ = Δ_{S}*a*, subject to the initial condition *a*(**x**,*t* = 0) = *z*/*R*, where *R* is the sphere radius, **x** is a location on the sphere surface, and *z* is the *z*-coordinate of the location. We solve for the evolution of *a* using the finite-element method and compare the results against the exact analytical solution *a*(**x**,*t*) = (*z*/*R*)*e*^{−2t} as shown in Fig. 2. The numerical solution agrees well with the analytical solution, as the species diffuses smoothly without any numerical oscillation. We find the *L*_{2} and *L*_{∞} errors to be on the order of 10^{−6} and 10^{−5}, respectively.

Next we validate the full reaction-diffusion equations used to model the activator-inhibitor system as described by Eqs. (4)–(6). The system is able to generate various types of patterns (Turing instabilities) on curved surfaces. In Fig. 3, we present different Turing patterns on a rigid sphere as obtained using our finite-element method. Overall, we observe stable patterns in the form of multiple (2 or more) patches of activator concentration and dynamic patterns in the form of single meandering patches, bifurcating patches, and stripes and ridges. Our method is capable of producing such known patterns of Turing instabilities and matches with multiple previous studies,^{5,6} thereby validating the numerical methodology for the activator-inhibitor system.

### H. Model parameters

The initial undeformed shape of the cell is assumed to be spherical. The simulation results are presented using dimensionless variables. Lengths are scaled by the cell radius *R*, time is scaled by *R*^{2}/*D*_{1}, and velocities are scaled by *D*_{1}/*R*. The major dimensionless parameters of interest are the ratio of protrusive force to membrane elastic force *α* = *ξ*/*RG*_{S}, the ratio of the cytoplasmic to extracellular fluid viscosity *λ*, and the ratio of inhibitor to activator diffusivities *β* = *D*_{3}/*D*_{1}. The force per actin filament is in the range 3–8 pN,^{64} cell radius *R* ∼ 10–50 *μ*m,^{1} *D*_{1} and *D*_{3} ∼ 1 *μ*m^{2}/s,^{64} and membrane shear modulus *G*_{S} ∼ 10^{−6} N/m.^{65} Membrane stiffness is altered in malignant and drug-treated cells.^{14–16,66} Our interest is in the role of membrane elasticity, cytoplasmic viscosity, and protein diffusivity on cellular motility. We therefore consider the range *α* = 1–10 and *λ* = 0.1–10. Note that cell deformation is controlled by both *α* and *λ*. Larger values of *α* and smaller values of *λ* result in “softer” cells, and vice versa. For the reaction-diffusion system Eqs. (4)–(6), we use dimensionless values as *r*_{1} = 100, *r*_{3} = 65, *k*_{1} = 500, *k*_{2} = 25, *s*_{1} = 10^{−4}, and *s*_{2} = 0.2, following Ref. 32.

The range of diffusivity ratio *β* is chosen such that a Turing pattern relevant for the swimming amoeba is generated. The key to pseudopod-based swimming is the bifurcation of an existing pseudopod into two “offspring” pseudopods. Thus the Turing pattern that needs to be generated on the surface of the cell should be a dynamic bifurcating pattern as shown in Fig. 3(c). In this dynamic pattern, an existing activator patch bifurcates into two offspring patches. After sometime, one of the patches decays while the other survives. Subsequently, the surviving patch also bifurcates, thereby repeating the cycle. Also note that the patches usually move over the sphere surface. Our numerical experiments suggest that a bifurcating instability can be generated in the range *β* = 1.5–3.5; thus this range is chosen to generate pseudopods.

## III. RESULTS

### A. General swimming behavior

Time-progression images of a swimming amoeba as obtained in our simulations are shown in Fig. 4 (Multimedia view). Similarity is observed between the predicted cell shapes and experimentally observed shapes in Refs. 1, 8, and 13, although the latter provide 2D images only. The sequence shows *de novo* pseudopod generation at the beginning, followed by pseudopod bifurcation, and cell direction changes caused by pseudopod meandering. Consider the swimming sequence shown in Fig. 4 (top) (Multimedia view). Starting from an initial spherical cell shape subject to random noise for a sufficient time, a Turing instability develops, leading to a region of high activator concentration. As the activator is coupled directly to the protrusive force in our model, the local cell membrane begins to protrude outward resulting in a *de novo* pseudopod [a in Fig. 4 (top) (Multimedia view)]. Soon after, a bifurcation event takes place, resulting in two distinct instabilities (b). The activator concentration at these new instabilities grows resulting in two distinct pseudopods diverging from one another (c). The pseudopods meander about the origin and propel the cell forward (d). (e) and (f) show the bifurcation of an existing pseudopod and decay of another. In (g) and (h), the cell changes direction due to the growth and backward movement of a pseudopod along its side. (h)–(o) show a persistent, nearly unidirectional forward movement of the cell, followed by another direction change in (o)–(q). Significant 3D movement of the pseudopods is predicted by our simulations as evident in the figure. These novel simulations also predict significant deformation of the whole cell due to dynamic growth and decay events of the pseudopods.

Cell shape is observed to strongly depend on membrane stiffness, protrusive force, protein diffusivity, and fluid viscosity. Figure 5 (Multimedia view) [and Fig. 4 (Multimedia view) also] presents a comparison of cell shapes for varying parameters. For relatively larger values of *α* (typically 3 and higher), which correspond to reduced membrane stiffness, cells assume elongated shapes as already shown in Fig. 4 (Multimedia view). Similar elongated shapes are also observed at higher values of diffusivity ratio *β* (usually 2.5 and higher) and smaller values of viscosity ratio *λ* (<1) [Fig. 5(b) (Multimedia view)]. More rounded shapes are observed for smaller values of *α* and *β*, but higher values of *λ* [Figs. 5(a) and 5(c) (Multimedia view)]. Pseudopods are observed to be slender and longer for softer cells at a larger diffusivity ratio, while they are shorter and thicker for stiffer cells at a smaller diffusivity ratio. For smaller *α*, protrusive forces are not strong enough to significantly overcome the membrane tension resulting in nearly rounded shapes. In the same vein, for larger *λ*, dissipation in the cytoplasm acts against protrusive forces, which, in turn, results in more rounded shapes. One interesting observation at large *λ* is the ridges forming over the cells which are the “footprints” of retracted pseudopods [Fig. 5(c) (Multimedia view)]. It is because of the increase in shape, recovery time at larger *λ*, which in turn, has resulted in slower retraction of the pseudopods, even though the activator concentration has already decayed.

### B. Instantaneous velocity and cell trajectory

In the absence of external cues, as is the case in our simulations, the pseudopod-driven swimming of an amoeba should resemble a random walk. The Turing instabilities on the cell surface used to generate the pseudopods bifurcate and meander around the cell facilitating random generation of pseudopods. Indeed, random swimming is observed in our simulations for a certain parameter range. However, we find that cell deformation and inhibitor diffusion alter the nature of this randomness and may lead to persistent behavior in swimming. This can be observed from sample cell trajectories that are shown in Fig. 6(a). They show that cells with smaller values of *β* and larger values of *λ* travel much shorter distances and exhibit frequent and sharp turns. As noted earlier in Fig. 4 (Multimedia view), turning events occur either due to pseudopod bifurcation along the sides of the cell or due to pseudopod meandering over the cell. In contrast, cells with larger values of *α* and *β* and smaller values of *λ* swim longer distance without changing their directions. These cells also exhibit much slower and less frequent turns, resulting in nearly persistent unidirectional swimming.

The instantaneous velocity *V*_{C} of the cell center is shown in Fig. 6(b) for several cases. The random nature of cell swimming is evident here from the fluctuations of *V*_{C}. Fluctuations increase with increasing *α* but decreasing *β*. Note that fluctuations in *V*_{C} can arise due to change in either overall cell velocity or pseudopod dynamics. In the time series plots, pronounced maxima and minima indicate large changes in the direction, while smaller maxima and minima are likely due to the pseudopod dynamics. As evident from the figure, cells at larger *α* and *β* exhibit large maxima and minima less frequently. This implies that these cells maintain a nearly persistent unidirectional motion over a larger distance. Similar results are obtained at smaller values of *λ*.

The question whether the individual cell motion represents a diffusive process is often an exclusive subject in many articles. In particular, Refs. 47 and 72 showed that at short times, the motion is ballistic due to the persistence in cell motion, while at longer times, it is diffusive. We have computed mean-squared displacement (MSD) for a few cell trajectories based on the location of the cell centroid (see supplementary material). The MSDs are nearly linear at larger times, thus suggesting a diffusive type motion, in agreement with previous findings. We caution, however, that our MSD calculations are based on a relatively small number of simulations. To get a clear behavior, MSDs must be calculated using many simulations for each cell type. Given the 3D nature of the simulations and their complexity, running multiple simulations for the same parameter values is computationally expensive and beyond the scope of the present work.

### C. Flow field

Figure 7 presents instantaneous fluid velocity vectors for a few representative cases. Our simulations reveal transient fluid velocity fields both in the cytoplasm and in the external fluid. Inside the cells, fluid streaming from the cell body towards the growing pseudopods is observed. A growing pseudopod is also observed to draw the surrounding fluid with it. In contrast, fluid flows from a receding pseudopod towards the cell body. Higher velocity is observed near the growing pseudopods, but smaller velocity is observed near the receding end. In the exterior, the flow pattern is complex and highly unsteady. In particular, vortex cores are observed to shed during the growth of a pseudopod and during rapid turns.

Two paradigms, pushers and pullers, have been identified in the literature for many swimming micro-organisms. Pushers are described by fluid repulsion along the directional axis of swimming and fluid attraction along the sides of the cells. Pullers are the opposite; they draw the fluid towards the cell in the swimming axis and push the fluid outwards along the sides. In our simulations, pusher motion is rare and can only be seen during the initial transience when the spherical cell is beginning to deform under the first Turing instabilities [Fig. 8(a)]. Pusher motion also requires that two activator patches initially grow in opposite directions. Eventually, one of the patches decays while the other one survives and in time exhibits tip splitting. Puller motion, on the other hand, is seen more often in our simulations [Figs. 8(b) and 8(c)]. When pseudopods grow along the sides of a cell, fluid repulsion is observed along the sides and fluid flows towards the cell in the swimming direction. Puller type motion is observed across all values of *α* for smaller values of *β*, when cells are less elongated. In general, however, our simulations show that the cell behavior cannot be uniquely classified as a pusher or puller [Figs. 8(d)–8(f)]. This is particularly due to the focusing of pseudopods observed while using higher values of *α* and *β* as described before. In such cases, fluid repulsion along the sides cannot occur as the pseudopods are generated in the front only.

Note that the conclusion about pullers or pushers is *not* based on just one instant of flow vectors. We have carefully analyzed the flow field over the entire simulation time for each case and made the conclusion. A quantitative way to characterize a puller-type or pusher-type motion is to compute the stresslet tensor **S** as derived in Ref. 69 as

For a highly complex cell shape with dynamically extending, bifurcating, and retracting pseudopods as considered in our study, the tensor **S** changes continually, has all non-zero components, and is not easy to interpret. In contrast, for an axisymmetric swimmer moving along, say the direction *e*_{z}, Eq. (22) can be simplified in terms of a single scalar *Ψ* as

which lends to a clear physical interpretation, as often is the case for simple shaped swimmers.^{70} In our simulations, *nearly* axisymmetric shapes may occur rarely during the initial transience. Thus, we can compute *Ψ* for only a few cases and interpret it to make any conclusion about pullers or pushers. Specifically for the pusher shown in Fig. 8(a), *Ψ* is negative (−15, 500), and for the puller in Fig. 8(b) (initial axisymmetric phase only), it is positive (190). This trend of *Ψ* (positive for a puller and negative for a pusher) is in agreement with Ref. 70. However, we re-emphasize that the Dicty shapes as predicted in our simulations [Fig. 4 (Multimedia view), for example] are highly complex, and some information about a puller or pusher, at best qualitatively, can be obtained for these complex shapes using the flow vectors as done in Ref. 33.

### D. Pseudopod dynamics

As evident from Fig. 4 (Multimedia view) and Fig. 5 (Multimedia view), the transient cell shape is dictated by the pseudopod dynamics, which in turn are strongly dependent on membrane stiffness, fluid viscosity, and protein diffusivity. One interesting observation in our simulations [Fig. 4 (Multimedia view) and Fig. 5(b) (Multimedia view)] is that the pseudopods are more directed and confined near the anterior of the cell for larger values of *α* and *β* and for smaller values of *λ*. In contrast, for smaller values of *α* and *β* and for larger values of *λ*, pseudopods appear all over the cell without any directional preference. The directed orientation of pseudopods at larger values of *α* and smaller values of *λ* is one reason why the cell shape becomes elongated. Also noteworthy is the appearance of multiple pseudopods, 3-4 in general, which are produced at larger values of *α* and *β* and for smaller values of *λ*, while a maximum of two pseudopods are observed for stiffer cells and at reduced inhibitor diffusivity.

The directionality of pseudopods is quantified by computing the angular position of pseudopod tips using the spherical coordinates *θ* and *ϕ* as shown in Fig. 9. The role of the three parameters (*α*, *β*, and *λ*) on the pseudopod directionality is evident here. For relatively stiff cells as in the case of smaller values of *α* [Fig. 9(a)] and larger values of *λ* [Fig. 9(d)], a uniformly scattered distribution is noted which indicates that pseudopods are generated around the entire cell in all directions. A scattered distribution is also observed for smaller values of *β* [Fig. 9(a)]. In contrast, for softer cells as in the case of larger *α* [Fig. 9(c)] and smaller *λ* (<1) [Fig. 9(e)], a tight grouping of data points in a small range of *θ* is observed, which indicates that pseudopods are focused and generated around the same direction. Cells with a larger diffusivity ratio also show tight groupings of data points and, hence, a persistent focused directionality of pseudopods.

A closer inspection of the instantaneous cell shapes, as shown in Fig. 9(f), illustrates a mechanism for the focused directionality of pseudopods. After one bifurcation, the resulting two pseudopods initially move away from each other. Later, they move back again towards the front of the cell before one of them retracts. Thus, the next bifurcation occurs in nearly the same location as the preceding one, and hence the same direction of cell motion is maintained.

Our simulations also reveal three-dimensionality in pseudopod generation: We observe that a new pair of pseudopods often migrates in a direction normal to the previous bifurcation path [Fig. 9(g)]. This characteristic is not inferred from experimental images which are two-dimensional. Essentially, this is a three-dimensional extension of the dynamics seen in Fig. 9(f), which occurs repeatedly until its rhythm is broken by a third errant pseudopod.

Further quantification of pseudopod dynamics is presented in Fig. 10 where the number of active pseudopods and their lifetime are shown. For the parameter range considered, the average number of pseudopods ranges from about 1 to 2.5, and the maximum number of pseudopods observed at any time ranges from 2 to 4. Both the average and maximum numbers show similar trends with respect to *α*, *β*, and *λ*. The number of pseudopods increases as the cells become more deformable (i.e., *α* increases or *λ* decreases) or as the diffusivity ratio increases. For relatively stiffer cells, a maximum of two pseudopods are observed at any time instant, while for more deformable cells, three to four pseudopods are observed. Likewise, cells with smaller values of the diffusivity ratio exhibit at most two pseudopods at any time, while those with a larger diffusivity ratio are characterized by the presence of three or, in few instances, four pseudopods.

Figure 10 shows that the average lifetime of a pseudopod is strongly dependent on *α* and *β* but not on *λ*. The pseudopod lifetime decreases with increasing *α*. This also means that the frequency of tip-splitting, which occurs through bifurcation of a Turing instability, increases with increasing *α*. We attribute this to the greater tendency of pseudopods to bifurcate when exposed to a higher degree of deformation. Turing instabilities become more severe with increasing deformation. Then, an existing pseudopod can bifurcate before the other retracts, resulting in multiple pseudopods existing simultaneously. In addition, cells exhibiting multiple bifurcating pseudopods at any given instant in time will result in a higher bifurcation frequency as well.

Pseudopod lifetimes are observed to increase with increasing *β* (Fig. 10). This is because at larger *β*, the activator concentration increases at a faster rate than the inhibitor concentration. Therefore, the pseudopods can survive longer.

### E. Swimming speed

Pseudopod dynamics, which are directly linked to the mechanical properties of cells and to protein diffusivity as shown above, strongly influence the cell’s speed and trajectory. Figure 11 shows the time-averaged cell speed $V\xafC$ (scaled by *D*_{1}/*R*) as a function of *α* for two values of diffusivity ratio. The cell speed is observed to increase with increasing *α*, which is expected due to an increase in the protrusive force. Earlier, Figs. 9 and 10 suggested, however, that the change in pseudopod dynamics which occurs as *α* increases from smaller to higher values also contributes to the increase in cell speed. Figure 9 has shown that the pseudopods are always confined near the front of the cells at higher *α*. As a result, the protrusive forces from different pseudopods act in a similar direction, increasing the cell speed. The presence of multiple (3–4) pseudopods at higher *α* as shown in Fig. 10(a) also results in higher protrusive forces and hence a higher cell speed. At smaller values of *α*, on the other hand, pseudopods are generated all around the cell as noted in Fig. 4 (Multimedia view), and at any given time, at most two pseudopods co-exist as noted in Fig. 10, resulting in a reduced cell speed.

Interestingly, Fig. 11 shows that the slope of $V\xafC$ versus *α* curves increases significantly with increasing diffusivity ratio. Therefore, cell deformation has a greater influence on cell speed at higher inhibitor diffusivity. Similarly, inhibitor diffusivity has a greater influence as the cell becomes more deformable. This result implies a strong coupling between cell deformation and surface protein diffusion, which affects the pseudopod dynamics and in time the cell speed. The effect of the diffusivity ratio *β* on $V\xafC$ is plotted in Fig. 12(a). Swimming speeds are relatively constant for *β* < 2.25, but they start increasing dramatically as *β* is further increased. This trend is again entirely due to the alteration in the pseudopod dynamics. As noted previously in Figs. 9 and 10, pseudopods at higher *β* become confined near the front of the cell. The protrusive force vectors then add to increase the cell speed. Furthermore, Fig. 10 showed that at higher *β*, multiple pseudopods appear and their lifespan increases. These two factors also lead to a higher cell speed.

The influence of cytoplasmic to external fluid viscosity ratio *λ* on the time-averaged cell speed is shown in Fig. 12(b). Two distinct trends, separated by a minimum, are observed here. For *λ* < 1, the swimming speed and viscosity ratio are inversely correlated: A decrease in $V\xafC$ with increasing *λ* is observed in this range. After the minimum is reached at *λ* ≈ 3, $V\xafC$ increases with increasing *λ*. This trend is also directly related to the pseudopod dynamics. In the range *λ* < 1, multiple pseudopods were observed to be focused near the front of the cells (Figs. 9 and 10) resulting in persistent unidirectional motion and higher cell speeds. Simulations with *λ* > 1, on the other hand, exhibit pseudopods randomly generated around the entire cell, resulting in a loss in the directionality of cell motion and a reduced cell speed. For *λ* > 3, however, the pseudopod lifespan increases [Fig. 10(e)]; as a result, the cells are able to maintain their direction for a longer time and achieve higher speeds.

## IV. DISCUSSION

We presented a computational modeling study of a swimming amoeba freely suspended in a fluid. The cell swims without any adhesion or any external cues. The model is based on the immersed boundary/front tracking methods, and it combines cell deformation with coarse-grain biochemical reactions and fluid motion in the cytoplasm and extracellular fluid. Because of the coupling of fluid flow with the rest of the model, the effect of fluid drag on the cell motion is directly resolved without the need for an *ad hoc* modeling as done in some earlier studies.^{12,13,25,26,28,32} Furthermore, the nonlinear membrane model used here is a significant improvement over the Hookean models that have been used in prior models. These latter models often experience unwanted mesh distortion and require re-meshing of the cell surface or other *ad hoc* corrections. For example, in Ref. 32, a tangential velocity along the membrane was introduced to avoid surface mesh distortion. Our method does not suffer from such numerical issues and is therefore robust. Our simulations can be continued for a long time even for cells that generate very complex shapes with multiple elongated pseudopods. As evident from the cell shapes, the model predicts extremely deformed and complex cell shapes without any numerical instability, suggesting the robustness of the method. The predicted cell shapes are similar to those observed in experiments investigating the motility of *Dicty*.^{1,8–10,13,67,68,72}

We used the dynamic nature of Turing instabilities to evolve membrane protein concentration, which is then correlated to the protrusive force. Experimental studies have shown a large density of proteins like Ras or F-actin clustered at the site of a pseudopod.^{25,71} We model this effect using Turing instabilities, which result in the generation, migration, and bifurcation of activator patches on the cell surface. Consequently, pseudopods are shown to bifurcate in the front, move backwards, and eventually retract. The pseudopod dynamics result in a paddle-like non-reciprocal motion of the cell, which is in agreement with Purcell’s scallop theorem.^{11} Furthermore, pseudopod bifurcation, often referred to as “tip-splitting,” has been perceived as the key event in amoeboid swimming motility.^{71,72} Indeed, except for *de novo* pseudopods seen during the initial transience, our model generates pseudopods exclusively from splitting of a pre-existing pseudopod.

Our numerical data are quantitatively comparable to experimental measurements. It may be noted that experimental data on freely suspended swimming amoeba is relatively scarce as most studies have considered cells crawling on substrates. Reference 10 reported data on both crawling and swimming amoeba. In converting the numerical data to dimensional values, we use a cell size of 20 *μ*m and the activator diffusivity of 1 *μ*m^{2}/s. Then, a typical swimming speed of $V\xafC\u2248$ 1–3 as obtained in our simulations (Figs. 11 and 12) for moderate values of the parameters yields a dimensional speed of 3–9 *μ*m/min which is in agreement with the experimental measurements of 3 *μ*m/min for swimming cells and 7–12 *μ*m/min for crawling cells.^{10,47} Additionally, the simulations yield an average pseudopod tip velocity 15–21 *μ*m/min which is also in excellent agreement with the experimental measurement of about 18 *μ*m/min.^{10}

A prediction from our simulations is the persistent unidirectional swimming of cells observed even in the absence of external cues. The cells in our simulations are only exposed to and react solely from random noise introduced in the reaction-diffusion model. In the absence of an external cue, it is expected that an amoeba would migrate in a random manner because of the random nature of pseudopod generation. Indeed our simulations predict such random migration. In addition, our simulations predict that a persistent unidirectional swimming may also occur in the absence of an external cue. The persistent unidirectional motion of amoeboid cells in the absence of any external chemoattractant gradient was observed in experiments.^{10,47,72} In Refs. 10 and 72, it was found that cells moved persistently because pseudopods maintained the same direction. In our simulations, we show that such persistent motion is a direct result of pseudopod dynamics, which in turn depend on membrane elasticity, cytoplasmic viscosity, and protein diffusivity. Pseudopods were observed to exhibit cycles of migration away from and then towards the point of bifurcation. These characteristics were also observed in 3D where pseudopod paths were nearly perpendicular to one another. We show that pseudopods have a definite polarity for reduced membrane elasticity and cytoplasmic viscosity, but higher inhibitor diffusivity. We further show that the persistent motion is not only because of a focused orientation of the pseudopods but also because of an increase in their numbers. This alteration in pseudopod behavior is shown to occur also with decreasing membrane elasticity and cytoplasmic viscosity, but increasing inhibitor diffusivity.

The simulations allow us to study cell swimming by independently varying membrane elasticity, cytoplasmic viscosity, and protein diffusivity. We find that softer cells swim faster. This observation also directly relates to the experimental finding that metastatic cells are significantly softer than benign cells and further corroborates the proposition that cell stiffness is a potential biomarker. In fact, nearly an order of magnitude variation of swimming speed is found in our simulations for the parameter range considered. We further show that the increased swimming speed is a result of focused orientation of pseudopods and persistent unidirectionality of cell motion.

Another finding is the coupling between membrane elasticity and protein diffusivity. Swimming speed is nearly independent of inhibitor diffusivity *β* for stiffer cells, but it is strongly dependent on *β* as the cell becomes more deformable. The swimming speed is shown to increase with decreasing membrane elasticity at a faster rate as the inhibitor diffusivity is increased. In fact, pseudopod dynamics and motility are observed to be highly sensitive to inhibitor diffusivity for softer cells. Not only do pseudopods become more focused, their number and lifespan also increase with increasing inhibitor diffusivity.

Several experimental studies have used wild-type and mutant cells to characterize amoeboid movement. Reference 71 found that wild-type amoeba assumed an elongated shape and migrated with higher speeds (∼22 *μ*m/min), while mutant cells maintained a nearly spherical shape and moved at much slower speeds (∼3 *μ*m/min). Our simulation results agree with this observation and relate it to membrane stiffness, cytoplasmic viscosity, and protein diffusivity. We found that the cell shape becomes rounder and the cell speed decreases with increasing membrane stiffness and cytoplasmic viscosity, but decreasing inhibitor diffusivity. The cell speed in our simulations changes from about 1.5 *μ*m/min to 18 *μ*m/min as *α* varies from 1 to 10 (for *β* = 3). The predicted variation in cell speed is therefore in agreement with the experimental data for the wild-type and mutant cells. Time series of instantaneous cell velocity also show similarities with those in Ref. 71. Fluctuations in cell velocity are observed for all cells as swimming directions change, but for more deformable cells using a higher inhibitor diffusivity, these changes occur due to a form of continuous turning where pseudopods meander over the cell surface and is in agreement with the experimental study in Ref. 72. At higher *α*, cells display an elongated shape with longer pseudopods as observed in several experiments.^{1,10,13,71,72} For *α* < 1, cells barely move since the protrusive force is not strong enough to deform the membrane.

Inclusion of the cytoplasm and extracellular fluid in the present model enables us to address the role of fluid viscosity on cell motility. The ratio of cytoplasmic to extracellular fluid viscosity is an important parameter that controls cell deformation: Cells become stiffer with increasing viscosity ratio. Our simulations show a novel result that the swimming speed varies non-monotonically with the changing viscosity ratio, with the minimum speed occurring when the cytoplasm viscosity is close to the extracellular viscosity. This trend is shown to be a direct consequence of changing pseudopod dynamics in terms of their directional orientation, numbers, and lifespan. Furthermore, swimming speed increases significantly as the cytoplasm viscosity is lowered below the extracellular fluid viscosity. This may seem surprising at first since a less viscous cell is able to penetrate through a more viscous surrounding. However, this is similar to the so-called viscous fingering phenomenon observed in interfacial fluid flows, in which a less viscous fluid displaces a more viscous one forming dynamic patterns at the interface.

Inclusion of fluid interaction also enables us to extract instantaneous 3D velocity fields both inside and outside a swimming amoeba and address whether a pusher-type or puller-type motion occurs. References 33 and 34 considered vesicles as a model of deformable cells and spherical harmonics for the protrusive force. Out of the different strokes in each displacement cycle, they found that one showed puller behavior and another showed pusher behavior, while the remaining strokes showed a hybrid state which encompassed aspects of both pullers and pushers. So in one cycle, both behaviors existed. In our full simulations where pseudopods are auto-generated, we find that the pusher motion exists only during initial transience while the puller motion may occur thereafter. Predominantly, however, the motility is more complex and often shows both hybrid modes and motion lacking any classification.

## V. CONCLUSION

In conclusion, we presented a 3D modeling study of swimming of amoeboid cells using immersed boundary/front tracking methods. Our model combines large deformation of the cell with protein reaction-diffusion and cytoplasmic and extracellular fluid motion. We predict that cell swimming changes from random-like to persistent unidirectional motion, and the swimming speed increases with an increase in cell deformability and inhibitor protein diffusivity. This persistent swimming is observed without any external cues, and we show that this change is a direct result of the changes in pseudopod dynamics, which in turn depend on membrane stiffness, protein diffusivity, and fluid viscosity. The influence of these parameters on pseudopod dynamics is quantified in terms of pseudopod directionality, quantity, and lifespan. We find a non-monotonic behavior of the swimming speed with respect to the cytoplasmic to extracellular fluid viscosity ratio. We further find that the speed increases significantly as the cytoplasm becomes less viscous compared with the extracellular fluid, resembling the viscous fingering phenomenon observed in interfacial flows. While these results support the experimental observations that softer cells migrate more aggressively, they also suggest a strong coupling between cell deformation and membrane protein diffusion in dictating cell motility.

## SUPPLEMENTARY MATERIAL

See supplementary material for a plot of computed mean-squared displacements.

## ACKNOWLEDGMENTS

The work is supported by Grant No. CBET 1438255 from the National Science Foundation and Peter B. Cherasia fund from Rutgers. Computational resources from NSF’s XSEDE resources at TACC are acknowledged.