Methodology article

Modeling cellular deformations using the level set formalism

Liu Yang^{* 1} , Janet C Effler^{* 1} , Brett L Kutscher¹ , Sarah E Sullivan¹ , Douglas N Robinson² and Pablo A Iglesias¹

¹Electrical & Computer Engineering, Johns Hopkins University, Baltimore, MD 21218, USA

²Cell Biology, Johns Hopkins School of Medicine, Baltimore, MD 21205, USA

author email corresponding author email* Contributed equally

BMC Systems Biology 2008, 2:68doi:10.1186/1752-0509-2-68

The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1752-0509/2/68

Received:	16 January 2008
Accepted:	24 July 2008
Published:	24 July 2008

© 2008 Yang et al; licensee BioMed Central Ltd.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 cell-signaling 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 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 [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 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 [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 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 [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 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 [14]. 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.

Biological background

Cells derive their mechanical properties from actin, actin-associated proteins, and motor proteins such as myosin-II [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 [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 [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].

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:

Γ(t) = {x|ϕ(x, t) = 0}.

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 [24], defined by:

(1)

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:

(2)

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 [24]. This can be done using the re-initialization 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³) and a viscous element τ_c (nNs/μm³). The latter describes the association and dissociation dynamics of the cross-linkers. We model the cytoplasm by a purely viscous element, τ_a (nNs/μm³), 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² found in the denominators of the parameters in our model.

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) [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²) 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². 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 [31], 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 R_c is

P_eq = 2γ_ten/R_c.(3)

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:

P_round(L_P) = 2γ_ten/R_cap - P_eq,

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:

P_crit = P_round(R_p) = 2γ_ten(1/R_p - 1/R_c).(4)

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². 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 cross-linkers, just beneath the cell membrane. The corresponding mathematical model is:

where w₀ represents the initial position of the cell cortex when no force is applied to the system. We define ℓ such that:

x_c = x_m - ℓ - w₀.

With this variable change, the transformed system can be written as:

(5)

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, L_p, is given by:

(6)

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³, τ_c = 0.064 ± 0.018 nNs/μm³, and τ_a = 6.09 ± 1.44 nNs/μm³ (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. [35] and Hochmuth [30], is equivalent to an elastic Young's modulus of 70 pN/μm², which is similar to the value of 95 pN/μm² 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 (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 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:

P_vol = K_vol(V_resting - V_actual)n(7)

where n is the outward normal. The cell's volume is evaluated by assuming the cell is radially symmetric:

V_actual = ∫_{cell length}πr(x)²dx.(8)

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⁵, which was sufficiently high to ensure that volume changes were minimal (Fig. 3G) while maintaining stability of the simulations.

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:

P_ten(x) = 2γ_tenκ(x)n.(9)

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_cK_vol).

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:

P_total = P_ext + P_vol - P_round.(10)

The formulation of in Eqn. 5 provides us with the pressure-velocity relationship:

(11)

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 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 [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:

(12)

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, 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 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

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:

(13)

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 desk-top computer.

We simulated the micropipette aspiration experiment under several different aspiration pressures. Using an aspiration pressure of 0.65 nN/μm² (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 (R²) of 0.79; the simulation has 0.74 μm and 0.78 respectively. Using different parameter values: k_c = 0.1 nN/μm³, τ_c = 0.08 nNs/μm³, and τ_a = 4.6 nNs/μm³, 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² value of 0.79).

Using 0.35 nN/μm² 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 [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₃) [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₃ correlate temporally with increased levels of actin polymerization [44].

Rather than implementing the complete reaction-diffusion equations describing the PI(3,4,5)P₃ model, we simplify it by using a steady-state distribution of PI(3,4,5)P₃ along the cellular membrane. It was shown that the membrane concentration of PI(3,4,5)P₃ is an amplified response of the relative cAMP concentration observed on the membrane [41,45]:

PI(3,4,5)P₃ ∝ [cAMP/mean(cAMP)]³.(14)

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₃ concentration is above its mean level:

(15)

Similarly, we assume myosin-II retraction occurs wherever PI(3,4,5)P₃ concentration is below its mean level:

(16)

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² [46]. We chose K_prot = 0.5 nN/μm² and K_retr = 1 nN/μm².

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:

P_area = K_area(A₀ - A)n.

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:

P_ten = K_ten(κ(x) - 1/R_c)n.(17)

Values of K_area = 0.2 nN/μm⁴ and K_ten = 1 nN/μm were used in these simulations.

Summing all these components yields the total force normal to the cell membrane:

P_total = P_prot + P_retr + P_area - P_ten.

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.

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).

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 [48] 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.

Given a stable cell shape Γ₀ traveling at velocity u, we let Γ_u be the displaced cell at time Δt, and ϕ₀ and ϕ_u be the potential functions representing Γ₀ and Γ_u respectively. The effective velocity field necessary for this displacement is:

(18)

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:

(19)

where P_ten 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)P₃ 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.

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 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 [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 R01-71920) and the NSF (0621740).

References

1. Devreotes P, Janetopoulos C: Eukaryotic chemotaxis: distinctions between directional sensing and polarization.

   J Biol Chem 2003, 278(23):20445-20448. PubMed Abstract | Publisher Full Text

2. Reichl EM, Effler JC, Robinson DN: The stress and strain of cytokinesis.

   Trends Cell Biol 2005, 15(4):200-206. PubMed Abstract | Publisher Full Text

3. Iglesias P, Devreotes P: Navigating through models of chemotaxis.

   Curr Opin Cell Biol 2008, 20:35-40. PubMed Abstract | Publisher Full Text

4. Meyers J, Craig J, Odde D: Potential for control of signaling pathways via cell size and shape.

   Curr Biol 2006, 16:1685-1693. PubMed Abstract | Publisher Full Text

5. Jean RP, Chen CS, Spector AA: Finite-element analysis of the adhesion-cytoskeleton-nucleus mechanotransduction pathway during endothelial cell rounding: axisymmetric model.

   J Biomech Eng 2005, 127(4):594-600. PubMed Abstract | Publisher Full Text

6. 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):494-501. PubMed Abstract | Publisher Full Text

7. Hartmann C, Delgado A: Numerical simulation of the mechanics of a yeast cell under high hydrostatic pressure.

   J Biomech 2004, 37(7):977-987. PubMed Abstract | Publisher Full Text

8. Peskin C, McQueen D: A general method for the computer simulation of biological systems interacting with fluids.

   Symp Soc Exp Biol 1995, 49:265-276. PubMed Abstract

9. Bottino DC: Modeling viscoelastic networks and cell deformation in the context of the immersed boundary method.

   J Comput Phys 1998, 147:86-113. Publisher Full Text

10. Jadhav S, Eggleton CD, Konstantopoulos K: A 3-D computational model predicts that cell deformation affects selectin-mediated leukocyte rolling.

    Biophys J 2005, 88:96-104. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

11. Graner F, Glazier J: Simulation of biological cell sorting using a two-dimensional extended Potts model.

    Phys Rev Lett 1992, 69:2013-2016. PubMed Abstract | Publisher Full Text

12. Marée AFM, Jilkine A, Dawes A, Grieneisen VA, Edelstein-Keshet L: Polarization and movement of keratocytes: a multiscale modelling approach.

    Bull Math Biol 2006, 68(5):1169-1211. PubMed Abstract | Publisher Full Text

13. 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):230-237. PubMed Abstract | Publisher Full Text

14. Rubinstein B, Jacobson K, Mogilner A: Multiscale two-dimensional modeling of a motile simple-shaped cell.

    SIAM J Multiscale Modeling and Simulation 2005, 3(2):413-439. Publisher Full Text

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

16. Gardel ML, Shin JH, MacKintosh FC, Mahadevan L, Matsudaira P, Weitz DA: Elastic behavior of cross-linked and bundled actin networks.

    Science 2004, 304(5675):1301-1305. PubMed Abstract | Publisher Full Text

17. Wachsstock DH, Schwarz WH, Pollard TD: Cross-linker dynamics determine the mechanical properties of actin gels.

    Biophys J 1994, 66:801-809. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

18. 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:471-480. PubMed Abstract | Publisher Full Text

19. Osher S, Sethian JA: Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations.

    J Comput Phys 1988, 79:12-49. Publisher Full Text

20. Fablet R, Pujolle S, Chessel A, Benzinou A, Cao F: Variational level-set reconstruction of accretionary morphogenesis from images.

    IEEE Intern Conf Image Proc 2006, 221-224. Publisher Full Text

21. Machacek M, Danuser G: Morphodynamic profiling of protrusion phenotypes.

    Biophys J 2006, 90(4):1439-52. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

22. Almomani T, Udaykumar H, Marshall J, Chandran K: Micro-scale dynamic simulation of erythrocyte-platelet interaction in blood flow.

    Ann Biomed Eng 2008.PubMed Abstract | Publisher Full Text

23. 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):857-861. Publisher Full Text

24. Mulder W, Osher S, Sethian JA: Computing interface motion in compressible gas dynamics.

    J Comput Phys 1992, 100(2):209-228. Publisher Full Text

25. Sussman M, Smereka P, Osher S: A level set approach for computing solutions to incompressible two-phase flow.

    J Comp Phys 1994, 114:146-159. Publisher Full Text

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

27. Trickey WR, Lee GM, Guilak F: Viscoelastic properties of chondrocytes from normal and osteoarthritic human cartilage.

    J Orthop Res 2000, 18(6):891-898. PubMed Abstract | Publisher Full Text

28. Leipzig ND, Athanasiou KA: Unconfined creep compression of chondrocytes.

    J Biomech 2005, 38:77-85. PubMed Abstract | Publisher Full Text

29. Effler JC, Kee YS, Berk JM, Tran MN, Iglesias PA, Robinson DN: Mitosis-specific mechanosensing and contractile-protein redistribution control cell shape.

    Curr Biol 2006, 16(19):1962-1967. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

30. Hochmuth RM: Micropipette aspiration of living cells.

    J Biomech 2000, 33:15-22. PubMed Abstract | Publisher Full Text

31. Dai J, Ting-Beall HP, Hochmuth RM, Sheetz MP, Titus MA: Myosin I contributes to the generation of resting cortical tension.

    Biophys J 1999, 77(2):1168-1176. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

32. 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):5275-5286. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

33. 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):9088-9095. PubMed Abstract | Publisher Full Text

34. Diez S, Gerisch G, Anderson K, Müller-Taubenberger A, Bretschneider T: Subsecond reorganization of the actin network in cell motility and chemotaxis.

    Proc Natl Acad Sci USA 2005, 102(21):7601-7606. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

35. Theret D, Levesque M, Sato M, Nerem R, Wheeler L: The application of a homogeneous half-space model in the analysis of endothelial cell micropipette measurements.

    J Biomech Eng 1988, 110:190-199. PubMed Abstract

36. Herant M, Heinrich V, Dembo M: Mechanics of neutrophil phagocytosis: behavior of the cortical tension.

    J Cell Sci 2005, 118(Pt 9):1789-1797. PubMed Abstract | Publisher Full Text

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

38. Zhao H, Chan T, Merriman B, Osher S: A variational level set approach to multiphase motion.

    J Comput Phys 1996, 127:179-195. Publisher Full Text

39. Mitchell IM: A toolbox of level set methods.

    UBC Department of Computer Science Technical Report TR-2007-11 2007.

40. Manahan CL, Iglesias PA, Long Y, Devreotes PN: Chemoattractant signaling in Dictyostelium discoideum.

    Annu Rev Cell Dev Biol 2004, 20:223-53. PubMed Abstract | Publisher Full Text

41. 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)P₃ response in Dictyostelium cells.

    Biophys J 2004, 87(6):3764-3774. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

42. 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):603-14. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

43. van Haastert PJ, Keizer-Gunnink I, Kortholt A: Essential role of PI3-kinase and phospholipase A2 in Dictyostelium discoideum chemotaxis.

    J Cell Biol 2007, 177(5):809-16. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

44. Iijima M, Devreotes P: Tumor suppressor PTEN mediates sensing of chemoattractant gradients.

    Cell 2002, 109(5):599-610. PubMed Abstract | Publisher Full Text

45. 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.

    Proc Natl Acad Sci USA 2004, 101(24):8951-6. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

46. Ananthakrishnan R, Ehrlicher A: The forces behind cell movement.

    Int J Biol Sci 2007, 3:303-317. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

47. 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.

    Eukaryot Cell 2004, 3(3):646-662. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

48. Asano Y, Mizuno T, Kon T, Nagasaki A, Sutoh K, Uyeda TQP: Keratocyte-like locomotion in amiB-null Dictyostelium cells.

    Cell Motil Cytoskeleton 2004, 59:17-27. PubMed Abstract | Publisher Full Text

49. Levine H, Kessler D, Rappel W: Directional sensing in eukaryotic chemotaxis: a balanced inactivation model.

    Proc Natl Acad Sci USA 2006, 103:9761-9766. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

50. Robinson D, Cavet G, Warrick H, Spudich J: Quantitation of the distribution and flux of myosin-II during cytokinesis.

    BMC Cell Biol 2002, 3:4. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text

Have something to say? Post a comment on this article!