Computational estimates of mechanical constraints on cell migration through the extracellular matrix

Cell migration through a three-dimensional (3D) extracellular matrix (ECM) underlies important physiological phenomena and is based on a variety of mechanical strategies depending on the cell type and the properties of the ECM. By using computer simulations of the cell’s mid-plane, we investigate two such migration mechanisms—‘push-pull’ (forming a finger-like protrusion, adhering to an ECM node, and pulling the cell body forward) and ‘rear-squeezing’ (pushing the cell body through the ECM by contracting the cell cortex and ECM at the cell rear). We present a computational model that accounts for both elastic deformation and forces of the ECM, an active cell cortex and nucleus, and for hydrodynamic forces and flow of the extracellular fluid, cytoplasm, and nucleoplasm. We find that relations between three mechanical parameters—the cortex’s contractile force, nuclear elasticity, and ECM rigidity—determine the effectiveness of cell migration through the dense ECM. The cell can migrate persistently even if its cortical contraction cannot deform a near-rigid ECM, but then the contraction of the cortex has to be able to sufficiently deform the nucleus. The cell can also migrate even if it fails to deform a stiff nucleus, but then it has to be able to sufficiently deform the ECM. Simulation results show that nuclear stiffness limits the cell migration more than the ECM rigidity. Simulations show the rear-squeezing mechanism of motility results in more robust migration with larger cell displacements than those with the push-pull mechanism over a range of parameter values. Additionally, results show that the rear-squeezing mechanism is aided by hydrodynamics through a pressure gradient.


S1.1. Solution method: regularized Stokeslets
Because the Stokes equations are linear, boundary integral methods can be used to solve the system in Eqs. (S1.1) and (S1.2). In particular, the solution to the Stokes equations in free space can be determined by convolving the external force function f against the free space Green's function for Stokes flow over the interface or immersed structure where the force is applied [2]. However, the singularity of the Green's function at the point where the force is applied must be treated with accurate quadrature rules, or the singular integral must be regularized. We choose to regularize the problem with the method of regularized Stokeslets [3]. In particular, we suppose that the external force of strengthF 0 is distributed over a small ball centered on a point x 0 (rather than confined exclusively to the point x 0 ), so that (S1. 3) The solution to the Stokes equations with external forcing given by Eq. (S1.3) can be derived from the form of the "blob" or "cutoff" function φ . For example, if φ (x) = 3 3 2π( x 2 + 2 ) 5/2 , (S1.4) then p (x, x 0 ) = 1 2π (F 0 · (x − x 0 )) r 2 0 + 2 2 + r 2 0 + 2 ( r 2 0 + 2 + )(r 2 0 + 2 ) 3/2 , (S1. 5) and u (x, x 0 ) = −F 0 4πµ ln r 2 0 + 2 + − ( r 2 0 + 2 + 2 ) ( r 2 0 + 2 + ) r 2 0 + 2 + 1 4πµ (F 0 · (x − x 0 ))(x − x 0 ) r 2 0 + 2 + 2 ( r 2 0 + 2 + ) 2 r 2 0 + 2 (S1. 6) are the pressure and velocity that result from the force in Eq. (S1.3), where r 0 = x − x 0 . The derivation of these expressions can be found in [3]. An equivalent formulation is to write the result of Eq. (S1.6) as a matrix-vector product, so that u(x) = M(x, x 0 )F 0 , (S1. 7) where the matrix M is determined from the locations of the points x and x 0 .
The key point here is that given the external forcing term(s)F 0 , the velocity at any point x in the domain can be determined by Eq. (S1.7). The units in Eq. (S1.6) imply that the "force"F 0 has units pN/µm.

S1.2. Stokes' paradox and force imbalances in 2D
A complication arises when modeling objects that have nonzero net forcing in 2D. Here we wish to consider an extracellular matrix (ECM) composed of fibers that are tethered together by springs. In that case, it is virtually impossible for the net force on the ECM to be zero during a dynamic simulation, where by definition the fibers are moving and the resulting forces are changing. A net force on the domain inevitably results from our model of the ECM. The fundamental solutions of Stokes flow in free space (regularized in Eq. (S1.6)) contain a logarithmic term so that net force imbalances result in velocities as ||x|| → ∞. Thus the "usual" free space boundary condition of zero velocity at ∞ is incompatible with the net forcing, and a new boundary condition must be imposed [4]. One option is to specify a rigid ECM, which leads to a no-flow boundary condition on the surface of the ECM fibers. Here we study the effect of ECM elasticity and want to account for fiber flexibility. We therefore follow our previous work [5] and impose the boundary condition of zero mean flow on a large circle of radius R that surrounds the domain. As R → ∞, this boundary condition is equivalent to the usual free space boundary condition. We use the value of R = 1000 (cell diameters) as in [5].
The boundary condition on the large circle of radius R allows for net nonzero forcing within the domain and can be enforced by adding a constant velocity throughout the domain of flow. In [5], we derived this velocity for a single point force as If there are N points where force is applied, each with forceF k , the constant velocity added throughout the domain is This ensures that the additional boundary condition is satisfied and yields a well-posed problem. The full numerical scheme is therefore as follows: given a collection of N points with applied forcesF k , Eq. (S1.7) is used to compute an intermediate velocity at each of the N points. Then the boundary condition is enforced through the addition of a constant velocity in Eq. (S1.9), so that the total fluid velocity at each point is the sum of the two. Explicit first-order time-stepping is then used to update the positions of the points.
One of the numerical parameters in the model is , the width of the blob function in hydrodynamic Eq. (S1.4). As discussed in [3] and [6], should be on the order of the spacing of the grid in the system. We set = 0.075 cell diameters.

S2. Elastic forces
In this section, we describe how forces due to stretching and bending on the cortex and nucleus are implemented in our numerical algorithm.
Defining notation first, we denote by X m the mth point on a fiber X; i.e. X m = X(s m ), where s is a parameterization in reference arclength coordinates and m = 1, . . . N . Any function φ is defined at the mth point by φ m . The reference point corresponding to X m is defined as X 0 m . We begin by defining a central difference operator to approximate the (reference arclength) derivatives along the periodic fiber [7] All indices are interpreted in the periodic sense, so that if m = N , the number of fiber points, the index m + 1 maps to index 1. In the case of an equispaced reference configuration, the denominator of Eq. (S2.1) becomes ∆s = 2πr/N , where r is the radius of the contour in its circular configuration. For the case of elastic energy due to stretching, we define the tangent vectors as and the tension is Notice that k m+1/2 and r m+1/2 are the tension constant and rest length at s m+1/2 . In mechanism 2, these change as a function of the position on the fiber. The unit tangent vector τ is dimensionless and k has units of pN/µm. The elastic force density at point m in pN/µm 2 is then defined as [7] F el m = To obtain a force density in pN/µm from force density in pN/µm 2 , we multiply by X 0 The calculation for equispaced reference nodes is the same, but the denominator of τ m+1/2 can be simplified to ∆s in that case. We model bending energy with a preferred curvature. Let the cortex/nucleus radius in the reference state be given by r. We model the preferred curvature κ = 1/r. The bending energy (units of pN·µm) is then defined as [7], In 2D simulations, we define the bending constant to have dimensions pN·µm. The notation K b,m indicates the bending constant depends on the node index. For mechanism 2, the bending constant can be zero in regions where the cortex is loose and relaxed, and its value can be nonzero where the cortex is stiff. Carrying out the numerical differentiation in Eq. (S2.1) twice, we have Differentiating the (negative of the) energy with respect to X , we have the force density (in pN/µm), Unlike for elastic forces, the denominator X 0 m+1/2 − X 0 m−1/2 , which is unknown, no longer vanishes. We therefore approximate so that the force in pN/µm at node is given by For the case of equispaced nodes in the reference configuration, Eq. (S2.10) simplifies to (S2.11) which is the form of [7, Eq. (17)], modified to include a preferred curvature.

S3. Remeshing algorithm
Here we test our remeshing algorithm outside of the cell motility simulations. We consider a closed elastic fiber (representing the cell cortex or nuclear envelope) of radius 1 discretized with N = 50 points, sitting by itself in a domain of fluid of viscosity µ = 1 Pa·s inside and outside. The back half of the cell (points N/4 to 3N/4 ) is stiff with k = 100 pN/µm and resting spring length r = 0.5, while the front of the cell (blue circles) is loose, with k = 1 pN/µm and r = 1 (these constants refer to the tension equation, (4) in the main text). These parameters are similar as those used in the mechanism 2 simulations, with the exception that we choose here to keep r = 0.5 so that the dynamics of the back of the cell are similar to those shown in mechanism 2 simulations. Elastic forces using a remeshed configuration are computed as described in Section S2. The key ingredient is to resample the reference configuration of the nodes, so that deviations from the original circular configuration (with all nodes equispaced) are penalized by the elastic forcing.
We simulate the dynamics of the circle within the method of regularized Stokeslets with equal to the initial point spacing. We begin by simulating the cell to steady state (t = 7 with ∆t = 0.01) without remeshing. Fig. S1(a) shows the deformed shape at this point. The original back half of the cell (points N/4 to 3N/4 marked with red diamonds) has collapsed down, leaving the front-half potentially unresolved. Fig. S1(b) shows the corresponding arclength function s(j) at this point in the algorithm as a blue line. From indices 13 to 37 ( N/4 to 3N/4 ), the slope of the arclength function changes drastically because the points are more tightly packed. In Fig. S1(b), the new sample locations are shown as red circles. These locations are more sparse in the region of decreased slope. This indicates that the arclength function s(j) is being sampled at equally spaced node indices, S (j eq ). Indeed, as shown in Fig. S1(c), the points X(S (j eq )) are equally spaced on the cell after they are remeshed (we always keep the 2 points bound to the ECM nodes). Finally, Fig. S1(d) shows the new reference configurations, X 0 (S (j eq )), which become more sparse at the back of the cell.
The remesher always includes the two points bound to ECM nodes. Because of this, an ideal test for the remesher is to consider the displacement of one of these points with and without remeshing. That is, we vary the number N of points on the cortex and simulate to t = 7 (empirically determined to be steady state) with ∆t = 0.01 as before. Next, we simulate to t = 8 with and without remeshing (only once at t = 7) and determine how much the first red point in the second quadrant (moving counterclockwise) has moved relative to not remeshing (i.e., the error is the norm of the position of the point at t = 8 with remeshing minus the position without remeshing). Table S1 shows the results of this test. It appears that our remesher makes an error on the order O(10 −4 ) for ≈ 100 points, which is what we use in our simulations. This error could be decreased by increasing the number of points. The procedure for remeshing the nucleus is exactly the same as that described for the cortex, with the exception that there are no "bound points" that must be preserved by the remesher. Thus the sampling of the arclength function is truly uniform in this case.    Figure S4: T > E > N : The cell migrates by squeezing the nucleus, which is relatively more deformable than the cortex and ECM.   Figure S6: N > T > E: The cell moves by a combination of nuclear and ECM deformation. However, in this parameter regime, the cell must deform the ECM to migrate because the nucleus is relatively stiffer (compare to the data in Fig. S7). Black ×'s are used to show the initial position of the key ECM nodes, which achieve their maximum displacement at t = 4.53.   Figure S8: N > T > E for a small nucleus: We use the same parameters as in Fig. S6, but halve the diameter of the nucleus. The behavior is qualitatively more similar to the T > N > E regime shown in Fig. S5. Thus halving the nuclear size simply decreases the characteristic force N required to deform the nucleus.

Quantifying migration distance
We quantify the distance moved by the nucleus and the amount that the nucleus and cortex penetrate the ECM. The nuclear displacement is measured by tracking the mean of the coordinates on the nuclear contour (the nucleus is centered at the origin at the start of the simulation). In Figs. S4−S9, the cell begins migrating by moving through a gap between two ECM nodes. We quantify penetration into the network by measuring the fraction of nucleus and cortex points that pass through this gap (past the two nodes). Table S2 shows these statistics at the end of the first motility cycle.
The data in Table S2 show that halving the diameter size of the nucleus with the N > T > E parameters, i.e., keeping all other parameters fixed atk c,r = 50,k n = 100, results in a larger distance moved and more penetration into the network. While the nuclear penetration is larger (since the nucleus is smaller), the total distance moved by the nucleus is similar to the T > N > E regime with the original larger nucleus. We conclude that halving the nuclear diameter for mechanism 1 is equivalent to reducing the force required to move the nucleus.  Penetration is calculated using the fraction of points on the nucleus and cortex that move past the line dividing the two ECM nodes that the cell moves past, approximately located at the points (0.5, ±0.5). Data indicated by * are simulated with a small nucleus (rnuc = 0.225).

Bending rigidity
We determine the effect of bending rigidity by performing simulations using the same parameters as in Table S2, but including a finite bending rigidity of K b = 0.05 (and preferred curvature of a circle) for both the nucleus and cortex. For the cortex, we introduce the bending energy only after the cell binds to an ECM node. Table S3 shows the percentage change for distance and penetration when bending energy is included. We see that the measurements are generally within ±10% of the data from Table S2, with the exception of E > T > N . Surprisingly, the distance increases when bending rigidity is included for this regime. Simulation results in Fig. S9 show the cortex becomes rounder at the front and back. Increased migration distance results when the round cortex at the leading edge fills out the empty space and pulls the nucleus forward. Rounding at the cell rear also pushes the nucleus through the gap between the ECM nodes.

Parameters
Distance Nucleus Penetration Cortex penetration Table S3: Distance and penetration data for mechanism 1 with bending energy included on both the nucleus and cortex (K b = 0.05, preferred curvature of a circle). Bending energy on the cortex is only included after the cell binds to an ECM node, while it is always included on the nucleus. Percentages are the percentage change from the data in Table S2.  Figure S10: E > T and N > T : When the contraction of the cell cortex is too weak to deform either the nucleus or ECM, the cell is not able to squeeze through a gap in the ECM.  Figure S12: T > N > E: The nucleus begins to buckle (bottom row) because of high cortical tension and increased cortical stiffness. After releasing the nodes, the nucleus partially regains its rounded shape (compare the position of the nucleus at t = 5.18 to t = 4.14).  Figure S13: N > T > E: In spite of the nucleus' relatively large stiffness, the cell is still able to migrate through the ECM. However, the cell is unable to entirely pass through the gap after one motility cycle.  Figure S14: E > T > N : Although the cell makes some progress, there is not enough force to drive it completely through the gap. After one motility cycle, the cell retracts slightly after releasing the nodes.  Figure S15: N > T > E for a small nucleus: We use the same parameters as in Fig. S13, but halve the diameter of the nucleus. The behavior is qualitatively more similar to the T > E > N regime shown in Fig. S11.

Quantifying migration distance
As in the previous section, we quantify the distance traveled and the penetration of the nucleus and cortex into the ECM network. We measure distance (nuclear displacement) by computing the center of mass of the nuclear contour (since the nucleus is centered at the origin at the start of the simulation). In Figs. S11−S16, the cell begins migrating by binding to two ECM nodes and moving through the gap between them. We quantify penetration into the network by measuring the fraction of nucleus and cortex points that pass through this gap (past the nodes). Table S4 shows these data at the end of the first motility cycle.
The nucleus travels a larger distance when its size is reduced. Data in Table S4 show that halving the nuclear diameter while keeping the parameters from the N > T > E regime in mechanism 2 (k c,r = 100,k n = 1000) results in behavior more characteristic of the T > E > N regime in simulations with the original larger nucleus. In mechanism 2, decreasing the diameter allows the nucleus to move more easily because the fluid from the back of the cell is able to push the nucleus forward.

Parameters
Distance Nucleus penetration Cortex penetration T > E > N  Table S4: Distance traveled after the first cycle in Figs. S11−S15 computed by tracking the nucleus' center of mass. Penetration is calculated using the fraction of points on the nucleus and cortex that move past the line dividing the two ECM nodes approximately located at the points (0.5, ±0.5). Data indicated by * are simulated with a small nucleus (rnuc = 0.225).

Bending rigidity
We analyze the importance of bending rigidity for mechanism 2 by including a finite bending rigidity of K b = 0.05 on both the nucleus and cortex (with preferred curvature of a circle). Bending energy is included in the cortex in the stiff cell rear after it binds to the two nodes. Other simulation parameters are the same as used to generate data in Table S4. Table S5 shows the relative distances travelled with bending stiffness are approximately the same as those in Table S4, with the exception of the regime E > T > N . We note this is the same parameter regime where bending forces cause increased distance traveled in mechanism 1. In contrast to mechanism 1, motility decreases when the nucleus has bending energy in mechanism 2. Data in Fig. S16 show the nucleus is unable to deform enough to squeeze through the gap between the ECM nodes. The qualitative behavior of the cell with bending rigidity is similar to simulations without bending rigidity in the parameter regime E > N > T where migration fails.

Parameters
Distance Nucleus penetration Cortex penetration Table S5: Distance and penetration data for mechanism 2 with bending energy included on both the nucleus and cortex (K b = 0.05, preferred curvature of a circle). Bending energy on the stiff part of cortex is only included after the cell binds to an ECM node, while it is always included on the nucleus. Percentages are the percentage change from the data in Table S4.