A theoretical model of efficient phagocytosis driven by curved membrane proteins and active cytoskeleton forces

Abstract

Phagocytosis is the process of engulfment and internalization of comparatively large particles by cells, and plays a central role in the functioning of our immune system. We study the process of phagocytosis by considering a simplified coarse grained model of a three-dimensional vesicle, having a uniform adhesion interaction with a rigid particle, and containing curved membrane-bound protein complexes or curved membrane nano-domains, which in turn recruit active cytoskeletal forces. Complete engulfment is achieved when the bending energy cost of the vesicle is balanced by the gain in the adhesion energy. The presence of curved (convex) proteins reduces the bending energy cost by self-organizing with a higher density at the highly curved leading edge of the engulfing membrane, which forms the circular rim of the phagocytic cup that wraps around the particle. This allows the engulfment to occur at much smaller adhesion strength. When the curved membrane-bound protein complexes locally recruit actin polymerization machinery, which leads to outward forces being exerted on the membrane, we found that engulfment is achieved more quickly and at a lower protein density. We consider spherical and non-spherical particles and found that non-spherical particles are more difficult to engulf in comparison to the spherical particles of the same surface area. For non-spherical particles, the engulfment time crucially depends on the initial orientation of the particles with respect to the vesicle. Our model offers a mechanism for the spontaneous self-organization of the actin cytoskeleton at the phagocytic cup, in good agreement with recent high-resolution experimental observations.

---

Phagocytosis, also termed ‘cell eating’,^1–6 is a cellular process by which cells internalize extracellular particles that are relatively large (0.5–20 μm). Phagocytosis plays an important role in the immune system, enabling immune cells to destroy foreign elements and clear dead cells and other debris.⁷ This process is also important for diagnostic and therapeutic approaches to cancer, for example, using specialized engulfed particles.^8,9 Despite its biological importance, it is still not understood how the actin cytoskeleton is coordinated spatio-temporally during the phagocytosis process.

Phagocytosis (and endocytosis) in general requires receptors on the cell membrane that bind to the surface of the target particles.¹⁰ Binding-driven engulfment, which does not involve the forces of the cytoskeleton, was modelled in ref. 10–15 with purely diffusive dynamics of the receptors^10,12 and with an active drift.^11,14 The engulfment of charged particles was previously modelled and shown to depend on the lateral segregation of mobile charged membrane components (lipids and proteins), leading to a discontinuous wrapping transition.¹⁶ The effect of signalling molecules was also considered in ref. 11, where using drift and diffusion the model successfully recovered the two stages of the process: an initial slow engulfment followed by a fast engulfment of the second half of the spherical object, compared to the experiments of the neutrophil phagocytosis of polystyrene beads.¹¹ Phagocytosis is, however, an active process which is known to involve the forces of actin polymerization that push the engulfing membrane forward.^17,18 In subsequent models, a detailed description of active forces^18–21 and also the bending and elastic energies of the membrane and cortex are considered.^22–26 The cytoskeletal components of phagocytosis were investigated during the engulfment of artificially produced particles^14,27–29 or during bacterial engulfment by the immune cells.^30–32 In ref. 14, it was found that the engulfment process either stalls before the half engulfment of a spherical object or leads to complete engulfment. In ref. 27 and 33, the experiment was performed using micro-pipettes to hold both the immune cell and the target particle, allowing the measurement of the cortical tension. The recent experiments provide high-resolution details of the organization of the actin cytoskeleton at the rim of the phagocytic cup,^34–36 challenging the theoretical modelling of this process.

In reality, however, the shape of the engulfed particles, such as bacteria or viruses, may be highly non-spherical, which motivated several experiments with the artificial particles of various geometries.^28,37–39 In ref. 28, it was found that the engulfment time for non-spherical shapes can be five times larger than that for spherical particles.³⁷ In ref. 38, it was observed that an oblate spheroid was engulfed more easily compared to spheres, while spheres were easier to engulf compared to prolate spheroids, rod-shaped particles, or needles.⁴⁰ Note that the surface area was not kept constant for the particles with different geometries. In theoretical studies, several non-spherical shapes were modelled, such as ellipsoids, rod-like particles, and capped cylinders;^41–43 however, in these studies, the effect of active forces was not considered.

From these studies, we know that the phagocytosis process involves highly complex dynamics of cytoskeleton rearrangements, membrane shape deformations and protein aggregations.⁷ At present, there is no complete theoretical understanding of the dynamics of the self-organization of the membrane and the actin cytoskeleton, including the active forces it exerts, during the engulfment process. Here, we address this problem using a coarse-grained theoretical model.^44–46 In this model, a three-dimensional vesicle that contains curved membrane complexes,^47–49 which recruit the protrusive forces of the cytoskeleton,⁵⁰ was shown to describe spreading over a flat substrate.⁴⁶ The spreading of the vesicle in this model is driven by the self-organization of the curved proteins in the form of a leading edge aggregate, stabilized by the protrusive forces that describe the actin polymerization pressure. We therefore propose to use this model to describe the dynamics of the spontaneous recruitment of the actin cytoskeleton and the membrane spreading during phagocytosis. This modelling approach can be further motivated by the observations of curved membrane proteins that are also related to the actin recruitment at the leading edge of the cell membrane protrusions.^51–54

Using this model, we aim to expose the role of curved membrane proteins, coupled with adhesion, in self-organizing and regulating the spatio-temporal dynamics of the active cytoskeletal forces, during phagocytosis. This minimal model allows us to explore the basic physical mechanisms during this process, which exist in addition to the more complex biochemical signaling, which plays a role in regulating phagocytosis.^6,55

Materials and methods

Theoretical modelling

We consider a three-dimensional vesicle (Fig. 1), which is described by a closed surface having N vertices, each of them connected to its neighbours with bonds, to form a dynamically triangulated, self-avoiding network, with the topology of a sphere.^{44–46,56,57} For the self-avoiding membrane network, the constraint of fixed topology also influences the vesicle shape determination.⁵⁸ The vesicle is kept in contact with a particle, with which the vesicle has an attractive interaction, as shown in Fig. 1. The vesicle energy has four following contributions: the bending energy is given bywhere κ is the bending rigidity, C₁ and C₂ are the principal curvatures and C₀ is the spontaneous curvature. The discrete version of the bending energy is described in the ESI† Section S1. The direct binding energy between protein complexes on nearest neighbour nodes is given bywhere [script letter H] is the Heaviside step function, is the distance between proteins, are the position vectors of i- and j-th proteins, r₀ is the range of attraction, and w is the strength of attraction. The range of attraction is such that only the proteins that are in the neighbouring vertices can attract each other. The active energy cost for displacing a vertex i containing a protein complex is given bywhere F is the magnitude of the active force, representing the protrusive force due to actin polymerization that is acting in the direction of the outward normal vector of the local membrane surface (along n̂_i) and is the displacement vector of the protein complex. The “active” forces in our simulations are implemented as external forces that act on specific nodes of the system. This is performed by giving a negative energy contribution when the points on which these forces act move in the direction of the force. This means that the work done by these forces effectively enters the energy calculation.

Fig. 1 Schematic representation of our model. The vesicle is formed by a closed triangulated surface, having N vertices connected to neighbours by bonds. The red vertices on the surface of the vesicle represent the curved membrane protein complexes, while the blue vertices represent the bare membrane. A zoomed version of a small section of the vesicle surface is shown in the inset. The vesicle is kept in contact with a particle, having an attractive interaction between them. We consider spherical and non-spherical particles, such as spheroid, spherocylinder, dumb-bell, etc.

The above equation indicates that when the proteins are distributed inhomogeneously, there will be a net force on the vesicle in a particular direction.⁴⁶ However, in the present work, we only simulate vesicles that are adhered to particles that are fixed in their location, which effectively links our vesicle to the lab frame thereby restoring momentum conservation (fixing the engulfed particle acts as an infinite momentum reservoir).

Finally, the adhesion energy is given by

where E_ad is the adhesion strength, and the sum runs over all the vertices that are adhered to the phagocytic particle.^45,46,59 By ‘adhered vertices’, we mean all such vertices, whose perpendicular distance from the surface of the particle is less than ε. We choose ε to be equal to the length l_min, which is the unit of length in our model, and defines a minimal length allowed for a bond. Here, we assume that the particle is perfectly rigid such that no vertices can penetrate the particle. Note that in our model we measure E_ad in energy units, since it is defined as the energy per adhered vertex or a node. In the continuum description, however, the adhesion strength is defined in units of energy per area.

We update the vesicle with mainly two moves: (1) vertex movement and (2) bond flip. In a vertex movement, a vertex is randomly chosen and attempted to move by a random length and direction, with a maximum possible distance restricted by 0.15 l_min. In a bond flip move, a single bond is chosen, which is a common side of two neighbouring triangles, and this bond is cut and reestablished between the other two unconnected vertices.^45,46,59 While the vertex movement allows the vesicle to change its shape, a bond flip mechanism provides lateral fluidity to the membrane, by allowing vertices (and the vertices that are occupied by protein complexes) to diffuse on the membrane surface. This lateral motion allows the protein complexes to move on the membrane surface and aggregate with other protein complexes to form larger aggregates. The maximum bond length is restricted to l_max = 1.7l_min. The particle is assumed to be static and does not have any independent dynamics. We update the system using the Metropolis algorithm, where any movement that increases the energy of the system by an amount ΔW = ΔW_b + ΔW_d + ΔW_F + ΔW_A occurs with the rate exp(−ΔW/k_BT); otherwise, it occurs with rate unity. Note that our model only captures the wrapping and engulfment events, while the other upcoming processes, such as the membrane fission,⁶⁰ and the movement of the engulfed particle inside the vesicle are not captured in the present model.

The forces exerted at the locations of the curved proteins are “active” since they give an effective energy (work) term that is unbounded from below (3). This property defines them as out-of-equilibrium, as they can not arise in an energy conserving system. These active forces represent the pushing forces exerted on the membrane by polymerizing actin filaments,^61–64 which we do not model explicitly to reduce the complexity of the model.⁶⁵

Unless specified, we use a vesicle of the total number of vertices, N = 3127 (radius ∼ 20 l_min), where l_min is the unit of length in our model and defines the minimum bond length. The bending rigidity is κ = 20 k_BT, the protein–protein attraction strength is w = 1 k_BT, and ρ = N_c/N is the protein density, with N_c vertices occupied by curved (convex) membrane proteins having spontaneous curvature: C₀ = 1.0 l_min⁻¹. Note that we do not conserve the vesicle volume, although this could be maintained using an osmotic pressure term as the volume constraint.⁴⁵ The membrane area is very well conserved (ΔA/A < 1%).

Note that our MC simulations do not describe the correct dynamics of the fluid motion around the membrane and the engulfed object. However, it gives the energetically most favorable trajectory that controls the dynamics, without the correct time-scales, which depend on hydrodynamic dissipation. In addition, the actin filaments that exert the pushing forces at the leading edge of the phagocytic cup also transmit traction forces to the engulfed particle, which can therefore move and rotate it. These traction forces are not explicitly described, as we maintain the fixed position of the engulfed particle, similar to some experimental set-ups.⁶⁶

In this numerical study, we chose a fixed vesicle and spherical particle size, for which we have a reasonable balance between accuracy and simulation time. The vesicle radius was chosen to be R_vesicle ∼ 20 l_min, while the particle radius was taken to be R = 10 l_min (smaller and larger particles are explored in the ESI† Section S3 and Fig. S1). The mesh size for these dimensions allows us to describe the membrane shapes during the engulfment process with good accuracy, even along the sharp (highly curved) leading edge. The leading edge curvature is determined by the spontaneous curvature of the curved proteins, which we kept from our previous work on vesicle spreading.⁴⁶ We quantified the accuracy by comparing the calculated mean curvature of the membrane with the discrete version and found that, even for spherical vesicles with a radius of ∼1 l_min, the error is only ∼10% (Fig. S21, ESI†). Clearly, membrane features that are on smaller lengthscales than l_min are not captured by this calculation.

Note that our continuum description of the membrane is valid for lengthscales that are larger than that of the membrane width or the single-protein size. Our simulations are therefore valid for particles that are larger than this lengthscale. Furthermore, we use the term “curved proteins” to describe any membrane-bound or membrane embedded protein complexes, or larger nano-domains, that have the properties of the spontaneous curvature and ability to recruit the nucleation of actin polymerization. Let us emphasize that our continuum model of the membrane is not a microscopic model, where each node of the mesh represents a single protein or a lipid.

The length scale l_min that we use in our simulation does not necessarily correspond to a particular length in reality. By assigning a different real length to the basic length-scale of the simulation (l_min), we do not change the resulting dynamics for the following reason: the dynamics of the engulfment process depends on the energy gain/loss per additional adhered area increment (vertex). The adhesion energy per adhered vertex rescales as 1/l_min², as does the bending energy per vertex (because the total bending energy of a specific shape is scale-invariant, provided we also rescale the spontaneous curvature (C₀) accordingly). The active work is also independent of the absolute scale of l_min, as we define the force exerted by each vertex to be in units of k_BT/l_min, and therefore, its contribution to the total energy increment per unit area, also scales as 1/l_min². The absolute scaling of l_min does not therefore change the qualitative nature of our simulation results: both the adhesion energy and the bending energy per unit area, scale as 1/l_min².

There are important length scales in this system which we wish to discuss. The first is the lengthscale given by the ratio between the adhesion and bending energy of the adhered membrane. This ratio determines the minimal radius of the spherical particle: , below which the membrane does not lower its energy when adhering and therefore adhesion vanishes completely.⁶⁷ Above the critical radius, the engulfment is inhibited by the interplay between the adhesion energy of the adhered membrane and the bending energy of the non-adhered membrane (see the ESI†, Fig. S2). The bending energy in this regime depends on the ratio of the vesicle to the particle size,⁶⁷ such that large particles are not engulfed due to the finite size of the vesicle. We demonstrate this using our simulations for the protein-free vesicle in the ESI† (Fig. S4). This analysis and the dynamics of the engulfment in this regime are modified when the vesicle contains curved proteins, as well as active forces, and is the main focus of this work. We leave the detailed exploration of the effects of some of the parameters of our model for future studies, such as C₀.

Imaging of macrophage phagocytosis

The culturing of RAW264.7 macrophages (ATCC) and generation of murine bone-marrow derived macrophages is described in a previous study.³⁴ Macrophages were transfected with mEmerald-Lifeact (Addgene, #54148) to label F-actin using the Neon electroporation system, according to manufacturer's instructions, and allowed to recover for 24 h. Flash Red polystyrene beads (Bangs Laboratories Inc., 7 μm diameter) were washed three times in sterile PBS and opsonized overnight at 4 °C in 3 mg mL⁻¹ mouse IgG (Invitrogen). To remove excess antibody, beads were washed three times with PBS and resuspended in sterile PBS. Beads were applied to macrophages, plated on 5 mm round coverslips, in a 37 °C-heated, water-coupled bath in an imaging medium (FluoroBrite (Thermo Scientific), 0–5% FBS, Pen/Strep) prior to imaging acquisition. Imaging was performed using a lattice light-sheet microscope operated and maintained by the Advanced Imaging Center at the Howard Hughes Medical Institute Janelia Research Campus (Ashburn, VA). 488, 560, or 642 nm diode lasers (MPB Communications) were operated between 40 and 60 mW initial powers, with a 20–50% acousto-optic tunable filter transmittance. The microscope was equipped with a Special Optics 0.65 NA/3.75 mm water dipping lens, excitation objective, and a Nikon CFI Apo LWD 25 × 1.1 NA water dipping collection objective, which used a 500 mm focal length tube lens. Images were acquired using a Hamamatsu Orca Flash 4.0 V2 sCMOS camera using the custom-written LabView Software. Post-image deskewing and deconvolution were performed using the HHMI Janelia custom software and 10 iterations of the Richardson–Lucy algorithm.

Results

We simulate the dynamics of a three-dimensional vesicle, which is described by a closed surface having N vertices, each of them connected to its neighbours with bonds, forming a dynamically triangulated, self-avoiding network, with the topology of a sphere.^45,46 On this network, we describe a population of curved proteins, which can diffuse in the membrane, and exert active forces. Within our model, the actin polymerization of individual filaments is not explicitly described,⁶⁵ but is mapped to out-of-equilibrium (active) forces that locally push outwards on the membrane. We first study the case where the phagocytic particle is spherical and then generalize it for non-spherical particles as well.

We chose the surface area of the vesicle to be much larger (at least four times larger) than the particles, so that the engulfment is not inhibited by the membrane surface area limitation. During phagocytosis, the cell is able to increase the local membrane area,⁶ to allow for the phagocytosis of larger particles. Note that our vesicle does not represent the entire shape of the cell and does not include the nucleus and other constraints. Furthermore, we do not explicitly take into account the lateral⁶⁸ and surface⁶⁹ membrane tension, but our calculation conserves the membrane area (total area fluctuations of <1%). In addition, we do not fix the vesicle volume, corresponding to the vanishing osmotic pressure gradient across the membrane, as maintained in cells.⁷⁰ This allows our vesicle to undergo both large volume and shape changes during the engulfment process. Introducing a high osmotic pressure inside our simulated vesicle would act to conserve the volume,⁴⁵ and in this way hinder the engulfment process by inhibiting large vesicle shape deformations. The effects of osmotic pressure on the phagocytosis process can be further investigated in the future.

Spherical particles – passive curved proteins

As a first validation of our model, we simulated the engulfment process of a spherical particle by a protein-free vesicle (see the ESI† Sections S3, S4, Fig. S1–S4 and Supplementary Movies S1, S2). In agreement with previous studies,⁷¹ we found an optimal size for the full engulfment of the particle in this system; very small particles cannot be engulfed due to the high bending energy, while very large particles require larger vesicles in order to be fully engulfed. For our further simulations, we choose the particle radius to be R ∼ 10 l_min (where l_min is the minimal allowed length of an edge on the triangulated surface), which is a suitable size for the parameter ranges we simulate here. The protein-free case can be compared to a simplified analytical model (see the ESI† Section S4) and used to calibrate the relevant parameter regime. Note that a fully engulfed state in our simulations is when the vesicle encloses the particle, such that only a small hole (of size 1–3 triangles) remains exposed, as shown in Fig. 2d (leftmost image). Since we do not allow for membrane fusion, this hole remains.

Fig. 2 Sphere engulfment by a vesicle with passive proteins. (a) Phase diagram in the E_ad–ρ plane, where E_ad is the adhesion strength and ρ is the density of curved membrane proteins. The background color shows the adhered area fraction. The red dashed line is approximately separating the fully engulfed state and the partially engulfed state. The snapshots are shown for E_ad = 0.60, 0.80, 1.0, and 1.20 (in units of k_BT) and ρ = 1.6%, 4.8%, 9.6%, 12.8%, 16%, and 19.2%. (b) Comparison of the adhesion energy gain (dashed horizontal line) and the bending energy cost (ΔW_b, points) per adhered node (both having units of k_BT) as a function of the engulfment fraction, for a fixed adhesion strength E_ad = 1.0 k_BT and various values of ρ. Blue symbols denote ρ = 0.8%, green circles denote ρ = 2.4% and red triangles denote ρ = 4.0%. The origin of the large error bars is due to the large fluctuations in the bending energy at the leading edge of the membrane and having a relatively small ensemble of simulations for each case. (c) Mean cluster size (the mean number of proteins in a cluster) as a function of the engulfed area fraction for both partial and complete engulfment cases. The color codes are same as shown in Fig. 2b. The inset shows the increase in the density of curved proteins at the rim of the phagocytic cup during the engulfment process for ρ = 4.0%. (d) Examples of the cross-sections of the vesicle membrane around the particle, at the end of the engulfment process, for different protein densities. We use E_ad = 1.0 k_BT here. the other parameters are the total number of vertices (nodes) in the vesicle is N = 3127, the radius of the spherical particle is R = 10 l_min, the bending rigidity is κ = 20 k_BT, the protein–protein attraction strength is w = 1 k_BT, and the spontaneous curvature of curved (convex) membrane proteins is C₀ = 1.0 l_min⁻¹.

Next, we consider vesicles with passive (F = 0, where F is the active protrusive force exerted by the proteins) curved proteins. In Fig. 2a, we show the phase diagram in the ρ–E_ad plane (where E_ad is the adhesion energy per adhered node and ρ is the areal density of the proteins on the vesicle membrane). For small E_ad, below ∼0.60 k_BT, we only obtain partial engulfment, with the engulfment area increasing with ρ but not reaching the full engulfment. For intermediate E_ad (0.80–1.0 k_BT), there is an incomplete engulfment for small ρ (Movie S3, ESI†), but the engulfment is complete for larger ρ (Movie S4, ESI†). This is the interesting regime, where engulfment is facilitated by the presence of the curved proteins.

Note that we calculate the adhered fraction of the spherical particle (A_ad/A_max), using the maximum adhered area for A_max. This may be larger than the true surface area of the particle, since we defined a finite width of the adhesive interactions, and within this width the adhered membrane may fluctuate, therefore containing a larger membrane area than the area of the rigid bead surface (see for example the cross-sections in Fig. 2d, and for more details see the ESI,† Section S2).

Increasing ρ reduces the bending energy cost of the engulfment by the spontaneous aggregation of the curved proteins at the sharp leading edge of the engulfing membrane and allows the vesicle to adhere more and engulf the particle. This is the same mechanism that allows the passive curved proteins to drive vesicle spreading on flat adhesive substrates.⁴⁶

The engulfment process can be understood by considering the competition between the adhesion and the bending energy cost of the vesicle as the engulfment proceeds. In Fig. 2b, we compare the adhesion energy gain (E_ad) and the bending energy cost (ΔW_b), per additional adhered membrane node, as a function of the adhered fraction of the spherical particle. The bending energy cost increases with the adhered fraction, and for small ρ, the engulfment process stops where the bending energy cost becomes higher than the adhesion energy gain per node. In Fig. 2b, we show two cases of partial engulfment, for ρ = 0.8% and 2.4% (blue and green points), where the bending energy cost becomes larger than the adhesion energy gain, and the system reaches steady-state at A_ad/A_max values of ∼0.2 and 0.35, respectively. At these low densities, the curved proteins fail to form large aggregates at the leading-edge of the engulfing membrane, as shown by the mean cluster size (Fig. 2c), and therefore do not lower the bending energy cost sufficiently (see the insets of Fig. 2c for the adhered fractions that are smaller than 0.5).

The complete engulfment is achieved for higher densities, where the adhesion energy gain is larger than the bending energy cost. For ρ = 4.0%, the bending energy cost is always smaller than the adhesion energy gain, and the system proceeds to complete engulfment (red points in Fig. 2b). This is driven by the formation of a large protein cluster along the leading-edge of the engulfing membrane (Fig. 2c and insets for adhered fractions that are larger than 0.5). However, the aggregation at the leading edge is transient, since as the radius of the hole gets too small, the mean curvature of the rim of the hole decreases and is less favorable for the curved proteins (ESI† Section S5 and Fig. S5), which leads to the dispersal of the curved proteins on the membrane (Fig. 2d, ρ = 4.8%). This is the origin for the appearance of a maximum value of the mean cluster size during complete engulfment (Fig. 2c).

Note that in addition to the adhesion energy gain, there is also the gain in the protein–protein binding energy, which further helps drive the membrane over the bending energy barrier, but this is a small contribution due to the small size of the ring-like clusters of proteins (Fig. 2c). Following the completion of the engulfment, this ring-like cluster of curved proteins disperses spontaneously. Since the topology of the membrane is fixed, we do not allow for membrane fusion, and there remains a tiny open membrane neck that connects the engulfed sphere with the space outside the vesicle (Fig. 2d).

For very large ρ (Fig. 2a), however, the wrapping is again incomplete. The large number of curved proteins stabilizes the cluster at the membrane rim and prevent the hole from closing up (Movie S5, ESI†). In Fig. 2d, we show an example of the cross-sectional view of the vesicle at the end of the engulfment process. For high densities (ρ = 9.6% or larger), the protein cluster remains stable at the rim of the phagocytic cup. This cluster maintains the curvature of the rim to fit better the spontaneous curvature of the proteins, and this prevents the rim from shrinking, resulting in incomplete engulfment (ESI† Section S6 and Fig. S6). Note that at large densities the curved proteins form pearled clusters, which are not directly related to the rim cluster at the edge of the phagocytic cup.

For large E_ad(≥1.2 k_BT), the adhesion interaction between the vesicle and the particle is strong enough to drive complete engulfment even in the absence of curved proteins. The red dashed line in Fig. 2a encloses the region within which we find complete engulfment. The shape of this complete engulfment transition line can be qualitatively reproduced by a simplified analytical model (ESI† Section S7 and Fig. S8, S9).

Finally, we also calculated the engulfment time as a function of ρ (see the ESI† Section S8 and Fig. S10), where we found that it is decreasing with increasing ρ. The engulfment time diverges for ρ ∼ 2.25%, below which there is no engulfment (for E_ad = 1.0 k_BT).

Spherical particles – active curved proteins

Next, we study the engulfment of a rigid particle by a vesicle containing active curved proteins (Fig. 3). In order to expose most clearly the role of activity during the engulfment process, we start with a low concentration of curved proteins (ρ = 1.6%), such that their passive effect is weak (Fig. 2a). In Fig. 3a, we show the engulfment phase diagram in the F–E_ad plane. We found that the active forces enable full engulfment even at low adhesion strengths, where the passive curved proteins (F = 0) are not sufficient to drive full engulfment. For small E_ad, the engulfment is incomplete even for large F, while for intermediate E_ad, the engulfment becomes complete as F increases (Movie S6–S8, ESI†). For large values of E_ad, the engulfment is achieved even for the passive system.

Fig. 3 Sphere engulfment by a vesicle with active proteins. (a) Phase diagram in the E_ad–F plane for a low concentration of curved proteins ρ = 1.6%, where F represents the protrusive force due to actin polymerization that is coupled with the curved membrane proteins. The background color shows the adhered area fraction. The dashed yellow line approximately separates the fully engulfed (above the line) states from the partially engulfed states. The snapshots are shown for E_ad = 0.60, 0.80, 1.0, and 1.20 (in units of k_BT) and F = 0.20, 0.40, 0.60, 0.80, and 1.0 (in units of k_BT/l_min). (b) Comparison of the bending energy cost with the adhesion energy gain (both having units of k_BT per adhered node) for different values of F. The blue boxes, green circles and red triangles are for F = 0.1, 0.4 and 2 (in units of k_BT/l_min), respectively. The red dashed line shows the combined effective energy gain per adhered node, including the adhesion and the work performed by the active forces (for F = 2 k_BT/l_min). (c) The mean cluster size as a function of the engulfed area fraction for the same cases shown in (b) (same color code). Lower snap shots correspond to F = 0.1 and 0.4 and top snap shots correspond to F = 2. (d) Snap shots of the vesicle (F = 2 k_BT/l_min), indicating the contribution of the active forces to the the engulfment. The arrows show the tangential component of the individual forces (small arrows), and the direction of the total force (large arrow), for proteins that are close to the spherical surface. For (b–d), we use E_ad = 1 k_BT and ρ = 1.6%. (e) Mean cluster size for large E_ad and ρ. Green circles are for ρ = 6.4% and F = 1.0 k_BT/l_min, and blue boxes are for ρ = 4.8% and F = 2.0 k_BT/l_min. We use E_ad = 1.5 k_BT for both the cases.

For the passive case, we analyze the engulfment process by plotting the adhesion energy gain and the bending energy cost per adhered node, as a function of A_ad/A_max. In Fig. 3b, we compare both the case of the complete and incomplete engulfment, as the force is increased. As in the passive case, incomplete engulfment occurs when the bending energy cost per node becomes larger than E_ad (blue boxes in Fig. 3b). For a higher force (F = 0.40 k_BT/l_min), the active force enhances the aggregation of the proteins around the rim and thereby lowers the bending energy cost, resulting in complete engulfment (green circles in Fig. 3b).

Even higher forces (F = 2 k_BT/l_min) produce large shape fluctuations of the vesicle, and the bending energy cost increases. However, a complete engulfment is achieved nevertheless (red triangles in Fig. 3b and top insets in Fig. 3c). In this case, the engulfment is enabled by the active forces that directly push the membrane over the spherical surface, and complete engulfment is achieved even when the bending energy cost per node goes above E_ad. In order to extract the contribution due to active forces in this case, we define a region close to the particle (phagocytic cup or rim) and assume that only the forces from these proteins are helping the vesicle to engulf the particle. We then calculate the magnitude of active forces in the tangential direction to the surface of the particle (Fig. 3d). We then divide this total force F_total by the length of the leading rim of the membrane to extract the effective active work done per adhered node (see the ESI† Section S9 and Fig. S11 for details). We add this active work to the adhesion energy gain per adhered node (the red dashed line in Fig. 3b), which is shown to be sufficient to offset the bending energy cost per node. The magnitude of the total active force (without scaling by the rim length) as a function of the adhered fraction is shown in the ESI,† Fig. S11(c). It is clear from Fig. S11(c) (ESI†) that the total active force is very small during the early stages when the curved protein complexes are still dispersed over the membrane surface, so that their active forces are pointing in random directions, before they form the organized aggregate at the leading edge of the phagocytic cup.

Note that in the active case the complete engulfment can proceed at much lower densities of proteins, compared to the passive case, such that there are not necessarily enough of them to complete a ring-like cluster at the leading edge. An arc-like cluster of active proteins, which forms at an earlier stage of the engulfment process (red triangles in Fig. 3c), can be sufficient to drive a flat protrusion that spreads over the spherical particle to complete the engulfment (Fig. 3c and d). This is further shown in Fig. 3e for larger protein densities (Movie S9 and S10, ESI†). The spherical particle is transparent to ensure the visibility of the protein clusters along the rim of the phagocytic cup. Due to the higher density, the rim clusters form at a lower adhered fraction (compared to Fig. 3c), but they are still highly fragmented and do not merge to form a complete circular ring until very late stages of the engulfment (Fig. 3e).

Many of the qualitative spatio-temporal features that we found for the curved active proteins in our simulations match very well with the recent high-resolution experimental images of the actin organization along the leading edge of the engulfing membrane protrusion during phagocytosis (Fig. 4 and Movie S11, S12, ESI†).³⁵ In both the experiments and the simulations, we found that the actin ring⁵⁰ is initially discontinuous, and highly fragmented, in the form of dispersed “teeth” and “arcs” along the leading edge. The actin aggregate becomes more cohesive in the final stages of the engulfment, before it disperses after the engulfment is complete (also see the ESI† Section S10 and Movie S13). Note that in our model the curved membrane complexes are also the sites of the recruitment of actin polymerization, which therefore ends up exerting a protrusive force on the membrane at these sites. We therefore identify the membrane accumulation of the curved proteins in our simulations with the accumulation of cortical actin in the experimental observations.

Fig. 4 Comparison of the actin organization at the rim of the phagocytic cup observed in simulations (a) and experiments (b). In both the in vitro experiments and the simulations, we found that the actin ring is highly fragmented, in the form of dispersed “teeth” (yellow arrows) or in more continuous “arcs” (red arrows) along the leading edge. The “arcs” often connect to form a complete ring around the leading-edge rim during the final stage of the engulfment process, both in the simulations and the experiments. The actin aggregate becomes more cohesive in the final stages of the engulfment, before it disperses after the engulfment is complete. (a) Snapshots from a simulation, using E_ad = 1.50 k_BT, ρ = 6.4% and F = 1.0 k_BT/l_min. (b) Time-lapse sequences of maximum intensity projection images illustrating the engulfment of immunoglobulin-G-coated polystyrene beads (7 μm diameter) by (upper) RAW264.7 macrophage-like cell lines and by (lower) murine bone-marrow derived macrophages, imaged by lattice-light sheet microscopy. In (b), cells were transfected with mEmerald-Lifeact to label F-actin. Scale bar 5 μm, the time is indicated in min:sec.

When comparing to the experiments, we have to remember the limitations of our numerical simulation, which can only describe membrane features that are on lengthscales larger than l_min (for a particle of radius 3.5 μm, l_min = 350 nm). However, the main features of the actin organization that we highlight in Fig. 4a, namely the fragmented ring and its continuous closure, are on larger lengthscales and therefore reliably emerge in our simulations. We therefore found this choice of mesh resolution to be sufficient in order to explain the outstanding features of the membrane and actin organization during engulfment, as observed in the experiments, while maintaining a reasonable computation time. Future studies using a finer mesh may allow smaller scale features of this process to be studied.

Spherical particles – engulfment dynamics

The engulfment dynamics of the different systems are shown in Fig. 5a. We note that MC simulations do not provide us with the correct real time dynamics, as the hydrodynamic dissipation processes are not described. Nevertheless, we can obtain the qualitative form of the dynamics and obtain relative engulfment times for comparison between different cases.

Fig. 5 Engulfment dynamics for spherical particles. (a) Different growth behaviours of the adhered fraction with time for protein-free, passive and active cases. We scale the area axis by the maximum adhered area (A_max) and time axis by the engulfment time (t_eng). The blue cross symbols are for the protein free case with E_ad = 1.30 k_BT; the green star symbols are for the passive proteins with E_ad = 1.0 k_BT and ρ = 4.0%; the red boxes are for E_ad = 1.0 k_BT, ρ = 1.6% and F = 0.40 k_BT/l_min; the yellow circles are for E_ad = 1.0 k_BT, ρ = 1.6% and F = 2.0 k_BT/l_min; the black triangles are for E_ad = 1.50 k_BT, ρ = 6.4% and F = 1.0 k_BT/l_min. (b) Engulfment time as a function of F, showing the non-monotonic behavior: below and above a critical force (green and blue shaded areas, respectively), there is partial engulfment. (c and d) Distributions of the engulfment time for F = 1,2 k_BT/l_min, respectively.

For the protein-free and passive proteins systems, driven by adhesion alone, we found in our simulations the typical “two-stage” behavior: an initial slow growth up to an adhered area fraction of ∼0.5, and then a very fast stage to complete engulfment. A similar behavior is also observed for the active system in the regime of low protein densities. This dynamics was observed experimentally, and for the passive engulfment system it was also motivated by theoretical modelling.^11,16

However, for the active system at high protein densities, where the rim clusters form at an early stage of the engulfment process (Fig. 3e), we found a much more uniform dynamics with steady growth of the adhered fraction (Fig. 5a). This agrees qualitatively with the observations in,⁶⁶ as well as,^34,35 which indeed correspond to this regime (Fig. 3e and 4).

In Fig. 5a, the engulfment dynamics are shown as a function of the normalized time, so the actual engulfment duration is not presented. The dependence of the mean engulfment time on the active force is found to be non-monotonic (Fig. 5b). Below a critical force, in the regime of low E_ad and low ρ, where the passive system does not fully engulf (Fig. 3a), there is no complete engulfment (the green shaded area). Above this regime, the increasing force induces complete engulfment with faster and more robust engulfment as F increases. The engulfment time distribution in this regime is relatively narrow (Fig. 5c).

However, increasing F beyond some values, which are set by the relationship of the adhesion energy and the active work, leads to increasing both the mean value and the variability of the engulfment time. This occurs due to active forces stretching the membrane away from the adhered particle (top insets in Fig. 3c), which gives rise to long “waiting times” during which the engulfment does not progress. Following this “waiting time”, once the protein aggregate forms at the membrane rim (Fig. 3c), the engulfment proceeds rapidly as shown in Fig. 5a (yellow circles). This variability is manifested in the very broad distribution of engulfment times (Fig. 5d).

At even larger forces, the stretching of the membrane sideways by the active forces pulls the membrane away from the adhesive particle and prevents engulfment (the blue shaded area in Fig. 5b, and ESI† Section S11, Fig. S13(a), Movie S14).

In this regime of large F, if we further increase ρ, we reach the regime where a free vesicle will form a flat pancake-like shape, where all the proteins are aggregated at the rim.⁴⁵ In this regime, the vesicle can engulf the particle if this process is faster than the timescale of the formation of the pancake-like shape (ESI† Section S11, Fig. S13(b) and Movie S15). For very large forces, we found again a non-engulfed regime (similar to the blue shaded area in Fig. 5b and ESI† Fig. S13(c), Movie S16).

Engulfment of non-spherical particles

Immune cells often engulf highly non-spherical objects, such as bacteria. Furthermore, studies with different particle shapes indicated a strong shape-dependence of the engulfment time and the success rate.^12,14,28 We therefore study the engulfment of non-spherical particles by our vesicles that contain passive or active curved proteins.

We first consider spheroids. We keep the surface area constant, i.e. same as for the sphere of radius 10 l_min. We use prolate shapes (R_x = R_y < R_z) and oblate shapes (R_x = R_y > R_z), such that the aspect ratios are R_x/R_z < 1 for prolate and R_x/R_z > 1 for oblate.

In Fig. 6, we show the mean engulfment time for oblate and prolate spheroids, for both the cases when the vesicle is initially in contact from the top (the poles of the shape) or from the side (along the equator), as shown in the insets. For both passive and active curved proteins, we found that the engulfment times increase as the shape deviates from a sphere. This is in agreement with previous experimental observations²⁸ and can be attributed to the higher local curvature on the non-spherical shapes, which increase the bending energy barrier that opposes the progression of the leading edge membrane during their engulfment. Furthermore, we found that the engulfment time increases more rapidly for the oblate shapes, compared to the prolate shapes, for both passive and active, as was observed in experiments.²⁸

Fig. 6 Mean engulfment time for non-spherical particles. The vertical dashed line is separating oblate (right) and prolate (left) shapes. Here, for the passive case, we use E_ad = 1.0 k_BT and ρ = 4.8% and for the active case E_ad = 1.0 k_BT, ρ = 1.6%, and F = 1.0 k_BT/l_min. The surface area of all the shapes is constant and is equal to that of a sphere of radius 10 l_min. The insets demonstrate the initial conditions of either starting from the top (poles) of the shapes (along their axis of rotational symmetry) or from the side (equator).

For the particles with a small aspect ratio (R_x/R_z close to unity), the engulfment time of all the cases is comparable. However, for the particles with large aspect ratios, there are significant differences in the engulfment times, between the passive and active systems. For the passive curved proteins (blue triangles and red boxes), we found that the engulfment times increase sharply with the deviation from the spherical shape. Furthermore, we found that initialising from the side (red boxes) results in relatively shorter engulfment times, in particular for the oblate shapes. For an oblate particle with a large aspect ratio (for R_x = R_y > 11 l_min), the particle is engulfed from side, while it is not engulfed from the top (Fig. 6). For prolate shapes with large aspect ratios, however, the engulfment times seem comparable for both initial conditions. The differences between the engulfment dynamics on the oblate and prolate shapes arise from the different local curvatures that the edge membrane encounters during engulfment: when adhering to the side (rim) of the oblate shape, the leading edge moves mostly on the flat surfaces, which allow the balance of adhesion versus bending to be favorable for the smooth progression of the forward motion. When spreading from the top of the prolate shape, the membrane edge encounters the highly bent surface all along the leading edge, and therefore faces a higher bending energy barrier for the forward movement of the membrane.⁷² When spreading from the side of the prolate shape, it is similar to the spreading of the vesicle on a cylinder, where one direction has low curvature, and therefore faster engulfment (see these dynamics in terms of the energy of adhesion and bending in Fig. S15, ESI†).

By comparison, for active curved proteins, the engulfment time is lower than that for the passive case, even with a lower protein density, and the active curved proteins enable complete engulfment over a larger range of aspect ratios. We also note that for the active curved proteins it is the prolate shape that exhibits a dependence of the engulfment time on the initial conditions, where starting from the side (black circles) leads to faster engulfment. For the oblate shapes, however, the engulfment time is largely independent of the initial conditions, opposite to the behavior of the passive curved proteins.

For a non-spherical particle with high aspect ratios, due to the highly curved poles (in the prolate) or equator (oblate) regions, the engulfment process is slower and the variability in engulfment times increases. The slower engulfment allows the vesicle to reorient on the surface of the particle, driven by both passive energy minimization and active work. We show the details of this reorientation process in the ESI† Section S12 (Movies S17–S22).

When the vesicle is initially at the top of a prolate shape with a high aspect ratio, the vesicle takes a longer time to form the initial adhesion since the particle's surface facing the vesicle has high curvature in all directions. Following this slow initial adhesion time, the vesicle reorients itself always to the side and engulfs the particle, along the easier direction where the long axis of the prolate shape presents a direction of the low curvature (Fig. S15e and f, ESI†). For passive curved proteins, this reorientation process is faster and thereby the engulfment time when started from the top or from the side are comparable (blue triangles and red boxes in Fig. 6 for R_x/R_z < 1). For active curved proteins, this reorientation process is slower because of large shape fluctuations of the vesicle, and thus the engulfment time when started from the side is smaller than that from the top (green inverted-triangles and black circles in Fig. 6 for R_x/R_z < 1).

For the oblate shapes, we found a similar reorientation process when the vesicle starts the adhesion from the highly curved sides. In this geometry, the passive vesicle can engulf from side without any reorientation, which is very fast, or it reorients to the top of the oblate shape and ends up with no engulfment. Therefore, when we calculate the engulfment time, we only average over those realizations where the particle was engulfed without any reorientation (red boxes for R_x > 11). Otherwise, the mean engulfment time of oblate shapes by passive vesicles diverges for both the top and side orientations around an aspect ratio of ∼1.25 (Fig. 6). This value can be explained by the bending energy barrier along the oblate rim and has the local radius of the curvature similar to the smallest spherical particle that can still be engulfed (ESI† Section S13 and Fig. S17).

For the active case, the vesicle reorients in every realization and engulfs from the top. This reorientation process is very fast and the engulfment time is similar for starting either from the side or from the top (green inverted-triangles and black circles in Fig. 6 for R_x/R_z > 1).

We further explore more extreme non-spherical shapes, such as elongated sphero-cylinders and dumb-bell shapes in the ESI† (Section S14 and Movies S23–S26). The elongated particles resemble the shapes of different bacteria or fungi such as E. coli, Bacillus subtilis and budding yeast.^12,50,73

Discussion

In this study, we demonstrated that curved (convex) membrane-bound protein complexes, or curved membrane nano-domains,⁴⁹ both passive and active, can strongly affect the engulfment process of adhesive and rigid particles by a vesicle. The engulfment process that we calculate, in particular using active curved proteins, exhibits many features that are observed during the phagocytosis of rigid particles and pathogens by living cells. We therefore propose that our model exposes the physical mechanisms for the self-organization of the membrane and the actin cytoskeleton in the phagocytic cup. These physical mechanisms operate in the cell, in addition to the layers of signaling that allows the cell to optimize the activity of the cytoskeleton near the engulfed surface. However, the curvature-based feedback is shown here to be sufficient to drive the basic underlying dynamics.

We found that the curved proteins spontaneously aggregate around the highly curved rim of the phagocytic cup and enable the engulfment to proceed efficiently by reducing the bending energy costs. The active forces in our model, which represent the protrusive forces due to the recruited actin polymerization, further facilitate the engulfment, due to several effects that we have identified: the protrusive force at the membrane leading edge promotes the aggregation of the curved proteins (reducing the bending energy cost) and directly pushing the membrane over the engulfed particle. However, our model predicts that if the active forces are too large, they induce membrane ruffles that are inefficient for smooth engulfment and may hinder the phagocytosis process.

In addition, the active forces can stretch the membrane to wrap around sharp corners, which enables the engulfment of non-spherical particles. Compared to the passive curved proteins, the curved active proteins drive more efficient engulfment at lower adhesion energy, at lower curved protein density, and of highly non-spherical particles. For oblate and prolate shapes, we recover the experimentally observed increase in engulfment time with the aspect ratio. We expose the mechanisms that allow the active forces to more robustly overcome the local curvature barriers during the engulfment of such non-speherical objects. These results of our model highlight the central role that the recruitment of actin plays during phagocytosis, which is usually abolished by the inhibition of actin polymerization (see for example⁷⁴).

Our model predicts that the actin aggregation at the rim of the phagocytic cup is typically fragmented and does not form a complete circular ring. Nevertheless, it is still highly effective in driving the engulfment. Furthermore, the actin aggregate forms a more cohesive ring in the final stages of the engulfment, before it spontaneously disperses as the engulfment is complete. These features of our model are verified and observed in recent high-resolution imaging of the actin during phagocytosis³⁵ (Fig. 4).

The model we presented currently lacks the contractile forces that arise at the phagocytic cup due to myosin-II activity.³⁵ The contractile forces play an important role during the engulfment of soft objects, such as cells,^75–77 but they are dispensable for the phagocytosis of more rigid particles.^35,74 The contribution of myosin-II-mediated contractility may be more or less pronounced depending on the particle size and stiffness and, in addition, could influence mechanosensitive signaling within phagocytes rather than directly drive phagocytic cup extension and closure. Our results should therefore be important even without this feature, which we plan to add to our model in the future, in particular when modelling the engulfment of flexible particles.

To conclude, we demonstrate that curved membrane proteins can play a key role during phagocytosis. Convex-shaped proteins or a curved complex of proteins, which furthermore recruit actin polymerization, spontaneously aggregate at the leading edge of the membrane and drive robust engulfment. This is similar to the role of these proteins that was recently found at the leading edge of the lamellipodia, in experiments^53,54 and theory.⁴⁶ We also show that concave-shaped proteins can further enhance the engulfment process by reducing the bending energy cost of the membrane that adheres to the particle's surface^78,79 (ESI† Section S15, Fig. S20 and Movie S27). Our results should motivate future experimental studies aimed at the identification and characterization of the different curved membrane components that participate in phagocytosis.^3,52,80

Author contributions

RKS, SP, AI and NG developed the theoretical simulation model. RKS carried out the numerical simulations. SP developed the original MC simulation program under the supervision of AI. RKS and NG did the theoretical analysis. SRB, MK and NCG did the experiments and their analysis. RKS and NG wrote the manuscript and edited by all the authors.

Conflicts of interest

The authors declare no competing interest.

Acknowledgements

We acknowledge useful comments by Orion Weiner. N. S. G. is the incumbent of the Lee and William Abramowitz Professorial Chair of Biophysics and acknowledges support by the Ben May Center for Theory and Computation and the Israel Science Foundation (Grant No. 207/22). This research was made possible in part by the historic generosity of the Harold Perlman Family. A. I. and S. P. were supported by the Slovenian Research Agency (ARRS) through the Grant No. J3-3066 and J2-4447 and Programme No. P2-0232. LLSM imaging was performed at the Advanced Imaging Center (AIC)-Howard Hughes Medical Institute (HHMI) Janelia Research Campus. The AIC is jointly funded by the Gordon and Betty Moore Foundation and the Howard Hughes Medical Institute. MK was supported by the National Institute of General Medical Sciences of the NIH under the Award R01GM138652.

References

1. G. Cannon and J. Swanson, J. Cell Sci., 1992, 101, 907–913.
2. S. Kumari, S. MG and S. Mayor, Cell Res., 2010, 20, 256–275.
3. R. S. Flannagan, V. Jaumouillé and S. Grinstein, Annu. Rev. Pathol.: Mech. Dis., 2012, 7, 61–98.
4. I. Ellinger and P. Pietschmann, Wien. Med. Wochenschr., 2016, 166, 193–195.
5. S. Mylvaganam, S. A. Freeman and S. Grinstein, Curr. Biol., 2021, 31, R619–R632.
6. F. Niedergang and P. Chavrier, Curr. Opin. Cell Biol., 2004, 16, 422–428.
7. D. M. Richards and R. G. Endres, Rep. Prog. Phys., 2017, 80, 126601.
8. R. Imani, R. Dillert, D. W. Bahnemann, M. Pazoki, T. Apih, V. Kononenko, N. Repar, V. Kralj-Iglič, G. Boschloo, D. Drobne, T. Edvinsson and A. Iglič, Small, 2017, 13, 1700349.
9. W. Xu and X. Liu, J. Appl. Phys., 2022, 132, 174702.
10. H. Gao, W. Shi and L. B. Freund, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 9469–9474.
11. D. M. Richards and R. G. Endres, Biophys. J., 2014, 107, 1542–1553.
12. D. M. Richards and R. G. Endres, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 6113–6118.
13. F. Frey, F. Ziebert and U. S. Schwarz, Phys. Rev. Lett., 2019, 122, 088102.
14. J. S. van Zon, G. Tzircotis, E. Caron and M. Howard, Mol. Syst. Biol., 2009, 5, 298.
15. A. Khosravanizadeh, P. Sens and F. Mohammad-Rafiee, Soft Matter, 2019, 15, 7490–7500.
16. M. Fošnarič, A. Iglič, D. M. Kroll and S. May, J. Chem. Phys., 2009, 131, 105103.
17. M. Akamatsu, R. Vasan, D. Serwas, M. A. Ferrin, P. Rangamani and D. G. Drubin, eLife, 2020, 9, e49840.
18. M. Herant, V. Heinrich and M. Dembo, J. Cell Sci., 2006, 119, 1903–1913.
19. M. Herant, C.-Y. Lee, M. Dembo and V. Heinrich, PLoS Comput. Biol., 2011, 7, 1–6.
20. V. Heinrich, Biophys. J., 2015, 109, 469–476.
21. E. A. Francis and V. Heinrich, PLoS Comput. Biol., 2022, 18, 1–26.
22. A. H. Bahrami, Soft Matter, 2013, 9, 8642–8646.
23. A. H. Bahrami, M. Raatz, J. Agudo-Canalejo, R. Michel, E. M. Curtis, C. K. Hall, M. Gradzielski, R. Lipowsky and T. R. Weikl, Adv. Colloid Interface Sci., 2014, 208, 214–224.
24. S. Dmitrieff and F. Nédélec, PLoS Comput. Biol., 2015, 11, 1–15.
25. J. Liu, M. Kaksonen, D. G. Drubin and G. Oster, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 10277–10282.
26. S. Tollis, A. E. Dart, G. Tzircotis and R. G. Endres, BMC Syst. Biol., 2010, 4, 149.
27. M. Herant, V. Heinrich and M. Dembo, J. Cell Sci., 2005, 118, 1789–1797.
28. J. A. Champion and S. Mitragotri, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 4930–4934.
29. H. T. Spanke, R. W. Style, C. Fran çois Martin, M. Feofilova, M. Eisentraut, H. Kress, J. Agudo-Canalejo and E. R. Dufresne, Phys. Rev. Lett., 2020, 125, 198102.
30. T. Eierhoff, B. Bastian, R. Thuenauer, J. Madl, A. Audfray, S. Aigal, S. Juillot, G. E. Rydell, S. Müller, S. de Bentzmann, A. Imberty, C. Fleck and W. Römer, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 12895–12900.
31. J. C. Patel and J. E. Galán, Curr. Opin. Microbiol., 2005, 8, 10–15.
32. K. Rottner, T. E. Stradal and J. Wehland, Dev. Cell, 2005, 9, 3–17.
33. V. Heinrich and W. Rawicz, Langmuir, 2005, 21, 1962–1971.
34. S. R. Barger, N. S. Reilly, M. S. Shutova, Q. Li, P. Maiuri, J. M. Heddleston, M. S. Mooseker, R. A. Flavell, T. Svitkina and P. W. Oakes, et al. , Nat. Commun., 2019, 10, 1–18.
35. D. Vorselen, S. R. Barger, Y. Wang, W. Cai, J. A. Theriot, N. C. Gauthier and M. Krendel, eLife, 2021, 10, e68627.
36. P. P. Ostrowski, S. A. Freeman, G. Fairn and S. Grinstein, Dev. Cell, 2019, 50, 397–410.e3.
37. D. Paul, S. Achouri, Y.-Z. Yoon, J. Herre, C. E. Bryant and P. Cicuta, Biophys. J., 2013, 105, 1143–1150.
38. G. Sharma, D. T. Valenta, Y. Altman, S. Harvey, H. Xie, S. Mitragotri and J. W. Smith, J. Controlled Release, 2010, 147, 408–412.
39. N. Doshi and S. Mitragotri, PLoS One, 2010, 5, 1–6.
40. Z. Lu, Y. Qiao, X. T. Zheng, M. B. Chan-Park and C. M. Li, Med. Chem. Commun., 2010, 1, 84–86.
41. S. Dasgupta, T. Auth and G. Gompper, Soft Matter, 2013, 9, 5473–5482.
42. S. Dasgupta, T. Auth and G. Gompper, Nano Lett., 2014, 14, 687–693.
43. R. Vácha, F. J. Martinez-Veracoechea and D. Frenkel, Nano Lett., 2011, 11, 5391–5395.
44. N. Ramakrishnan, P. Sunil Kumar and R. Radhakrishnan, Phys. Rep., 2014, 543, 1–60.
45. M. Fošnarič, S. Penič, A. Iglič, V. Kralj-Iglič, M. Drab and N. S. Gov, Soft Matter, 2019, 15, 5319–5330.
46. R. K. Sadhu, S. Penič, A. Iglič and N. S. Gov, Eur. Phys. J. Plus, 2021, 136, 495.
47. N. Ramakrishnan, P. Sunil Kumar and J. H. Ipsen, Biophys. J., 2013, 104, 1018–1028.
48. H. Alimohamadi and P. Rangamani, Biomolecules, 2018, 8, 1–25.
49. L. Mesarec, M. Drab, S. Penič, V. Kralj-Iglič and A. Iglič, Int. J. Mol. Sci., 2021, 22, 2348.
50. M. Clarke, U. Engel, J. Giorgione, A. Müller-Taubenberger, J. Prassler, D. Veltman and G. Gerisch, BMC Biol., 2010, 8, 154.
51. M. Chabanon, J. C. Stachowiak and P. Rangamani, Wiley Interdiscip. Rev.: Syst. Biol. Med., 2017, 9, e1386.
52. J. Linkner, G. Witte, H. Zhao, A. Junemann, B. Nordholz, P. Runge-Wollmann, P. Lappalainen and J. Faix, J. Cell Sci., 2014, 127, 1279–1292.
53. I. Begemann, T. Saha, L. Lamparter, I. Rathmann, D. Grill, L. Golbach, C. Rasch, U. Keller, B. Trappmann and M. Matis, et al. , Nat. Phys., 2019, 15, 848–857.
54. A. Pipathsouk, R. M. Brunetti, J. P. Town, B. R. Graziano, A. Breuer, P. A. Pellett, K. Marchuk, N.-H. T. Tran, M. F. Krummel and D. Stamou, et al. , J. Cell Biol., 2021, 220, e202003086.
55. T. A. Masters, B. Pontes, V. Viasnoff, Y. Li and N. C. Gauthier, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 11875–11880.
56. M. Drab, R. K. Sadhu, Y. Ravid, A. Iglič, V. Kralj-Iglič and N. S. Gov, Plasma Membrane Shaping, Academic Press, 2023, pp. 415–429.
57. N. S. Gov, V. Kralj-Iglič, R. K. Sadhu, L. Mesarec and A. Iglič, Plasma Membrane Shaping, Academic Press, 2023, pp. 393–413.
58. N. Ramakrishnan, J. H. Ipsen and P. B. S. Kumar, Soft Matter, 2012, 8, 3058–3061.
59. S. Penič, A. Iglič, I. Bivas and M. Fošnarič, Soft Matter, 2015, 11, 5004–5009.
60. E. Gongadze, L. Mesarec, S. Kralj, V. Kralj-Iglič and A. Iglič, Membranes, 2021, 11, 1–20.
61. R. K. Sadhu and S. Chatterjee, Phys. Rev. E, 2016, 93, 062414.
62. R. K. Sadhu and S. Chatterjee, Phys. Rev. E, 2018, 97, 032408.
63. R. K. Sadhu and S. Chatterjee, Eur. Phys. J. E: Soft Matter Biol. Phys., 2019, 42, 15.
64. R. K. Sadhu and S. Chatterjee, Phys. Rev. E, 2019, 100, 020401.
65. J. Weichsel and P. L. Geissler, PLoS Comput. Biol., 2016, 12, 1–13.
66. A. Zak, S. Dupré-Crochet, E. Hudik, A. Babataheri, A. I. Barakat, O. Nüsse and J. Husson, Biophys. J., 2022, 121, 1381–1394.
67. M. Deserno, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 69, 031903.
68. E. A. Evans and R. Skalak, Mechanics and Thermodynamics of Biomembranes, CRC Press, West Palm Beach, United States, 1980.
69. S. Ehrig, B. Schamberger, C. M. Bidan, A. West, C. Jacobi, K. Lam, P. Kollmannsberger, A. Petersen, P. Tomancak, K. Kommareddy, F. D. Fischer, P. Fratzl and J. W. C. Dunlop, Sci. Adv., 2019, 5, eaav9394.
70. P. Sens and J. Plastino, J. Phys.: Condens. Matter, 2015, 27, 273103.
71. R. Lipowsky and H.-G. Döbereiner, Europhys. Lett., 1998, 43, 219–225.
72. A. Khosravanizadeh, P. Sens and F. Mohammad-Rafiee, J. R. Soc., Interface, 2022, 19, 20220462.
73. J. Mölller, T. Luehmann, H. Hall and V. Vogel, Nano Lett., 2012, 12, 2901–2905.
74. A. Prashar, S. Bhatia, D. Gigliozzi, T. Martin, C. Duncan, C. Guyard and M. R. Terebiznik, J. Cell Biol., 2013, 203, 1081–1097.
75. R. K. Tsai and D. E. Discher, J. Cell Biol., 2008, 180, 989–1003.
76. N. Araki, T. Hatae, A. Furukawa and J. A. Swanson, J. Cell Sci., 2003, 116, 247–257.
77. S. R. Barger, N. C. Gauthier and M. Krendel, Trends Cell Biol., 2020, 30, 157–167.
78. M. J. Sánchez-Barrena, Y. Vallis, M. R. Clatworthy, G. J. Doherty, D. B. Veprintsev, P. R. Evans and H. T. McMahon, PLoS One, 2012, 7, 1–15.
79. K. Hanawa-Suetsugu, Y. Itoh, M. Ab Fatah, T. Nishimura, K. Takemura, K. Takeshita, S. Kubota, N. Miyazaki, W. N. I. Wan Mohamad Noor, T. Inaba, N. T. H. Nguyen, S. Hamada-Nakahara, K. Oono-Yakura, M. Tachikawa, K. Iwasaki, D. Kohda, M. Yamamoto, A. Kitao, A. Shimada and S. Suetsugu, Nat. Commun., 2019, 10, 4763.
80. M. Bohdanowicz and S. Grinstein, Physiol. Rev., 2013, 93, 69–106.

---

Footnote

† Electronic supplementary information (ESI) available: High resolution supplementary movies. See DOI: https://doi.org/10.1039/d2sm01152b

---