|Home | About | Journals | Submit | Contact Us | Français|
Mathematical and computational modeling is rapidly becoming an essential research technique complementing traditional experimental biological methods. However, lack of standard modeling methods, difficulties of translating biological phenomena into mathematical language, and differences in biological and mathematical mentalities continue to hinder the scientific progress. Here we focus on one area—cell motility—characterized by an unusually high modeling activity, largely due to a vast amount of quantitative, biophysical data, ‘modular’ character of motility, and pioneering vision of the area’s experimental leaders. In this review, after brief introduction to biology of cell movements, we discuss quantitative models of actin dynamics, protrusion, adhesion, contraction, and cell shape and movement that made an impact on the process of biological discovery. We also comment on modeling approaches and open questions.
Cell motility has occupied minds of thousands of scientists ever since 1675 when van Leeuwenhoek discovered ‘pleasing and nimble’ motions of the tiny creatures in rainwater that put forth little ‘horns’, extended and contracted . Three centuries of research did not diminish the significance of this discovery: these extensions and contractions constitute main steps of the cell crawling [1, Fig. 1A]: first, the cell pushes out the front, then it assembles tight adhesions to the surface at the leading edge and weakens such adhesions at the rear, and finally the cell develops contractions that pull up the weakly adherent rear toward the strongly adherent front, completing the motility cycle. This process is an important part of wound healing, morphogenesis and cancer, among many other biological and medical phenomena [13,36], but it is the elegance of the seemingly simple, yet underlined by layers of complexity, motile cycle that inspired thousands research papers in the last four decades.
The devil, of course, is in the details, and it is these increasingly meticulous and numerical molecular details  of the motile machinery of the cell that occupied biologists’ minds. Here, we will discuss the role, methods and impact of mathematical modeling of cell motility that started roughly 25 years ago and continues to grow. The modeling emerged from the need to understand actin treadmill, which is at the heart of the cell movement [88, Fig. 1A].
Actin is a globular protein that polymerizes into a linear two-stranded right-handed double helix twisting around itself every 37 nm. Each monomer (G-actin) is ~5.4 nm in size, so that assembly of a monomer imparts a 2.7-nm axial rise to the filament (F-actin). Both strands of the helix have the same polarity, so the two ends of the filament are structurally different, termed barbed and pointed ends. The barbed end assembles and disassembles monomers two orders of magnitude faster than the pointed end. In the absence of nucleotide hydrolysis, the critical concentration of G-actin—when the rates of polymerization and depolymerization balance—is the same at both ends. At this concentration, the average length and position of a filament do not change. However, actin monomers bind ATP, and the filaments they polymerize into become dynamically asymmetric: the barbed end is ATP-capped; ATP inside the filament is hydrolyzed, and the phosphate group released, so there is ADP-actin, with a tendency to disassemble, at the filament’s pointed end. As a result of this asymmetry, the effective affinity for new monomers at the barbed end is high and the critical concentration is low, ~0.2 μM. On the pointed end, the monomer affinity is low so that the corresponding critical concentration is higher, ~0.7 μM. The consequence of this asymmetry is the non-equilibrium process of treadmilling: net depolymerization from the pointed end balanced by net polymerization onto the barbed end without changing its average length. In the cytoplasm, the monomers release ADP and bind ATP and recycle by diffusion. If such filament adheres to the surface and is enveloped by the plasma membrane, within which ADP-ATP re-charge is taking place, then the actin treadmill can serve as the basic cell motility mechanism (Fig. 1a).
First modeling efforts [46,79,118] were directed at quantifying this actin treadmill and using thermodynamics to understand the nature and magnitude of the polymerization force. These early works introduced fundamental ideas that are still used in developing increasingly complex models. Another crucial step in quantifying motility was the measurement of the rates of actin polymerization and disassembly  that provided parameter values for the previously conceptual models. The modelers realized very early that cell motility is the process that involves a combination of chemical kinetics, transport and physical forces and movements of a polymer network interacting with a vast number of other proteins, and so the problem can be formulated and treated mathematically by combining models of continuum mechanics and biochemical kinetics [34,80,123].
Here, we review the modern state of the field of mathematical and computational cell motility. Despite it being specialized, it is already so vast that we had to limit ourselves to relatively few examples, the choice of which is no doubt skewed by the author’s tastes. Readers who would like wider and more balanced picture will have to read other reviews, from early ones [25,67,105] to more recent [18,59,88,93], as well as textbooks [13,47,89]. We start with introduction into the relatively well understood lamellipodial motility and review modeling of three main processes—protrusion, contraction and adhesion. We then discuss integrative models of whole motile cells, and after a survey of related models of cell polarization and other motile systems, conclude with discussion of open questions and modeling methodology. For the benefit of the reader who might be interested in getting involved into cell motility modeling, we use the bold question sign (?) to mark the open questions.
Actin treadmill described above illustrates the principle of locomotion but does not explain how thousands of actin filaments self-organize into motile networks. This self-organization was unraveled by a number of studies that culminated in dendritic-nucleation model  and myosin-powered network contraction model  that together explained lamellipodial motility. Not all cells use lamellipodia as motile appendages, but many do, and it is the best understood motile system that attracted most of the modeling effort. In this review, we focus on this mode of locomotion and discuss two other types in the last section. Reader interested in amoeboid mode of locomotion would have to look up the reviews [36,37,45] and references therein.
Lamellipod, the motile engine of the cell, is its front, advancing part. It is a flat leaf-like extension filled with a dense rectangular actin network  It is only a few tenths of a micron thick but is tens of microns long and wide (Fig. 1a, b). The cell body, containing the cell nucleus and other organelles, appears to be a mechanically passive structure pulled forward by the lamellipodial action. According to the dendritic-nucleation hypothesis, nascent actin filaments branch from the sides of existing filaments in an angularly precise way, so that there is ~70° angle between the ‘mother’ and ‘daughter’ filaments, and all leading edge filaments’ barbed ends are oriented forward at roughly ±35° to the direction of protrusion (Fig. 1b). The growing barbed ends at the front elongate ~0.1 μm/sec and push out the leading edge membrane, but are capped within seconds, so individual filaments grow only to 0.1–1 μm. As the filaments adhere to the surface, they lag behind the leading after capping edge being replaced at the front by the next generation of filaments, so the growing barbed ends abut the extreme leading edge of the lamellipod while the disassembling pointed ends dominate far from the leading edge.
Thousands of myosin molecules are distributed throughout the cell. These molecules are motor proteins that use ATP hydrolysis energy to glide toward actin barbed ends, and if stalled, develop a pN-range force sliding a filament with its pointed end forward (Fig. 5). This sliding action is the key to muscle contraction, however, in the lamellipod, there is no crystalline telescopic arrangement of actin-myosin sarcomeres designed for macroscopic contractions. The branching adhesive network of short filaments constitutes a relatively rigid gel (the network is further reinforced by crosslinking proteins), which myosin motors are unable to perturb. However, actin network disassembles throughout the lamellipod and weakens mechanically toward the rear (Fig. 1b). At the same time, in the coordinate system of the moving cell, myosin molecules attached to the F-actin network are effectively ‘swept’ to the rear and concentrate there. This allows dense myosin clusters at the cell rear to generate contractile stresses and to collapse the largely isotropic lamellipodial actin network into a bipolar actin-myosin bundle at the very rear of the lamellipod (Fig. 1b). Subsequent muscle-like sliding contraction of the bundle pulls lamellipodial actin filaments into the bundle, advancing the rear boundary of the lamellipod forward. At the same time, adhesion complexes assembled at the front ‘age’ and largely disassemble at the rear, allowing the myosin-powered contraction to pull the bulk of the cell forward.
A big part of the dendritic-nucleation model is a beautiful molecular machinery of the lamellipodial treadmill: first, active GTPases and phosphoinositides activate WASP/Scar proteins, which in turn activate Arp2/3 protein complexes nucleating nascent actin filaments at the side of the existent filaments. Proteins of the ADF/cofilin family accelerate F-actin disassembly across the lamellipod. G-actin produced in the process of the disassembly is sequestered by ADF/cofilin, but then rapid exchange with two other important actin binding proteins, profilin and thymosin, ensues, as well as ADP-ATP exchange. Actin-thymosin complexes cannot polymerize and are maintained by the cell as a ‘backup’ pool, while actin-profilin complexes are able to assemble onto the uncapped barbed ends. G-actin is delivered from the depolymerization sites to the leading edge by diffusion. The important point of the model is that the lamellipodial network advances by array treadmilling, rather than by treadmilling of individual filaments: barbed ends of a few filaments grow at the leading edge, while their pointed ends are stable; other filaments have both ends capped; yet others have the barbed ends capped and the pointed ends depolymerizing.
This combined dendritic-nucleation and network contraction hypotheses presented an unusually complete qualitative scheme begging for physical and mathematical scrutiny. Another unusual feature of this model was availability of many parameter values such as reaction rates, concentrations, diffusion coefficients, etc. , largely due to the leadership of T. Pollard who tirelessly promoted quantitative and modeling approaches in cell biology for decades. These hypotheses pose a number of questions that required combined modeling/experimental approaches. These are:
In what follows, we review the modeling efforts accompanying experimental studies to answer these questions.
Protrusion is the best researched and understood of all locomotion steps. Especially, force of protrusion was modeled vigorously following pioneering ideas of T. Hill. Modeling protrusion was recently reviewed; the reader is referred to  for more detail. ‘Polymerization Brownian ratchet’ model  accounts for the force generated by polymerization when the polymers are rigid. According to this model, an object in front of the filament diffuses away from the tip, creating a gap sufficient for monomers to intercalate and assemble onto the tip, thereby inhibiting the object from diffusing backward (Fig. 2). Even when a resisting load force is applied to the object, Brownian motion can still create a sufficient gap, and so the object’s diffusion is biased forward. If the polymer is not rigid, its own thermal undulations can create a gap between its tip and the load . If the gap is large enough and persists long enough, a monomer can intercalate into the gap and bind onto the tip of the growing actin filament. This increases the filament’s length so that when the tip contacts the load the polymer is bent and the resulting elastic force pushes on the load. In general, both filament and load (cell membrane) are fluctuating, and the combined effect can be sufficient to allow intercalation of monomers and consequent force generation.
These ‘ratchet’ models used differential-difference equations of the following type (Fig. 2):
where p (x, t) is the probability density function for the gap of width x, D is the effective diffusion coefficient, f is the load force opposing polymerization, kB T ≈ 4 pN·nm is the thermal energy, δ is the half-monomer size, koff ≈ 1/sec is the monomer disassembly rate, kon ≈ 10/(μM sec) is the assembly rate, a ~ 10 μM is the G-actin concentration. In this equation, the first term describes effective thermal diffusion of the bending filament and membrane in front of it, the second one is responsible for the force (membrane resistance)-generated drift closing the gap, while the difference terms account for the ‘jumps’ of the gap size after the assembly/disassembly acts. Characteristic parameter values δ ≈ 3 nm, D ~ 105 nm2/sec are such that the thermal diffusion is much faster than the relatively slow actin kinetics (D/δ2 ~ 104 sec kona ~ 100/sec), so singular perturbation theory allows finding the quasi-steady gap probability distribution p ~ exp [−f δ/kB T], and then to derive the force–velocity relation:
where V is the net filament growth rate. One important conclusion from this relation is the estimate of the force that stalls the filament growth: if V = 0, then
This formula predicts that one filament can generate force in the pN range: kB T/δ ~ 1pN, while the logarithmic factor ~ 5 is insensitive to exact parameter values. Experimental measurements [35,50], indeed, report the pN range forces. However, one has to be cautious in claiming that these measurements support the ratchet theories: Kovar and Pollard  used filaments with protein formin at the growing barbed end that alters the polymerization process, while Footer et al.  used a filament bundle, in which the details of the force distribution among dynamic filaments is still unclear.
These theories made two other simple predictions: to grow against membrane resistance, the filaments should be neither too short (filaments shorter than ~70 nm are too rigid to bend enough), nor too long (filaments longer than 500 nm are too soft, so the load would simply buckle them). Also, the optimal growth of the filaments is achieved when they subtend the leading edge at an angle: filaments parallel to the leading edge obviously can grow fast, albeit without causing protrusion, while those normal to the leading edge are effectively too rigid to bend enough. These two predictions were instrumental for formulating the dendritic-nucleation model, when the experimentalists realized that the capping limits filaments’ lengths to tenths of a micron, and Arp2/3-mediated branching keeps the filaments at ±35° to the direction of protrusion. Recently, ratchet models underwent further development when it was suggested, following earlier theoretical work , that some of actin filaments are attached to the surface they push . The ‘tethered ratchet’ model explains the mechanism in this case by assuming that the filaments attach to the surface transiently, dissociating fast and growing freely until getting capped and losing contact with the surface altogether. During this process, the attached filaments are in tension generating effective resistance and friction, while the dissociated filaments are in compression, and generate the pushing force.
One alternative model of protrusion proposes that all filaments are attached to the surface: all pushing barbed ends are clamped to an end-tracking protein (peculiar stepping motor coupling protrusion to ATP hydrolysis) associated with the surface . This model would be able to explain forces larger and a variety of force-velocity relations wider than those predicted by the ratchet models. Another influential model draws attention to elastic shear stress developed by recoil of the actin meshwork generated by growth and pushing force at the boundary . This model was developed for geometry different from that of the lamellipod, and it is based on clever scaling estimates. Derivation and mathematical and numerical analysis of a meso-scale model of viscoelastic (see below) growing actin network with microscopic force generators at the boundary (?) is an urgent and potentially very rewarding task.
An important aspect of the protrusion process is the self-organization of the dendritic actin network. Its angular organization was modeled in a very elegant paper . The authors assumed that the rate of capping, c (), depends on angle between the growing filament and direction of protrusion (Fig. 3a) because filaments more normal to the leading edge bend away from it less and are more protected from capping by some molecular complexes at the leading edge. Furthermore, as discussed above, if a ‘mother’ filament growing at angle ′ sprouts a ‘daughter’ filament, then the latter grows on average at either angle = ′ + ψ, or = ′ − ψ, where ψ = 70°. Introducing the angle-dependent rate of branching K (θ) with a peak at ψ = 70°, the integral equation describing the filament angular distribution n (, t) was derived:
Solution of this equation revealed two peaks at ±35° (also confirmed experimentally) important for effective protrusion, as predicted by the ratchet models. Arp2/3 mediated branching imposes a 70° angle between the peaks, obviously, but this does not explain the symmetric ±35° orientation of the filaments relative to the leading edge. The model explains this symmetry mathematically demonstrating that the angularly symmetric mother-daughter filament pairs ‘survive’, while the asymmetric (relative to the leading edge) pairs do not (Fig. 3a). At the same time, this model provides a plausible explanation for the actin polarization: barbed ends growing away from the leading edge get capped fast, while those growing forward do not.
Actin filaments self-organize not only in angular space, but also in physical space : filament density along the leading edge looks like an inverted parabola (Fig. 3b), with density high at the center and low at the sides. A simple explanation of this phenomenon  is the balance of capping, branching and lateral flow . The latter is effective sliding of the barbed ends laterally along the advancing leading edge (Fig. 3b). These three processes are described by respective terms in the equation for the filament spatial distribution n (x, t):
Here index ± corresponds to right-/left-oriented filaments. Solution of this equation predicts the data very well . A very promising open modeling project is to derive a combined spatial-angular model (using integro-PDEs) of the leading edge actin dynamics (?) without some of severe simplifications that were used earlier.
All models described above rely on solving differential equations for deterministic continuous densities. This classical method of applied mathematics can go only that far in cell biology, where systems often consist of thousands of molecules and so are not microscopic, but they are not macroscopic either: fluctuations of chemical or physical quantities in the cell are comparable in magnitude to the average values of these quantities. Therefore, the discrete nature of the molecular components cannot be neglected, so computational methods taking into account stochastic effects, albeit missing some of the insight and mathematical elegance, must be used. Such Monte Carlo-like methods are usually simulations of agent-based models, in which each molecule is represented explicitly and interactions between them obey simple rules. Excellent recent example of such approach, pioneered in , is the recent studies [99,100] that simulated the lamellipodial dendritic network growing against resistive elastic membrane. This model confirmed many predictions of models described above without simplifications necessary to use continuous math tools.
Recently, all protrusion models are facing challenges from advanced biophysical experiments measuring the force-velocity relations for growing actin networks. The results [35,84,91] are not easy to interpret in terms of the simple models (with some notable exceptions when results are explained by calculations published in the same paper ), and often with other experiments, most likely because the existing models are not advanced enough and do not take into account properly the force-dependent network self-organization. One notable, albeit limited, exception is the autocatalytic branching theory  that assumes that the rate of the filament branching is proportional to the density of the existent leading edge filaments. An unexpected conclusion of this theory is that the protrusion rate does not depend on the load: effectively, greater load force causes faster branching, and therefore greater actin density, so the load per filament remains constant leaving the growth rate unchanged. A very important open problem is to develop a comprehensive model of coupled actin network self-organization and force generation (?) that would explain fascinating history-dependent effects .
Besides the force generation, other important factor of protrusion is the rate of recycling of G-actin to the leading edge where it assembles onto the growing barbed ends. Note that Eq. (2) for the actin growth rate shows that the protrusion rate is limited not only by the membrane resistance, f, but also by the G-actin concentration at the front, a. Note also that the rate of the single actin filament treadmill is equal to koff δ~ 3 nm/sec, which is two orders of magnitude less than the observed protrusion rate of rapid cells. The answer to this paradox was suggested by the pioneering funneling model , which suggested that in the steadily treadmilling lamellipodial array, the barbed-end capping produces a shortage of the number of growing barbed ends, which is N ~ 100 times less than that of the shortening pointed ends (Fig. 4). This keeps the supply of monomers plentiful, and accelerates growth of any barbed ends that are temporarily uncapped to ~ N koff δ ~ 0.3 μm/sec, as observed. Another issue, however, is that the G-actin disassembling away from the leading edge has to diffuse to the front, and also go through the chain of nucleotide- and accessory protein-exchanging reactions. Respective diffusion–reaction equations were solved in [69,72,95], and it was shown that diffusion is fast enough to move actin monomers across the lamellipod. These models also made an interesting, supported indirectly but still unconfirmed directly, prediction that the number of the leading barbed ends has to be fine-tuned for faster protrusion: too few filaments are too weak to push, while too many deplete the G-actin pool.
Besides diffusion, convective D’Arcy flow of the cytoplasm that can be caused by myosin-driven pressure gradients [20,95] can also transport G-actin. This flow direction and magnitude would drastically depend on cell membrane permeability to water, which is unknown. Thorough model of the porous flow of cytoplasm in the lamellipod is suggested in , yet further work in the direction of developing and testing 3D models is needed (?).
Another piece of the actin recycling puzzle has to do with coupling ATP hydrolysis with the actin treadmill. Recent theories of such coupling in a single filament [9,111,106] modeled interactions of each individual monomer as they exchange on and off the ends and undergo ATPase reactions By analyzing nucleotide profiles within the filament, Bindschadler et al.  predicted that a combination of enhancing phosphate release, increasing the off rate of ADP-bound subunits at pointed ends, and optimal level of capping and are necessary to accelerate the treadmilling to rates observed in vivo. Extension of such model from the filament treadmill to the actin array treadmill (?) is necessary.
Finally, let us mention that in some cells, the lamellipodial leading edge is interspersed with filopodia—bundles of actin filaments that are packed tightly together and protrude forward. Recent evidence points out that the lamellipodial filaments can bend together and ‘zipper’ into such parallel bundles . Some preliminary modeling of both emergence and maintenance of filopodia is done [6,43,72,99], but comprehensive model of this phenomenon (?) is lacking. Such model can assist answering interesting question whether filopodia are ‘guiding’ devices probing space ahead of the lamellipod, or mechanical devices ‘penetrating’ the environment and serving as a robust scaffold for the lamellipodial protrusion.
It is not completely clear what is the physical mechanism that pulls forward the cell body after the leading edge of the cell advances. In fact, the nature of dynamic mechanical attachment between the nucleus and organelles and actin network in the motile cell is still unknown. It is not completely out of question that no contraction is necessary at all—in principle, it is possible that the actin array treadmill inside an inextensible bag of plasma membrane is the minimal system for locomotion. Significant evidence, however, points at the actin-myosin contraction as the dominant mechanism for contraction of the cytoskeleton and forward translocation of the cell.
The qualitative myosin-powered network contraction model  is one of few studies describing respective mechanism in molecular detail, yet many mechanistic details remain murky. One of the fascinating open questions is how does contraction take place in a relatively disordered actin cytoskeleton: A pair of free actin filaments would be effectively pulled by motor heads of a myosin cluster into the configuration, in which the barbed ends of the filaments end up co-localized, while the lengths of the filaments extend outward (Fig. 5), so compaction of the filaments cannot be achieved. It is easy to imagine scenario, in which actin binding proteins exist that preferentially associate with the barbed ends and are able to bind to more than one end. Then, ‘mini-sarcomeres’, as in muscle, would appear, and interaction of those with each other mediated by myosin could lead to contraction (Fig. 5). The problem is that such binding proteins were not observed. Relevant questions were investigated mathematically with significant elegance in 1D in  where the assumptions of pair-wise actin filament interactions and additive character of myosin sliding force were made and integro-differential equations amenable to perturbation analysis and numerical solutions were derived.
Usual modeling approach to contraction in continuous models is to assume that there is an isotropic contractile stress (negative pressure) in the cytoskeleton proportional to the myosin density that has to be added to passive stresses [4,42,54]. Some discussion of possible anisotropic theory, which, if biologically important, would be of great mathematical interest, can be found in . The magnitude of the contractile stress can be estimated from the cell size and measurements of the total traction force, presumably myosin-generated, exerted by the cell on the substrate, ~104–105 pN, that compares well with the total number (~104–105) of myosin molecules, each able to generate a few pN of force [33,78].
Finding the myosin contractile stress is, in fact, easier part of the problem. Harder task is to model another aspect of contraction—mechanical properties of the actin-myosin meshwork that determine the rates and character of its deformations in response to the myosin-generated stress. A good introduction into relevant mechanical behavior of the cytoskeleton can be found in . A great deal of experiments and modeling is done on elasticity of in vitro actin gels. Respective modeling involves formulation and minimization of a free energy functional and requires very delicate estimates and knowledge of statistical physics. There were very few numerical studies; one of the most interesting and clear ones is . Rheology of the cytoskeleton (?) is ripe for mathematical examination. Note, that length distribution of actin filaments, which is regulated in the cell by a host of accessory proteins, is an important factor for mechanics. The length distribution was modeled mathematically in .
Physical estimates posited that the mechanical properties of such gel have a complex nature determined, in different regimes, by either enthalpic (bending of the filaments), or entropic (thermal undulations of the filaments) effects. These physical theories estimated the magnitude of Young modulus, E (main characteristic of the actin gel elasticity), of the order of
where α ~ 2, β ~ 5. In this formula, kB T ≈ 4 pN · nm is the thermal energy, lp ~ 10 μm is the so called persistence length of actin filament, and ζ ~ 0.1 μm is the mesh size in the actin gel  so E ~ 104 pN/μm2. Elasticity of the keratocyte lamellipod of this order was actually measured using atomic force microscopy ; however, this estimate is very sensitive to the concentration and stiffness of actin filaments, so that a few fold fluctuations of these parameters could lead to orders of magnitude variance in mechanical properties. Recently, a number of high profile combined experimental and theoretical papers were published [22,38,68] that examined effects of crosslinking and myosin contractions on the actin network elasticity. Very complex nonlinear response of the elastic moduli to deformations were discovered including either softening or hardening regimes. Incorporation of corresponding constitutive relations into the continuous mechanical models of the motile cells (?) is an interesting mathematical problem. Interesting novel theory of dynamics of ‘active’ actin gel is suggested in 
Finally, let us mention that the actin network exhibit viscous, as well as elastic behavior . Respective viscoelastic properties that can be approximated by a combination of Maxwell and Kelvin–Voight models are very complex, nonlinear, and sensitive to gel parameters. Engineering ‘spring-and-dashpot’ modeling and continuous models of the viscoelastic actin cytoskeleton can be gleaned from  to [42,54], respectively. Numerical analysis approaches to modeling the cytoskeleton rheology, such as immersed boundary method , and reactive interpenetrating viscous flows  are a powerful start for future predictive models.
Adhesion is perhaps the hottest subfield of the cell motility these days, with tens of biological papers published on this topic every month. Mathematical biologists who would be well advised to get in on the ground floor should start with reviews [8,117]. The property of adhesion crucial for cell locomotion is being strong at the leading edge and weak at the trailing edge—without it, contraction and protrusion would just uselessly recycle cytoskeleton without forward translocation of the cell. This graded adhesion is achieved as a result of a very complex and incompletely understood dynamics of many molecules: first, integrins—transmembrane proteins that are key ingredients of adhesion molecular machinery—associate with the substrate at the leading edge, then a vast host of proteins—talin, vinculin, paxillin, FAK, alpha-actinin—to name but a few, form adhesion complexes ultimately connecting actin network with the substrate (Fig. 6). Some molecules in these complexes are signaling, other play mechanical roles. Altogether, they probably constitute a Velcro-like mesh of weak and transient links that are both malleable and resilient. The cell assembles tens of different types of adhesions in response to various environments ; it is possible that the reason for the great variety of proteins making up the adhesion complexes is the need to diversify the adhesion type.
How the graded adhesion is maintained in the moving cell is incompletely understood. This maintenance has something to do with rapid hierarchical assembly of integrin, talin, vinculin and other adhesion molecules at the leading edge and slow subsequent ‘aging’ of the adhesion complexes . Because of the cell translocation, this dynamics lead to weakening of the adhesions by the time the cell rear ‘catches up’ with them (Fig. 6). Spatial gradients of regulatory molecules in the cell, such as calpain ), are other important factors.
Thermodynamics and kinetics models of the adhesion complexes’ assembly and turnover started long time ago [27,82,116]. The most important achievement of these models was theoretical discovery confirmed by experiment that the cell speed has a biphasic dependence on the adhesion strength [30,81]: when the adhesion is too strong, the cell cannot retract its rear, when it is too weak, the front is pulled back by contraction negating the protrusion. A poignant modeling problem is unraveling signaling-dependent adhesion assembly on the molecular level. One notable nascent modeling effort is .
One of the most fascinating properties of adhesion is its complex response to force: in many circumstances, adhesion assembles and becomes stronger with increasing force . Three conceptual qualitative models of this phenomenon are discussed in ; quantitative physical models of some aspects of this phenomenon, including dependence of the adhesion dynamics on softness of the substratum can be gleaned from [76,101]. Pioneering thermodynamics model of the force-dependent adhesion growth  is based on the idea that application of the pulling force results in stretching of the adhesion aggregate, consequent accumulation of elastic stresses within the aggregate, and the assumption that insertion of new proteins into the aggregate results in stress relaxation.
Finally, this mechanosensitivity of the adhesion complexes can be the reason for the interesting self-organization of adhesions and actin-myosin stress fibers: adhesions inside the cell are pulled into all directions by the myosin contraction, so corresponding forces cancel and such adhesions remain weak. Meanwhile, adhesions at the cell boundary are pulled only inward and grow under the influence of the contractile force. Integro-differential mathematical model of this phenomenon  predicts interesting observed effect: adhesions concentrates along curvier parts of the cell boundary.
One of the fundamental questions in cell biology is how to connect the microscopic, subcellular cytoskeleton dynamics with geometry of the whole cell, in other words, how to predict the cell shape from the knowledge about the protrusion, contraction and adhesion. Mathematically, this question can be formulated in terms of usually very difficult free boundary problem: the cell is defined as a time dependent domain Ω(t) (in 3D, or, if approximations are made, in 2D or even in 1D). Evolution of its boundary, Ω (t), is governed by protrusion-retraction boundary conditions that are in turn determined by essential densities defined by dynamic equations on Ω(t). Formulating such self-consistent problem is difficult; even harder is to solve it numerically. Developing respective software and methodology (?) is one of the most worthy and potentially rewarding sub-area of cell motility modeling, because any reasonable experimental effort in this direction will need quantitative backing. One has to be careful with moving from lower to higher dimensional models, as some physical mechanisms do not work the same way in 1D, 2D and 3D and can give rise to different behaviors.
The simplest relevant mathematical problem is determining the 1D front-to-rear cell length, ignoring the lateral and ventral-dorsal cell dimensions. Computational approach to this problem is modeling the cell as a 1D chain of springs (contractile actin-myosin units) connecting nodes (material points of actin cytoskeleton) [30, Fig. 7]). The nodes, in turn, link with the substrate through dashpots (adhesion complexes). At the front and rear of the cell there are two possibilities—either advance the leading node and pull up the rear node preserving their identity, or to mimic actin polymerization at the front and depolymerization at the rear by adding nascent nodes at the front and deleting nodes at the rear, respectively (Fig. 7).
Continuous approach to the same problem can be considered as a limit of the previous problem, in which the distance between the nodes tend to zero, and all other variables and parameters are scaled accordingly. A few recent models [42,54,70] have the following generic structure. Equation for the displacements u (x, t) of the cytoskeletal material points with origin at x and at time t has the form:
where β(x, t) is the adhesion viscous drag, and σ(x, t) is the stress in the cytoskeleton. This stress in purely elastic model has the form:
where E is the Young modulus of the actin meshwork, and C is the contractile stress. In viscoelastic models, the term of the form , or other similar one (μ is effective vis-cosity of F-actin), is added to the stress. These mechanical equations are complemented by appropriate boundary conditions, equations for F-actin, myosin and adhesion transport. Usually, these equations (for appropriate density ρ) have reaction-drift-diffusion character:
with appropriate boundary conditions. Finally, algebraic constitutive relations are formulated, for example, the contractile stress is assumed to be proportional to the myosin density, etc. Most importantly, these equations are solved on the moving segment [r (t), f (t)], where r and f, coordinates of the rear and front, respectively, are governed by the boundary conditions determining the protrusion and retraction rates, P and R, as functions of stress, distance, actin and myosin densities, etc.:
Analytical or numerical solutions of system (7–10) aim at finding asymptotically stable stationary cell length (f − r) and distributions of the stress, cytoskeleton drift and essential densities. Mathematically, for the steadily crawling cell, a solution in the form of the traveling pulse has to be found. This is usually done ad hoc, but the first rigorous study  demonstrated that such problem can be solved up to high mathematical standards. A very promising untouched problem is to investigate transient movements of the cell (?) using this 1D approach.
Ignoring the lateral dimensions of the cell is, of course, a gross simplification, so a 2D model would be much more satisfying. Most attention was paid to modeling characteristic fan-like shape of the lamellipod of the motile keratocyte. In this case, a 2D model is completely adequate, as the dorsal–ventral thickness of the lamellipod is minuscule. The first important achievement in this direction was the so-called graded radial extension model that shed light on kinematic principles underlying cell shape . This model, supported by the authors’ experiment, suggested that cell boundary extension/retraction is locally normal to the boundary, and that the corresponding local rate of extension/retraction is graded, decreasing from the center to the sides of the cell (Fig. 8). This principle results in a remarkably simple geometric formula:
where θ(x) is the orientation of the normal to the leading edge at position x, Vn (x) is the normal rate of protrusion at this position, and Vn (0) is the normal rate of protrusion at the cell center. The function Vn (x) can be calculated as a function of relevant cytoskeletal stresses and densities on the cell boundary based on modeling assumptions, and then the profile of the leading edge (and, similarly, of the rear edge) can be calculated by solving simple geometric equations.
This model did not explain what are the subcellular dynamic mechanisms underlying the proper graded extension/retraction rates. Attractive hypothesis is that at the leading edge the protrusion rate is graded because of the graded F-actin density: greater number of actin filaments at the center pushes the membrane faster at the center . Myosin in the steadily crawling cell is ‘swept’ to the rear of the cell (see explanation in the Sect. 7), where it contracts, pulling the rear and the sides of the cell inward (Fig. 8). Mathematically, the effective protrusion/retraction rate at the cell boundary is the vector sum of the actin polymerization rate and myosin powered rate of F-actin centripetal flow [108, Fig. 8]. These data and qualitative explanations were incorporated into the 2D models of the leading edge  and of the whole lamellipod  shapes. Notable other efforts include the model of the keratocyte shape that includes spatial distributions of the actin-regulatory proteins , the ‘dendritic actin network enveloped by elastic membrane’ keratocyte shape model , and the spring-and-dashpot 2D model of the nematode sperm lamellipod (; discussed more below). All these quantitative models are based on severe simplifications, and much work is needed to model the cell shape adequately.
Part of the question about the lamellipodial shape is the question about what limits its area. There is a great number of theoretical possibilities—limiting total area of plasma membrane, limiting number of important molecular complexes, such as Arp2/3, balance of differentially size-dependent protrusion and contraction rates, balance between adhesion and membrane mechanical energies, to name but a few. Data to select one of the answers is lacking.
Some general principles governing the keratocyte shape-graded protrusion/retraction determining the shape, and the balance of actin polymerization and myosin powered retraction underlining the graded protrusion/retraction—are very likely applicable to other cells. However, specific physics, biology and biochemistry vary, so different cells would have to be approached differently. Also, the roles of the cell membrane recycling and mechanics and of the adhesion dynamics in determining the cell shape remain unclear. ‘Top-down’ mathematical models that do not presume detailed biophysical mechanisms, but instead use quantifiable rules of the cell boundary protrusion and retraction, can complement ‘bottom-up’ models in exploring the principles of cell migration. For example, Satulovsky et al.  modeled the cell as a shape machine regulated by a circuit of local stimulation—global inhibition. By varying parameters, this model mimicked the shapes of motile Dictyostelium, keratocyte, neuron, and fibroblast cells with shocking faithfulness.
We discussed the cell motility based on the actin treadmill implying persistent front-back asymmetry. Natural related questions are: How is symmetry of a stationary cell broken? How is the asymmetry of the motile cell maintained? How does the cell choose the direction of locomotion?
The remarkable experiment  with lamellipodial fragments not containing nucleus or any organelles (but containing actin, myosin, ATP and all accessory proteins in the cytoplasm) is perhaps the minimal system allowing addressing the first two questions. When excised from the whole cell, the fragments spread into discoid shapes and remain stably non-motile (Fig. 9). Microscopy shows that myosin is distributed almost evenly throughout the disc, and that there is a very slow centripetal F-actin flow caused by ineffective myosin contraction. When water is blown at the side of the fragment, the fragment is deformed by the hydrodynamic load. If the flow is weak, and the deformation is not great, the fragment remains stationary and recovers its discoid shape (Fig. 9). However, if the flow is faster and is sustained for tens of seconds, the fragment deforms to an almost half-disc shape, and then, even if the flow stops, its symmetry is broken—the crescent-shaped fragment evolves, so that almost all myosin concentrates along the rear edge, actin growth takes place along the leading edge, and the fragment crawls steadily for hours, much like the whole keratocyte cell. On a rare occasion, when the motile fragment encounters an obstacle, it stalls and resumes the discoid non-motile shape. This duality in fragments’ behavior is similar to that of the whole cell with one important distinction: stationary cells become motile spontaneously, which is likely to involve complex RhoGTPase-dependent regulatory pathways.
The following qualitative scenario suggests purely mechanical explanation for the symmetry break and polarization maintenance : in the symmetric state, F-actin growth, stalled by the high membrane tension, is very slow at the edges. Myosin molecules turn over (dissociate from F-actin, diffuse and re-associate with F-actin) and spread across the whole fragment. Thinly spread myosin generates weak contraction causing the centripetal F-actin flow to be slow and ineffective in biasing the myosin distribution. When the fragment is pushed from the side, actin network at this side collapses into a bi-polar bundle, while myosin concentrates near this side creating conditions for actin-myosin contraction relieving tension in the whole fragment. F-actin starts to polymerize rapidly against the small tension along unperturbed boundary, so the fragment starts to move as a whole away from the deformed side. In the framework of the now motile fragment, myosin is swept to the rear, and the mechanical actin-myosin treadmill becomes stable. Additional adhesion dynamics could be important for this scenario that still awaits mathematical modeling.
Another possible mechanism is suggested by the ‘competition for resources’ model: if we start with two oppositely oriented F-actin arrays of equal density (the symmetric state), then it is possible that these two arrays ‘compete’ for a common pool of ‘resources’ diffusing in the cytoplasm, i.e. G-actin or Arp2/3, or other protein complexes responsible for the F-actin branching. If, as a result of a fluctuation, more than half of these resources attach to filaments from one array, then more than half daughter filaments that grow in the same direction as the mother filaments would belong to that array. In the next generation of filaments, the greater number of filaments from the randomly ‘winning’ array would get the proportional share of the resources, and sprout yet greater number of filaments. It is well known from mathematical ecology models  that such ‘competition for resources’ may cause ‘win-loose’ dynamics such that denser actin array ‘consumes’ more resources and grows further depleting ‘loosing’ array supplies and eventually taking over remaining the only stable polarized treadmill. For dendritic actin array, a 1D computational model of this scenario was developed in .
These two mechanical/transport scenarios are probably a small part of the whole picture, because extensive studies implicate complex and redundant spatial–temporal feedbacks between multiple signaling pathways and cytoskeleton in most of the motile cells [5,113], so it is likely that spatial-temporal biochemical signaling participates in initiation and maintenance of the cell treadmill. Mathematically, this more complete situation was modeled recently in [26,64], where reactions and diffusion of Rho GTP-ases and their interactions with actin/Arp2/3 system were considered. It was found that interactions of small G-proteins Cdc42, Rac and Rho are responsible for the polarity switch, and that maintenance and robustness of polarity was due to the rapid diffusion of the inactive forms of the small G-proteins in the cytoplasm. These papers are a good start for next generations of models.
Above, we discussed the self-polarization of the cells in the random direction. Most of the cells, however, polarize and move in response to an extracellular cue. This process attracted a lot of modeling recently, and we direct the reader to another review  for references to the original models. Most of these models are guided by experiments suggesting that the chemotactic response of the eukaryotic cells is based on signaling reaction-diffusion processes and respective well studied PDEs. One type of models explored Turing type mechanism [66,75] based on chemical reactions of local fast excitation, local slow inhibition and global slow inhibition triggered by a chemoattractant gradient. The most popular ‘local excitation-global inhibition’ (LEGI) model proposed that directional sensing depends on a balance between a rapid, local excitation and a slower global inhibition processes . Receptor (on the cell surface reacting with the spatially varying chemoattractant) occupancy controls the steady-state levels of each process, and the difference between the two regulates the response. As inhibition is assumed to depend on average receptor occupancy, its steady-state level is less than that of local excitation at the front of the cell. At the back, the situation is reversed. In the third type of models , the cell responds temporally to the ‘first hit’: a second messenger, generated by local activation of the membrane, diffuses rapidly (faster than the chemoattractant) through the interior of the cell, suppresses the activation of the back of the cell, and converts the temporal gradient into an initial cellular asymmetry. All these models are so far conceptual, rather than detailed, but rapidly evolving experimental techniques continue to provide necessary data, so the task of developing predictive models will soon be realistic.
Quite a number of other motile systems provide fascinating data and food for modeling. Here, we mention only two examples. Many intracellular pathogens are propelled by the polymerization of actin filaments serving as simplified model systems for eukaryotic cell motility: a purely protrusive system without adhesion, contraction, and cell membrane. In particular, Listeria propels itself by assembling the host cell actin into a comet-like tail of cross-linked filaments, with their polymerizing barbed ends oriented toward the bacterial surface [14, Fig. 10]. Only one protein on the surface of Listeria, ActA, is required for motility. In addition to actin, ActA and ATP, only a handful of proteins in cytoplasmic extracts are essential for bacterial propulsion, including Arp2/3, capping proteins and the depolymerizing factor, ADF/cofilin .
Listeria propulsion attracted a great deal of modeling; one of the earliest notable efforts was the elastic propulsion model  suggesting that the curved rear surface of the cell is not merely pushed by growing filaments, but rather ‘squeezed’ forward by an elastic stress. This model treats the actin network as an isotropic elastic continuum not considering explicitly the microscopic mechanism of the force generation at the surface. The squeezing stress develops when the growth of actin at the surface pushes the actin gel outward stretching it and generating tangential tension balanced by radial compression at the surface (Fig. 10). Existence of this squeezing stress was experimentally observed (see  for review).
One of the latest ‘state-of-art’ modeling studies is the ‘nano-propulsion’ model —the first in silico reconstruction of Listeria’s movement. In this model, the filaments propel the virtual bacterium by the tethered ratchet mechanism, but also the realistic geometry and actin network architecture are simulated stochastically. The model also takes into account the reaction–diffusion process of actin recycling and hydrolysis. The simulations result in a vivid and realistic mimicking of the Listeria’s propulsion. A great open problem is combining macroscopic elastic propulsion model and microscopic force generation model at the surface into a mesoscopic theory (?) for Listeria.
Other phenomena attracting modeling efforts are: (i) Listeria appears to rotate about its long axis as it moves, indicating that, in addition to developing an axial thrust, the actin tail generates a net torque about the bacterial long axis . (ii) Symmetry breaking : a ‘cloud’ of actin growing around beads or flexible vesicles loses its symmetry, the actin comet develops, and the bead’s motility ensues. The ratchet models could explain this phenomenon as cooperative stochastic acts of filament growth at one side and actin/crosslinks disassembly at the other side . Elastic models are more successful in describing the sequence of events for large beads: growth of actin at the bead’s surface pushes the outer actin layer outward stretching it and generating growing tangential stress. When critical tangential stress is reached, a crack at the gel outer surface develops and propagates to the bead’s surface . Most of relevant modeling is based on non-rigorous estimates; consistent mathematical theories (?) are pending.
Another example is crawling of the nematode Ascaris suum, the locomotion machinery of which is dramatically simplified . Instead of actin, movement is powered by a cytoskeleton built from filaments of major sperm protein (MSP). MSP monomers associate into symmetrical dimers that polymerize into helical filaments. Hydrophobic and electrostatic interfaces allow these filaments to bundle together, forming thick rope-like fibers that extend from the front to the rear of the lamellipod. The sperm maintain a pH gradient that decreases from the front to the rear of the cell and probably has an important regulatory function. In the more alkaline environment at the leading edge, MSP assembles into filaments. A few models [12,48,70] suggest that these filaments then bundle together into fibers, which forces flexible filaments to assume an end-to-end distance that is longer than it would be in solution, and drive finger-like protrusions from the cell surface (Fig. 10). These bundles are stiffer than, and contain the stored elastic energy of, their constituent filaments. In the acidic cytoplasm at the rear, the fibers lose their adhesion to the substratum, unbundle, and disassemble releasing stored elastic energy and contracting.
Serious mathematical work on this system is just starting, pioneered by a rigorous treatment . One important aspect of the Ascaris motility, and probably of actin-based locomotion, almost untouched theoretically with the notable exception of , stems from the fact that both actin and MSP are highly charged. Because of the counterions to the fixed charges, the filaments of a cross-linked polyelectrolyte cytoskeletal networks are always in a state of elastic tension. Transient changes in mechanical properties of the network are accompanied by changes in local osmotic pressure that could be very significant. Relevant physical models are very difficult, yet mathematically beautiful, and beg for accurate mathematical work.
Three corner stones of cell biology—biochemistry, microscopy, and genetics—have been recently complemented by experimental biophysical and theoretical modeling methods. There are two important roles of mathematical modeling in biology. The first one is a part of the reductionist approach : after the molecular inventory is complete, and concentrations and rates are measured, both biochemical in vitro and mathematical in silico reconstitutions serve as quantitative test for hypotheses about working of biological modules. The second role is to provide a framework for the next step—synthetic, systems approach—understanding how these modules work together.
Significant success of reductionist modeling of motility modules, such as protrusion, adhesion, contraction, and polarity does not mean that the effort in this direction can relax. Quite the contrary—with rare exceptions, the existent models are conceptual, rather than detailed and predictive. Four things are to be done: (1) New conceptual models have to be suggested, but only if they are essentially different from the existent ones, and if there is compelling qualitative data supporting such models. (2) Much more importantly, conceptual models have to be made detailed and predictive by quantifying actual biochemical pathways and physical mechanisms suggested by data and inferring values of model parameters from these data. (3) There are usually a number of alternative explanations for the same phenomenon. When enough data is available, it becomes possible to limit the number of possible explanations using modeling, and ultimately end up with the only true one. The caveat here is that very often in biology redundant mechanisms work in the same cells to make the cell behavior more robust and resilient. (4) Most of the models have ‘theoretical physics’ level of rigorousness; neither modelers deeply involved with experiment, nor, needless to say, experimentalists, are interested in actually proving mathematical results of the model. Meanwhile, rigorous mathematical proofs, fallen out of favor in the past decade, sometime bear unexpected fruit overlooked by hasty ‘front-line’ modelers, besides being intellectually satisfying and beautiful. Therefore, traditional mathematical biologists by no means should shy away from cell motility, but rather be on the lookout for interesting mathematical structures.
As far as systems-level modeling is concerned, in cell motility it is in its infancy. However, valuable lessons can already be learned from other areas where systems approaches are used. One of these lessons is that it will be important to bridge micro- and macro-scales, so the problems will become multi-scale . Two related methodological problems that will arise are: some software for numerical analysis of multi-scale models already exists, but normally it requires formidable computer power and consumes a lot of time. However, models in cell biology are short-lived and malleable, and for them to be useful, one computer simulation cannot be longer than minutes, in order to scan parameter space, compare results to data, change the model many times and screen again. Another problem is that intuition fails to deal with complex multi-scale models, so errors cannot be easily captured, and results interpreted. The only suggestion we have is always to combine small-scale ‘caricature’ models with detailed large-scale ones; then, going back-and-forth between different levels both build intuition and ensures correctness.
We hope that both mathematically oriented and biologically inspired modelers will be inspired by this article to try their hand in cell motility, because there is a number of both mathematical problems and biological questions in the area that require urgent attention. Many interesting mathematical projects do not require following hundreds of papers and memorizing tens of pathways. These projects, to name but a few, are: developing easy to use numerical methods for (i) solving PDE systems (both reaction-drift-diffusion, and continuum mechanics equations) on free boundary domains, (ii) simulating stochastic models of multiple molecules engaged in both chemical and mechanical interactions; (iii) developing adequate theory of the actin-myosin active gel.
Biologically motivated modeling is very rewarding: a natural high of recognizing that your model helped to understand reality is incomparable. There are a number of practical matters that a beginner has to keep in mind when starting on this path. Such work is a long haul–it takes years of reading biological papers, hanging around labs and talking to experimentalists to know a biological subfield intimately enough to be able to formulate right questions and answer them in a way interesting for biologists. One has to be prepared for short life of models and for not being afraid to make specific testable predictions, some of which will be proven wrong. The most crucial is the right choice of the biological questions. So, what are some of these questions in cell motility field that require modeling, in addition to the ones already discussed?
One of the most poignant and important questions: why is the lamellipod flat (?) Why do actin filaments not grow upward from the surface ‘puffing the lamellipod up’? Some models argued that myosin based contraction, or, in general, mechanical effects are responsible. Another possibility is that for some reason actin filaments can grow only in a close contact with adhesion complexes at the surface. Yet another option is some peculiar property of the dorsal membrane inhibiting actin polymerization. Similar more general question is: what regulates the height of more three-dimensional protrusions of cells like leukocytes or neutrophils?
Dendritic actin array in the lamellipod has a very peculiar geometry revealed by electron microscopy studies: there is relatively irregular ‘brushwork’ of short actin filaments near the leading edge, and more smooth and regular, nearly rectangular lattice of longer filaments farther behind the front. What are the dynamics underlining this structure (?)
There are many experimental indications that the leading edge of the motile cells is structurally and biochemically very special: WASP/Scar proteins responsible for activating Arp2/3 complexes are there, as well as adhesion assembly, specific membrane structure, etc. What makes and maintains this distinction (?) One attractive possibility is a positive feedback loop: if there is an elevated concentration of adhesion/polymerization-enhancing complexes at the boundary, then filaments’ growing barbed ends would follow this boundary’s extension; barbed-end-directed molecular motors of myosin family could transport more protein complexes to the barbed ends, in turn maintaining high growing end concentration, etc. Another possible loop is: polymerizing barbed ends push and curve the leading edge membrane; hypothetic membrane proteins have a curved shape making it preferential for them to concentrate in the curved part of the membrane; these proteins may be polymerization-enhancing. No doubt, both of these explanations are too simplistic.
Each of the protrusion models described above is supported by a number of observations, and contradicts some other observations. Multiple explanations of this state of affairs are possible, including especially redundancy of simultaneously acting protrusion mechanisms, and interactions of the protrusive system with membrane mechanics [6,41] and adhesion dynamics . Comprehensive model of redundant protrusion mechanisms (?) interacting with other subcellular systems is one of the most urgent challenges in cell motility. There are other physical mechanisms, such as gel swelling  and osmotic pressure , that also have to be considered.
Recently, the dendritic-nucleation hypothesis of the lamellipodial motility met a very serious challenge, when quantitative experiments  hinted at existence of two lamellipodial networks—frontal lamellipodial actin array and ‘lamellar’ meshwork at the rear, which is thicker, has more chaotic actin geometry, and turns over differently. Theoretical explanation of this duality, which seems to be general (?), is pending. Similarly, protrusion–adhesion–contraction cycles are more complex than what we described for the simple shaped steadily crawling cells . Detailed models of such realistic cycles (?), going much farther than first schematic attempts , are needed.
The role of the cell membrane, other than few phenomena discussed earlier, remains largely unexplored. Pathways of its recycling from the rear to the front during locomotion (endocytosis at the rear and exocytosis at the front, or simple flow of lipids in the plane of the membrane are possible, as well as more complex processes) are not elucidated. How rapidly the plasma membrane area can change is unclear; if this change is slow, whether or not total membrane area is a limiting factor (?) for the cell size is not clear either. Modeling could be a great help in addressing these questions.
Last, but not least, is the problem of modeling three-dimensional cell motility. Essentially 2D lamellipodial motility attracted much attention in the past decades because it was convenient for using experimental methods. However, in most biological and medical situations, the cells move not on the flat surfaces, but rather through 3D extra-cellular matrices. One of the most important recently discussed phenomena is cancer cells’ ability to exhibit both 2D and 3D motility depending on the environment, as well as a switch between them . Scarce data are available about 3D cell crawling through the extra-cellular matrices, but it is already clear that character of the 3D protrusions is very different from that in 3D. Also, adhesion and active remodeling of the environment play much greater roles in 3D. Modeling of the 3D cell motility barely started : this is perhaps the most promising avenue (?) for future modelers. Needless to say, mathematical and numerical challenges in this direction will be formidable.
Answering these questions will require suggesting physical-chemical mechanisms for modules of motility, and adding layers of complexity of spatial–temporal regulation, as well as analyzing systems-level features of sensitivity, robustness and redundancy of the motile networks. Ultimately, the modeling will be complete when not only relatively simple qualitative diagrams, like dendritic-nucleation model, but also comprehensive mechanochemical ‘interactomes’ involving tens of essential molecular players (see the mind-boggling figures in  and ), are quantified and simulated. So, to go back to the title of this article, do we have cell motility’s number? Not quite yet, but we are getting there …
This work was supported by NIH GLUE Grant ‘Cell Migration Consortium’ (NIGMS U54 GM64346) and by NSF Grant DMS-0315782. We apologize for not being able to cite many excellent modeling papers.