Mechanical plasticity of cell membranes enhances epithelial wound closure

During epithelial wound healing, cell morphology near the healed wound and the healing rate vary strongly among different developmental stages even for a single species like Drosophila . We develop deformable particle (DP) model simulations to understand how variations in cell mechanics give rise to distinct wound closure phenotypes in the Drosophila embryonic ectoderm and larval wing disc epithelium. We ﬁnd that plastic deformation of the cell membrane can generate large changes in cell shape consistent with wound closure in the embryonic ectoderm. Our results show that the embryonic ectoderm is best described by cell membranes with an elasto-plastic response, whereas the larval wing disc is best described by cell membranes with an exclusively elastic response. By varying the mechanical response of cell membranes in DP simulations, we recapitulate the wound closure behavior of both the embryonic ectoderm and the larval wing disc.

The supplemental material includes seven sections that provide additional details concerning the deformable particle simulations and analyses of the experimental data of wound closure in epithelial monolayers.In Section I, we describe the intercellular forces, initialization of the cell monolayer, generation of the wound, equations of motion, and tracking of the wound boundary over time in the numerical simulations of wound closure.In Section II, we provide a novel method for calculating smooth intercellular forces in the deformable particle model simulations.In Section III, we describe the cell segmentation, cell shape parameter calculations, and error estimation from the analyses of confocal microscopy images of wound healing in epithelial tissues in Drosophila [1].In Sections IV and V, we relate the cell bulk modulus to the deformable particle area stiffness, and express the bulk modulus and other simulation quantities in physical units.In Section VI, we study the deformation of a cell experiencing an extensile force dipole to give further evidence that a larger bulk modulus leads to enhanced shape changes in plastic cells.In Section VII, we discuss multiple contributing factors to the overall process of cell shape plasticity.

I. WOUND CLOSURE SIMULATION PROTOCOL
We model cell monolayers as nearly confluent packings of deformable particles.Each deformable particle (cell) i obeys the shape energy in Eq. 1 in the main text and interacts with other cells through vertex-segment forces between vertex α on cell i and segment ⃗ l βj on cell j as shown in Fig. 2 in the main text.The intercellular forces can be derived using the following pair potential: where ϵ is the cell membrane interaction strength, d ≡ d αi,βj is the vertex-segment distance, σ is half the membrane width, l 2 = 0.3σ sets the attractive range, l 1 controls the membrane attraction strength, and l 1 < l 2 .The vertex-segment pair potential u int and the force f int = −∂u int /∂d are plotted against d/σ in Fig. S1.See Table S1 for a list of default parameters and other quantities used in the simulations.
To generate unwounded cell monolayers, we first place cells with random initial positions within a circular boundary at an initial packing fraction ϕ 0 = 0.78, and then compress the system in small packing fraction increments ∆ϕ = 0.005 until ϕ = 0.92.After each compression step, we use the FIRE algorithm [2] to minimize the total potential energy of the cells within the circular boundary: where U rep is equal to U int with l 1 = 0 and Here, r αi is the distance between a vertex α on cell i and the center of the circular boundary of radius R. We set ϵ cb = ϵ as the strength of the interaction energy between the circular boundary and the vertices.The resulting packings of deformable particles with A = 1.2 and ϕ = 0.92 serve as initial conditions for the numerical simulations of wound closure.First, we add cell-cell adhesion to the packings by setting l 1 /σ = 0.1 and carry out constant energy dynamics without the circular boundary for a total time ∼ √ a 0 /ω to allow the cells to explore intercellular gaps.We then minimize the total potential energy U DP (Eq. 4 in the main text) using an overdamped equation of motion, where each vertex α on cell i obeys Here, ⃗ f αi is the total force on the vertex i, and b is the damping coefficient.Overdamped dynamics are frequently employed in active particle models [3] to balance the input of energy due to activity, which corresponds to the purse-string contraction in this model.The effects of the cells' environment on cell dynamics are commonly modeled using a single damping coefficient [3,4], as we have done in Equation S5.
We then introduce a wound by removing 5 central cells in the packing, define the purse-string (PS) on the wound boundary, and integrate Eq.S5 using a modified velocity-Verlet algorithm with timestep dt * = 0.05/ k * a .We record cell properties, such as the cell shape parameter A, and wound area A throughout the course of the simulation.The simulated wound boundary is tracked over time by monitoring vertex-vertex contacts and determining the largest cluster of vertices near the center of the tissue using the Newman-Ziff union-find algorithm [5].We terminate the wound closure simulation when the wound area satisfies A < 10 −2 a 0 .In both the numerical simulations and image analyses, a cell is considered wound-adjacent if it is within a 0 /π of the wound boundary.

II. SMOOTH SLIDING INTERCELLULAR FORCES
Intercellular forces are bumpy in the deformable particle model when the intercellular potential is a function of the vertex-vertex distances.We develop a novel method (within the deformable particle model) to compute smooth (frictionless) sliding intercellular forces.To do this, we assume that the intercellular pair potential U int is only a function of the closest distance d αi,βj between vertex α on cell i and line segment ⃗ l βj on cell j (Eq. 3 in the main text), which ensures that there is no component of the intercellular force tangential to the membrane at ⃗ l βj as shown in Fig. S2a.This force law mimics smooth surfaces that consist of connected rectangles (cyan) and wedges (blue) in Fig. S2b.Because d αi,βj is a function of ⃗ r αi , ⃗ r βj , and ⃗ r (β−1)j , forces computed using the pair potential u int (d αi,βj ) will affect the dynamics of all three vertices α, β, and β − 1.
For interaction potentials that depend on the closest distance between a point and a line segment, one must determine whether there are discontinuities in the force that can occur when the closest distance changes discontinuously even though a vertex or line segment moves by an infinitesimal amount.We consider two cases, concave and convex sections of the cell membrane, classified by the interior angle θ defined by three successive vertices.In the convex case, θ < π (Fig. S2c), a vertex that overlaps with the membrane surface at ⃗ r contact can slide along the surface with a continuous vertex-segment distance d αi,βj to the line segment ⃗ l βj .In the concave case, θ > π (Fig. S2d), a vertex on a similar trajectory will experience a discontinuity in the vertex-segment dis-tance, since the concave surface (Fig. S2b) lacks a wedge where there is an overlap of the two rectangles.One method to remove this discontinuity is to add a wedgeshaped patch, shown as a red grid in Fig. S2e.Within the wedge-shaped region, a vertex α overlapping with the membrane surface l βj at ⃗ r contact incurs an additional force such that in the patch region, f patch provides an equal and opposite force to offset the discontinuity in f int (d αi,βj ) (Fig. S1b) that occurs when the vertex α enters the patch region.

III. IMAGE ANALYSIS
To measure the cell shape parameters and wound area over time of the wounded embryo and wing disc epithelia, we analyze segmented images of the wound closure process from Ref. [1].We use the segmented cell boundaries in the 5 wing disc wounds found in Ref. [1].For the 2 embryo wounds in Ref. [1], we perform our own segmentation procedure by first taking a maximum intensity projection (Fig. S3a) along the z-axis of confocal microscope z-stack time-series images.Next, we use Tissue Analyzer [6], an ImageJ plugin for segmentation of singlelayered epithelia, to obtain a first pass segmentation of the cells and wound (Fig. S3b).Then we make manual corrections to account for under-and over-segmented regions near the wound, which yields a binary image of the cell boundaries (Fig. S3c).We employ regionprops in MATLAB R2022a, Update 5 (Fig. S3d-e) to calculate the area and perimeter of each unique segmented region, which we use to report the wound areas and cell shape parameters A for each time point.
To generate error estimates for the cell shape parameter measurements, we carry out a similar process on synthetic data.We use 26 test shapes: ellipses with eccentricities 0.996, 0.987, 0.968, 0.933, 0.872, 0.768, 0.586, and 0 (circle), polygons with 3 to 12 sides, and 8 different 7-sided concave shapes.A subset of the test shapes is shown in Fig. S4a.We generate images of these shapes with a range of resolutions (Fig. S4b) and compare the measured A to the true value A t for each shape (Fig. S4c).We define the fractional error of A as For a cell with area a, we report an estimate of δA using the bounding fractional error E(a) as shown in Fig. S4c.E(a) is calculated by taking the maximum δA over all synthetic test shapes at each area, and storing the running maximum as a function of decreasing area.For each measurement of the cell shape parameter A for a cell with area a (px 2 ) in the embryo and wing disc wounding experiments, we associate E(a) with the measurement error in A.
We plot A(t) taking the variance-weighted mean over cells adjacent to the wound boundary (Fig. S3d), with error bars given by the standard error of this weighted mean.Onto the variance-weighted mean and error bars, we overlay the simulation results as in Fig. 4b in the main text to show that A(t), from simulations of the embryo and the wing disc using cells with shape memory, falls within the margin of error for the A(t) from the experimental measurements of the embryo and wing disc wounds.

IV. DERIVATION OF BULK MODULUS
The bulk modulus B of a single cell is related to the cell area stiffness k a through Eq. 9 in the main text.To derive this relation, we start with the definition where P is the pressure and V is the cell volume.Rearranging, we obtain and P = B log(V /V 0 ), where V = V 0 at zero pressure.The energy under isothermal compression is given by To second order in V /V 0 − 1, we obtain Comparing Eq.S8 to the shape-energy function (Eq. 1 in the main text), the energy due to compression in two dimensions is given by Assuming that the energy scale of volume changes is equal to that of area changes, and that

V. CONVERTING SIMULATION UNITS TO PHYSICAL UNITS
Simulation units can be converted into physical units using three physical quantities that set the mass, length, and time scales of the simulation.We are able to determine these scales using a choice of force, velocity, and area (See Table S1.).Atomic force microscopy can determine single cell forces [7], which allows us to estimate the unsticking force f adh = 1 nN based on cohesion between Zebrafish embryo ectodermal cells [8].We are unaware of any measurements of f adh on Drosophila embryo ectoderm, and we assume that the measurements of f adh on Zebrafish embryo ectoderm give a reasonable estimate.We choose ω = 0.3 µm/sec based on a typical actin ring constriction rate [9,10].We find that typical Drosophila cell areas are a 0 ∼ 16 µm 2 for the late-stage larval wing disc epithelium and ∼ 25 µm 2 for the late-stage embryo ectoderm.
We define the unsticking force in our simulations by f adh = l 1 ϵN v /3σ 2 with units of ϵ/σ.A factor of N v /6 comes from the assumption that two cells are adhered to each other through 1/6 of their membranes on average, given that a cell has approximately 6 neighbors.The maximum vertex-vertex adhesion force in Eq.S1 is 2l 1 ϵ/σ 2 , again with units of ϵ/σ.

VI. STIFFNESS ENHANCES SHAPE CHANGE IN PLASTIC CELLS
To understand how greater stiffness can lead to enhanced shape change in plastic cells, we conduct simulations varying the plastic relaxation timescale τ and the cell area stiffness k a of a single cell experiencing an extensile force dipole (Fig. S5a), i.e. the cell experiences two equal and opposite forces that generate net zero force.For elastic-like cells with large τ , we find that increasing k a leads to a reduction in the final A (Fig. S5b), which matches the expectation that stiffer cells are less deformable.However, the trend reverses for plastic-like cells with small τ , as increasing k a leads to an increase in the final A. We find that stiff, plastic cells are able to lengthen their membranes without changing their area.In contrast, soft plastic cells increase their area when lengthening their membranes, which results in more modest shape changes.These results show that plastic cells become more deformable as they become stiffer, and corroborates a similar trend in Fig. 3 in the main text.

VII. ALTERNATIVE MECHANISMS INFLUENCING CELL SHAPE PLASTICITY
Cell shape plasticity, the property of cells allowing them to retain their new shapes after deformation, is the result of several mechanisms, which include actin cortex remodeling, caveolae acting as membrane reservoirs, and vesicle trafficking processes like exocytosis and endocytosis.These processes are involved in mechanoprotection, and allow the cell to modulate the cell membrane surface area in response to stress.In our model, a natural way to incorporate changes in cell membrane surface area due to membrane reservoirs and vesicle trafficking is to add plasticity in the membrane rest length.We describe cell shape plasticity as membrane plasticity using Equation 5. A different model for cell shape plasticity could add plasticity in the equilibrium bending angle of each membrane segment, which would account for how actin cortex remodeling contributes to cell shape plasticity independently of membrane surface area relaxation processes.Since curvature is dependent on both membrane segment lengths and the angles between the segments, we note that membrane rest length plasticity also includes relaxation of membrane curvature.S1.Default parameters and other quantities used in the deformable particle simulations of wound closure.An asterisk denotes non-dimensional simulation units.Bolded quantities have physical units and are based on experimental measurements.If only one parameter value is listed, then the parameter is not varied in our simulation studies.In the right image, we also show the exterior-facing half of the cell membrane of a "smooth" deformable particle, which consists of circulo-line segments (cyan) and vertex sectors (blue).Within a smooth deformable particle, the section around each vertex is considered concave when the angle the vertex makes with its two neighboring vertices has an interior angle θ > π, and convex otherwise.(c) For locally convex geometries, when a vertex that overlaps the membrane surface at ⃗ rcontact moves along the indicated trajectory (red line), the vertex-segment distance relative to ⃗ l βj is continuous.(d) For locally concave geometries, the vertex-segment distance (relative to ⃗ l βj ) for a vertex overlapping with the membrane along the trajectory ⃗ rcontact (red line) undergoes a discontinuous jump when crossing the point indicated by the black star.(e) The wedge-shaped patch (red grid) indicates the region over which the force in Eq.S6 acts to remove the discontinuity in force that occurs in ∂uint/∂d αi,βj in the case of locally concave membrane geometries.The wedge also removes a similar discontinuity relative to l (β−1)j .
(c data is plotted versus the number of pixels contained inside the shape over an experimentally relevant range.Estimate of the error in A as a function of the area of the shape in the synthetic images (blue line).The red and black solid lines represent the mean areas in the experimental data sets of embryo and wing disc cells, respectively, and the spacings between the dotted lines represent the standard deviations.(d) We plot A(t) for the embryo (red) and wing disc (black) experiments using the variance-weighted mean (open circles) over cells adjacent to the wound, where the variance is given by E 2 (a).We include error bars using the standard error of the weighted mean.The solid lines are moving averages with the same window size as in Fig. 1b.We overlay the simulation results for the embryo (yellow triangles), the wing disc (cyan triangles), and the wing disc using cells with shape memory (blue squares), all of which are the same data as in Fig. 4 in the main text.
FIG. S1. (a)Vertex-segment pair potential energy uint and (b) vertex-segment force fint is plotted against the vertex-segment distance d for several values of the attraction strength l1/σ.uint, fint, and l1 are all nondimensionalized using the half membrane width σ and membrane interaction strength ϵ. l2 is fixed at 0.3σ for all numerical simulations of wound closure.
FIG. S2. (a)The vertex-segment distance d αi,βj is the shortest distance from vertex α on cell i to line segment ⃗ l βj on j.If the projection P of ⃗ d αi,βj to the segment ⃗ l βj (gray arrow) falls outside of l βj , then d αi,βj is the distance between vertex α on cell i and vertex β or (β − 1) on cell j.(b) Each deformable particle (left) is modeled as a collection of vertices (dark blue circles) connected by line segments (cyan region).In the right image, we also show the exterior-facing half of the cell membrane of a "smooth" deformable particle, which consists of circulo-line segments (cyan) and vertex sectors (blue).Within a smooth deformable particle, the section around each vertex is considered concave when the angle the vertex makes with its two neighboring vertices has an interior angle θ > π, and convex otherwise.(c) For locally convex geometries, when a vertex that overlaps the membrane surface at ⃗ rcontact moves along the indicated trajectory (red line), the vertex-segment distance relative to ⃗ l βj is continuous.(d) For locally concave geometries, the vertex-segment distance (relative to ⃗ l βj ) for a vertex overlapping with the membrane along the trajectory ⃗ rcontact (red line) undergoes a discontinuous jump when crossing the point indicated by the black star.(e) The wedge-shaped patch (red grid) indicates the region over which the force in Eq.S6 acts to remove the discontinuity in force that occurs in ∂uint/∂d αi,βj in the case of locally concave membrane geometries.The wedge also removes a similar discontinuity relative to l (β−1)j .
FIG. S3. (a)The first step in the image analysis pipeline is to transform raw confocal microscopy images into maximum intensity projections along the z-axis, shown here for a Drosophila embryo ectoderm immediately after wounding from Ref.[1].The scale bar has width 10 µm.(b) We next perform automated cell boundary segmentation using Tissue Analyzer[6], an ImageJ plugin for segmentation of single-layered epithelia.(c) We make manual corrections to cell boundaries near the wound in each frame.From the segmented cell boundaries, we calculate (d) the cell area a and the wound area, (e) the cell perimeter p, and (f) the shape parameter A of cells adjacent to the wound.In this example, the wound area is 9048 px 2 or 157 µm 2 .
FIG. S5.(a) A deformable particle is shown with an extensile force dipole with magnitude f that generates elongation over time t.(b) Final cell shape A(T ) of the extensile deformable particle as a function of the cell area stiffness k * a and the normalized plastic relaxation timescale τ † = τ /T , where T is the duration of the stretching simulation.