Mechanics of neutrophil phagocytosis: experiments and quantitative models.

To quantitatively characterize the mechanical processes that drive phagocytosis, we observed the FcγR-driven engulfment of antibody-coated beads of diameters 3 μm to 11 μm by initially spherical neutrophils. In particular, the time course of cell morphology, of bead motion and of cortical tension were determined. Here, we introduce a number of mechanistic models for phagocytosis and test their validity by comparing the experimental data with finite element computations for multiple bead sizes. We find that the optimal models involve two key mechanical interactions: a repulsion or pressure between cytoskeleton and free membrane that drives protrusion, and an attraction between cytoskeleton and membrane newly adherent to the bead that flattens the cell into a thin lamella. Other models such as cytoskeletal expansion or swelling appear to be ruled out as main drivers of phagocytosis because of the characteristics of bead motion during engulfment. We finally show that the protrusive force necessary for the engulfment of large beads points towards storage of strain energy in the cytoskeleton over a large distance from the leading edge (∼0.5 μm), and that the flattening force can plausibly be generated by the known concentrations of unconventional myosins at the leading edge.


Introduction
Although many components of the molecular machinery involved in phagocytosis have been identified (for reviews, see May and Machesky, 2001;Swanson and Hoppe, 2004), mechanistic explanations of cup formation and subsequent engulfment of a solid particle have lagged behind. Amidst the transformations witnessed by cell biology, it is remarkable that decade-old interpretations of the mechanics of phagocytosis such as the zipper model of Griffin et al. (Griffin et al., 1975), the contractile shear-rigidity gradient model of Hartwig et al. (Hartwig et al., 1980), or the cortical cytoskeleton expansion model of Southwick and Stossel (Southwick and Stossel, 1983) have barely been readdressed since their inception. Textbooks tend to present the process of phagocytosis in a cartoon manner with rods, wires and anchor points forming a phagosomal pocket, which finally closes around the target. Whereas such cartoons are useful implements to convey complicated chemical mechanisms in an intelligible way, the fact that they intuitively seem to imply the right kind of behavior does not guarantee that it is so. For that, a rigorous quantitative analysis of deterministic models is needed along with a comparison with quantitative experimental data.
Regardless of the cell type and phagocytic target, there are four main ingredients that need to be specified to construct reasonably meaningful theoretical models of phagocytic mechanics. (1) The rapid change of cell morphology that accompanies the formation of a phagosome necessarily deforms the cytoplasm, and the forces driving phagocytosis must therefore overcome the viscoelastic resistance of the cytoskeleton. Specification of cytoplasmic rheological properties is thus essential for quantitative models of phagocytosis.
(2) Similarly, the formation of a phagosome is accompanied by a significant increase in macroscopic cell surface area, for which plasma-membrane-folds and -villi are unfolded, and internal vesicular membrane stores are recruited (Herant et al., 2005). The mechanical cost of this process is embodied by the cortical tension against which the forces driving phagocytosis must work to overcome the natural tendency of the cell to minimize surface area by rounding up. Specification of the cortical tension is thus also essential for a quantitative model of phagocytosis. (3) Specific membrane receptor interactions with ligands on the phagocytic target are of course central to the initiation of the signaling cascade that triggers the phagocytic program of a cell. But beyond this, adhesion to the phagocytic target can provide mechanical anchorage for the engulfment process (in contrast to macropinocytosis, where free pseudopods surge around a spatial region without anchorage). Moreover, the cell interface with the phagocytic target can be the locus of special forces that differ from the free boundary of the phagocyte. A quantitative model of phagocytosis thus requires a detailed specification of the mechanical boundary conditions in the contact region. (4) Finally, the driving force of phagocytosis probably lies in protrusive force developed at the advancing To quantitatively characterize the mechanical processes that drive phagocytosis, we observed the Fc␥ ␥R-driven engulfment of antibody-coated beads of diameters 3 m to 11 m by initially spherical neutrophils. In particular, the time course of cell morphology, of bead motion and of cortical tension were determined. Here, we introduce a number of mechanistic models for phagocytosis and test their validity by comparing the experimental data with finite element computations for multiple bead sizes. We find that the optimal models involve two key mechanical interactions: a repulsion or pressure between cytoskeleton and free membrane that drives protrusion, and an attraction between cytoskeleton and membrane newly adherent to the bead that flattens the cell into a thin lamella. Other models such as cytoskeletal expansion or swelling appear to be ruled out as main drivers of phagocytosis because of the characteristics of bead motion during engulfment. We finally show that the protrusive force necessary for the engulfment of large beads points towards storage of strain energy in the cytoskeleton over a large distance from the leading edge (~0.5 m), and that the flattening force can plausibly be generated by the known concentrations of unconventional myosins at the leading edge.
rim (the leading edge) of the phagocytic cup that hugs the phagocytic target. The so-called 'polymerization force' has been invoked for this process, but although single-filament analyses have been performed, calculations of the active behavior of whole cells have been few. It is clear, though, that such driving forces need to be specified for any quantitative model of phagocytic mechanics.
Although the four ingredients listed above make abstraction of the prodigiously complex biomolecular machinery that operates in phagocytosis, their specification already represents a tall order, even for a primitive study of the mechanics of phagocytosis. Among available experimental systems, the phagocytosis of spherical beads by initially round neutrophils is closest to minimizing these mechanical unknowns by providing empirical constraints that are as tractable as possible. First and foremost, the geometry is simple in that it consists of 'a sphere that eats a sphere', with concomitant formal twodimensional cylindrical symmetry along the axis joining the center of the two spheres. Next (ingredient 1), the rheological properties of the neutrophil have been studied for decades and are fairly well defined (e.g. Evans and Kukan, 1984;Herant et al., 2003). The behavior of the cortical tension with neutrophil deformation during phagocytosis has recently also been characterized (ingredient 2) (Herant et al., 2005). Although they are at the heart of the phagocytic process, the last two ingredients (3 and 4) remain, however, poorly defined. It is the principal aim of this paper to constrain them better by comparing simulations of a variety of mechanistic models of phagocytosis with experimental data.
In the present work, we study the mechanics of Fc␥R-driven phagocytosis of antibody-coated beads (diameters ranging from 3 m to 11 m) by initially spherical neutrophils. We present side by side the results of microscopic observation and numerical modeling. We comment on the compatibility of various theoretical models with empirical data that encompass the time-dependent shape of the neutrophil (especially the duration of the phagocytosis), the time course of bead position and the time course of the cortical tension during phagocytosis. We show that repulsive interactions between cytoskeleton and plasma membrane are most likely to be the main mechanical drivers for phagocytosis. We further show that successful theoretical models require transverse contractile forces at the neutrophil-bead interface to yield the appropriate lamellar geometry observed during phagocytic engulfment. Finally, we show that the data appear to rule out phagocytosis models in which protrusion is caused by cytoskeletal expansion (gel swelling).

Experimental findings
Some of the experimental observations used in this study have already been presented in Herant et al. (Herant et al., 2005), especially with respect to the behavior of the cortical tension. Briefly, when an antibody-coated bead is presented to a neutrophil, near instant adhesion occurs and is followed by a latency period during which no cellular activity is apparent under light microscopy. The duration of the latency period is variable (seconds, minutes) but is typically shorter for larger beads. The transition to active phagocytosis is sudden with extension of the cell-bead contact area followed by the emergence of pseudopods surging around the bead to form a Journal of Cell Science 119 (9) phagocytic cup that progresses without interruption to complete engulfment of the bead (unless the bead is too big). It is this phase of active engulfment that is examined here.
At first inspection of the data, a pair of aggregate proxy parameters emerge as those easiest to determine. The first is the duration of phagocytosis from beginning to end of engulfment. The second is the maximum cortical tension observed during phagocytosis. Those parameters are listed in Table 1 and, unsurprisingly, they both monotonically increase with bead diameter. Of particular interest is the maximum tension of ~1 mN m -1 , which is observed during frustrated phagocytosis of 11 m beads and whereby neutrophils manage to engulf only about 50% of the circumference. This tension is essentially the stall force (per unit length of leading edge), which halts further protrusion of the leading edge of the phagocytic cup. Thus, for the neutrophil, the stall force is 1000 pN m -1 for a total force ~30,000 pN. Interestingly, this is consistent with data on traction forces exerted by spreadmotile neutrophils that yield an average stress of 140 pN m -2 (Daniel Hammer, University of Pennsylvania, PA, personal communication). For the area of the half-sphere (radius 5.6 m, surface 200 m 2 ) this gives a total force of 28,000 pN, remarkably close to our estimate derived from the tension. This is also consistent with locomotive forces measured in single neutrophils with magnetic tweezers (10,000 pN) (Guilford et al., 1995) or with micropipette counter-pressure (30,000 pN) (Usami et al., 1992).
Despite the difficulty in quantifying and rigorously delineating a sequence of shapes imperfectly observed by microscopy, some basic qualitative features of the morphologic changes that take place from the beginning to the end of phagocytic engulfment can be described (see Fig. 1). After attachment, the contact area between neutrophil and bead grows, but the cell retains a blunt, nearly spherical shape. We call this the 'pedestal-phase'. Eventually, lamellar pseudopods emerge from the cell body on either side of the bead and protrude forward around the bead. We call this the 'cup-phase'. During the cup-phase, there is a distinctive joining angle between the cell surfaces belonging to the cell body and to the lamellae. We refer to this angle as the 'body-lamella contact angle'. Finally, the bead is pulled into the cell while the cell rounds up and the body-lamella contact angle vanishes. We call this the 'rounding phase'.
The transition from pedestal-to cup-phase seems to occur fairly consistently at the moment the diameter of the contact area is approximately equal to the cell radius, regardless of the size of the bead. However, the transition between the cup-phase and the rounding-phase (relative to beginning and end of engulfment) depends strongly on the size of the bead. For large beads it takes place early, whereas for smaller beads it takes place late, and for the smallest beads it may not occur at all before the end of engulfment. The motion of the bead with respect to the cell is affected by the three phases. During the pedestal-and cup-phases the bead moves inward slightly, but during the initial part of the rounding-phase it is drawn in rapidly. Notice that, we never observe the bead moving away from the cell to any significant degree, as might be expected under the outward push of growing pseudopods. As discussed below, this has important consequences for theoretical mechanical models.
The surface tension remains low during the pedestal-and cup-phase but increases strongly during the rounding-phase to eventually decay back down after engulfment is completed. This pattern has been interpreted (see Herant et al., 2005 for details) to be the result of a first stage, during which slack in the membrane is sufficient to accommodate the initial increase in macroscopic cell area without an increase in surface tension. If the bead is large enough, the slack runs out and further surface area has to be recruited by breaking bonds that stabilize surface-folds and -villi. This leads to a rise in surface tension that partially subsides when phagocytosis is complete.

Theoretical models
Much of our study centers on the successes and failures of various theoretical models of neutrophil phagocytosis. We deem the minimum requirements for a successful model to be the following.
(1) The model should incorporate self-consistent mass and momentum transport. This is, de facto, the case in this study because our modeling is based on mass and momentum conservation laws.
(2) The model should yield the correct time development of neutrophil morphology during phagocytosis. This includes the succession of the pedestal, cup and rounding-phases as well as the proper duration of engulfment.
(3) The model should be compatible with the observed time-dependent motion of the bead with respect to the neutrophil. (4) The model should yield the correct time course for the cortical tension during phagocytosis. We should remark here that, as discussed in the Materials and Methods, because the theoretical behavior of the surface tension with surface area was derived from observations of phagocytosis, a kinematically accurate model is expected to automatically provide an approximately correct tension. (5) The constitutive relations (e.g. magnitude of protrusive force along the rim of the phagocytic cup) should be biologically plausible. This, of course, is the crux of model believability but in many cases cannot be fully ascertained because the data on the biophysical details are simply not available. (6) Since the cell does not have any 'external' means of assessing the bead size, there cannot be explicit dependence of the model on bead dimensions. In other words, once the set of free parameters of a model are fixed, they should remain valid for different bead diameters. We point out, however, that there exists implicit dependence of the model on bead size as, for instance, the total area of contact and its associated total 'flattening force' increase with bead diameter.
In the process of testing mechanical models of phagocytosis, we have tried many different ideas, most of which proved wrong in some place or another, by yielding incorrect shapes, incorrect time courses, etc. Altough it is not possible to recapitulate them here, we would like to make the point that what does not work is in many ways just as interesting as what does or, to follow Popper, much of the worth of what we do is about invalidating models that might have at first blush seemed reasonable. In what follows, we first outline the optimal model that worked best. We then move on to show that a flattening contractile force is an important component of that model (and probably any potential model), by demonstrating what happens when it is removed. We later present a gel-swelling model that fails. We end with a more detailed view of the importance of cortical tension in the time-dependent development of cell morphology.
A model that succeeds: membrane-cytoskeleton repulsion with flattening force We begin with a cartoon that summarizes the key features of the model (Fig. 2). As contact occurs between bead and cell, Fc␥R receptors are sequentially engaged like in the zipper model of phagocytosis (Griffin et al., 1975). Newly formed adhesion sites anchor the internal cytoskeleton relative to the bead. We then postulate two active mechanisms leading to engulfment. The first is that polymerization of actin at the leading edge (the cell-bead contact line) drives protrusion through repulsive interactions between the cytoskeleton and the free membrane. The second is that contractile activity at the cell-bead interface drives flattening through attractive interactions between the cytoskeleton and the bound membrane.
The implementation of this model in the reactive interpenetrating flow formalism is done as follows. A small region of membrane (~1 m) is taken to emit polymerization messenger at high level, thus driving up the network concentration in that area by a factor of about 20. Simultaneously, during engulfment, every patch of cell surface newly adherent to the bead is taken to exert an attractive force on the internal network that is initially in the order of 1000 pN m -2 and that decays over a time scale of 1 minute. Since there is new surface adhering as engulfment progresses, the attractive force remains principally localized near the leading edge.
At the beginning of a simulation, the virtual spherical cell, which is filled with a low-baseline network concentration, is placed in tangential contact with the bead. Emission of polymerization messenger is turned on and, after a timescale of ~ n =20 seconds, the network reaches an equilibrium value that drives protrusion around the bead. Simultaneously, the attractive force of newly formed cell-bead contact areas is ramped-up linearly to its regular value (~1000 pN m -2 ) over Journal of Cell Science 119 (9) 20 seconds. This tends to pull the bead into the cell (or the cell onto the bead).
The main physical parameters that determine the dynamical behavior of the model are the cytoplasmic viscosity, the cortical tension, the membrane-cytoskeleton interaction strength and the flattening force. These four factors are obviously coupled with the cytoskeletal density that is strongly influenced by chemistry. In neutrophils, viscosity and cortical tension are well constrained by prior experimental and numerical results (Herant et al., 2003;Herant et al., 2005), so that they are not amenable to further tuning. Moreover, the maximum cytoskeletal (network) density can at the most reach a few percent of the volume fraction. This leaves only the membrane-cytoskeleton interaction strength as a significant tunable parameter with which to obtain the right dynamical timescales (the flattening force mostly impacts morphology). This parameter was therefore varied until agreement with the aggregate parameters (phagocytosis duration and maximum surface tension) of Table 1 was good.
Figs 3 and 4 show a direct comparison between simulations and experimental data (see also supplementary material Movies 1, 2). Notice that, for each bead radius, we used images from our best dataset -namely the one in which the cell and the bead were in the clearest focus, not the one that agreed best with our simulations. Three images were picked for each bead size, depicting the three phases of engulfment (Fig. 1). These images were then matched with simulation data that corresponded to the same time intervals, e.g. the time span between first and second microscopic image is the same as the time span between first and second simulation image. We also provide plots of the cortical tension and bead position for the displayed dataset versus that calculated by simulation. Finally, to illustrate the natural variability of neutrophil behavior, we also include a case where measured surface tension differs significantly from the one originally selected. Inspection of the data presented in Figs 3 and 4 shows that the overall match between experiment and simulation is fair -about the best one could expect with such a complex . Tension plots show experimental data in grey (filled squares correspond to shown micrographs, empty squares to other phagocytoses of same diameter beads purposefully selected for significant difference). For graphs on the left, the heavy black line shows total surface tension from simulation (red is the elastic contribution ␥ el and blue is the viscous contribution ␥ vis ). On the right, bead center position is shown with respect to the back of the cell. A magnification of the boxed area shown in the middle image of the second half of this figure is shown in Fig. 5. biological system. We comment on a few selected points that we consider of special interest. (1) Generally consistent with experimental data, the simulation of the 3-m-bead phagocytosis does not show a rounding-phase (at least before the end of engulfment), as witnessed by the fact that the cellbody to lamella contact-angle remains sharply defined throughout. This is because the increase of surface area of the cell is less than the threshold that drives surface tension up significantly (see section on surface tension and cell morphology).
(2) The simulation of the 4.6 m bead phagocytosis yields a tension curve in poor agreement with that observed experimentally, especially in that the latter stays elevated after engulfment is complete. This is evidence that in contrast to our essentially passive model of elastic and viscous cortical tension, there exists a contribution from active processes operating in the absence of surface expansion. (3) In the phagocytosis of 4.6 m, 6.2 m and 11.2 m beads, rapid inward motion of the bead is observed as soon as the tension rises significantly (when the change of cell surface area goes above a threshold). As the energetic cost of recruitment of new surface area increases, pseudopod length is minimized by the cell drawing the bead into itself.
The focus of this work is on whole-cell mechanics rather than lamellar dynamics. Nevertheless, we should like to comment briefly on the findings of our model regarding the lamella that spontaneously develops at the leading edge of our virtual neutrophils. Fig. 5 shows an expanded view of the lamella for the 4.6 m phagocytosis displayed in Fig. 3. Whereas it is apparent that the polymerization messenger is horizontally rather than vertically stratified, the cytoskeleton is mostly vertically stratified. This is owing to two reasons. The first is that the contractile force exerted by the area attached to the bead pulls the cytoskeleton down. The second is that, since the boundary condition at the bead interface is 'no-slip' in addition to 'stick' (zero velocity), the cytoskeleton there cannot easily be dissipated by advection flows. The resulting highdensity layer of network adjacent to the bead provides excellent viscous bracing against which the leading edge can thrust out.
One may also notice the region of somewhat lower network concentration next to the free membrane at the leading edge, which is simply the result of network-membrane repulsion.
In the frame of reference of the substratum (i.e. the bead), forward progression of the lamella is evident and is accompanied by a slow rearward and downward motion of the interior cytoskeleton. In the frame of reference of the cell, centripetal flow of cytoskeleton is obviously more noticeable. This is in line with observations of treadmilling cytoskeletal flow in the lamellae of migrating cells (e.g. Wang, 1985;Vallotton et al., 2004). Finally, we point out that lamellar geometry and flow fields were obtained without any assumptions regarding cytoskeletal anisotropy.
A contractile flattening force is necessary In our exploration of mechanical models of phagocytosis, we found it necessary to postulate that molecular motors (presumably, but not necessarily myosins) anchored to the substratum (the bead) pull on the interior cytoskeleton (Fig. 2). Although prior studies have hinted at contractile phenomena during phagocytosis (Evans et al., 1993;Vonna et al., 2003), this type of contractile activity has not been proposed before. We postpone a review of the experimental evidence on the role of myosins in phagocytosis until the discussion section. However, it is worthwhile to examine what happens when this component of the model is removed. Fig. 6 shows simulations of 3-m-bead and 6-m-bead phagocytoses performed with and without contractile flattening force. The following differences are observed in the absence of force. Not unexpectedly, the lamellae are significantly thicker (by >50%). The 3 m bead is pulled in less (apparent length pulled in l=0.6 m without flattening vs l=1.4 m with flattening compared to observed l~1.3 m). At the close of engulfment of a 6 m bead, the neutrophil is much rounder (apparent eccentricity of cell outline at the end of engulfment e=0.17 without flattening vs e=0.45 with flattening compared to observed e~0.5). All these findings go against the experimental evidence and, in our view, make the membrane-cytoskeleton-repulsion model inadequate in the absence of a flattening force. Whereas it is certainly possible that alternative mechanisms might be invoked to its rescue, this approach seemed the simplest and possibly the most plausible.
Journal of Cell Science 119 (9) A model that fails: gel swelling An alternative to the polymerization force as a driver of cellular protrusion is gel swelling -or, in other words, an alternative to membrane-network repulsion is network-network repulsion (Herant et al., 2003). The idea is that repulsive forces (entropic, electrostatic, etc.) between filaments lead to the expansion of overdense regions created by locally enhanced polymerization. This mechanism has previously been shown to be able to reproduce observations of multiple neutrophil phenomena such as pseudopod extension (Herant et al., 2003). However, we will argue here that it is extremely difficult for this type of model to account for what is observed in phagocytosis.
There are two main problems inherent to a gel-swelling model of phagocytosis. The first lies with the key experimental observation that the bead never moves away from the cell. This immediately places strong constraints on swelling models since a build-up of network just interior to the bead is expected to expand and at least initially, push the bead out. The second problem is that pseudopods have a tendency to swell and become thick. Both these problems can be partially remedied by the introduction of a contractile flattening force of the type already described above. Such a force helps to prevent both the outward motion of the bead and the thickening of the pseudopods, but even then, a number of difficulties remain.
For a swelling stress capable of generating protruding force corresponding to the stall surface tension (~1 mN m -1 ) observed in the incomplete phagocytosis of an 11 m bead, the flattening force per surface area has to be large (~6000 pN m -2 ) at the leading edge: this is improbable (see Discussion). Worse, the delicate balance between large swelling and flattening forces leads to a strong sensitivity of the simulations on the precise choice of force parameters. If the flattening force is slightly too low, the bead is pushed out and the phagocytic cup is much too thick. If the flattening force is slightly too high, the bead moves inward faster than is observed and pseudopods may pinch off. Finally, even with optimally chosen parameters, transient shapes are observed that do not match observations (Fig. 7).
It would be possible to damp-out this oversensitivity to parameter choice by increasing the cytoplasmic viscosity but the latter is already well constrained in the case of neutrophils (Herant et al., 2003) and cannot be adjusted. We believe that the requirement for fine-tuning of parameters to obtain a reasonable match to the observed dynamics of phagocytosis invalidates cytoskeletal expansion as the principal mechanical driver. Biological behaviors should be robust to the expected high level of perturbations and variations that occur in living processes, and a putative mechanism that displays an exquisite parameter dependence has to be rejected. This does not mean, however, that some degree of cytoskeletal swelling does not occur -there are reasons to believe that it is necessary to prevent separation of cytosol and cytoskeleton in externally imposed deformations such as pipette aspiration, although such a role can also be carried out by cytoskeletal elastic stiffness (Herant et al., 2003).

Surface tension determines transient cellular morphology
The synchronization of the transition from the cup-phase to the rounding-phase (see Fig. 1) with the increase in cortical tension is not a coincidence; it is explainable in simple The white patch of membrane corresponds to the region of enhanced messenger emissivity. Notice that, the velocity of the cytoskeleton vanishes in the substratum frame at the cell-bead interface due to the no-slip boundary condition. Notice also the high density of cytoskeleton at that same interface due to the contractile force directed down toward the substratum. physical terms, and such an understanding can only buttress confidence in the validity of the models developed in this study. Let us start by considering a viscous, initially hemispherical, droplet whose base begins to spread for unspecified reasons. Let the outward velocity of the dropletsubstrate contact line be v c . Does the droplet send out a precursor film or does it spread as a whole spherical cap (Fig.  8)? The answer depends on the relative importance of the surface tension ␥ that tends to minimize surface area by maintaining a spherical cap shape, and the interior viscosity that resists deformation of the droplet body. By simple dimensional analysis, the velocity of deformation of the droplet body is b~␥ / (this upper limit is approximately correct as long as the body is not too thin) (see de Gennes, 1985 for details). Thus, if v b ӷv c the droplet spreads as a cap but if v b Ӷv c , a precursor film travels ahead of the body.
Turning to neutrophil phagocytosis, it is found that the typical velocity of the leading edge moving around the circumference of the bead is v c~0 .1 m second -1 . The overall viscosity of the cell body is ~ 0 0 =6ϫ10 2 Pa seconds. The cortical tension on the other hand is highly variable due to a threshold effect in area expansion. For area increases (over spherical area) of less than 26%, the cortex remains fairly slack and ␥ slack~0 .03 mN m -1 . For larger area increases, the cortex becomes taut and ␥ taut~0 .3 mN m -1 . Using those numbers to compute the velocity of the body-lamella contact line (Fig. 1) one gets v slack~0 .05 m second -1 <v c and v taut~0 .5m second -1 >v c . At first the lamella extends ahead of the cell body, but when the area expansion reaches the threshold that causes a sharp increase in tension, the body of the cell subsumes the lamella by rounding up.
To illustrate this, we have simulated the phagocytosis of an 11 m bead by a cell with double than usual membrane slack (52% vs 26%). Fig. 8 shows that, indeed, the model with the lower surface tension produces longer lamellae when measured from the leading edge to the body-lamella contact angle. Finally, we point out that this analysis only pertains to transient shapes observed during engulfment. Once the leading edge stops moving (i.e. after engulfment is complete), the cell always evolves towards a spherical shape.

Discussion
Whereas failure of a mechanistic model to reproduce experimental data is usually a fairly strong argument for rejection, success in reproducing observations is unfortunately a weaker argument for acceptance. To strengthen the case, it is necessary to establish, if not the validity, at least the plausibility of the model's key assumptions. We therefore consider the biological underpinnings of the constitutive laws governing protrusion and flattening as implemented in our optimal model of phagocytosis.

Leading edge protrusion
There is good evidence that phagocytic engulfment and migratory cell spreading involve the same cytoskeletal machinery (May and Machesky, 2001) and it is thus in this generalized setting that we discuss our models for lamellar progression. Our optimal model of phagocytosis (Fig. 2) postulates that the driving force for protrusion is provided by a repulsion between cytoskeleton and plasma membrane. This repulsion is, for now, uncharacterized but the question naturally arises whether its magnitude -as determined by our computational experiments -can plausibly be achieved.
The disjoining force between cytoskeleton and membrane derives from a stored energy per surface area, much in the way that a bunch of parallel coiled springs can store energy to push a plate. The constitutive law that we have used in our calculations implies that the energy available for work is 0 nM d M m n , where d M is the effective range over which the cytoskeleton membrane interaction contributes the energy density ~ 0 nM m n (see Table 2). In our calculation 0 nM d M =1.3ϫ10 -3 pN m -1 , and m~60 so that the stored energy is (in the natural unit of energy k B T~4.3ϫ10 -21 Joules) 2.3ϫ10 7 n k B T m -2 . Since the turnover of the actin network is rapid, it is unlikely that the free energy associated with a monomer addition is more than a few k B T, but let us assume a generous 5 k B T. The energy density associated with this is 5 k B T n /V n where V n is the volume of a single actin monomer. Taking, V n =r -(2.7 nm) 3 , the energy density becomes nM n = 6.1ϫ10 7 n k B T m -3 . This yields d M~0 .4 m as a lower limit of the range of the stress field. This is a large distance compared with the monomeric and other solvated molecular scales, but it is smaller than the F-actin persistence length (e.g. Kovar and Pollard, 2004). This means that the stress field has to be stored in the large-scale strain energy from deformation of the cytoskeleton. For this to be the case, the timescale for filament lengthening over a distance d M must not be longer than the time for filament-strain relaxation. Or, in other words, because strained filaments tend to relax by flowing backward without producing forward motion, this 'loss of strain' must be compensated by the 'production of strain' with the polymerization of new filaments. The relaxation time is given by the ratio of the viscosity to the stress energy at 0 n / nM n~2 seconds. However, the monomeric addition rate to free filaments is of the order of at least a hundred per second (Pollard, 1986) leading to a lengthening that is indeed at least d M over the relaxation time. The magnitude of the protrusion is therefore plausible because it does not lead to a rate of strain relaxation by viscous dissipation faster than the rate of strain creation by polymerisation.

Contractile flattening force and the role of myosins in phagocytosis
Our optimal model of phagocytosis (Fig. 2) postulates the existence of molecular motors that are anchored to the cell-bead interface and pull down on the actin cytoskeleton. Although this mechanism has not been proposed before, a number of observations give it plausibility. One cannot exclude a role for actin-cross-linking proteins or microtubules as agents of the flattening force, but the most likely motor(s) have to be myosin(s). It is well established that myosin is present in the phagocytic cup, especially at the leading edge of pseudopods (Valerius et al., 1981). However, in the case of Fc␥R-mediated phagocytosis by macrophages, it appears that (conventional) myosin II activity is not necessary for engulfment (Olazabal et al., 2002), and we can confirm this for our own experimental conditions because addition of the myosin II inhibitors ML-7 or blebbistatin did not have a major impact on neutrophil phagocytic activity (data not shown). The case of unconventional myosins is more complicated because there are many. Nevertheless, multiple unconventional myosins have been shown to localize at phagocytic cups in various model systems such as Acanthamoeba (Ostap et al., 2003), Entamoeba (Voigt et al., 1999), Dictyostelium (Titus, 2000), and macrophages (Diakonova et al., 2002). Despite these findings, the functional dependence of phagocytosis on unconventional myosins has been difficult to delineate clearly because they might have overlapping functions, though an impact on phagocytosis has been demonstrated by knockout experiments in Dictyostelium (Jung et al., 1996) and macrophages (Cox et al., 2002). In addition, because our simulations without a flattening force (see Fig. 6) show only a moderate change in morphology during phagocytosis rather than outright phagocytic failure, the effects of successfully turning off the 'flattening force' in experiments are probably subtle.
Interpretations of the role of myosins in phagocytosis have Journal of Cell Science 119 (9) varied. For instance, the suggestion that they act to tighten the purse-string at the stoma of the phagosome has been made (Swanson et al., 1999), but this would only be effective past midpoint of engulfment. More often heard is the notion of shearing activity by which a sliding interaction between myosin and actin produces a longitudinal (with respect to the pseudopod) rather than a transverse force -as in our model (e.g. Chavrier, 2002). However, a potential difficulty of models based on shearing motors is that they require a high degree of cytoskeletal organization that clearly distinguishes between front and rear. That this is the case is not clear from the tangles of filaments seen in electron microscopy (Ryder et al., 1984). On the other hand, a flattening force offers a number of advantages. Orientation is natural because the substratum is de facto 'down'. Efficiency of forward thrust tends to be improved because wasteful growth in the 'up' direction is restrained, and also because contraction lays down, adjacent to the nascent phagosome, a thick mat of actin that provides viscous bracing for forward pseudopodial motion against cortical tension. Whereas the force exerted by a single, active myosin molecule is well determined to be around 1 pN (Finer et al., 1994;Molloy et al., 1995) the number of myosins active in a cell remains uncertain. Quantitative estimates of a few hundred to a few thousand molecules per m 2 of cell area have been published for Acanthamoeba (Baines et al., 1992;Baines et al., 1995). The magnitude of the attractive force posited in the optimal model of 1000 pN m -2 in the area of the leading edge is thus reasonable. A force of 6000 pN m -2 as necessitated by a swelling model is, though not impossible, less likely.

Experiments
The experimental procedure has been described in detail in Herant et al. (Herant et al., 2005). Briefly, polystyrene beads (Polysciences, Inc., Warrington, PA) were incubated with bovine serum albumin (BSA) overnight, washed, incubated for 1 hour at room temperature with polyclonal rabbit anti-BSA antibody (ICN Biomedicals, Aurora, OH) and then saved for future use. On days of experiments, a drop of whole blood was obtained by finger prick and microliter quantities of both blood and antibody-coated beads were separately introduced into an experimental chamber filled with heparinized Hanks balanced solution (with Ca 2+ and Mg 2+ ; Sigma, St Louis, MO) with 10% autologous plasma. Individual beads and quiescent spherical neutrophils were then selected visually among the background erythrocytes, picked up with opposing micropipettes and brought into brief contact. Once adhesion occurred (usually immediately), beads were released and the ensuing phagocytosis was observed using a 63ϫ objective under brightfield illumination (Axiovert inverted microscope, Carl Zeiss Microimaging, Inc., Thornwood, NY) and a cooled 12-bit video camera (SensiCam, The Cooke Corp., Auburn Hills, MI). Application of the same experimental protocol in control experiments using albumin-coated beads without antibody did not elicit adhesion or phagocytosis from neutrophils, providing strong evidence that the phagocytosis we observed was indeed antibody-mediated.
The aspiration pressure required to maintain a cytoplasmic projection of stationary length inside the micropipette holding the neutrophil was monitored during the experiments using a setup previously described (Heinrich and Rawicz, 2005). Pressure data were then transformed into cortical tension data by using the law of Laplace relating the tension ␥ to aspiration pressure ⌬P and membrane curvatures, which yields ␥=G⌬P/(1/R p -1/R) where R is the cell radius, and R p is the pipette radius.

Theory
The methodology used here to model phagocytosis has been described and validated in the context of other neutrophil phenomena in Herant et al. (Herant et al., 2003). Briefly, it assumes that at the mesoscopic scale (a scale that is small compared with the whole cell, but large compared to individual molecules), the principal actors of neutrophil mechanics are the membrane, the cytosol and the cytoskeleton, and that their roles can be captured by the methods of continuum mechanics. The membrane defines the boundary of the cell and, to the extent that cytoplasmic volume is conserved (as observed in our experiments), it can be thought to be impermeable to bulk volume fluxes. It contributes to the mechanical behavior of the neutrophil through the cortical tension at free surfaces and boundary conditions at contact surfaces (e.g. at the interface with the phagocytic target). The aqueous cytosol is a passive player but its incompressibility allows it to transmit pressure stresses. It is also a medium through which chemical signals can diffuse. The cytoskeleton is endowed with elastic and/or viscous properties that oppose deformation. It also can be stimulated to produce active forces that result in spontaneous movement and, thus, changes of shape. These active forces can take many forms, among which cytoskeletal-cytoskeletal interactions (e.g. gel swelling), and cytoskeletalmembrane interactions (e.g. cytoskeleton-membrane disjoining force) represent generalized formulations (Herant et al., 2003). Finally, the distinction between cytosol and cytoskeleton, although well-defined, is not permanent. For example, a solvated G-actin monomer belongs to the cytosolic phase but can convert to the cytoskeletal phase when polymerized into F-actin. Notice that, our models do not consider distinct mechanical properties for the nucleus like some others have done (Kan et al., 1999). In the case of neutrophils the nucleus is segmented in smaller lobules that are easily movable. Furthermore, our observations show that even in extreme deformations, such as during the phagocytosis of large beads, there is no evidence of stiffer intracellular organelles of significant size getting in the way.
The reactive interpenetrating flow formalism originally formulated by Dembo and Harlow (Dembo and Harlow, 1986) provides an appropriate framework in which to cast the ideas described above in a quantitative manner. In all that follows, the subscript ' s ' denotes 'solvent' (cytosol) related quantities and ' n ' denotes 'network' (cytoskeletal) quantities. When subscripts are omitted, network is implied. The basic idea is that at every point of the cellular interior, one can define volume fractions of solvent and network s and n , and the velocity fields of each of the two phases v s and v n . Notice that, these two velocities usually differ (as is the case when water flows through a squeezed sponge, for instance). Within this picture, problems of cell mechanics are reduced to determining the change in time of the volume fraction and velocity field for the network and the solvent. To do so, the laws of mass and momentum conservation can be used to define evolution equations.

Network-conservation equation
where J is the net rate of conversion of solvent to network (cytosol to cytoskeleton) by polymerization at a given location.

Solvent-momentum equation
where P is the pressure, and Ᏼ is the solvent-network drag coefficient (this assumes that the cytosol is essentially inviscid except for interactions with the cytoskeleton via Darcy's law).

Network-momentum equation
In this equation, ⌿ nn is the stress (tensor) of the network due to interfilament forces, ⌿ nM is the network-membrane interaction term, is the viscosity, F el is the elastic force due to deformations, and F ext includes applied body forces.

Constitutive relations
Prescriptions are necessary to provide closure to the mass and momentum equations: these are the constitutive equations that link physical laws and biological behaviors. The constitutive parameters that are used for our calculations are listed in Table 2.

Network polymerization
For simplicity, the net rate of polymerization of network is set by first order kinetics: where eq is the equilibrium network concentration, and n =20 seconds is the timescale to reach that equilibrium (e.g. see Cassimeris et al., 1990). In our simulations, eq is set by the (dimensionless) local concentration m of a polymerization messenger: (6) eqn n J = , where 0 is a baseline (quiescent) network concentration, 0 = 0.1% as per estimates of Watts and Howard (Watts and Howard, 1993). The polymerization messenger is taken to be produced by the plasma membrane with emissivity ⑀ m (m second -1 ), lifetime m (1 second), and diffusion coefficient D m (m 2 second -1 ) so that with appropriate von Neumann boundary conditions: where n is the unit normal vector to the membrane. D m and m are chosen such that the penetration depth ͱහ D m හ m is small (0.3 m); as a result n = 0 = 0.1% for most of the interior of the cell. However, when the neutrophil is stimulated, the emissivity rises sharply in certain areas resulting in localized polymerization and increase in equilibrium network concentration up to about 40 times the baseline concentration ( eq max ~40 0 ). Notice that, the 'messenger' is, of course, a catch-all for a complicated collection of biochemical intermediates, and is not intended to represent a real single signaling molecule.

Network-network interaction
In simulations of the cytoskeletal 'swelling model', we have implemented an isotropic stress term proportional to the network concentration, i.e.: where ⌱ ⌱ is the unit tensor. If 0 nn is positive, one has a repulsive pressure-like term that causes overdense regions of network to expand.

Network-membrane interaction
In simulations of the cytoskeletal 'membrane repulsion model', we implement an additional stress term present only in a cortical layer of thicknes d M ; this term has the form: where nn is the dyadic of the unit vector outward normal to the membrane. For nM positive the effect is to cause a normal stress or pressure that acts to push the membrane out. Naturally, the network experiences an equal and opposite reaction so that this stress corresponds to a disjoining force that tends to expand the cortical network layer. In the case of a polymerization force model, nM naturally depends  The dynamically relevant term is the network-membrane potential energy times its range d M . § Specific network swelling pressure was set to 3ϫ10 5 Pa for the gel swelling calculations. ¶ This is an approximate 3D-value derived from the actual 2D cylindrical geometry implementation. on the local polymerization rate (driven by the messenger m) and network density: in the calculations we have taken nM = 0 nM m n .

Network viscoelasticity
The viscosity is taken to be approximately linearly dependent on network concentration: = 0 n exp ( n / g ) .
For strain rates of order 0.1 seconds -1 in activated neutrophils, it was previously found that 0 =6ϫ10 5 Pa second for baseline viscosity 0 0 = 6ϫ10 2 Pa second was appropriate (Herant et al., 2003). This baseline viscosity is higher than that found in aspiration experiments of passive neutrophils (~100 Pa second) (Hochmuth, 2000), because neutrophil activation and lower strain rates both appear to increase the cytoplasmic viscosity, presumably through enhanced cross-linking of the cytoskeleton. A network gelation concentration threshold g =0.1 is also added; it helps stabilize lamellae against collapse and increases protrusive forces in stall conditions by improving bracing characteristics. On physical grounds, positing such a gelation process is reasonable to the extent that high-density networks can become highly crosslinked, but we also note that almost everywhere in our simulated cell n Ӷ g . The elastic stress term F el (used in Herant et al., 2003) was found to have little impact on the phagocytic models (results not shown) and was therefore set to zero for simplicity.

Body forces
As explained in Results under 'A contractile flattening force is necessary', an attractive force between the cell-bead interface and the cytoskeleton was found to be necessary. We set the force between a unit volume element of cytoskeleton and a unit area of bead neutrophil interface to be where n is the unit normal to the interface, d is the vector joining the surface to the volume element, d 0 =0.3 m is a force decay length scale, and f 0 b is the 'flattening force' coefficient. This functional form was chosen as the simplest that gave a smooth force with a well-defined range. Other, more complicated forms were tested as well; they tended to give similar results as long as the range and magnitude were comparable. Integrating the force exerted by an elementary surface patch over the covered half-space volume beneath it, one can show that the total force exerted by a unit surface area on the interior cytoskeleton is of order f 0 b d 0 3 .

Boundary conditions, surface tension
Boundary conditions have to be specified for both solvent-and network-phases. If we assume impermeability of the plasma membrane to both water and cytoskeleton, we have: where n is the unit normal vector to the membrane and v M is the velocity of the membrane. For convenience, the frame of reference of the simulations is such that the bead is stationary while the neutrophil moves around it. At the solid neutrophil-bead interface normal velocities vanish, i.e. v s ·n=v n ·n=0. We further assume that due to anchoring of the cytoskeleton through transmembrane adhesion sites, the tangential network velocity must vanish (no-slip) at the neutrophil-bead interface, i.e. v n ·t=0 (where t is the membrane tangent vector). Finally, we assume strong adhesion: whenever a patch of free boundary comes into contact with the bead surface it is irreversibly immobilized and we do not allow detachment, regardless of boundary stresses.
For free boundaries, stress balance across the membrane requires equating the internal stress tensor with external stresses while taking into account the contribution of surface tension: [ٌv n + (ٌv n ) T ] · n -⌿ · n -Pn = -2␥n , where ⌿ is the full interior network stress tensor, ␥ is the cortical tension, and is the mean curvature of the surface. Specification of the membrane tension demands care. Let A 0 be the macroscopic area of the cell in its spherical state, and A the area of the cell in its deformed state. During phagocytosis, the area A increases, and an area ⌬A is recruited by unfurling folds and villi: ⌬A = (A+0.05A 0 ) -A 0 -A exo .
The 0.05A 0 correction corresponds to the extra cellular area due to the projection in the micropipette used to probe the tension (which, for simplicity we do not model in our simulations). The correction A exo =A 0 [1+0.5 min (t/300 seconds,1)] corresponds to the new area delivered to the plasma membrane by exocytosis, which we take to be linearly dependent on t over 300 seconds up to a maximum delivery of membrane of 50% of A 0 (Herant et al., 2005).
In prior work (Herant et al., 2005), it was found that threshold effects determine tension during phagocytosis-driven membrane expansion, probably in part because folds and villi are enzymatically 'unzipped' to increase spare membrane availability. The critical threshold for area expansion above which surface tension increases substantially is 26% (with a a standard deviation of 11%, reflecting significant cellto-cell variability). We distinguish two contributions to the surface tension such that ␥=␥ el +␥ vis . The first comes from an elastic term that is a function of area expansion ⌬A only. We had previously found that The second term is the viscous contribution to the surface tension, deriving from the dynamical disruption of bonds that occurs when area expansion leads to the unraveling of membrane stored in folds and villi. We had previously found that We have used ␥ =75 mN m -1 second (Herant et al., 2003;Herant et al., 2005).

Numerical methods
The simulations presented in this paper were obtained by solving the model equations through a Galerkin finite element method using a mesh of quadrilaterals as described previously (Dembo, 1994;Herant et al., 2003). Briefly, the calculation is advanced over a time step ⌬t (determined by the Courant condition or other fast time scale of the dynamics) by means of five sequential operations: (1) We advect the mesh boundary according to the network flow and then reposition mesh nodes for optimal resolution while preserving mesh topology, boundaries and interfacial surfaces (Knupp and Steiberg, 1994). (2) We advect mass from the old mesh positions to the new mesh using a general Eulerian-Lagrangian scheme with upwind interpolation (Rash and Williamson, 1990). (3) We conduct diffusive mass transport and simultaneously carry out any chemical reactions. This is done according to a backward Euler (implicit) scheme coupled with a Galerkin finite element treatment of spatial derivatives and boundary conditions. (4) We use constitutive laws to compute necessary quantities such as viscosities and surface tensions. (5) Finally the momentum equations and the incompressibility condition together with the applicable boundary conditions are discretized using the Galerkin approach and the resulting system is solved for the pressure, network velocity and solvent velocity using an Uzawa style iteration (Temam, 1979). The computational cycle is repeated until engulfment is complete (the leading edges meet around the phagocytic target) or 300 seconds have elapsed (incomplete phagocytosis of 11 m beads).
Cylindrical symmetry allows the use of a two-dimensional mesh to solve the axisymmetric version of the equations described above. Numerical convergence was confirmed by checking that the results were not sensitive to variations of the tolerance of the different iterations performed by the code as well as to variations of the spatial resolution. Calculations were conducted using 64-bit arithmetic on a Linux PC workstation. The code was compiled with the Absoft Fortran 90 compiler. Post-processing was performed with a variety of publicly available software packages (SuperMongo, DISLIN, GMV, Origin) as well as with customized code.