3D Computational Modeling of Bleb Initiation Dynamics

Blebbing occurs in cells under high cortical tension when the membrane locally detaches from the actin cortex, resulting in pressure-driven flow of the cytosol and membrane expansion. Some cells use blebs as leading edge protrusions during cell migration, particularly in 3D environments such as a collagen matrix. Blebs can be initiated through either a localized loss of membrane-cortex adhesion or ablation of the cortex in a region. Bleb morphologies resulting from different initiation mechanisms have not been studied in detail, either experimentally or with theoretical models. Additionally, material properties of the cytoplasm, such as elasticity, have been shown to be important for limiting bleb size. A 3D dynamic computational model of the cell is presented that includes mechanics and the interactions of the cytoplasm, the actin cortex, the cell membrane, and the cytoskeleton. The model is used to quantify bleb expansion dynamics and shapes that result from simulations using different initiation mechanisms. The cytoplasm is modeled as a both viscous fluid and as a poroelastic material. Results from model simulations with a viscous fluid cytoplasm model show much broader blebs that expand faster when they are initiated via cortical ablation than when they are initiated by removing only membrane-cortex adhesion. Simulation results using the poroelastic model of the cytoplasm provide qualitatively similar bleb morphologies regardless of the initiation mechanism. Parameter studies on bleb expansion time, cytoplasmic stiffness, and permeability reveal different scaling properties, namely a smaller power-law exponent, in 3D simulations compared to 2D ones.


INTRODUCTION
Blebs are round membrane protrusions that are used in important cellular processes, such as migration [1] and cytokinesis [2,3]. The hallmark of a bleb is a separation of the cell membrane from the cortex, a thin layer of actin cytoskeleton that is normally attached to the cell membrane by linker proteins. Cells that bleb typically have high intracellular pressure compared to outside the cell. The source of this pressure is attributed to cortical tension due to the myosin molecular motors that slide actin filaments with respect to each other [4]. A bleb is initiated in a localized region by either a loss of membrane-cortex adhesion or by a defect in the actin cortex caused by laser ablation [5] or myosindriven contractility [1]. After a bleb is nucleated, contractile stresses within the cortex are no longer transmitted to the membrane in this localized region, resulting in a pressure gradient and cytoplasmic flow that expands the membrane to create a bleb. A bleb is fully expanded after about 30 s [1], but the timescale is considerably shorter for cells such as Dictyostelium discoideum, where bleb expansion can occur in as little as 0.2 s [6]. Cortical components such as actin and myosin then diffuse into the bleb, the cortex reforms under the naked membrane, and the bleb retracts. During cell migration, the bleb may not retract because adhesion to the substrate (or ECM) may stabilize the bleb.
Theoretical modeling of blebbing has addressed various aspects of the process including expansion [7][8][9][10], retraction [11,12], and migration [13][14][15]. One class of models is based on fluid-structure interaction, where the elastic membrane and cortex are immersed in fluid. Deformations of the elastic structures affect the fluid flow, and the fluid flow in turn influences the motion of the structures. First we summarize results from this approach using two-dimensional models. Using the framework of the immersed boundary (IB) method [16], simulation results from one model showed that cytoskeletal (cortical) drag, and not the viscosity of the cytoplasm, determined the timescale of bleb expansion [7]. When the cytoplasm was modeled as a poroelastic material, two bleb experiments from [5] were simulated. Results from these experiments found that the second bleb was approximately 30% smaller than the first bleb regardless of the location of the second bleb with respect to the first one. The authors found these experimental results could be explained by an initial fast time scale for pressure propagation across the cell combined with a slow timescale for pressure equilibration [9]. The models from [13,17] follow a similar approach from [7,9], but used a boundary integral method instead of the immersed boundary method. The authors simulated cell migration using blebs in confined and unconfined environments (swimming). A 2D agent-based model from [14,18] was used to simulate blebbingbased migration. This model focused on interactions of the cell with different environments, and intracellular pressure was constant in simulations.
Several models consider membrane dynamics in blebbing cells [19][20][21]. These models assume a constant or specified intracellular pressure in the membrane energy. The 1D membrane model from [19] was used to find minimum conditions on pressure and membrane length for bleb nucleation, while the model from [20] found that membranecortex adhesion was critical in determining bleb initiation. The model from [19] was also used to quantify conditions for "circus" blebs that travel around the periphery of the cell. The 2D membrane model from [21] found conditions for circus bleb velocity in terms of biophysical parameters and was used to hypothesize that heterogeneity within the cell surface was necessary to maintain compact circus blebs.
Blebbing dynamics have been modeled in 3D for some simplified cases. The models in [8,12,15,22] assume the cell is axisymmetric and intracellular pressure is constant. In the continuum mechanics model from [8], the membrane reference configuration was dynamically updated to maintain small increases in the area of the cell membrane. A similar approach that involved updating the membrane and cortical reference configuration to model bleb retraction was presented in [12]. In order to produce "smallnecked" blebs observed in experiments, the authors in [22] needed to include either localized membrane growth or global cortical contraction. In a different approach, molecular dynamics simulations of a surface particle-based model were used to simulate bleb expansion in [23]. The authors determined that bleb formation is energetically favorable when the membrane area is larger than its attached cortical area. A membrane model of bleb expansion was simulated with a boundary integral method in [10]. Simulations showed conditions where either no bleb was nucleated, a bleb was nucleated without membrane peeling, and a bleb was nucleated with additional membrane peeling. Although this model included membrane-cortex adhesion dynamics, the membrane was assumed to be axisymmetric and the cortex was fixed in space. A model of bleb expansion that included reaction-diffusion of membrane-cortex adhesion proteins with limited membrane deformation was presented in [24]. Lastly, a full 3D model of bleb expansion with a viscous fluid cytoplasm (an extension of the model in [7]) was presented in [25] to illustrate a numerical method for computing forces on an elastic shell using the immersed boundary method [26].
All of the aforementioned studies highlight the importance of intracellular pressure and mechanics of different cellular components in blebbing dynamics. I consider blebs similar to those in experiments from [5], namely a single cell with one bleb without interactions to an extracellular matrix. These blebs expand on the order of 10 s. The contribution of this paper is a 3D dynamic model of bleb expansion that includes the mechanics of the membrane, cortex, and cytoplasm. In particular, the model presented here includes a dynamic cortex, different cytoplasmic models (viscous fluid and poroelastic material), and makes no a priori axisymmetric assumptions. I then use the model to simulate different bleb initiation using two different mechanisms: cortical ablation and loss of membrane-cortex adhesion. The 3D model with a poroelastic cytoplasm is then simulated to determine whether pressure dynamics follow the same behavior and scaling as described models in 1D [27] and 2D [9].

Overview of Model
The mathematical model of the cell consists of the membrane, cortex, membrane-cortex adhesion, and cytoplasm (see Figure 1). In particular, I consider two models of the cytoplasm: a viscous fluid model and a poroelastic model. In the poroelastic model, the cytoplasm includes a permeable elastic cytoskeleton. Detailed descriptions of the 2D models can be found in [7,9]. The models are summarized here with detailed descriptions of notable differences between the 2D and 3D models. The cell membrane is modeled as an impermeable elastic structure that moves with the velocity of the cytosol (fluid portion of the cytoplasm) while the cortex and the cytoskeleton (in the case of the poroelastic model) are modeled as permeable elastic materials. The model equations consist of force balances on the liquid cytosol, cell cortex, and cytoskeleton, constitutive equations, and equations of motion for the structures. The membrane and cortex are represented by continuous 2D infinitely thin shells immersed in a 3D fluid domain. The cytoskeleton is represented by a 3D structure immersed in the fluid.
The immersed boundary (IB) method is used to account for the interactions between the components of the cell and cytosol/ cytoplasm [16]. In the IB method, structures such as the membrane are represented in a moving, Lagrangian coordinate system, while fluid variables such as cytosolic velocity and pressure are located on a fixed, Eulerian coordinate system. A surface force density on an immersed structure is communicated to the fluid coordinates as follows, where s ∈ Γ is the material coordinate and X (s, t) denotes the physical position of material point s at time t. The interpolation operator is given by where Ω represents the fluid domain. In this paper, capitalized letters represent Lagrangian variables and lower case letters indicate Eulerian variables.

Viscous Fluid Cytoplasm Model
The model of bleb expansion in [7] consisted of the cell membrane modeled as an impermeable elastic stucture, the actin cortex treated as a one-dimensional poroelastic structure attached to the membrane, and the cytoplasm modeled as a viscous fluid. A 3D extension of this model was published in [25] to demonstrate different methods for computing forces on a deforming surface. In this paper, I use the approach from [28] to compute forces due to elasticity on the cortex instead of the method in [25]. The fluid equation includes terms from the elastic membrane, membrane-cortex adhesion, and drag with cortex: where u represents fluid velocity, p pressure, and f i denotes force density on the fluid grid. Force densities due to membrane elasticity, membrane-cortex adhesion, and cortical drag are computed on their respective Lagrangian structures, then are spread onto the fluid grid using Eq. 1.
The drag force on the cortex is balanced by elastic forces within the cortex and adhesion to the membrane: The drag force density from the cortex moving through the fluid is explicitly given by and the drag force density on the fluid is related to the cortex drag force density by The membrane and cortex are each modeled as a hyperelastic shell that experiences forces due to surface tension and stretching. This model is consistent with [29], where the surface area of membrane on a bleb was shown to increase during expansion. A hyperelastic material is characterized by an energy functional E Γ0 W dq, where W is a given strain energy density and Γ 0 represents the reference configuration of the surface. In our model, the total membrane energy density is EdE NH + E ST , where E NH represents the neo-Hookean strain energy density and E ST represents the surface tension strain energy density. Bending forces are neglected because several studies have shown that they do not contribute significantly to blebbing mechanics [10,30]. Following the formulation from [28], the neo-Hookean surface strain energy density is given by where k E is the bulk modulus, μ E is the shear modulus and s (s 1 , s 2 ) represents the surface material coordinates. The following tensors G and G 0 are defined as where X(s) denotes the current surface position, Z(s) denotes the reference membrane position, and the scalar J (det(GG −1 0 )) −1 . Energy due to surface tension is where c i represents the parameter for surface tension on the membrane or cortex. The force density per unit reference configuration is then computed by where i indicates either the membrane or the cortex. Membrane-cortex adhesion is modeled by elastic springs attaching the membrane to the cortex with force density Frontiers in Physics | www.frontiersin.org November 2021 | Volume 9 | Article 775465 The membrane-cortex adhesion stiffness coefficient k mem/cortex adh was chosen so that the cortex and membrane are within a computational grid cube of each other and velocity of the cortex is close to zero if no bleb is initiated. The adhesion force density on the membrane is the opposite of the corresponding force density on the cortex. The adhesion force densities satisfy the equation Given the stiffness coefficient k mem/cortex adh , the corresponding stiffness coefficient for the cortex is obtained by k cortex/mem adh k mem/cortex adh dA mem j /dA cortex j , where dA i j represents the surface area differential of the jth triangle of the membrane or cortex in reference coordinates.
Given a configuration of the membrane and cortex, the forces are computed, and the velocities of the fluid and cortex are obtained by solving Eq. 3 and Eq. 5 as described in Section 2.2. Then, the positions of the membrane and cortex are updated with their respective velocities.
dX cortex dt Table 1 lists the default parameters for blebbing simulations. The values for membrane and cortex surface tension were decreased by 50% compared to the values in [9] so that the initial pressure difference across the membrane matches the difference in the 2D model. The bulk modulus of the cortex was taken to be several orders of magnitude higher than the membrane to reflect stiffness due to the cross-linked cortical actin network.

Poroelastic Cytoplasm Model
The model formulation is the same as above with the addition of a poroelastic cytoplasm throughout the cell interior (a 2D version of the model is described in [9,32]). The cytoskeleton is represented by a porous elastic network in Figure 1). The fluid equations have an additional term for cytoskeletal drag.
The force density balance on the cortex includes an additional adhesion term to link the cortex to the cytoskeleton, Similarly, the force density balance on the cytoskeleton is where cytoskeletal drag is defined as where κ is the permeability of the cytoskeleton. In this formulation of poroelasticity, the volume fraction of the network (cytoskeleton) is negligible [32]. The cytoskeleton is modeled as a porous, neo-Hookean elastic structure. Elastic forces are computed using the energy functional-based version of the IB method proposed in [33]. The neo-Hookean strain energy is the same as given in Eq. 8, but the 2 × 2 tensors G and G 0 now have dimensions 3 × 3 for the solid cytoskeleton. The equivalence of the solid and surface neo-Hookean energy formulas is described in [28].
The force density for cortex-cytoskeleton adhesion is computed similarly to cortex-membrane adhesion in Eq. 12 with the appropriate scaling of the force densities. Given the stiffness coefficient for cortex-cytoskeleton adhesion k cortex/cyto adh represents the reference volume differential of the cytoskeleton at the ith point of a tetrahedron.
The structures are updated with their respective velocities. Table 1 lists the values of parameters for blebbing simulations. Most simulations in the Results section use a value of 20 pN/μm for membrane stiffness, 500 Pa for k cyto , and 1 · 10 -2 µm 2 for permeabilty.  35] μm and spacing h 30/32. The cytoskeleton is represented by an adaptive unstructured tetrahedral mesh with radius 9.99 μm. The mesh is more refined near the cortex with approximately two Lagrangian points per Eulerian grid cell and one Lagrangian point per Eulerian grid cell in the interior. The mesh consists of 28,153 points and 153,202 tetrahedra. The cortex is taken to be the boundary of the cytoskeleton, and membrane points are initialized to be the same as cortical points, except adjusted to have a radius of 10 μm. The membrane and cortex each consisted of 5,018 points and 10,032 triangles. The unstructured meshes were generated using distmesh [34].
To compute elastic forces on the membrane, cortex, and cytoskeleton, methods from [28,33] are implemented. Elastic forces are computed directly from an energy functional without the use of stress tensors by taking the variational derivative of the energy (see Eq. 11). I assume the deformation map X (s, t) is a piecewise linear function on each triangle of a discretized surface (tetrahedron in the cytoskeleton) so that the deformation gradient tensor and strain energy are constant on each triangle (or tetrahedron). This simplifies the integral of the strain energy density in Eq. 8 and Eq. 10. The variational derivative in Eq. 11 can then be computed analytically on each discretized element, and the force at each vertex i is the sum of all elements that contain i. Details for computing forces due to shell elasticity can be found in [28], and details for computing the forces due to the elasticity of a solid are located in [32,33].
The time update follows [9,32]. Given the current position of the structures (membrane, cortex, and cytoskeleton): 1. Compute elastic forces based on the current membrane, cortex, and cytoskeleton configuration (X n i X i (s, t n ), where i denotes the structure: membrane, cortex, or cytoskeleton) using the constitutive laws described in Section 2.1. 2. Spread the force densities onto nearby Eulerian points using Eq. 1. 3. Solve the forced Stokes equations to obtain the fluid velocity u. 4. Interpolate the fluid velocity to the structure using Eq. 2 to obtain U. 5. Compute the porous structure velocities by Eq. 15 and Eq. 23, and update the structure by where ζ i indicates the drag coefficient of the cortex or cytoskeleton, and the forces acting on the respective structure denoted by F j i are in Eq. 15 and Eq. 23. The membrane is updated by the fluid velocity.
The time step for simulations using a pure fluid cytoplasm is Δt 1 · 10 -4 (seconds). Cytoplasmic elasticity introduces significant stiffness in the model, and the time step is reduced to Δt 7.5 · 10 -6 − 3 · 10 -5 . Small time steps are required for numerical stability when cytoplasmic permeability and/or elastic moduli are relatively large.
Velocity and pressure satisfy periodic boundary conditions on the Eulerian (fluid) domain. A Fourier-spectral method is used to solve the Stokes equations (Eq. 3 and Eq. 4) for the viscous fluid model, Eq. 16 and Eq. 17 for the poroelastic fluid model).
The IB method involves approximate δ functions for the spreading and interpolating operators in Eq. 1 and Eq. 2 . Here, I use spectral delta functions as described in the Supporting Material of [9] to avoid unphysical velocities that occur in regions with a large pressure jump, such as across the cell membrane. The discretized spreading and interpolation integrals are approximated in Fourier space by a nonuniform fast Fourier transform (NUFFT) described in [35]. The result yields obtain a computationally efficient approximation of the Fourier transform of the spread force density. I use an mth order cardinal B-spline with compact support over m + 1 mesh points as a smoothing kernel [36]. In this paper, the oversampling factor is r 1.5 and m 6. The Fourier transform of the spread forces are filtered with a second-order raised cosine filter [37], to remove the Gibbs phenomenon in the numerical solution to the pressure field.

RESULTS
At the beginning of a simulation, the cell is pressurized due to membrane and cortical tension. According to Laplace's law, the change in pressure ΔP across a spherical membrane satisfies ΔP c/(2R), where c represents surface tension and R is the radius of the cell. Here, ΔP ≈ (c mem + c cortex )/(2 r mem ) 44 Pa. A bleb is initiated by either removing membrane-cortex adhesion in a circular region at the top of the cell or by cortical ablation (see Figure 2). The initial radius of the circle is approximately 2.5 μm and centered at the point on the membrane with the highest z-coordinate. Numerically, k mem/cortex adh k cortex/mem adh 0 on surface triangles inside the region. The boundary between adjacent triangles on the surface where the membrane-cortex adhesion parameters change from nonzero to zero forms a round circular ring around the cell, referred to as the bleb ring. Figure 2C) shows a top down view of the cell with a dashed line indicating the bleb ring. When the cortex is ablated, the additional parameters k cortex 0, k cortex/cyto adh , and k cyto/cortex adh are set to zero on triangles within the bleb ring.
For both initiation mechanisms, when cortical tension is no longer transmitted to the membrane, the resulting pressure gradient leads to fluid flow and membrane expansion in the localized region. Bleb expansion ceases when membrane tension balances intracellular pressure.
Parameter values for simulations are located in Table 1. cytoplasm use a value of 500 Pa for the bulk modulus of the cytoplasm and 1 · 10 -3 for the cytoskeletal permeability unless otherwise stated.

Bleb Initiation Dynamics
I begin by simulating bleb expansion by removing only membrane-cortex adhesion within the bleb ring region.   Frontiers in Physics | www.frontiersin.org November 2021 | Volume 9 | Article 775465 Figure 3 shows the cell membrane position and speed (magnitude of the velocity on the cell membrane) at several time values during bleb expansion. The speed is highest after removing membrane-cortex adhesion and decreases over time. Similar to 2D simulations from [9], the bleb appears smaller in simulations using a poroelastic cytoplasm model compared to the pure fluid model. Intracellular pressure dynamics are shown in Figure 4 for bleb simulations using fluid and poroelastic cytoplasm models. Intracellular pressure appears spatially uniform inside the cell body (outside of the bleb) for the simulation with a fluid cytoplasm (Figures 4A,C). The maximum pressure inside the cell decreases by about 4 Pa over 40 s of simulation time. For the poroelastic cytoplasm, a pressure gradient extends across the entire cell (Figures 4B,D). Pressure also equilibrates to a smaller value: 34 Pa for the poroelastic cytoplasm model compared to 41 Pa for the fluid cytoplasm model. Maximum pressure inside the cell decreases by approximately 10 Pa over 40 s of simulation time. The decrease in intracellular pressure results in decreased fluid speed and bleb size.
Decreased intracellular pressure in the poroelastic cytoplasm is a result of compression of the cytoskeleton as the bleb expands. The total volume of the cell is conserved so that as the bleb expands, the main cell body is compressed. In [9], the cytoskeletal pressure that acts against compressive stresses was given as − k cyto (J − 1), where J − 1 is the strain in Eq. 11. Cytoskeletal pressure at several time values is shown in Supplementary Figure S1 when k cyto 500 Pa. Data show the cytoskeletal network is compressed near the nucleation site, and cytoskeletal pressure approaches a spatially nonuniform profile. Cytoskeletal pressure over time when k cyto 250, 500, and 1,000 Pa at the center of the cell is shown in Supplementary Figure S2. The pressure approaches a steady state value and increases with k cyto .
Bleb size is quantified two ways, depending on how the bleb is initiated. When a bleb is initiated through a loss of membranecortex adhesion, I measure the relative bleb volume over time. The volume of the cortex is subtracted from the volume of the membrane, then divided by the initial volume of the membrane. Since there is an initial small volume between the membrane and cortex, this small initial volume is subtracted from the relative bleb volume equation, given by Note that the volume of the membrane is approximately the same value over time due to incompressibility of the fluid. When a bleb is initiated through cortical ablation, the cortex is no longer a closed surface. I use bleb height as a measurement of bleb size and is defined as follows. The point on the membrane with the largest z-coordinate (see Figure 3A) and the point with the smallest z coordinate are identified. Initially, the difference between these z coordinates is the diameter of the cell, 20 μm. The difference between these z coordinates is measured over time. The initial distance is then subtracted from the difference between the z values, Bleb height, relative bleb volume, and maximum intracellular pressure over time are shown in Figure 5 for blebbing simulations using the viscous fluid and poroelastic cytoplasm models. Bleb height and volume initially increase after bleb nucleation, then approach a steady state value. In this paper, steady state is defined as a quantity having less than 5% relative change over 10 s of simulation time. Following [9], bleb expansion time is defined as 90% of the steady state value. White circles in Figure 5A indicate bleb expansion time using bleb height, and black circles in Figure 5B denote bleb expansion time defined as 90% of steady state relative bleb volume. Bleb height and volume are larger for the viscous fluid cytoplasm model compared to the poroelastic one. Bleb expansion time is also faster in simulations using the viscous fluid cytoplasm model. Maximum intracellular pressure decreases over time as the bleb expands for simulations with fluid and poroelastic cytoplasm models, but significantly more pressure is relieved when bleb expansion is simulated using a poroelastic cytoplasm model (Figures 4, 5C). To determine the time scale of pressure equilibration, the change in pressure from its initial value over 40 s is measured (approximately 4 Pa for a simulation with a fluid cytoplasm and 10 Pa for a simulation with a poroelastic cytoplasm). Pressure equilibration time is computed as the time when the maximum intracellular pressure equals the initial value minus 90% of the change in pressure over 40 s. Figure 5C) shows maximum intracellular pressure evaluated at pressure equilibration time, bleb expansion time using bleb height, and bleb expansion time using relative bleb volume. The data show close agreement between pressure equilibration time and bleb expansion time using relative bleb volume. Although bleb expansion time using bleb height is about 5 s faster than pressure equilibration time, 80% of the change in pressure is achieved by this time for simulations with both the fluid cytoplasm and poroelastic cytoplasm models. Therefore, bleb expansion time measured by using height or relative volume approximately correspond to pressure equilibration time.

Cortical Ablation
Bleb initiation by cortical ablation is simulated for several values of cortical elastic modulus, k cortex 100, 500, and 1,000 with both the fluid and poroelastic cytoplasm models. Figure 6B shows steady state membrane shape and mean curvature for simulations with cortical ablation. Membrane shape and curvature for bleb expansion by a loss of membrane-cortex adhesion with cortical elastic modulus k cortex 1,000 is shown in Figure 6A for comparison. The magnitude of mean curvature is computed along edges of triangles and assigned to vertices as described in [28,38]. Following [25], the convex hull of the membrane is computed to determine the sign of mean curvature; points on the convex part of the surface are assigned negative mean curvature values, and other points on the membrane maintain their positive sign. Figure 6 I show that steady state bleb shape is very sensitive to changes in cortical elastic modulus in the fluid cytoplasm model. As the values of cortical elastic modulus decrease, the bleb becomes more broad with a decrease in positive mean curvature near the bleb neck as compared to data when the bleb is initiated by only a loss in membranecortex adhesion ( Figure 6IA). When blebbing is simulated using the poroelastic model, membrane shape and curvature do not change significantly as the cortical elastic modulus decreases. Data in Figure 6IIB show slightly broader blebs with a small   Table 1. The scale bar has dimensions 5 × 1 × 1 μm. The colorbar indicates the value of mean curvature H over the surface of the membrane.

Data in
Frontiers in Physics | www.frontiersin.org November 2021 | Volume 9 | Article 775465 decrease in positive mean curvature near the bleb neck as compared to membrane shape and mean curvature in Figure 6IIA.
To quantity the broadness of the bleb, bleb ring diameter, defined as the distance between the point on the bleb ring closest to origin (in Euclidean distance) and the point furthest away from the origin as illustrated in Figure 2C, is computed. The change in bleb ring diameter for cortical ablation simulations is graphed in Figure 7 for values of cortical elastic modulus, k cortex 100, 500, and 1,000 with both the fluid and poroelastic cytoplasm models. The change in bleb ring diameter is defined as steady state bleb ring diameter minus the steady state bleb ring diameter for a simulation when a bleb is initiated by a loss of membrane-cortex adhesion with cortical elastic modulus k cortex 1,000. Data in Figure 7 show the diameter of the bleb ring increases as the cortical elastic modulus decreases, but the increase is relatively much larger for simulations with the fluid cytoplasm model than for simulations with the poroelastic model.
Bleb expansion time decreases when a bleb is initiated by cortical ablation compared to a loss of membrane-cortex adhesion. Expansion time is calculated by computing the time value when 90% of the steady state value of bleb height is reached. When blebbing is simulated with a fluid cytoplasm and a loss in membrane-cortex adhesion ( Figure 6IA), bleb expansion time is 1.95 s. In comparison, bleb expansion time from cortical ablation simulations ( Figure 6IB) is 0.016, 0.023, and 0.028 s for k cortex 100, 500, and 1,000 pN/μm, respectively. Bleb expansion time for the simulation with a poroelastic cytoplasm model initiated with loss of membrane-cortex adhesion ( Figure 6IIA) is 10 s. When a bleb is initiated with cortical ablation, bleb expansion time decreases to 7.37, 6.96, and 6.60 s for k cortex 100, 500, and 1,000 pN/μm, respectively. Data from bleb simulations with a fluid cytoplasm show expansion time decreases by two orders of magnitude compared to simulations when blebbing is initiated by a loss in membrane-cortex adhesion. Simulations from the poroelastic cytoplasm model with cortical ablation show a decrease in bleb expansion time compared to initiation by a loss in membrane-cortex adhesion, but the relative decrease is approximately 30%.

Effect of Poroelasticity on Bleb Expansion Time
This section focuses on bleb expansion dynamics when the cytoplasm is modeled by poroelastic material, and a bleb is initiated by a loss in membrane-cortex adhesion. Several studies have used a poroelastic model of the cytoplasm to interpret experimental results. For example, when blebbing was locally inhibited by myosin-II-inhibiting drugs, blebbing was unaffected in other parts of the cell, suggesting that intracellular pressure is not equilibrated on the time scale of bleb expansion [27]. Using a simple 1D model, the authors in [27] showed that pressure diffuses over a length x ∼ Dt √ , where x is a characteristic length (10 μm, the radius of the cell), D is a diffusion coefficient proportional to both cytoplasmic permeability and stiffness. This scaling law can be expressed equivalently as bleb expansion time t ∼ x 2 /D. Since D is proportional to κ G, bleb expansion time is expected to scale like κ −1 . Simulations from the 2D version of the blebbing model with a poroelastic cytoplasm presented are in approximate agreement with the scaling t ∼ κ −1 [9]. Figure 8 shows bleb expansion time as a function of cytoplasmic permeability κ and bulk elastic modulus G. Expansion time is computed by both relative bleb volume using Eq. 26 ( Figure 8A) and bleb height Eq. 27 ( Figure 8B). The data show qualitative agreement with results from the 2D model in [9]. Expansion time decreases as cytoplasmic permeability and stiffness increase. Bleb height (and volume) decrease as cytoplasmic stiffness increase. Note that bleb height and volume each maintain the same value for different values of permeability, but change as a function of elastic modulus. One difference from 2D model results is that data in Figure 8 show bleb expansion time scales slower than κ −1 . When relative bleb volume is the metric, expansion time scales like ∼ κ −0.65 . If bleb height is used, the scaling is slightly larger, ∼ κ −0.69 . These results suggest the cytoskeletal network does not obey a simple diffusion equation in 3D as described in [9,27].

DISCUSSION
In this work, fully 3D simulations of bleb expansion with two different rheological descriptions of the cytoplasm, purely viscous fluid and poroelastic material, are presented. The model formulation follows [9], with significant extensions to model the elasticity of the membrane and cytoskeleton. It should be noted that no geometric simplifications, such as axisymmetry, were made. 3D models allow quantitative comparision to experimental data. For example, values for relative bleb volume reported in Figure 5 match those computed from experimental data in [5].
3D simulations are necessary to simulate the dynamics of cortical ablation, rupturing the cortical actin shell, which has been performed in experiments and has been hypothesized to occur due to localized regions of myosin-driven contractility. In a 2D model, the cortex can be thought of as a rubber band that quickly recoils after being cut. In the full 3D model, an additional line tension helps to maintain the structure of the elastic cortical shell. Simulations with a viscous fluid cytoplasm show much broader blebs that expand several orders of magnitude faster than expansion times from experiments. Bleb shape and expansion times in simulations with a poroelastic cytoplasm are qualitatively similar to those when blebbing was initiated with only removing membrane-cortex adhesion. It is energetically unfavorable for the membrane to maintain regions of high curvature, such as those observed in blebs, and these results suggest that the internal cytoskeleton along with attachments from the membrane to the cytoskeleton help maintain the shape of bleb-like protrusions.
In cortical ablation simulations with a viscous fluid cytoplasm, the timescale of bleb expansion is determined primarily by fluid viscosity in the absence of cortical drag [7]. Cytoplasm drag is still present during simulations with a poroelastic cytoplasm when blebs are initiated with cortical ablation. Although bleb expansion is slightly faster compared to the case when blebs are initiated by a loss of membrane-cortex adhesion, cytoplasmic drag dominates the timescale of bleb expansion, and biologically relevant values are obtained from these simulations. The decrease in bleb expansion time can be explained by a loss of cortical drag within the bleb ring region.
Pressure dynamics are quantified for both models of the cytoplasm with qualitatively similar behavior to 2D model simulations results from [9]. Pressure in the main cell body, i.e. the cell except the bleb, is approximately uniform for simulations with a viscous fluid cytoplasm, and diffuses across the cell body in simulations with a poroelastic cytoplasm. Instances of uncontrolled bleb growth for soft membranes were reported in [9] for simulations with a purely viscous cytoplasm. Although I observed large blebs in 3D simulations with a soft membrane (data not shown), I was unable to confirm the growth was actually uncontrolled. As the bleb expands, the spacing between the discretized nodes became too large to obtain reliable results from the simulation. Extending the numerical method to include mesh refinement within the bleb would be necessary to ascertain whether uncontrolled bleb growth can occur in 3D.
In a parameter study where bleb expansion is simulated with different values of permeability and cytoplasmic stiffness, I found that bleb expansion time follows a power law, where t ∼ κ −p , where p 0.65, 0.69, depending on whether bleb volume or height was used to compute expansion time. These results are in contrast to the scaling law calculated from analyzing 2D simulations in [9] and 1D model analysis in [27]. If poroelasticity alone determines the time scale of bleb expansion, it can be shown that cytoskeletal displacement follows a diffusion equation after some simplifying model assumptions, i.e., deformation only in the radial direction and applying a small displacement in the radial direction (a reduced model in polar coordinates is located in [32]). Because the scaling of the diffusion equation remains the same in any dimension, my results point to the influence of other important contributing factors to limiting bleb expansion, such as stiffness of the cortex and membrane, and cortical permeability. The influence of these factors is a subject of future work.
The model presented here is limited in that I consider one cell and focus on the expansion of one bleb to focus on intracellular pressure dynamics with different models of the cytoplasm. Recent work on modeling blebbing-based migration involves simulations with multiple blebs and cycles of bleb expansion and retraction, but assume uniform cytoplasmic pressure [14,15]. I have also neglected to include membrane-cortex adhesion dynamics [10,19,21]. Model extensions are currently underway to include these important effects into the current modeling framework.
Finally, I note that simulations of fully 3D models are computationally expensive; some simulations in this paper took weeks to run using parallelized C++ code on a cluster. Numerical methods for stiff fluid-structure interaction systems will need to be investigated and implemented in order to simulate multiple blebs and cell migration in 3D environments, particularly with an internal elastic cytoskeletal network.