Many cellular processes involve substantial shape changes. Traditional simulations of these cell shape changes require that grids and boundaries be moved as the cell's shape evolves. Here we demonstrate that accurate cell shape changes can be recreated using level set methods (LSM), in which the cellular shape is defined implicitly, thereby eschewing the need for updating boundaries.
We obtain a viscoelastic model of Dictyostelium cells using micropipette aspiration and show how this viscoelastic model can be incorporated into LSM simulations to recreate the observed protrusion of cells into the micropipette faithfully. We also demonstrate the use of our techniques by simulating the cell shape changes elicited by the chemotactic response to an external chemoattractant gradient.
Our results provide a simple but effective means of incorporating cellular deformations into mathematical simulations of cell signaling. Such methods will be useful for simulating important cellular events such as chemotaxis and cytokinesis.
Many cellular processes are characterized by substantial shape changes. For example, chemotaxing cells become polarized, assuming a highly elongated form, and crawl across solid substrates in the direction of increasing concentrations of chemoattractant . During cytokinesis, a single cell undergoes significant cytoskeletal deformation, reforming into two daughter cells . These cellular processes are fundamentally mechanical, utilizing force generation at the molecular scale to generate shape changes. Properly simulating cellular shape change requires that we have a description of the underlying mechanical properties of the cell.
To understand fully the mechanisms that regulate these cell shape changes requires knowledge of the signaling pathways as well as their effect on the mechanical properties of cells. For example, a complete model of chemotaxis would require a description of the gradient sensing capability of cells together with a physical model for the cellular migration . Few such models exist, even though it is now appreciated that the response of cell-signaling pathways can be regulated in response to alterations in cell size and shape . The traditional method of simulating cellular deformations is by specifying the boundary of the cell explicitly through a finite-element model (FEM) [5-7]. One problem is that simulation of biological shape deformations – which invariably involves solving partial differential equations on moving boundaries – can be computationally expensive particularly when the cellular deformations are not small. During many processes including cytokinesis and chemotaxis, cellular shape deformations tend to be large and occur rapidly. Here, we demonstrate how the Level Set Method (LSM) can be used to couple mechanical models of the cell with biochemical models of signaling pathways to simulate large cellular deformations.
We briefly contrast the LSM approach to other methods that have been used to account for cellular deformations.
The immersed boundary method (IBM), introduced by Peskin  was developed to simulate the interaction of flexible tissues with the surrounding incompressible fluid. It has been used to simulate cell shape changes during motility . In the IBM, the Navier-Stokes equation describing the fluid flow can be solved on a fixed grid, simplifying this computationally expensive step. The membrane and cytoskeleton is discretized by assigning a series of nodes that are connected by viscoelastic elements. As the cell deforms, nodes and their corresponding links have to be inserted or deleted. This book-keeping comes at a considerable computational complexity. For this reason, the IBM may best be used in situations where the cell shape does not change considerable .
More recently, the cellular Potts model (CPM) has become a popular vehicle to simulate cell shape changes . In the CPM, a cell is described by a connected domain of pixels on a regular grid. The shape of the cell is evolved by updating each pixel based on a set of probabilistic rules. This method does not use an explicit viscoelastic description of the cell. Instead, cell shape is constrained by minimizing an energy function that penalizes size-deformations as well as membrane bending. Cellular Potts models have been used to simulate two-dimensional (2-D) models of cell motility in fish keratocytes  and amoebae . Unlike FEM or IBM, modeling large changes in the shape of the cell is no more computationally expensive than small changes. One drawback, however, is that the mechanical description of cells in the CPM framework is not as tightly integrated with experimentally-based measurements as the method presented here.
Models of cellular shape changes have all been derived based on explicit descriptions of the cell morphology that are updated based on the simulated behavior of the underlying cytoskeleton. For example, Rubinstein et al. provide a detailed 2-D computational model of the lamellipodium keratocyte motility . In this model, the cellular domain is updated at each time step based on the protrusive and retractive forces (actin polymerization and acto-myosin contraction) and re-gridded. This avoids the necessity for nodes and keeping track of the mechanical state of the system. However, the model relies on an elastic (rather than viscoelastic) network which may be appropriate for the thin keratocyte, but is not likely to be applicable to thicker cells.
The rest of the paper is organized as follows. We first provide some necessary background. We then develop a mechanical model within the LSM that accounts for the viscoelastic nature of the cell. We fit this model to experimental data obtained through micropipette aspiration experiments on Dictyostelium cells. We incorporate this viscoelastic model into a level set framework and illustrate how large-scale shape deformations can be accounted by the model. This is done through simulations by showing that the model accurately captures the behavior of the aspirated cell. Finally, using a simple gradient sensing model to generate internal force profiles, we simulate the changing morphology of a cell chemotaxing in response to an externally applied chemoattractant gradient. Using the framework developed here, we obtain the force profiles needed to achieve stable migrating cell morphologies observed for several strains. The methods developed here allow us to link forces acting on the cell and mechanical properties of the cytoskeleton to cell shape deformation explicitly, and will prove useful in studying cellular processes undergoing large-scale shape changes.
Cells derive their mechanical properties from actin, actin-associated proteins, and motor proteins such as myosin-II , which are components of the cytoskeleton. Though distributed throughout the cell, the actin cytoskeleton is concentrated along the periphery of the cell underneath the membrane, particularly in Dictyostelium, and is the molecular machinery that generates cellular shape changes during cell division and chemotaxis.
Cytoskeletal networks exhibit viscoelastic behavior, having both viscous and elastic properties [15-17]. Actin filaments alone do not create significant mechanical resistance; instead, cross-linking of actin filaments by various actin binding proteins imparts mechanical rigidity to the cell. Under applied load, cross-linked actin networks behave similarly to an elastic solid and can be described using Hooke's law. However, because cross-linking proteins bind to and dissociate from actin filaments, actin-based networks may also exhibit viscous flow. Myosin-II, in filament form, also binds to actin filaments and provides mechanical resistance of the cell, as well as influencing the binding kinetics of various actin crosslinkers [2,18]. The interior of the cell also contains cytoskeletal polymers, as well as organelles, a nucleus, and cytoplasmic fluid. Thus, owing to their viscoelastic nature, cells exhibit a time-dependent deformation in response to mechanical force.
Introduction to level set methods
Cell motion has been traditionally simulated by: discretizing the cell boundary, computing the displacement of each of the points according to the local velocity, and forming a new boundary with the displaced points (Fig. 1A). This method may run into difficulties when the spatial or temporal resolution of the simulations is not sufficiently fine, or when changes in topology occur (Fig. 1B). The Level Set Method (LSM) can be used to overcome these difficulties . LSM is a numerical technique for tracking interfaces and shapes which has been widely used in various fields including computer graphics , image processing , computational fluid dynamics  and material science .
Figure 1. Introduction to level set methods. A. The traditional method of tracking moving boundaries involves discretization of the boundary (dotted red) into a set of points, moving each point x = (x, y) according to the local velocity (v(x, t)), leading to a new boundary at the new locations (solid red). B. Difficulties can arise, however, when the geometry of the boundary becomes irregular. In this case, the point tracking method often fails to preserve the boundary topology. Special attention is required to resolve these errors, increasing computational costs. C. In the Level Set Method (LSM), the boundary Γ(t) is embedded into a higher dimensional potential function (ϕ(x, t)) as the zero-contour. Γ(t)moves as ϕ(x, t)evolves in time. D. Because the boundary is defined implicitly, the LSM framework overcomes some of the difficulties of boundary point tracking. E. This example illustrates how an arbitrary cell shape (black contour) can be embedded into a signed distance function to form the level set potential function ϕ(x, t). In this case, the potential function is given by the Euclidean distance to the cell boundary, with positive (resp. negative) sign when outside (resp. inside) the cell.
Suppose that the cell boundary at time t is described by the closed-contour Γ(t). The LSM requires a potential function (Fig. 1C.), denoted as ϕ(x, t), that is related to Γ(t) according to:
Thus, Γ(t) is the zero-level set of ϕ(x, t). It follows that, in the LSM, the cell membrane is represented implicitly through the potential function which is defined on a fixed Cartesian grid, thus eliminating the need to parameterize the boundary. This allows the LSM to handle complex boundary geometries efficiently (Fig. 1D).
One candidate for the potential function is the signed distance function , defined by:
where S identifies the area occupied by the cell and d(x, Γ) is the distance of position x to the curve Γ; see Fig. 1E for an example of a cell shape embedded in a potential function derived from the signed distance function.
We now manipulate Γ(t) implicitly through the function ϕ(x, t) according to the equation:
The vector v(x, t) is the velocity of the level set moving in the outward normal direction. In our case, v(x, t) intrinsically describes the cell's membrane protrusion and retraction velocities. These velocities can be driven by externally applied forces on the cell membrane (e.g. from a micropipette aspirator), or internally generated mechanical forces (e.g. actin polymerization or myosin-II retraction), or both. To determine how these forces translate to membrane velocity, however, first requires a mechanical model of the cell.
As the potential function corresponding to the cell shape is evolved, it can become quite steep or flat. To reduce the numerical errors caused by these effects, we re-initialize the potential function periodically . This can be done using the re-initialization equation :
where S(ϕ(x, 0)) is taken as +1 inside the cell, -1 outside the cell and zero on the cell membrane.
Results and Discussion
Viscoelastic model of cell deformation
The LSM relies on a continuum description of the material properties of the cell [2,26]. We use mechanical models to describe the viscoelastic behavior of the cell . Our mechanical model is based on a representation of cells that assumes a viscoelastic cortex surrounding a viscous core. For cells where intracellular components, such as the nucleus, take a considerable fraction of the cellular volume and play an active role in determining cell shape, the method described here will not be applicable without explicitly modeling these internal structures.
We model the cortex connecting the cell membrane and the cytoplasm with a Voigt model, which consists of the parallel connection of an elastic element kc (nN/μm3) and a viscous element τc (nNs/μm3). The latter describes the association and dissociation dynamics of the cross-linkers. We model the cytoplasm by a purely viscous element, τa (nNs/μm3), which is placed in series with the Voigt model (Fig. 2A). The element τa includes contributions from the interior of the cell as well as adhesion, friction and cytoskeletal reorganization. Strains of the cortex and cytoplasm are described by the variables xm and xc, respectively. Note that, in our simulations, we use pressure rather than force to induce the cellular deformations; this accounts for the extra μm2 found in the denominators of the parameters in our model.
Figure 2. Viscoelastic model of cell. A. Representation of the viscoelastic model of the cell. xc and xm denote the location of the cell cytoplasm and membrane, respectively; τc and kc define the mechanical model of the cell cortex; τa includes the viscous deformation of the cytoplasm as well as other components including adhesion. B. To validate our model and determine model parameters, we utilized micropipette aspiration technique. Relevant parameters include the radius of the micropipette (Rp), the radius of the cell (Rc), and the length of protrusion into the micropipette (Lp). C. Example of a Dictyostelium cell being aspirated into the micropipette at 0.65 nN/μm2. Time stamps are in seconds, scale bar shows 10 μm. D. Protrusion into the pipette was measured and was accounted for by the model. Different colored circles represent data from 22 individual experiments. Solid line represents the deformation defined by Eqn. 6 with parameters kc = 0.098 nN/μm3, τc = 0.064 nNs/μm3 and τa = 6.09 nNs/μm3.
As shown below, this combined Voigt-dashpot viscoelastic model reasonably approximates the mechanical properties of Dictyostelium where cross-linking proteins are predominantly enriched in the cortex. Extending our framework to other cell types may require different viscoelastic models to describe the cell of interest. For example, aspiration of chondrocytes suggests that these cells obey a Kelvin model (similar to the Voigt element, but includes an elastic component in series with the viscous element) . Once the appropriate viscoelastic model is developed, the implementation in the LSM framework introduced here is straightforward.
Experimental determination of model parameters
To determine appropriate parameters for the viscoelastic model, we used micropipette aspiration to apply step pressures (rapid increase of pressure from 0 to 0.65 nN/μm2) to individual cells [29,30]. In this technique, a small negative hydrostatic pressure is created at the tip of a micropipette. By bringing the micropipette into close proximity of the cellular surface, the cell is aspirated into the micropipette.
We applied step pressures to wild type interphase cells and measured cellular deformation as a function of time (Fig. 2B). Deformation is quantified by the length of cellular protrusion into the pipette tip, denoted Lp (Fig. 2C). We aspirated 22 cells with a radius of 4.3–6.1 μm, a pipette radius of 3.1 μm and a pressure of 0.65 nN/μm2. The cells deformed in two distinct phases (Fig. 2D). Within the first several seconds after application of the aspirator, the cellular deformation increased sharply, with the length of the aspirated cortex increasing to an average value of 4 μm. The deformation during this phase can be interpreted as being dominated by the elastic characteristics of the cytoskeletal network. Thereafter, the trajectory was dominated by slow cellular flow into the micropipette, increasing, on average, about 2 μm over the next 25 s.
The pressure applied by the micropipette aspirator is not the only pressure experienced by the cell. At rest, the cell is also under pressure from cortical tension (γten), which maintains the spherical shape of the cell. Under the cortical shell-liquid drop model , we assume that the cortical tension arises as surface tension (ignoring tangential stress). Following the Young-Laplace equation for liquid interfaces, the equilibrium pressure experienced by a spherical cell of radius Rc is
The cell's protrusion into the micropipette is driven by the aspiration pressure. As the cell is aspirated, the portion of the cell inside the micropipette will be a spherical cap of radius Rcap <Rc. Given a measured length of protrusion Lp, the radius of the spherical cap Rcap can be obtained (see Additional file 1). The cap's smaller radius gives rise to higher local curvature, creating a rounding pressure:
to oppose the aspiration
At the critical aspiration pressure, Pcrit, the cell extends a perfect hemispherical projection with radius Rp into the micropipette and does not protrude any further under this constant pressure. Thus, the critical pressure is:
The cortical tension has been measured to be 1–1.5 nN/μm in passive, wild type Dictyostelium cells [18,31,32]. Here, assuming cortical tension of γten = 1 nN/μm, pipette radius of Rp = 3.1 μm and cell radius of Rc = 5.1 μm, we can compute Pcrit to be approximately 0.25 nN/μm2. Because the applied pressure was greater than the critical pressure, the cell was continuously aspirated into the pipette. Cells were only tracked for 30 s, as longer timescales are dominated by cortical remodeling and turnover [33,34].
Pascal's law dictates that the hydrostatic pressure, Pext, in the micropipette is normal to the cell membrane inside the micropipette and has the same magnitude in all directions. Similarly, the cell's equilibrium pressure is normal to the cell membrane everywhere with the same magnitude. We used the total pressure, Ptotal = Pext - Pround, as the input to the cell's mechanical model. This pressure is applied to the cell membrane region around xm and is transferred directly to the cell's cortex, formed of cytoskeleton and its cross-linkers, just beneath the cell membrane. The corresponding mathematical model is:
where w0 represents the initial position of the cell cortex when no force is applied to the system. We define ℓ such that:
With this variable change, the transformed system can be written as:
Using the Voigt-dashpot model of Eqn. 5 to account for the viscoelastic response of the cell to an applied step pressure, the aspirated cellular length into the pipette, Lp, is given by:
Data from all 22 cells were combined. The following parameters in the viscoelastic model were obtained using a least squares fit (using Matlab's curve fitting toolbox): kc = 0.098 ± 0.007 nN/μm3, τc = 0.064 ± 0.018 nNs/μm3, and τa = 6.09 ± 1.44 nNs/μm3 (the ± value refer to a 95% confidence interval). With these parameter values, the 1-D model was able to capture the deformation trends observed in the experimental data (Fig. 2D). Note that the elastic constant obtained, when applying the methods of Theret et al.  and Hochmuth , is equivalent to an elastic Young's modulus of 70 pN/μm2, which is similar to the value of 95 pN/μm2 measured for Dictyostelium using different techniques .
Implementation of micropipette aspiration simulation
During micropipette aspiration, the cell's velocity is generated by externally applied pressure, as well as internally generated cellular pressures such as cortical tension. We now outline how the contribution of each pressure is computed and applied to the cell's potential function ϕ(x, t).
We choose to do the simulations in two dimensions. The level set method is directly applicable to three dimensions (3-D), and all of the level set equations either carry over without change into 3-D, or have natural extensions. In practice, however, the computational burden of 3-D simulations is significant and hence we restrict ourselves to two dimensions. To differentiate the forces (and hence pressures) which are 2-D, from the scalar pressures used above, we use bold characters.
Evolving the cell membrane
The simulation accounts for the effects of three pressures: those generated externally by the micropipette; those generated internally to maintain constant volume; and rounding pressure corresponding to the cell's cortical tension. Together, these pressures generate a velocity field that evolves the cell's membrane.
Externally applied pressure
To account for the force generated by the hydrostatic pressure in the micropipette, the pressure Pext, uniform in magnitude and normal to the cell membrane, is used in the LSM simulation. This force exists only inside the inner boundaries of the pipette.
Pressure due to volume conservation
We assume that, during aspiration, the cellular volume (V) remains constant. To enforce this constant volume condition numerically, we implement a pressure, acting normal to the surface:
where n is the outward normal. The cell's volume is evaluated by assuming the cell is radially symmetric:
To ensure that the cell's volume does not change during the course of the aspiration requires that Kvol be large. In our simulations, we set Kvol = 0.1 nN/μm5, which was sufficiently high to ensure that volume changes were minimal (Fig. 3G) while maintaining stability of the simulations.
Figure 3. Simulation of micropipette aspiration. A. To account for the solid surface of the micropipette, we introduce a mask potential function (Ψ) defined by the micropipette walls (black line). B. A cross section illustrates how the masking potential function is used to clip the evolving level set potential function. Based on the driving equations, the potential function evolves from ϕ(t)(solid blue) to ϕ(t + Δt)(dashed blue). However, this new position makes the cell cross into the pipette (defined by the mask function Ψ – green line). The level set function is then clipped to to account for the solid surface. C. Parallel spring-dashpot units are used to represent the viscoelastic state of the cell as the boundary evolves from Γ(t) to Γ(t + Δt). Each component consists of a viscoelastic model as defined in Fig. 2A. D. Simulation of micropipette aspiration implementing the viscoelastic model of the cell in the LSM (using the adjusted parameters; see main text). Shown is an overlay of simulation frames at t = 0 s (spherical cell, light grey), 0.5 s, 1 s, 5 s and 20 s (farthest protrusion, black). E. Measurements from aspiration simulations. Assuming an aspiration pressure of 0.65 nN/μm2, the protrusion into the cell from the simulation (black line) can account for the experimental data (grey dots; mean square error, MSE, is 0.74 μm; coefficient of determination, R2, of 0.78) nearly as well as the data fit (dotted line) from Fig. 2D. (MSE: 0.73, R2: 0.79). With slightly different parameters (see main text) the simulation (red line) overlaps the fitted data better (MSE: 0.73, R2: 0.79). Aspiration forces near the critical pressure (0.35 nN/μm2) can deform the cell initially, but do not draw it in further (green line). F. After 20 s of micropipette aspiration, the pressure in the aspirator is dropped, leading to a relaxation in the protrusion distance; a typical example is shown here. Time stamp indicate seconds after release of aspiration pressure, scale bar corresponds to10 μm. Note that the cell does not fully retract the aspirated portion. G. In a LSM simulation, the cell's retraction can be shown as the decrease of length of protrusion. Simulation of cell relaxation accurately demonstrates the lack of complete retraction observed (red). Also shown is the cell's volume (blue) during the simulation demonstrating that any volume changes are minimal.
Rounding pressure due to cortical tension
Resting cells experience cortical tension  which generates pressure, Peq, as shown in Eqn. 3. When a spherical cell is aspirated, the cell's cortex resists deformation.
The pressure generated depends on the local surface curvature
and a material property of the cortex referred to as the cortical tension (γten) according to:
Therefore the rounding pressure produce by the cell is Pround = Pten - Peq. This acts inward normal to the membrane.
We have chosen to define Pround as the difference between the tension and an equilibrium pressure. This is accordance to experiments on neutrophils that found that cortical tension depends on surface area . However, the latter term can also be incorporated into the volume conservation term. In particular, combining Pvol and Pround leads to:
The coefficient 2 in the pressure equation is introduced to account for the fact that our curvature calculation is based on the 2-D representation of cell shape, as the curvature of a sphere of radius r is 2/r, but the curvature of a 2-D circle is only 1/r. In the computation of curvature, spline-based interpolation was used to smooth out discretization noise.
Total pressure and cell evolution
In the above formulations, the total pressure outward normal to the cell membrane is:
The velocity vector, v, is defined for points on the cell membrane. This needs to be extrapolated to a velocity field to evolve the potential function ϕ. It is only the velocity variations tangential to a given interface that dictate the interface motion . A velocity field that minimizes the normal component of the field variation is achieved by extrapolating the membrane velocity with the nearest neighbor method. In other words, the velocity v(x) at a point x can be set equal to the membrane velocity v() at the membrane location closest to the point x. It has been shown that a signed distance function tends to stay a signed distance function when the closest neighbor extrapolation method is used . We can now use this velocity field to evolve the cell membrane according to Eqn. 2.
Eqn. 11 points to a difference between the LSM model of cellular deformation and the one-dimensional (1-D), scalar model used to obtain the viscoelastic parameters (Eqn. 6). In the latter, the pressure is co-aligned with the direction of the viscoelastic components, implying that the direction of motion is also always inline with the direction of the applied pressure. In the LSM simulation, the pressure is applied normal to the cell membrane, but the viscoelastic component, l, does not have to have the same directionality, and the resultant velocity vector is not always normal to the cell membrane. While providing us with good starting point for the parameter estimation, the 1-D formulation therefore can not be expected to explain the 2-D simulation completely.
Restricting cell shape inside micropipette
As the cell's level set potential function moves into the micropipette, its shape is restricted to remain inside the micropipette. This is achieved by first defining a mask potential function , Ψ, for the micropipette (Fig. 3A). The effect of the mask is to correct for the cell's potential function by clipping it (Fig. 3B) according to:
This restriction guarantees that the cell boundary never moves outside of the inner walls of the virtual micropipette. After this restriction, the net change in ϕ is: - ϕ(t), which translates (see Additional file 1) to an effective velocity that is normal to the cell membrane:
Thus, wherever clipping by the micropipette mask occurs, we must use this effective velocity to evolve the potential functions in simulation.
Evolution of the viscoelastic state of the cell
In our simulations, the cell can be represented by a series of parallel viscoelastic systems with the same parameters (Fig. 3C). These sub-systems are not interconnected, and the applied pressure on each system, Ptotal as defined in Eqn. 10, is normal to the cell membrane. We argue that applying total pressure to the parallel unconnected spring damper systems used in this model closely approximates cellular behavior when the following conditions are met:
1. Membrane pressure profile is piecewise smooth. This is a reasonable assumption as, in practice, pressure profiles are piecewise smooth. Even when a point force is applied to a particular location of the cell membrane, membrane elasticity will diffuse this force and make the pressure smooth locally.
2. Simulation grid density is dense enough for simulation stability, but not much denser than the discretization of the membrane pressure profile. With this assumption, the interpolation nature of level set method acts like a low pass filter, where effects of artificial abrupt jumps in the pressure profile are smoothed.
Let l(x, t), x ∈ Γ(t) be the viscoelastic state of the cell at time t and at a position x on the membrane. That is, |l| represents the length of the numerous parallel unconnected spring-damper systems. At a given position, x, on the membrane, there is a vector with length given by |l(x)| = |ℓ| in Eqn. 5, representing the state of a single spring-damper system. Then
Testing of model: Micropipette aspiration simulation
To summarize, the flow chart of the simulation steps is shown in Fig. 4. The implementation is derived from the Level Set Toolbox  and is coded in Matlab (Mathworks, Natick, MA). The simulations were implemented on a fixed grid of 10 μm in size, with density of 20 points/μm and 4 ms time steps. Simulating 15 seconds of aspiration takes approximately 8 h on a desk-top computer.
Figure 4. Algorithm for LSM simulation of micropipette aspiration.
We simulated the micropipette aspiration experiment under several different aspiration pressures. Using an aspiration pressure of 0.65 nN/μm2 (the pressure used to obtain our viscoelastic model parameters), our simulation reproduced the trend observed in real cells (black line in Fig. 3E). The result of this simulation did not completely overlap the least-squares fitted data, though the fit to the experimental data is nearly as good. The fitted data has a mean square error (MSE) of 0.73 μm and a coefficient of determination (R2) of 0.79; the simulation has 0.74 μm and 0.78 respectively. Using different parameter values: kc = 0.1 nN/μm3, τc = 0.08 nNs/μm3, and τa = 4.6 nNs/μm3, we were able to reproduce the fitted data slightly more accurately (Fig. 3D and red line in Fig. 3E; MSE of 0.73 μm and an R2 value of 0.79).
Using 0.35 nN/μm2 of pressure, the cell was rapidly and partially aspirated into the pipette. Thereafter, it remained nearly immobile. This simulation recreates the observed behavior of Dictyostelium cells at aspiration pressures near the critical pressure.
To test our model further, we simulated the relaxation of an aspirated cell and compared this to experimental results in which a cell is aspirated into the micropipette for approximately 20 s at which point the applied pressure is released. The cell responds by rapidly retracting the aspirated portion (Fig. 3F). The retraction gradually slows to a near halt, with a significant portion of the cell remaining inside the micropipette. This behavior was reproduced in our simulations. The simulated cell retraction from the micropipette is measured in the reduction of length of protrusion (Fig. 3G), matching the retraction behavior seen in live cells. As shown in Fig. 3G, the variation in cell volume during these simulations was less than 1%.
Simulating Dictyostelium cell shape changes using a simplified chemotaxis model
Having established that we can recreate the cellular shape during micropipette aspiration, in which externally applied pressures are driving cell shape changes, we consider a situation in which the pressures arise as a response to external stimuli. To this end we simulated the cell shape behavior of chemotactic Dictyostelium cells.
Dictyostelium cells have the ability to detect spatial differences in the concentration of the extracellular chemoattractant cAMP. They interpret these spatial differences and respond by localizing signaling molecules. These signaling molecules in turn bias the locations of actin polymerization driven protrusions and myosin-II motor mediated retractions, generating internal mechanical forces to deform the cell as well as propel the cell towards the chemoattractant [1,40].
Our goal in these simulations is not to propose new chemotaxis signaling mechanisms, or even to analyze the large number of proposed mechanisms (reviewed in ). Rather, it is to illustrate how cellular signaling can be coupled to the LSM framework to drive cellular deformations. Thus, we purposely implement a simple model connecting chemoattractant gradients with intracellular markers.
Implementation and testing
We base our model for pressure generation on a previously published signaling model that accounts for receptor mediated localization of phosphatidylinositol (3,4,5)-trisphosphate (PI(3,4,5)P3) . Though recent experimental data suggests that cells employ multiple parallel pathways to regulate chemotaxis [42,43], localization of this membrane lipid has been correlated with the appearance of pseudopods . Moreover, elevated levels of PI(3,4,5)P3 correlate temporally with increased levels of actin polymerization .
Rather than implementing the complete reaction-diffusion equations describing the PI(3,4,5)P3 model, we simplify it by using a steady-state distribution of PI(3,4,5)P3 along the cellular membrane. It was shown that the membrane concentration of PI(3,4,5)P3 is an amplified response of the relative cAMP concentration observed on the membrane [41,45]:
Next, we compute the pressure components contributing to cell motion, which include protrusion, retraction, volume conservation, and cortical tension pressures. To compute protrusion pressure, we first assume that actin polymerization creates a pressure wherever the PI(3,4,5)P3 concentration is above its mean level:
Similarly, we assume myosin-II retraction occurs wherever PI(3,4,5)P3 concentration is below its mean level:
Both of these act normal to the cell membrane. We let the proportionality constant in Eqn. 14 be absorbed into constants Kprot and Kretr. Eukaryotic cells can generate actin mediated protrusion pressures of 0.5–5 nN/μm2 . We chose Kprot = 0.5 nN/μm2 and Kretr = 1 nN/μm2.
When computing the conservation of volume pressure, we assume that the cell is flat with uniform thickness. Thus, volume conservation is equivalent to conserving the 2-D area of the cell:
The flat cell assumption also implies that the pressure generated by cortical tension depends only on the 2-D local surface curvature and the 2-D equilibrium pressure. The rounding pressure due to cortical tension is therefore given by:
Values of Karea = 0.2 nN/μm4 and Kten = 1 nN/μm were used in these simulations.
Summing all these components yields the total force normal to the cell membrane:
Finally, the membrane velocity is computed using Eqn. 11, with the same viscoelastic parameters τa, kc and τc. The simulation algorithm is similar to the micropipette aspiration case, and is summarized in Fig. 5.
Figure 5. Algorithm for LSM simulation of cell shape changes in response to external chemotaxis gradients. This algorithm includes only general steps required to generate the pressure profile. To simulate chemotaxis also requires that the chemoattractant gradient be generated and that the protrusive (Pprot) and retractive (Pretr) pressures be computed. These would be determined by specific models of chemotactic response. In our simulations, these were generated by Eqn. 15 and Eqn. 16, respectively.
This simulation successfully generated chemotaxis behavior (Fig. 6). In response to a chemoattractant gradient, the cell, whose shape was initialized as a circle, changed shape and migrated in the direction of the chemoattractant gradient (Fig. 6A). The pressure profile (Fig. 6B) and displacement (Fig. 6C) are shown as functions of local cAMP concentration and time, respectively. The cell achieved a velocity of 11.7 μm/min, which is similar to published velocities of Dictyostelium cells (e.g. 11.8 μm/min). During the simulation, the cellular area (and hence volume) remained nearly constant (Fig. 6C).
Figure 6. LSM simulations of cell shape changes during chemotaxis. A. We simulated the change in cellular morphology of a Dictyostelium cell exposed to a point source of chemoattractant (1 μM of cAMP). Shown is the resultant chemoattractant field (computed by solving the diffusion equation) as well as the location of the cell at times 0, 1.5, 10, 40, 60, 80 and 100 s. Initially, the cell is assumed to be round (red circle). B. In this model, pressure was determined by the concentration of PI(3,4,5)P3 on the membrane as described in the text. The maximum and minimum refer to the concentrations experienced by the cell around the membrane. C. The position of the cell (blue) was plotted as a function of time, showing fairly constant velocity (11.7 μm/min). Also shown is the cell's area (red) during the simulation demonstrating that changes are also minimal.
Membrane pressure profile and cell shape
While our simulations of Dictyostelium recreate the motion of the cell in response to the chemoattractant gradient, the resultant cell shape change is small and the steady-state morphology does not resemble that observed experimentally in chemotaxing. Wild type chemotaxing Dictyostelium cells become elongated (Fig. 7A). Other strains, including the amiB- mutants  can move stably in fan-like shapes that are reminiscent of keratocytes (Fig. 7D). Without determining the underlining molecular methods, we hypothesized that the difference in cell shape can be accounted for by the way that the force generation is distributed along the cell membrane. Our LSM simulation framework allows us to determine how these forces are distributed along the cell to generate the resulting cell shapes, both for wild type and mutants. To this end, we set out to replace our initial model, described by Eqn. 15 and 16, by one based on the observed morphologies.
Figure 7. Pressure profile drives cell shape. A. During chemotaxis, wild type Dictyostelium cells acquire a polarized, elongated morphology. B. Eqn. 19 was used to compute the pressure profile (red dots) necessary to maintain the elongated cell shape (inset) along the cell membrane, and this is plotted as a function of the local chemoattractant (cAMP) concentration. The maximum and minimum refer to the concentrations experienced by the cell around the membrane. Pressure profile used in the simulations (blue line) was obtained by fitting the computed pressure profile, details are given in the Additional file 1. C. Chemotaxing cell using the pressure profile of panel B. The shapes of the cell are shown at times 0, 1.5, 10, 20, 40, 60, 80 and 100 s. Other details in the simulation are as in Fig. 6. D-F. Simulations of chemotaxis in Dictyostelium amiB- cells. These mutant cells acquire a fan-like morphology (panel D) and move along their broad axis. This form of movement was recreated using the pressure profile of panel E (colors as in panel B). F. Chemotaxing cell using the force profile of panel E. Times of the shapes are as in panel C.
Given a stable cell shape Γ0 traveling at velocity u, we let Γu be the displaced cell at time Δt, and ϕ0 and ϕu be the potential functions representing Γ0 and Γu respectively. The effective velocity field necessary for this displacement is:
If the cell shape is at steady state, we can assume that the internal viscoelastic network is also in steady state, that is, . Therefore, from Eqn. 5, we compute the viscoelastic steady state ℓ = Ptotal/kc.
Taking into account the effect of cortex tension, and assuming that there is no cell volume deviations, we can compute:
where Pten is the cortical tension-driven rounding pressure defined in Eqn. 17. Using this formula, and a cell velocity of 10 μm/min, we calculated the pressure profiles required to generate cell shapes seen in wild type cells as well as in amiB- cells.
Obtaining these pressure profiles is straight-forward computationally, taking less than one minute of CPU time on a desk-top computer. It does require, however, a smooth shape. Thus, a certain amount of image processing is needed when using segmented images from experiments. Moreover, the formula in Eqn. 18 is based on a steady-state shape. Handling transient cell shape changes, such as protrusions or retractions, needs a local description of the velocity v(x).
Our results indicate that to generate polarized cell morphologies observed in wild type Dictyostelium cells, the protrusive forces must be primarily concentrated along the anterior ≈ 25% portion of the cell; see Fig. 7B. This is reminiscent of the PI(3,4,5)P3 threshold observed in cells [45,49]. At the sides, a smaller and less localized retractive force gives the cell its elongated shape. When this pressure profile was used to simulate a chemotaxing cell (Fig. 7C), the resulting virtual cell achieved an elongated shape and chemotaxed successfully to the source of chemoattractant achieving a stable velocity of 11.1 μm/min.
Clearly, a different pressure profile is needed to generate a fan like shape as observed in amiB- cells (Fig. 7D). Here, the maximum protrusive force is spread out considerably more at the front, while large amount of retraction force is still needed to pull the tail region along. Using this pressure profile in the chemotaxis simulation led to a migrating cell with stable shape similar to that seen experimentally (Fig. 7F). The resultant fan-shaped cell achieved the stable velocity of 9.7 μm/min.
We have shown that the simulation framework we have developed can be used to model cell shape deformations as well as cell motility. The simulations can produce deformations seen during micropipette aspiration experiments. This requires parameter values for the viscoelastic model which can be obtained experimentally. It should be noted, however, that 2-D simulations using parameters based on a 1-D model may not reproduce the 1-D model simulation precisely.
In the simulations of cell shape changes during chemotaxis, we saw that our simple model for generating the cell's protrusive and retractive forces in response to a chemoattractant gradient does not produce experimentally observed cell shapes. However, our techniques allow us to work backwards from shape to obtain the required forces. We determined that generating the elongated cell shape requires a large protrusive force at the front (the pressure profile there is positive). At the sides, there is a large retractive force (the pressure profile there is negative). While measuring this pressure profile directly would be difficult, it is possible to image fluorescently-tagged myosin-II to infer a measure of the forces acting on the cell. Under the assumption that the retractive force is being generated by myosin-II, we expect that myosin-II would be greatly enriched at the sides. Quantitatively, the spatial distribution of myosin could be used to estimate how much force is being generated along the membrane (as has been done during cytokinesis ).
LY implemented the LSM simulations and drafted the manuscript. JCE performed experiments to measure the viscoelastic properties of cells, under the guidance of DNR. BLK and SES participated in the implementation of the LSM algorithm. PAI conceived of the study, and participated in its design and coordination. LY, JCE, DNR and PAI wrote the manuscript which was read and approved by all the authors.
The authors declare that they have no competing interests.
We thank Charles Wolgemuth (U. Connecticut) who first suggested to us the use of level set methods for simulating chemotaxis. We also thank Dr Taro Q. P. Uyeda (Tsukuba University, Japan) for images of the amiB- cells. This work was supported in part by grants from the NIH (NIGMS R01-71920) and the NSF (0621740).
Symp Soc Exp Biol 1995, 49:265-276. PubMed Abstract
J Comput Phys 1998, 147:86-113. Publisher Full Text
SIAM J Multiscale Modeling and Simulation 2005, 3(2):413-439. Publisher Full Text
Reichl E, Ren Y, Morphew M, Delannoy M, Effler J, Girard K, Divi S, Iglesias P, Kuo S, Robinson D: Interactions between myosin and actin crosslinkers control cytokinesis contractility dynamics and mechanics.
J Comput Phys 1988, 79:12-49. Publisher Full Text
IEEE Intern Conf Image Proc 2006, :221-224. Publisher Full Text
Comp Mat Sci 2007, 39(4):857-861. Publisher Full Text
J Comput Phys 1992, 100(2):209-228. Publisher Full Text
J Comp Phys 1994, 114:146-159. Publisher Full Text
J Biomech Eng 1988, 110:190-199. PubMed Abstract
J Comput Phys 1996, 127:179-195. Publisher Full Text
Ma L, Janetopoulos C, Yang L, Devreotes PN, Iglesias PA: Two complementary, local excitation, global inhibition mechanisms acting in parallel can explain the chemoattractant-induced regulation of PI(3,4,5)P3 response in Dictyostelium cells.
Janetopoulos C, Ma L, Devreotes PN, Iglesias PA: Chemoattractant-induced phosphatidylinositol 3,4,5-trisphosphate accumulation is spatially amplified and adapts, independent of the actin cytoskeleton.
Wessels D, Brincks R, Kuhl S, Stepanovic V, Daniels KJ, Weeks G, Lim CJ, Spiegelman G, Fuller D, Iranfar N, Loomis WF, Soll DR: RasC plays a role in transduction of temporal gradient information in the cyclic-AMP wave of Dictyostelium discoideum.