- Methodology article
- Open Access
- Published:

# Modeling cellular deformations using the level set formalism

*BMC Systems Biology*
**volume 2**, Article number: 68 (2008)

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

where *S* identifies the area occupied by the cell and *d*(**x**, Γ) is the distance of position **x** to the curve Γ; see Fig. 1E for an example of a cell shape embedded in a potential function derived from the signed distance function.

We now manipulate Γ(*t*) implicitly through the function *ϕ*(**x**, *t*) according to the equation:

The vector **v**(**x**, *t*) is the velocity of the level set moving in the outward normal direction. In our case, **v**(**x**, *t*) intrinsically describes the cell's membrane protrusion and retraction velocities. These velocities can be driven by externally applied forces on the cell membrane (e.g. from a micropipette aspirator), or internally generated mechanical forces (e.g. actin polymerization or myosin-II retraction), or both. To determine how these forces translate to membrane velocity, however, first requires a mechanical model of the cell.

As the potential function corresponding to the cell shape is evolved, it can become quite steep or flat. To reduce the numerical errors caused by these effects, we re-initialize the potential function periodically [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^{3}) and a viscous element *τ*_{
c
}(nNs/*μ* m^{3}). The latter describes the association and dissociation dynamics of the cross-linkers. 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.

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^{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 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
}.

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

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 cross-linkers, 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:

*x*_{
c
}= *x*_{
m
}- ℓ - *w*_{0}.

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

Using the Voigt-dashpot model of Eqn. 5 to account for the viscoelastic response of the cell to an applied step pressure, the aspirated cellular length into the pipette, *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 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^{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 (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**

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*)^{2}*dx*.

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.

##### 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**.

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 ${\widehat{V}}_{\text{actual}}$ = *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 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}.

The formulation of ${\dot{x}}_{m}$ in Eqn. 5 provides us with the pressure-velocity 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**($\overline{x}$) at the membrane location $\overline{x}$ 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: $\overline{\varphi (t+\Delta t)}$ - *ϕ*(*t*), which translates (see Additional file 1) to an effective velocity that is normal to the cell membrane:

Thus, wherever clipping by the micropipette mask occurs, we must use this effective velocity to evolve the potential functions in simulation.

#### Evolution of the viscoelastic state of the cell

In our simulations, the cell can be represented by a series of parallel viscoelastic systems with the same parameters (Fig. 3C). These sub-systems are not interconnected, and the applied pressure on each system, *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, [*D* **l**]**v** represents the displacement of the whole cell membrane, and $\frac{\partial l}{\partial t}=\dot{\ell}n$ 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 desk-top computer.

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 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*^{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 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_{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 reaction-diffusion equations describing the PI(3,4,5)P_{3} model, we simplify it by using a steady-state 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]:

PI(3,4,5)P_{3} ∝ [cAMP/mean(cAMP)]^{3}.

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 myosin-II 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 2-D area of the cell:

**P**_{area} = *K*_{area}(*A*_{0} - *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**.

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:

**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 Γ_{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, $\frac{dl}{dt}=0$. Therefore, from Eqn. 5, we compute the viscoelastic steady state ℓ = *P*_{total}/*k*_{
c
}.

Moreover, the membrane speed at steady state is expressed as ${\dot{x}}_{m}$ = *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 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_{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 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]).

## References

- 1.
Devreotes P, Janetopoulos C: Eukaryotic chemotaxis: distinctions between directional sensing and polarization. J Biol Chem. 2003, 278 (23): 20445-20448. 10.1074/jbc.R300010200

- 2.
Reichl EM, Effler JC, Robinson DN: The stress and strain of cytokinesis. Trends Cell Biol. 2005, 15 (4): 200-206. 10.1016/j.tcb.2005.02.004

- 3.
Iglesias P, Devreotes P: Navigating through models of chemotaxis. Curr Opin Cell Biol. 2008, 20: 35-40. 10.1016/j.ceb.2007.11.011

- 4.
Meyers J, Craig J, Odde D: Potential for control of signaling pathways via cell size and shape. Curr Biol. 2006, 16: 1685-1693. 10.1016/j.cub.2006.07.056

- 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. 10.1115/1.1933997

- 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. 10.1007/s10439-005-2506-3

- 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. 10.1016/j.jbiomech.2003.11.028

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

- 9.
Bottino DC: Modeling viscoelastic networks and cell deformation in the context of the immersed boundary method. J Comput Phys. 1998, 147: 86-113. 10.1006/jcph.1998.6074.

- 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. 10.1529/biophysj.104.051029

- 11.
Graner F, Glazier J: Simulation of biological cell sorting using a two-dimensional extended Potts model. Phys Rev Lett. 1992, 69: 2013-2016. 10.1103/PhysRevLett.69.2013

- 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. 10.1007/s11538-006-9131-7

- 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. 10.1016/j.jtbi.2006.09.027

- 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. 10.1137/04060370X.

- 15.
Bray D: Cell Movements: From Molecules to Motility. 2000, Garland, 2

- 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. 10.1126/science.1095087

- 17.
Wachsstock DH, Schwarz WH, Pollard TD: Cross-linker dynamics determine the mechanical properties of actin gels. Biophys J. 1994, 66: 801-809.

- 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. 10.1016/j.cub.2008.02.056

- 19.
Osher S, Sethian JA: Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations. J Comput Phys. 1988, 79: 12-49. 10.1016/0021-9991(88)90002-2.

- 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. 10.1109/ICIP.2006.312465.

- 21.
Machacek M, Danuser G: Morphodynamic profiling of protrusion phenotypes. Biophys J. 2006, 90 (4): 1439-52. 10.1529/biophysj.105.070383

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

- 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. 10.1016/j.commatsci.2006.10.018.

- 24.
Mulder W, Osher S, Sethian JA: Computing interface motion in compressible gas dynamics. J Comput Phys. 1992, 100 (2): 209-228. 10.1016/0021-9991(92)90229-R.

- 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. 10.1006/jcph.1994.1155.

- 26.
Fung YC: Biomechanics: Mechanical Properties of Living Tissues. 1993, Springer, 2

- 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. 10.1002/jor.1100180607

- 28.
Leipzig ND, Athanasiou KA: Unconfined creep compression of chondrocytes. J Biomech. 2005, 38: 77-85.

- 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. 10.1016/j.cub.2006.08.027

- 30.
Hochmuth RM: Micropipette aspiration of living cells. J Biomech. 2000, 33: 15-22. 10.1016/S0021-9290(99)00175-X

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

- 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. 10.1091/mbc.E06-08-0767

- 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. 10.1074/jbc.M112144200

- 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. 10.1073/pnas.0408546102

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

- 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. 10.1242/jcs.02275

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

- 38.
Zhao H, Chan T, Merriman B, Osher S: A variational level set approach to multiphase motion. J Comput Phys. 1996, 127: 179-195. 10.1006/jcph.1996.0167.

- 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. 10.1146/annurev.cellbio.20.011303.132633 - 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

_{3}response in*Dictyostelium*cells. Biophys J. 2004, 87 (6): 3764-3774. 10.1529/biophysj.104.045484 - 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. 10.1016/j.devcel.2007.03.005

- 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. 10.1083/jcb.200701134 - 44.
Iijima M, Devreotes P: Tumor suppressor PTEN mediates sensing of chemoattractant gradients. Cell. 2002, 109 (5): 599-610. 10.1016/S0092-8674(02)00745-6

- 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. 10.1073/pnas.0402152101

- 46.
Ananthakrishnan R, Ehrlicher A: The forces behind cell movement. Int J Biol Sci. 2007, 3: 303-317.

- 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. 10.1128/EC.3.3.646-662.2004 - 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. 10.1002/cm.20015 - 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. 10.1073/pnas.0601302103

- 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- 10.1186/1471-2121-3-4

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

## Author information

### Affiliations

### Corresponding author

## Additional information

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

Liu Yang, Janet C Effler contributed equally to this work.

## Electronic supplementary material

## Authors’ original submitted files for images

Below are the links to the authors’ original submitted files for images.

## Rights and permissions

**Open Access**
This article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License (
https://creativecommons.org/licenses/by/2.0
), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## About this article

### Cite this article

Yang, L., Effler, J.C., Kutscher, B.L. *et al.* Modeling cellular deformations using the level set formalism.
*BMC Syst Biol* **2, **68 (2008). https://doi.org/10.1186/1752-0509-2-68

Received:

Accepted:

Published:

DOI: https://doi.org/10.1186/1752-0509-2-68

### Keywords

- Pressure Profile
- Viscoelastic Model
- Immerse Boundary Method
- Signed Distance Function
- Cell Shape Change