Computational estimates of membrane flow and tension gradient in motile cells.

All parts of motile cells, including the plasma membrane, have to translocate in the direction of locomotion. Both directed intracellular membrane transport coupled with polarized endo- and exocytosis and fluid flow in the plane of the plasma membrane can contribute to this overall plasma membrane translocation. It remains unclear how strong a force is required to generate this flow. We numerically solve Stokes equations for the viscous membrane flow across a flat plasma membrane surface in the presence of transmembrane proteins attached to the cytoskeleton and find the membrane tension gradient associated with this flow. This gradient is sensitive to the size and density of the transmembrane proteins attached to the cytoskeleton and can become significant enough to slow down cell movement. We estimate the influence of intracellular membrane transport and actin growth and contraction on the tension gradient, and discuss possible 'tank tread' flow at ventral and dorsal surfaces.

Introduction

The plasma membrane plays several crucial roles in cell life: separating the inside of the cell from the environment; serving as a scaffold for regulatory and structural proteins; and organizing cytoskeletal dynamics [1]. The plasma membrane's mechanical characteristics, such as flow [2], tension [1], [3] and curvature [4] are also important for cellular phenomena, especially for cell motility. Here we mathematically and computationally examine the mechanical effect that plasma membrane flow and the associated membrane tension have on motile cell behavior.

Cell migration on surfaces is a fundamental phenomenon underlying many physiological processes [5]. When a cell migrates, its parts, including its cytoskeleton, organelles, fluid cytoplasm and plasma membrane, have to translocate forward (Fig 1A). In many types of migrating cells, this forward translocation is driven by the dynamic actomyosin network, one of the main parts of the cytoskeleton, in which nascent actin filaments appear and grow at the cell front and push the leading edge forward, while older parts of the network disassemble and contract to pull the rear forward [6]. The actomyosin network adheres to the substrate via molecular complexes which contain integrins, and which span from actin to the substrate through the plasma membrane (Fig 1A). These adhesions are crucial for transducing the effective treadmill of the actomyosin array into forward propulsion of the cell [7]. The mechanisms by which organelles and cytoskeletal components move forward are not all entirely clear, but actomyosin contractions [8], microtubule-based motors [9] and membrane tension at the cell rear [10] contribute to these processes in various cells.

Figure 1

Possible types of membrane flow.

A–E: View of the motile keratocyte cell's lamellipodium from the side. Shaded ellipsoid is the cell body. A: Growth of actin network (green) inside the lamellipodium pushes the leading edge forward, resisted by the membrane tension at the front (yellow arrows,T). At the rear, membrane tension (yellow arrows, T) pushes forward the disassembling actin networks. Besides the breaking actin network, breaking adhesions (orange rectangle) also resist rear retraction. Myosin (dark red dumbbells) powered contraction assists rear retraction. The membrane flows forward in the lab coordinate system (blue arrows) on the ventral and dorsal surfaces. Adhesions (blue rectangles) linked to the stationary actin network resist this flow at the ventral surface; transmembrane proteins (light blue ovals) resist this flow at the dorsal surface. B–E: blue (red) arrows show the membrane flow in the lab (moving cell) coordinate systems, respectively. B: One possibility is that the membrane flow is the same at the ventral and dorsal surfaces. In this case, these flows' rates are the same as the cell speed in the lab coordinate system, and the flows are zero in the cell frame. C: Example of tank-tread flow. D: In this case the membrane is transported from the rear to the front solely through the intracellular transport of membrane vesicles. The membrane flow is zero in the lab coordinate system and directed to the rear in the cell frame. E: Hypothesis about how the membrane flow can be the same on the ventral and dorsal surfaces for any different protein concentrations at these surfaces. This can be explained if the membrane flow across the leading edge (black) is obstructed. In this case, the membrane tension at the rear of the dorsal and ventral surfaces is the same, but rear-to-front gradients of tension are different along the ventral and dorsal surfaces because the same ventral and dorsal flows are resisted by different protein concentrations at these surfaces. Thus, tensions at the fronts of the ventral and dorsal surfaces are different. F: View of the motile keratocyte cell's lamellipodium from above. Shaded ellipsoid is the cell body. Insert: Cartoon of membrane flow around solid circular obstacles (proteins). Proteins attach to the cytoskeleton and/or the external environment, and so stay stationary (in the lab frame of reference). Thus, as the cell moves, the membrane is forced to flow around the proteins.

The plasma membrane enveloping the cell also has to translocate from the rear to the front. Because lipids and proteins in the mosaic membrane diffuse rapidly [11] in the membrane plane, the membrane can simply flow forward (Fig 1A). The flow of the plasma membrane can be supplemented, or even replaced, by directed intracellular movement of membrane vesicles mediated by motor-driven transport, so that endocytosis is responsible for removing plasma membrane at the rear and exocytosis for adding membrane at the front (Fig 1D). Indeed, in some cases there is evidence of polarized membrane trafficking [2], [12].

There is sometimes confusion in the literature that stems from the fact that the membrane flow looks different in the frame of the moving cell and in the lab coordinate system. In this paper, we will consider a cell steadily moving forward with the rate (Fig 1A), such as fish epithelial keratocyte [13]. We illustrate possible types of membrane flow in Fig 1B–D. The simplest possibility is if in the lab coordinate system both ventral and dorsal membranes flow forward (Fig 1B, blue arrows) at rates equal to the cell speed: . In the cell frame, there is no flow in this situation. Such a case was observed in a number of motile cells, including fibroblasts, fish keratocytes [14]–[16], leukocytes [17], and Dictyostelium amoebae [18]. A more complex possibility is a tank-tread flow in which, in the lab coordinate system, both ventral and dorsal membranes flow forward with different speeds; for example, the dorsal flow is faster (Fig 1C, blue arrows). Conservation of membrane material requires that in this case . Then in the cell frame the dorsal flow is directed forward, and the ventral flow is directed rearward (Fig 1C, red arrows; note that ). Finally, if intracellular traffic is solely responsible for forward membrane translocation, then in the lab coordinate system the membrane is stationary, while in the cell frame the membrane flows backward with equal rates at the ventral and dorsal surfaces equal to the cell speed (Fig 1D, red arrows). Interestingly, in neuronal growth cones it was observed that membrane flow is directed from the front to the rear in the cell frame [19].

The membrane flow is determined by a force that drives it. This force arises from the gradient of the in-plane membrane tension [20] (Fig. 1A), so that the tension at the front, , is higher than that at the rear, , and so the more tensed membrane at the front pulls the plasma membrane forward against weaker tension at the rear. A front-rear membrane tension difference on the order of ∼1 pN/µm was indeed measured between the cell body and the growth cone in neurons [19], where this tension gradient was accompanied by membrane flow. For comparison, average membrane tension in different cell types varies widely, from a few pN/µm in neuronal growth cones to tens of pN/µm in melanoma cells [21] to hundreds of pN/µm in rapidly moving fish keratocytes [22].

Another set of factors determining the membrane flow are the mechanical properties of the membrane. With respect to out-of-plane deformations, the membrane has complex, partially elastic properties, but in-plane, the membrane is an incompressible viscous fluid [23]. Effective resistance to the flow of the membrane relative to the substrate arises largely due to transmembrane proteins, especially integrins, and proteins with domains which insert into or bind the membrane [24]. According to the mosaic model of plasma membrane structure, up to 50% of the membrane surface is occupied by such proteins, with lipid molecules filling the rest [25]. In general, lipids and membrane proteins are free to diffuse within the membrane, however, a subset of transmembrane proteins are restricted by binding to cytoskeletal structures [26], [27]. For example, membrane-associated proteins such as fodrin, plectin, and ankyrin attach to actin [28]. ERM and Ena/VASP are other examples of the protein families that associate both with actin and membrane. Given the abundance and strength of cell-substrate adhesion sites, they also contribute extensively to membrane-cytoskeleton interactions.

After a brief surge of interest in the role of membrane transport in cell locomotion [2], the efforts of motility researchers were largely concentrated on actin and myosin dynamics. The implicit assumption in most of the literature is that the plasma membrane flows forward effortlessly in motile cells. However, this assumption was never critically examined. Indeed, in the large part of the lamellipodium/lamellum (dynamic motile appendage of the cell, see Fig. 1), the actomyosin network is either almost stationary relative to the substrate [13] or undergoes rearward flow [29] in the direction opposite to that of cell migration. If a fraction of the transmembrane proteins are associated with this stationary or rearward moving cytoskeleton, then these proteins are effective buoys that obstruct the forward membrane flow. How much resistance these buoys exert on the flow is an open problem (an initial foray into this problem was recently made [24]). Recently, a few studies have reported that inhibition of membrane trafficking in several cell types, for example fibroblasts, endothelial cells and Dictyostelium, reduces persistent migration [30], [31]. One possible explanation is that membrane transport could be rate-limiting for cell forward translocation. Here, we compute the geometry of the plasma membrane flow and associated membrane tension gradients, and we predict conditions for high and low membrane tension gradients and for equal versus unequal ventral and dorsal flows.

Results

1. Model

1.1 Estimate for the drag coefficient of a single cylindrical protein

Variables and parameters in the models are defined and explained in Table 1 and Table 2. Saffman and Delbruck [32] derived an analytic estimate for the drag coefficient of a single cylindrical protein embedded in a flat, infinite sheet of membrane:

Table 1

Parameters in the models.

Variable notation	Value	Meaning	Reference
		cell speed	[13]
		viscosity of aqueous medium	[32]
		two-dimensional membrane viscosity	[51]
		Adhesion complex diameter	[43]
		Transmembrane protein diameter	[44]
	0.577	Euler-Mascheroni constant	[42]
		Order of magnitude ide-to-side lamellipodial width	[13]
		front-to-rear lamellipodial length	[13]
	n/d, varies	number of transmembrane proteins linked to actin	NA
φ		area fraction of the transmembrane proteins linked to actin	Assumed
		drag coefficient of a single cylindrical protein embedded in a flat, infinite sheet of membrane	Computed in this paper
		parameters defining the lamellipodial shape	Chosen in this paper to qualitatively fit the observed shape
,,	depends on other parameters	effective combined viscous drags of the transmembrane buoys in the ventral and dorsal surfaces	Computed in this paper
	Defined in combination with ,	intracellular membrane vesicles per second per surface	[45], [46]
	Defined in combination with ,	average membrane area per vesicle	[45], [46]
		characteristic membrane tension that stalls protrusion	[22], [41]
	Defined in combination with other parameters	characteristic retraction speed of the cell rear,	[10]
	Defined in combination with other parameters	characteristic membrane tension at the rear,	[10]

Table 2

Variables in the models.

Variable notation	Meaning	Dimension
	Speed of the membrane flow at the ventral surface
	Speed of the membrane flow at the dorsal surface
	Membrane tension at the rear
	Membrane tension at the front
	Local membrane flow rate
	Coordinates describing lamellipodial boundary
	Angle between the direction of cell motion and the normal to the boundary	n/d

Here is the three-dimensional viscosity of the membrane, is the thickness of the membrane, is the viscosity of the fluid surrounding the membrane, is the diameter of the protein, and  = 0.577 is the Euler-Mascheroni constant. Note that the two-dimensional membrane viscosity can be related to from equation (1) by the expression , so:

1.2. Flow and tension distributions on the flat membrane surface

We model transmembrane proteins or protein complexes, such as adhesions, attached to the actin network as rigid cylindrical solid objects embedded in the membrane (Fig 1F). We consider steadily moving cells of a constant shape, like fish keratocytes (Fig 1F). In keratocytes, the motile appendage – the lamellipodium – is a flat actomyosin network with a characteristic canoe-like shape (Fig 1F). Most of this network, other than a small part at the very rear of the cell, is almost immobile relative to the surface [33], so we assume that in the lab coordinate system the transmembrane cylindrical ‘buoys’ are immobile. Following previous studies [24], [34]–[36], we represent the membrane as a thin layer of a continuous incompressible viscous fluid. Since the membrane is extremely thin, we treat it as a 2D fluid, and so we also treat the membrane domain of a protein as a circular obstruction in the 2D fluid domain. As the transmembrane buoys are held stationary during cell migration while the membrane advances forward in the lab coordinate system, the membrane and proteins interact mechanically to produce a non-uniform membrane flow, as sketched in Fig 1F. We model the flow as 2D viscous Stokes flow:

Here is the local membrane tension, and is the local membrane velocity in the lab coordinate system. Note that in two dimensions, the tension (which can be considered as pressure with opposite sign) has the dimensions Force/Length instead of the usual Force/Length2 in the 3D Stokes equations. Similarly, the membrane viscosity has the dimensions of Force × Time/Length, rather than the usual 3D viscosity dimensions of Force × Time/Length2. Another assumption we make is that we neglect attachment (between membrane and cytoskeleton) energy contribution to the membrane tension, so in equation (3) is the in-plane membrane tension. The length scale in this model is the characteristic length of the lamellipodium, , and the speed scale is the characteristic keratocyte speed [13]. The natural membrane tension scale that we use in numerical calculations is .

1.3. Cell shape, boundary conditions and protein distributions

We solve equation (3) for a flat membrane domain, which represents either the dorsal or ventral membrane surface. We define the shape of this domain similar to that of the keratocyte lamellipodium (Fig 2A) and represent the boundary of this domain by the equations: for the front half of the boundary, and for the rear half in the Cartesian coordinate system with the y-axis passing from the rear to the front through the middle of the cell, where . Here the length is non-dimensionalized using the lamellipodial length scale defined above.

Figure 2

Lamellipodial shape and protein and velocity distributions.

Shape of cell's lamellipodium given by equations (4–5). A: Sketch of velocity boundary condition (equation (6)) given by the graded radial extension model. Velocity at the boundary is normal to the boundary, and decreases in magnitude from a maximum when the normal points in the direction of cell motion to zero when the boundary is tangent to the direction of motion. B: Membrane boundary with leading and rear edges in red. We integrate tension over these red curves to compute the average tension at the front and rear of the cell, which allows us to compute the tension drop in the membrane. C–D: Randomly generated placements of transmembrane proteins distributed uniformly throughout the membrane (C) and distributed within 1 µm of the cell front (D).

We follow the Graded Radial Extension model for the steady lamellipodial locomotion [16] and require that at each point along the lamellipodial boundary the local membrane velocity is directed normal to the boundary, with magnitude where is the angle between the direction of cell motion and the normal to the boundary, as shown in Fig 2A. This boundary condition ensures that the plasma membrane translocates forward in a manner that preserves a steady cell shape. In non-dimensional units, .

We use the no-slip boundary condition on the 2D circular surface of each protein. After the Stokes flow equation (3) with these boundary conditions is solved numerically (Methods), we calculate the average membrane tension at the front and rear edges of the cell (shown in Fig. 2B) as: where and are the average membrane tensions at the front and rear edges, respectively, and and are the arclengths along the front and rear boundaries, and , respectively.

We investigate the effect of transmembrane proteins or protein complexes on membrane tension by randomly placing circular obstructions in the interior of the membrane. We generate protein positions according to one of two distributions (Fig 2C,D). Proteins are either uniformly placed throughout the interior of the membrane (Fig 2C), or are restricted to being near the front of the cell (Fig 2D). We chose this latter distribution to account for the observations that at the ventral surface the adhesion complexes and respective traction are concentrated near the leading edge of the cell [37], [38] and that the membrane is loaded with proteins at the leading edge [39], [40]. For these front-loaded proteins, we required that each protein be entirely in the front half of the cell (defined by the line y = 0), and that the protein be within 1 µm of the front half of the cell boundary. For both the uniform and the front-loaded distributions, protein placement is generated in several steps: 1. Points are sampled from a uniform random distribution covering the entire cell. These points will become the centers for our circular proteins. 2. For each generated point, we check to ensure that it is in the correct region of space (either the cell interior for the uniform distribution or the front 1 µmof the cell for the front-loaded distribution). Any points failing this check are replaced by new random points. For numerical accuracy, we also require that all points be at least 0.25 µmaway from the membrane boundary. 3. To ensure that none of the circular proteins overlap with one another, we impose a short range repulsive interaction on neighboring points. This repulsion moves the protein centers until all centers are a distance of at least one diameter apart from one another.

1.4. Tank-tread flow at ventral and dorsal surfaces

So far, we have described the membrane flow in the flat plane of either the ventral or the dorsal surface of the cell in the case when the ventral and dorsal flows are the same. The computational 2D model is useful for estimates of the membrane tension gradient and transmembrane protein resistance to the flow, and for visualizing the non-uniform membrane flow in the membrane plane. In the more general cases when either the flows are different or effects of actin mechanics and intracellular membrane transport are included, full simulation of the 2D flow would be difficult because of the complex actin mechanics and geometry of the curved plasma membrane at the rear and leading edges of the cell and unknown distribution of intracellular transport at the edges. However, simple estimates in the general cases can be made if we consider a simplified 1D model of the plasma membrane in the side view of the cell in the lab frame (Fig 1A). In this model, we approximate the flows in the ventral and dorsal planes as uniform (this approximation is in general very good, as shown in the results reported below) with rates and , respectively. Let and be the effective viscous resistance (due to the transmembrane proteins) of the ventral and dorsal surfaces, respectively. Then the membrane tension gradient between the front and rear of the cell is: for the dorsal and ventral surfaces. Note that coefficients and have different dimensions (Force × Time/Length3) than that of the protein viscous drag . This can be best understood in the example of the linear approximation, when the drag imposed by proteins is additive. In this case, if proteins, each characterized by the protein viscous drag , are distributed over the rectangular region of the membrane in the presence of the uniform flow with speed , with the leading and rear edge length and the rear-to-front length , then the total resistance force from the proteins is equal to , and the rear-to-front difference in membrane tension can be estimated is . Introducing average density of proteins , we can write: , and . From comparison with equation (8), it is clear that coefficients and have dimension Force × Time/Length3 because the protein viscous drag has dimensions Force × Time/Length, and the average density of proteins has dimensions 1/Length2.

We complement equations (8) with a third equation for conservation of membrane material for the plasma membrane:

We will solve the three equations (8–9) for the three unknowns () below.

1.5. Effects of actin pushing at the front, actin contraction at the rear and intracellular membrane transport

Another similar 1D model considers the effect of intracellular membrane transport and force balances at the front and rear of the motile cell on both the membrane flow and the tension gradient. For simplicity, in this model we consider equal viscous resistance for the ventral and dorsal surfaces, and therefore equal respective membrane flows equal to at these two surfaces. In this case, the force balance for the membrane flow and membrane tension gradient has the form:

Let intracellular membrane vesicles with average membrane area per vesicle per second per surface (dorsal or ventral) be transported from the rear to the front of the cell (Fig 1D) per one micron of the leading edge. Then the protrusion length per time is equal to . Therefore, the cell speed, /, is equal to:

The two equations (10–11) for the four unknown variables () have to be complemented by the force balances at the front and rear of the cell. At the front, we use the force-velocity relation for the actin network polymerizing against the membrane load:

Here is the free polymerization rate in the absence of load, and the expression in brackets is the force-dependent factor that accounts for the slowing down of protrusion by membrane tension. Parameter is the stall tension at which the protrusion halts. For simplicity we use the linearized force-velocity relation, which is an approximation to the non-linear relations measured and fitted for keratocytes [13], [41]. At the rear, we use the force-velocity relation for the disassembly-weakened actin network being crushed and hauled forward by the membrane tension, which was suggested and fitted to the data in [10]:

Here is the rate of the cell rear retraction generated when the membrane tension at the cell rear is equal to parameter . Below, we solve the four equations (10–13) analytically, allowing us to examine the effects of intracellular transport and actin mechanics on membrane flow and tension gradients.

2. 2D membrane flow and tension distributions in the plane of the membrane

We performed numerical simulations on flat membrane domains, which can represent either ventral or dorsal membrane surfaces. For both uniform and front-loaded protein distributions, we performed simulations with total numbers of proteins ranging from 0 to 350. Fig 3 shows representative simulation results for both uniform and front-loaded distributions. Figs 3A,C and B,D show maps of magnitude and directions of flow and tension distributions, respectively, in sample simulations. More results at other protein densities and distributions are shown in Figures S1–S6. Interestingly, the flow concentrates in narrow meandering effective channels on the membrane surface. In these channels, which are effectively paths devoid of the attached proteins, the flow rate could be a few-fold higher than the cell speed.

Figure 3

Computed membrane velocity and tension.

A–B: Computed membrane velocity (A and C, units of cell speed) and tension (B and D, units of pN/µm) for 350 uniformly distributed (A,B) and front loaded (C,D) proteins. Proteins are shown in white in A,C and black in B,D.

For each simulation, we calculated the average tension drop across the cell as described above. We performed 10 simulations at each number of proteins, and computed the average front to rear tension drop over all runs at each protein number. The membrane tension increases almost linearly from the rear to the front for the random uniform distribution of the attached transmembrane proteins. When the proteins are front-loaded, the same proteins number as in the uniform case produces a greater total average tension increase from the rear to the front. This suggests that as the proteins become more tightly packed, they restrict the flow more than an equal number of well-spaced proteins, as expected.

3. Dependence of the membrane tension gradient on transmembrane protein size and area fraction in the 2D model

No matter how much resistance the transmembrane proteins generate, the average membrane flow rate in the 2D model is of the order of the cell speed. It is the membrane tension gradient that must adjust to maintain this flow rate. The Saffman-Delbruck formula (2) allows us to estimate the drag due to one transmembrane protein. For the parameter values given in Table 1, the drag coefficient is . In the linear approximation, the drag imposed by proteins is additive, and the rear-to-front difference in membrane tension can be estimated as: where is the lateral width of the lamellipodium.

The computations of the 2D flow allow us to investigate how the rear-to-front difference in membrane tension depends on protein density. We plotted respective results in Fig 4A,B. At low protein densities, there is good agreement between the analytic prediction (14) and our numerical results (Fig 4A). Our numerical simulations show that this linear approximation, however, is valid only up to the density of about one ‘buoy’ per square micron (or area fraction φ occupied by proteins less than 1%), in which case , which is very low compared to the measurements in the neuronal growth cone and which is well below the experimental noise in the tether-pulling experiments.

Figure 4

Computed membrane tension gradient.

A: Computed average membrane tension drop as a function of the area fraction of proteins for uniform distributions of proteins (red) compared with the analytic linear prediction (blue). B: Computed average membrane tension drop as a function of the area fraction of proteins for front-loaded proteins (red) compared with the analytic linear prediction from equation (16) (blue). The data in (A–B) are for variable number of proteins 0.1 µm in diameter. C: Numerical results (red) and equation (16) fit (blue) for tension drop as a function of protein diameter. The data points include front loaded simulations and show individual simulation results from the runs with varied diameters and protein number equal to 100. Other parameters are listed in Table 1.

At higher protein densities, however, both the uniform and front-loaded protein distributions quickly produce a much larger tension drop than would be predicted by the linear theory because closely spaced proteins affect the flow in a non-additive way, which is a well-known hydrodynamic phenomenon. Indeed, there is a certain analogy between the 2D membrane flow through a ‘maze’ of disks and 3D groundwater flow through a maze of spheroid objects. In the latter case, the well-known Kozeny-Carman equation describes the pressure drop for a fluid of viscosity flowing a distance through a porous medium at an average velocity [42]:

The porous medium is assumed to be made up of randomly placed, immovable particles that can be described by their sphericity , which is a measure of how spherical the particles are ( = 1 for spheres), the diameter of a sphere having the same volume as one of the particles, and the volume fraction φ, which is the fraction of total volume that is occupied by the solid particles.

We fitted our numerical data to a similar equation of the form:

Here is a small correction that stems from the fact that there is a small viscous resistance to the membrane flow even in the absence of the protein obstacles, is the dimensionless parameter to be found from fitting, and is the area on which disc-like obstacles are randomly placed. We minimized the least squared error over all our simulation runs, and found that for equation (16) is an excellent fit for our simulations (Fig 4B). Note that for values of the area fraction φ below a few percent, the numerical data shows agreement with the linear tension-density relation, not with the Kozeny-Carman quadratic relation. The reason is that the Kozeny-Carman relation comes partially from considering flow through narrow, twisting channels, and the geometry of the flow becomes quite different at very low protein density. However, tension gradients of biological interest develop at higher protein density, so for practical purposes equation (16) is valid. To additionally test equation (16) we did computations with proteins varying the protein diameter. Note that when parameter varies, so does parameter φ, so we plotted the results as shown in Fig 4C. The fit generated by equation (16) to this numerical data is good.

It is unknown what fraction of the transmembrane proteins covering a few tens of per cent of the plasma membrane area is attached to the cytoskeleton. It would be reasonable to assume that about 10% of the whole membrane area is covered by such attached proteins. Assuming φ = 0.1, we use equation (16) to estimate the tension gradient for (respective protein density is ): .

This estimate is very similar to the measurement in the neuronal growth cones [17].

So far we reported the simulations using the value of characteristic for nascent adhesion protein complexes at the leading edge [43]. Single proteins are smaller; the characteristic size is [44]. With such a value of this parameter, the estimate for the tension gradient for φ = 0.1 becomes , which is a great tension gradient, comparable in fact to the average membrane tension in many cells. (Respective protein density is ). Here we report two simulations with the parameter value that we performed for front-loaded proteins to verify that the predictions made from our simulations with are consistent with simulations at a smaller protein diameter. Figures S7 and S8 show individual simulations with 2000 and 3000 transmembrane proteins of diameter , respectively. These two simulations had area fractions of φ = 0.023 and 0.036, respectively. These simulations gave membrane tension drops of 1.03 and 2.65 , respectively, on the same order as the values predicted by our analytical fit. Further work is necessary to determine whether our analytic fit is as good in this small diameter regime as it is for larger diameters.

Since the membrane is a very viscous fluid, the tension and velocity distributions in the plane of the membrane reach steady state for a given transmembrane protein distribution almost instantaneously [1], so we can predict the tension and flow for a given transmembrane protein distribution without regard to the timescale over which these adhesions move relative to the cell boundary. However, to verify that our steady-state computations of membrane tension were valid over the timescales relevant for cell motility, we simulated an individual membrane with 100 embedded transmembrane proteins for 200 seconds. As in our other simulations, we computed the membrane flow velocity and tension in the lab frame, but we allowed the membrane boundary to move forward relative to the adhesions at the speed of the cell. To remain consistent with our simulations of transmembrane proteins concentrated near the leading edge, as the cell moved forward we removed any transmembrane proteins that became more than 1 µm from the leading edge, and replaced those proteins with new, randomly generated ones within 1 µm of the leading edge. Movies S1 and S2 (online) show the results of this simulation for both flow and tension. In the movies, all results are given in the lab frame, but we show a field of view that keeps pace with the boundary of the membrane. Figure S9 shows in red the spatial average membrane tension drop over the cell as a function of time, and in blue shows the predicted tension drop from the analytical fit for this cell, computed from the average value of φ for this cell over the 200 seconds of simulated time.

4. Tank-tread flow versus no membrane flow in the framework of the moving cell

Examination of the simple 1D model of membrane flow in the side view of the cell raises an interesting issue that stems from the fact that in most motile cells, there is a very high density of attachments of adhesion molecular complexes that go through the ventral membrane to the actin network, while at the dorsal membrane, the molecular nature and number of the actin-membrane attachments is likely very different. This suggests that the membrane flows differently at the ventral and dorsal surfaces, as the effective viscous drag coefficients are scaled by protein densities. Indeed, the analytical solution for the three linear equations (8–9), implies that if there is the same density of transmembrane proteins attached to the cytoskeleton on the dorsal and ventral surfaces (so that ) then the dorsal and ventral flows would be the same in the lab coordinate system (meaning no flow in the cell frame). It would be very surprising if there is the same density of transmembrane proteins attached to the cytoskeleton on the dorsal and ventral surfaces and so in general for , the model predicts a tank-tread flow. How can we reconcile this prediction with the fact that in a number of cell types no flow in the cell frame is observed?

There are reports on abnormally low diffusion coefficients of membrane-associated proteins at the very leading edge of motile cells [39], [40], suggesting that the leading edges harbors especially high concentrations of proteins, many of which could be linked to the cytoskeleton. If this is the case, then the tank-tread membrane flow across the leading edge would be largely obstructed by this ‘protein crust’, and conservation of lipid number would require that both dorsal and ventral flow in this case be equal to the cell speed in the lab coordinate system: . Interestingly, this means that the tensions at the leading edge at ventral and dorsal membrane surfaces are different (Fig 1E):

An experiment testing this prediction – pulling membrane tethers from ventral and dorsal leading edge and measuring respective forces – would be very difficult, but not impossible.

5. Influence of membrane-cytoskeleton attachments on cell locomotion

To investigate the effects of intracellular membrane transport and force balances at the front and rear of the motile cell on the membrane flow and tension gradient, we solve four linear equations (10–13). The solution has the form:

This formula shows that the cell speed is accelerated by intracellular membrane transport and is slowed down by the mechanical resistance of the rear (the first term in the second bracket) and of the front (the second term in the second bracket). The influence of intracellular membrane transport is determined by the relative magnitude of the effective speed of this transport, , and of the membrane flow speed generated by the stall tension (membrane tension that stalls the actin protrusion) against the transmembrane viscous resistance, .

The overall rates of transport between internal membranes and the plasma membrane vary among different cell types. In Dictyostelium amoeba it takes a few minutes to replace an area equal to the cell surface area [45] whereas in fibroblasts it takes approximately one hour [46]. If it takes 500 s to replace the whole plasma membrane, while it takes about 100 s to move one body length, then equation (11) suggests that . From [41], , and our estimate of the transmembrane viscous resistance above is . Then, and so the mechanical influence of intracellular membrane transport on the cell speed is minuscule. However, if the stall tension is an order of magnitude lower in non-keratocyte cells, while the membrane resistance is higher, then it is possible that intracellular membrane transport becomes the limiting factor in cell locomotion.

Based on the analysis of the keratocyte fragments in [10], the ratio of the stall force to the free polymerization rate is of the same order of magnitude as the ratio of the characteristic tension and motion rate at the rear , so the factor . Thus, the influence of the membrane flow on the cell speed is determined by the ratio of the free polymerization rate to the membrane flow speed generated by the stall tension against the transmembrane viscous resistance, . The free polymerization rate could be a few-fold greater than the observed cell speed. Taking , we estimate the ratio , depending on the value of . Thus, the effect of the transmembrane proteins' viscous drag is small but could be noticeable in slowing down cell movement. Note that in the limit of extremely high viscous drag , , tension at the front is almost at stall, and tension at the rear is very small.

Membrane tension at the front can also be easily estimated from equations (10–13):

This tension is relieved by intracellular membrane transport (negative term in the first bracket), but this effect is quite small for keratocytes because according to our estimates . This transport, however, can significantly relieve tension for slow-moving cells. The greater the transmembrane proteins' viscous resistance , the greater the membrane tension at the front, and in the limit of high viscous resistance the membrane tension approaches the stall tension. Our estimates show that if , then the transmembrane proteins' viscous drag increases the membrane tension at the leading edge by at most ∼10%. The approximate value of the leading edge tension in this case is , and at the rear the tension is almost the same.

Discussion

In this paper, we estimated the gradient in membrane tension from the rear to the front of the cell necessary to maintain the membrane flow that keeps up with cell locomotion. This gradient is predicted to be of the order of a few pN/µm if the area fraction of the transmembrane protein complexes attached to actin is about 10% and if their diameter is of the order of 0.1 µm. However, if at the same area fraction the protein diameter is of the order of 20 nm, then the tension gradient becomes tens of pN/µm. The apparent average membrane tension extracted from tether-force measurement values ranges from a few pN/µm in neuronal growth cones [21] to hundreds of pN/µm in keratocytes [22]. Thus, the tension gradient could be comparable to the average membrane tension and detectable in principle if membrane tethers are pulled from the front and rear of the cell. The caveat to interpretation of such measurements though is that it is not easy to separate the in-plane membrane tension from the effective tension due to membrane-cytoskeleton attachment. Another interesting model prediction is that the membrane flow could concentrate in narrow meandering effective channels on the membrane surface.

Furthermore, this tension gradient is generated by the actin network pushing the leading edge forward, and this force is of the order of hundreds of pN/µm [22], [41]. Thus, in principle a non-negligible part of this pushing force has to support membrane transport, in addition to pulling up the cell rear and overcoming resistance of the extracellular environment. We estimate that if the area fraction of the actin-attached transmembrane proteins becomes higher than 10%, the force needed to generate the membrane flow becomes non-negligible, and at 20% it becomes rate-limiting. With such a high area fraction, the directed intracellular membrane traffic can also become the principal mechanism of membrane recycling.

Our model predicts that in general, if the plasma membrane is an interconnected fluid domain then a tank-tread mode of the membrane flow has to be observed. The fact that in the majority of cases where the membrane flow was measured the flow does not show such a pattern [14]–[18] may indicate that at the leading edge the membrane is not fluid but rather a non-flowing (or very slowly flowing) ‘protein crust’ interspersed with lipids trapped in protein ‘corrals’. We predict that in this case the membrane tensions at the ventral and dorsal leading edges are different.

Our model has a number of limitations that stem from simplifying assumptions that we made in order to increase its transparency and tractability. We considered a constant rate of intracellular membrane traffic from the rear to the front, implicitly assuming that the rates of exocytosis at the front and endocytosis at the rear are not rate limiting. This issue has to be researched further. In addition, measurements have revealed that rates of endocytosis decrease as membrane tension increases [47], and exocytosis is likely to also be tension-dependent [20]. In the future, the model can be extended to include a more detailed and realistic sub-model for intracellular membrane recycling. We also did not explicitly consider a possible contribution to the membrane tension from adhesion breaking at the very rear of the cell, however, mathematically, this contribution can be included in model parameters and .

The shear between the substrate and ventral membrane for the motile cell is negligible. Respective mechanical contributions can be estimated as water viscosity times cell speed divided by the characteristic distance between the ventral surface and the substrate and multiplied by the cell length: . On the other hand, we likely underestimate the membrane tension gradient because we did not consider the effect of hydrodynamic interactions with multiple actin filaments, some of which are proximal to the plasma membrane. This effect was recently investigated [24], and it is not negligible. The authors of [24] computed analytically the membrane flow and forces generated by a number of proteins anchoring the actin network to the dorsal and ventral membrane sandwiched between two solid surfaces and found that the effective friction between the proteins anchoring the actin network to the membrane and the surface can propel a cell in the absence of any direct adhesion mechanism with velocities comparable to velocities of adhesive cells. The geometry of the problem, methods and goals of that important study are different from those in our paper, yet [24] also demonstrates that the membrane flow and associated forces could be important in the cell motility process.

In addition, effects of inhomogeneous distributions of many types of lipids and proteins in the plane of the plasma membrane and of various diffusion coefficients of the lipids and proteins on the flow and tension gradient in many cells also remain unclear. Moreover, most cells migrate in a non-steady way, so time dependent non-steady solutions for the membrane flow equations have to be explored to investigate membrane recycling in these cells. Also, there is often a significant retrograde flow of the actin network near the leading edge of motile cells. Many transmembrane proteins at the dorsal surface attached to actin should then flow to the rear of the cell, while some adhesions at the ventral surface remain stationary in the lab frame, and some are also drifting to the rear [29]. Our model can be easily extended to simulate this situation, but qualitatively it is clear that the effect of the actin retrograde flow will lead to further increase of the resistance to the forward membrane flow and of the tension gradient.

Last, but not least, perhaps the most interesting open problem stems from the fact that in most cells the plasma membrane enveloping the cell is not pulled taut, but rather is folded, ruffled and invaginated [48] in multiple local membrane reservoirs, i.e. cup-shaped caveolae [49]. These folds enable motile cells to undergo fast changes in cell surface area [48], and it is likely that a significant effect of the membrane fold dynamics on the forward membrane flow and tension gradient will be to ease the flow and decrease the necessary membrane tension gradient. Future modeling will clarify how important this effect is.

Methods

The numerical method used to produce the solutions to Stokes equations is the method of regularized Stokeslets. This method was introduced by Cortez, and the details of it are given in [50]. The method uses a boundary integral approach to approximate the solution to Stokes equations driven by a discrete collection of prescribed point forces. The velocity and pressure at any point can be computed by summing the contributions to the solution from each force. Alternatively, the method can be used to find the solution to Stokes equations for a set of points with prescribed velocities. For this second approach, a 2N×2N linear system (where N is the number of points) is solved to compute the force at each point. Those forces can then be used to compute the solution at any point in the plane. We used this second approach to prescribe the velocity of the membrane at the boundary of the membrane and at the boundary of each transmembrane protein.

Discretizing the cell

In order to use the method of regularized Stokeslets to compute membrane flow, we first needed to discretize both the cell boundary and the boundaries of the transmembrane proteins into a collection of points. To determine the spacing between neighboring points in our discretization, we started by discretizing the transmembrane proteins. For most of our simulations, we used circular proteins with a diameter of d = 0.1 µm. At this radius, the numerical method worked well with 25 equally spaced proteins around the boundary of the circle. This gave a spacing of s  =  πd = 25 = 0.0126 µm between neighboring points. We used this fixed spacing for all our simulations, including on the cell boundary and on circular proteins of varying radii. As suggested in [50], we use a blob function φ with the spreading parameter ε chosen to scale with this spacing s. We found that simply using ε = s worked well. Discretizing the cell boundary required us to modify the shape of the cell. The original cell shape is given by equations (4,5). To generate equally spaced collections of points along each of these curves, it is necessary to compute each curve's arc length, which requires numerical evaluation of the integral . The derivatives of (4,5) are both singular near x = 0 and x = 1, which made this arc length computation highly unstable. To fix this, we modified both curves, defining the shape of the cell as: , for a small value of the parameter δ. In our simulations, we took δ = 0.00001, which was enough for us to compute the arc length without difficulty.

Computing tension drops

With the cell boundary and proteins discretized, it is easy to use the method of regularized Stokeslet. We simply prescribed velocities at each of our boundary points according to the

GRE model, and prescribed a zero velocity at each point on a protein boundary. Inverting the mobility matrix as described in [50] gives us the appropriate force at each boundary point, from which we can compute pressure and membrane flow everywhere within the cell. Our main use of this technique was to integrate pressure at the leading and rear edges of the cell to determine average lead and rear pressures, from which we computed the tension drop over the length of the cell. The only complication to this approach was that the numerical method causes a sharp jump in the pressure very near the cell boundaries, as shown in Figure S10. This jump causes evaluations of the pressure on the boundary itself to be inaccurate. To rectify the inaccuracies caused by the jump at the boundary, we simply used pressure values that were very close to, rather than on, the leading edge. To integrate the pressure at the leading edge, we instead integrated over a curve parallel to the leading edge, but offset a distance of 4s = 0.0503 µm towards the interior of the cell, and similarly for the trailing edge. As can be seen in Figure S10, computing pressure values offset by this small amount was sufficient to avoid error caused by the jumps at the boundaries.

Computed membrane flow in cell without transmembrane proteins. Black lines and arrows give streamlines of the membrane velocity, while the color plot gives the local speed. The speed in this and all next supplementary figures is reported in units of the steady cell speed.

(PDF)

Click here for additional data file.

Computed membrane velocity (top) and tension (bottom) for 30 uniformly and randomly distributed proteins of diameter 0.1 µm. Proteins are shown in white in the top figure and black in the bottom figure. The speed in this and all other figures is reported in units of µm/s. The tension in this and all next supplementary figures is reported in units of pN/µm.

(PDF)

Click here for additional data file.

Computed membrane velocity (top) and tension (bottom) for 30 front-loaded proteins of diameter 0.1 µm. Proteins are shown in white in the top figure and black in the bottom figure.

(PDF)

Click here for additional data file.

Computed membrane velocity (top) and tension (bottom) for 100 front-loaded proteins of diameter 0.04 µm. Proteins are shown in white in the top figure and black in the bottom figure.

(PDF)

Click here for additional data file.

Computed membrane velocity (top) and tension (bottom) for 100 front-loaded proteins of diameter 0.1 µm. Proteins are shown in white in the top figure and black in the bottom figure.

(PDF)

Click here for additional data file.

Computed membrane velocity (top) and tension (bottom) for 100 front-loaded proteins of diameter 0.2 µm. Proteins are shown in white in the top figure and black in the bottom figure.

(PDF)

Click here for additional data file.

Computed membrane velocity (top) and tension (bottom) for 2000 front-loaded proteins of diameter 0.02 µm.

(PDF)

Click here for additional data file.

Computed membrane velocity (top) and tension (bottom) for 3000 front-loaded proteins of diameter 0.02 µm.

(PDF)

Click here for additional data file.

The membrane tension drop as a function of time. The membrane tension drop over the cell as a function of time (red), and the predicted tension drop from the analytical fit for this cell (blue), computed from the average values in the simulation of the moving cell over 200 seconds of simulated time. About 100 front-loaded proteins of diameter 0.1 µm are stationary, appear at the leading edge and disappear at the rear edge of the protein-loaded band, while the cell moves forward.

(PDF)

Click here for additional data file.

Pressure near cell boundaries. A: Pressure (pN/µm) in a simulated cell, with blue line through middle of the cell. B: Cross section plot showing pressure in cell along the cell's midline. Note that the pressure jumps sharply near both the rear and leading edges. C–D: Close up of the pressure cross section near the rear (C) and front (D) of the cell. Note the sharp jumps in pressure over very small spatial scales. In red (dots and dashed line), we show the position 4s≈0.05 µm away from the edges used for actual computation of edge pressures. Note also that membrane tension is equivalent to the pressure with the minus sign.

(PDF)

Click here for additional data file.

Membrane flow speed distribution in steadily moving cell. Results of the simulations with about 100 embedded transmembrane proteins for 200 seconds. As in other simulations, the membrane flow velocity and tension are computed in the lab frame, but the membrane boundary is allowed to move forward at the speed of the cell. As the cell moves forward, any transmembrane proteins that are displaced more than 1 micron from the leading edge are removed and replaced with new, randomly generated ones within 1 micron of the leading edge. The results are shown in a field of view that keeps pace with the boundary of the membrane.

(MP4)

Click here for additional data file.

Membrane tension distribution in steadily moving cell. This movie shows the results for the tension of the simulations described in the legend of the Movie S1.

(MP4)

Click here for additional data file.