Abstract
Background
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.
Results
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.
Conclusion
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.
Background
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 [1]. During cytokinesis, a single cell undergoes significant cytoskeletal deformation, reforming into two daughter cells [2]. 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 [3]. Few such models exist, even though it is now appreciated that the response of cellsignaling pathways can be regulated in response to alterations in cell size and shape [4]. The traditional method of simulating cellular deformations is by specifying the boundary of the cell explicitly through a finiteelement model (FEM) [57]. 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 [8] 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 [9]. In the IBM, the NavierStokes 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 bookkeeping 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 [10].
More recently, the cellular Potts model (CPM) has become a popular vehicle to simulate cell shape changes [11]. 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 sizedeformations as well as membrane bending. Cellular Potts models have been used to simulate twodimensional (2D) models of cell motility in fish keratocytes [12] and amoebae [13]. 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 experimentallybased 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 2D computational model of the lamellipodium keratocyte motility [14]. In this model, the cellular domain is updated at each time step based on the protrusive and retractive forces (actin polymerization and actomyosin contraction) and regridded. 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 largescale 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 largescale shape changes.
Biological background
Cells derive their mechanical properties from actin, actinassociated proteins, and motor proteins such as myosinII [2], 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 [1517]. Actin filaments alone do not create significant mechanical resistance; instead, crosslinking of actin filaments by various actin binding proteins imparts mechanical rigidity to the cell. Under applied load, crosslinked actin networks behave similarly to an elastic solid and can be described using Hooke's law. However, because crosslinking proteins bind to and dissociate from actin filaments, actinbased networks may also exhibit viscous flow. MyosinII, 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 timedependent 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 [19]. LSM is a numerical technique for tracking interfaces and shapes which has been widely used in various fields including computer graphics [20], image processing [21], computational fluid dynamics [22] and material science [23].
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 zerocontour. Γ(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 closedcontour Γ(t). The LSM requires a potential function (Fig. 1C.), denoted as ϕ(x, t), that is related to Γ(t) according to:
Thus, Γ(t) is the zerolevel 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 [24], 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 myosinII 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 reinitialize the potential function periodically [24]. This can be done using the reinitialization equation [25]:
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 [27]. 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 k_{c }(nN/μm^{3}) and a viscous element τ_{c }(nNs/μm^{3}). The latter describes the association and dissociation dynamics of the crosslinkers. We model the cytoplasm by a purely viscous element, τ_{a }(nNs/μm^{3}), 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 x_{m }and x_{c}, respectively. Note that, in our simulations, we use pressure rather than force to induce the cellular deformations; this accounts for the extra μm^{2 }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. x_{c }and x_{m }denote the location of the cell cytoplasm and membrane, respectively; τ_{c }and k_{c }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 (R_{p}), the radius of the cell (R_{c}), and the length of protrusion into the micropipette (L_{p}). C. Example of a Dictyostelium cell being aspirated into the micropipette at 0.65 nN/μm^{2}. 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 k_{c }= 0.098 nN/μm^{3}, τ_{c }= 0.064 nNs/μm^{3 }and τ_{a }= 6.09 nNs/μm^{3}.
As shown below, this combined Voigtdashpot viscoelastic model reasonably approximates the mechanical properties of Dictyostelium where crosslinking 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) [28]. 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/μm^{2}) 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 L_{p }(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/μm^{2}. 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 shellliquid drop model [31], we assume that the cortical tension arises as surface tension (ignoring tangential stress). Following the YoungLaplace equation for liquid interfaces, the equilibrium pressure experienced by a spherical cell of radius R_{c }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 R_{cap }<R_{c}. Given a measured length of protrusion L_{p}, the radius of the spherical cap R_{cap }can be obtained (see Additional file 1). The cap's smaller radius gives rise to higher local curvature, creating a rounding pressure:
Additional file 1. This document presents detailed derivation of several of the formulae in the text.
Format: PDF Size: 78KB Download file
This file can be viewed with: Adobe Acrobat Reader
to oppose the aspiration
At the critical aspiration pressure, P_{crit}, the cell extends a perfect hemispherical projection with radius R_{p }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 R_{p }= 3.1 μm and cell radius of R_{c }= 5.1 μm, we can compute P_{crit }to be approximately 0.25 nN/μm^{2}. 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, P_{ext}, 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, P_{total }= P_{ext } P_{round}, as the input to the cell's mechanical model. This pressure is applied to the cell membrane region around x_{m }and is transferred directly to the cell's cortex, formed of cytoskeleton and its crosslinkers, just beneath the cell membrane. The corresponding mathematical model is:
where w_{0 }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 Voigtdashpot 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, L_{p}, 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): k_{c }= 0.098 ± 0.007 nN/μm^{3}, τ_{c }= 0.064 ± 0.018 nNs/μm^{3}, and τ_{a }= 6.09 ± 1.44 nNs/μm^{3 }(the ± value refer to a 95% confidence interval). With these parameter values, the 1D 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. [35] and Hochmuth [30], is equivalent to an elastic Young's modulus of 70 pN/μm^{2}, which is similar to the value of 95 pN/μm^{2 }measured for Dictyostelium using different techniques [18].
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 (3D), and all of the level set equations either carry over without change into 3D, or have natural extensions. In practice, however, the computational burden of 3D simulations is significant and hence we restrict ourselves to two dimensions. To differentiate the forces (and hence pressures) which are 2D, 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 P_{ext}, 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 K_{vol }be large. In our simulations, we set K_{vol }= 0.1 nN/μm^{5}, 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 springdashpot 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/μm^{2}, 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, R^{2}, of 0.78) nearly as well as the data fit (dotted line) from Fig. 2D. (MSE: 0.73, R^{2}: 0.79). With slightly different parameters (see main text) the simulation (red line) overlaps the fitted data better (MSE: 0.73, R^{2}: 0.79). Aspiration forces near the critical pressure (0.35 nN/μm^{2}) 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 [31] which generates pressure, P_{eq}, 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 P_{round }= P_{ten } P_{eq}. This acts inward normal to the membrane.
We have chosen to define P_{round }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 [36]. However, the latter term can also be incorporated into the volume conservation term. In particular, combining P_{vol }and P_{round }leads to:
where = V_{actual }+ 2γ_{ten}/(R_{c}K_{vol}).
The coefficient 2 in the pressure equation is introduced to account for the fact that our curvature calculation is based on the 2D representation of cell shape, as the curvature of a sphere of radius r is 2/r, but the curvature of a 2D circle is only 1/r. In the computation of curvature, splinebased 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 formulation of in Eqn. 5 provides us with the pressurevelocity relationship:
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 [37]. 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 [38]. 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 onedimensional (1D), scalar model used to obtain the viscoelastic parameters (Eqn. 6). In the latter, the pressure is coaligned 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 1D formulation therefore can not be expected to explain the 2D 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 [39], Ψ, 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 subsystems are not interconnected, and the applied pressure on each system, P_{total }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 springdamper 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 springdamper system. Then
where D is the Jacobian operator, [Dl]v represents the displacement of the whole cell membrane, and as defined in Eqn. 5. The equation describing the evolution of l is:
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 [39] 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 desktop 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/μm^{2 }(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 leastsquares 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 (R^{2}) of 0.79; the simulation has 0.74 μm and 0.78 respectively. Using different parameter values: k_{c }= 0.1 nN/μm^{3}, τ_{c }= 0.08 nNs/μm^{3}, and τ_{a }= 4.6 nNs/μm^{3}, 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 R^{2 }value of 0.79).
Using 0.35 nN/μm^{2 }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 myosinII 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 [3]). 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)P_{3}) [41]. 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 [40]. Moreover, elevated levels of PI(3,4,5)P_{3 }correlate temporally with increased levels of actin polymerization [44].
Rather than implementing the complete reactiondiffusion equations describing the PI(3,4,5)P_{3 }model, we simplify it by using a steadystate distribution of PI(3,4,5)P_{3 }along the cellular membrane. It was shown that the membrane concentration of PI(3,4,5)P_{3 }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)P_{3 }concentration is above its mean level:
Similarly, we assume myosinII retraction occurs wherever PI(3,4,5)P_{3 }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 K_{prot }and K_{retr}. Eukaryotic cells can generate actin mediated protrusion pressures of 0.5–5 nN/μm^{2 }[46]. We chose K_{prot }= 0.5 nN/μm^{2 }and K_{retr }= 1 nN/μm^{2}.
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 2D area of the cell:
The flat cell assumption also implies that the pressure generated by cortical tension depends only on the 2D local surface curvature and the 2D equilibrium pressure. The rounding pressure due to cortical tension is therefore given by:
Values of K_{area }= 0.2 nN/μm^{4 }and K_{ten }= 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}, k_{c }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 (P_{prot}) and retractive (P_{retr}) 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[47]). 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)P_{3 }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 steadystate morphology does not resemble that observed experimentally in chemotaxing. Wild type chemotaxing Dictyostelium cells become elongated (Fig. 7A). Other strains, including the amiB^{ }mutants [48] can move stably in fanlike 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. DF. Simulations of chemotaxis in Dictyostelium amiB^{ }cells. These mutant cells acquire a fanlike 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 ℓ = P_{total}/k_{c}.
Moreover, the membrane speed at steady state is expressed as = P_{total}/τ_{a}. Combined with Eqn. 18, we find P_{total}:
Taking into account the effect of cortex tension, and assuming that there is no cell volume deviations, we can compute:
where P_{ten }is the cortical tensiondriven 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 straightforward computationally, taking less than one minute of CPU time on a desktop 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 steadystate 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)P_{3 }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 fanshaped cell achieved the stable velocity of 9.7 μm/min.
Conclusion
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 2D simulations using parameters based on a 1D model may not reproduce the 1D 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 fluorescentlytagged myosinII to infer a measure of the forces acting on the cell. Under the assumption that the retractive force is being generated by myosinII, we expect that myosinII 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 [50]).
Authors' contributions
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.
Competing interests
The authors declare that they have no competing interests.
Acknowledgements
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 R0171920) and the NSF (0621740).
References

Devreotes P, Janetopoulos C: Eukaryotic chemotaxis: distinctions between directional sensing and polarization.
J Biol Chem 2003, 278(23):2044520448. PubMed Abstract  Publisher Full Text

Reichl EM, Effler JC, Robinson DN: The stress and strain of cytokinesis.
Trends Cell Biol 2005, 15(4):200206. PubMed Abstract  Publisher Full Text

Iglesias P, Devreotes P: Navigating through models of chemotaxis.
Curr Opin Cell Biol 2008, 20:3540. PubMed Abstract  Publisher Full Text

Meyers J, Craig J, Odde D: Potential for control of signaling pathways via cell size and shape.
Curr Biol 2006, 16:16851693. PubMed Abstract  Publisher Full Text

Jean RP, Chen CS, Spector AA: Finiteelement analysis of the adhesioncytoskeletonnucleus mechanotransduction pathway during endothelial cell rounding: axisymmetric model.
J Biomech Eng 2005, 127(4):594600. PubMed Abstract  Publisher Full Text

Baaijens FPT, Trickey WR, Laursen TA, Guilak F: Large deformation finite element analysis of micropipette aspiration to determine the mechanical properties of the chondrocyte.
Ann Biomed Eng 2005, 33(4):494501. PubMed Abstract  Publisher Full Text

Hartmann C, Delgado A: Numerical simulation of the mechanics of a yeast cell under high hydrostatic pressure.
J Biomech 2004, 37(7):977987. PubMed Abstract  Publisher Full Text

Peskin C, McQueen D: A general method for the computer simulation of biological systems interacting with fluids.
Symp Soc Exp Biol 1995, 49:265276. PubMed Abstract

Bottino DC: Modeling viscoelastic networks and cell deformation in the context of the immersed boundary method.
J Comput Phys 1998, 147:86113. Publisher Full Text

Jadhav S, Eggleton CD, Konstantopoulos K: A 3D computational model predicts that cell deformation affects selectinmediated leukocyte rolling.
Biophys J 2005, 88:96104. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Graner F, Glazier J: Simulation of biological cell sorting using a twodimensional extended Potts model.
Phys Rev Lett 1992, 69:20132016. PubMed Abstract  Publisher Full Text

Marée AFM, Jilkine A, Dawes A, Grieneisen VA, EdelsteinKeshet L: Polarization and movement of keratocytes: a multiscale modelling approach.
Bull Math Biol 2006, 68(5):11691211. PubMed Abstract  Publisher Full Text

Nishimura SI, Sasai M: Modulation of the reaction rate of regulating protein induces large morphological and motional change of amoebic cell.
J Theor Biol 2007, 245(2):230237. PubMed Abstract  Publisher Full Text

Rubinstein B, Jacobson K, Mogilner A: Multiscale twodimensional modeling of a motile simpleshaped cell.
SIAM J Multiscale Modeling and Simulation 2005, 3(2):413439. Publisher Full Text

Bray D: Cell Movements: From Molecules to Motility. 2nd edition. Garland; 2000.

Gardel ML, Shin JH, MacKintosh FC, Mahadevan L, Matsudaira P, Weitz DA: Elastic behavior of crosslinked and bundled actin networks.
Science 2004, 304(5675):13011305. PubMed Abstract  Publisher Full Text

Wachsstock DH, Schwarz WH, Pollard TD: Crosslinker dynamics determine the mechanical properties of actin gels.
Biophys J 1994, 66:801809. PubMed Abstract  Publisher Full Text  PubMed Central 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.
Curr Biol 2008, 18:471480. PubMed Abstract  Publisher Full Text

Osher S, Sethian JA: Fronts propagating with curvature dependent speed: algorithms based on HamiltonJacobi formulations.
J Comput Phys 1988, 79:1249. Publisher Full Text

Fablet R, Pujolle S, Chessel A, Benzinou A, Cao F: Variational levelset reconstruction of accretionary morphogenesis from images.
IEEE Intern Conf Image Proc 2006, 221224. Publisher Full Text

Machacek M, Danuser G: Morphodynamic profiling of protrusion phenotypes.
Biophys J 2006, 90(4):143952. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Almomani T, Udaykumar H, Marshall J, Chandran K: Microscale dynamic simulation of erythrocyteplatelet interaction in blood flow.
Ann Biomed Eng 2008. PubMed Abstract  Publisher Full Text

Zhao X, Yin Y, Yang B, Zhu B, Tian X: Level set and geodesic active contours based measurement of material removal between serial sections.
Comp Mat Sci 2007, 39(4):857861. Publisher Full Text

Mulder W, Osher S, Sethian JA: Computing interface motion in compressible gas dynamics.
J Comput Phys 1992, 100(2):209228. Publisher Full Text

Sussman M, Smereka P, Osher S: A level set approach for computing solutions to incompressible twophase flow.
J Comp Phys 1994, 114:146159. Publisher Full Text

Fung YC: Biomechanics: Mechanical Properties of Living Tissues. 2nd edition. Springer; 1993.

Trickey WR, Lee GM, Guilak F: Viscoelastic properties of chondrocytes from normal and osteoarthritic human cartilage.
J Orthop Res 2000, 18(6):891898. PubMed Abstract  Publisher Full Text

Leipzig ND, Athanasiou KA: Unconfined creep compression of chondrocytes.
J Biomech 2005, 38:7785. PubMed Abstract  Publisher Full Text

Effler JC, Kee YS, Berk JM, Tran MN, Iglesias PA, Robinson DN: Mitosisspecific mechanosensing and contractileprotein redistribution control cell shape.
Curr Biol 2006, 16(19):19621967. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hochmuth RM: Micropipette aspiration of living cells.
J Biomech 2000, 33:1522. PubMed Abstract  Publisher Full Text

Dai J, TingBeall HP, Hochmuth RM, Sheetz MP, Titus MA: Myosin I contributes to the generation of resting cortical tension.
Biophys J 1999, 77(2):11681176. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Octtaviani E, Effler JC, Robinson DN: Enlazin, a natural fusion of two classes of canonical cytoskeletal proteins, contributes to cytokinesis dynamics.
Mol Biol Cell 2006, 17(12):52755286. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Robinson DN, Ocon SS, Rock RS, Spudich JA: Dynacortin is a novel actin bundling protein that localizes to dynamic actin structures.
J Biol Chem 2002, 277(11):90889095. PubMed Abstract  Publisher Full Text

Diez S, Gerisch G, Anderson K, MüllerTaubenberger A, Bretschneider T: Subsecond reorganization of the actin network in cell motility and chemotaxis.
Proc Natl Acad Sci USA 2005, 102(21):76017606. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Theret D, Levesque M, Sato M, Nerem R, Wheeler L: The application of a homogeneous halfspace model in the analysis of endothelial cell micropipette measurements.
J Biomech Eng 1988, 110:190199. PubMed Abstract

Herant M, Heinrich V, Dembo M: Mechanics of neutrophil phagocytosis: behavior of the cortical tension.
J Cell Sci 2005, 118(Pt 9):17891797. PubMed Abstract  Publisher Full Text

Osher S, Fedkiw R: Level Set Methods and Dynamic Implicit Surfaces. Springer; 2003.

Zhao H, Chan T, Merriman B, Osher S: A variational level set approach to multiphase motion.
J Comput Phys 1996, 127:179195. Publisher Full Text

Mitchell IM: A toolbox of level set methods.
UBC Department of Computer Science Technical Report TR200711 2007.

Manahan CL, Iglesias PA, Long Y, Devreotes PN: Chemoattractant signaling in Dictyostelium discoideum.
Annu Rev Cell Dev Biol 2004, 20:22353. PubMed Abstract  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 chemoattractantinduced regulation of PI(3,4,5)P_{3 }response in Dictyostelium cells.
Biophys J 2004, 87(6):37643774. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Chen L, Iijima M, Tang M, Landree MA, Huang YE, Xiong Y, Iglesias PA, Devreotes PN: PLA2 and PI3K/PTEN pathways act in parallel to mediate chemotaxis.
Dev Cell 2007, 12(4):60314. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

van Haastert PJ, KeizerGunnink I, Kortholt A: Essential role of PI3kinase and phospholipase A2 in Dictyostelium discoideum chemotaxis.
J Cell Biol 2007, 177(5):80916. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Iijima M, Devreotes P: Tumor suppressor PTEN mediates sensing of chemoattractant gradients.
Cell 2002, 109(5):599610. PubMed Abstract  Publisher Full Text

Janetopoulos C, Ma L, Devreotes PN, Iglesias PA: Chemoattractantinduced phosphatidylinositol 3,4,5trisphosphate accumulation is spatially amplified and adapts, independent of the actin cytoskeleton.
Proc Natl Acad Sci USA 2004, 101(24):89516. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Ananthakrishnan R, Ehrlicher A: The forces behind cell movement.
Int J Biol Sci 2007, 3:303317. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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 cyclicAMP wave of Dictyostelium discoideum.
Eukaryot Cell 2004, 3(3):646662. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Asano Y, Mizuno T, Kon T, Nagasaki A, Sutoh K, Uyeda TQP: Keratocytelike locomotion in amiBnull Dictyostelium cells.
Cell Motil Cytoskeleton 2004, 59:1727. PubMed Abstract  Publisher Full Text

Levine H, Kessler D, Rappel W: Directional sensing in eukaryotic chemotaxis: a balanced inactivation model.
Proc Natl Acad Sci USA 2006, 103:97619766. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Robinson D, Cavet G, Warrick H, Spudich J: Quantitation of the distribution and flux of myosinII during cytokinesis.
BMC Cell Biol 2002, 3:4. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text