Cell-scale dynamic modeling of membrane interactions with arbitrarily shaped particles

Abstract

Modeling membrane interactions with arbitrarily shaped colloidal particles, such as environmental micro- and nanoplastics, at the cell scale remains particularly challenging, owing to the complexity of particle geometries and the need to resolve fully coupled translational and rotational dynamics. Here, we present a force-based computational framework capable of capturing dynamic interactions between deformable lipid vesicles and rigid particles of irregular shapes. Both vesicle and particle surfaces are represented using triangulated meshes, and Langevin dynamics resolves membrane deformation alongside rigid-body particle motion. Adhesive interactions between the particle and membrane surfaces are modeled using two numerical schemes: a vertex-to-vertex mapping and a vertex-to-surface projection. The latter yields more accurate wrapping energetics, as demonstrated by benchmark comparisons against ideal spheres. The dynamic simulations reveal that lower particle-to-vesicle mass ratios facilitate frequent particle reorientation and complete membrane wrapping, while higher mass ratios limit orientation changes and stabilize partial wrapping. To illustrate the framework's versatility, we simulate interactions involving cubical, rod-like, bowl-shaped, and tetrahedral particles with spherical, cigar-shaped, or biconcave vesicles. This generalizable modeling approach enables predictive, cell-scale studies of membrane–particle interactions across a wide range of geometries, with applications in environmental biophysics and nanomedicine.

---

1. Introduction

Marine plastic waste fragments into micro- and nanoplastics (MNPs), broadly defined as particles with diameters ≤ 5 mm. Once released into the environment, these particles can enter a variety of cell types through endocytosis, macropinocytosis, and phagocytosis, potentially triggering harmful biological responses.^1,2 Among these pathways, phagocytosis plays a central role in the uptake of larger particles, where the plasma membrane actively deforms to engulf extracellular materials into an intracellular vesicle.^1–6 Unlike engineered colloids that are predominantly spherical, MNPs often display highly irregular and anisotropic shapes.^7,8 These morphological features, along with surface chemistry,^9–12 size,^13,14 and deformability,^15–17 critically influence how particles interact with and are internalized by cellular membranes.^18,19

Anisotropic particle shapes introduce curvature heterogeneity, resulting in membrane wrapping behaviors that differ fundamentally from those of spherical particles.^13,19–28 Computational and theoretical studies have demonstrated that for non-spherical particles, the local variation in curvature leads to spatially heterogeneous bending energy costs. This heterogeneity can stabilize partially wrapped states and significantly influence vesicle morphology.^{20,21,23,24,29} For example, Dasgupta et al.²³ reported that partially wrapped ellipsoids can remain stable because the membrane resists bending around the high-curvature tip of the particle. Other simulations have reported a two-stage sequence for ellipsoid particles: an initial flat-side adhesion, followed by a tip-first reorientation that lowers the overall bending penalty and enables complete engulfment. In such cases, the particle's large aspect ratio amplifies the energy barrier to reaching complete wrapping.^20,22,25,30

Experimental observations support these computational findings. A recent study by van der Ham et al.²⁶ confirmed orientation-dependent wrapping behavior, showing that a long, rod-shaped particle initially adopts a “surfing” state along the surface of a giant unilamellar vesicle (GUV) before becoming engulfed. During this process, the particle undergoes further reorientations prior to complete uptake. Similarly, Azadbakht et al.³¹ showed that variations in the neck curvature of dumbbell-shaped particles modulate both the speed and completeness of engulfment by GUVs. In particular, the constricted neck region imposes a kinetic barrier, rendering wrapping highly sensitive to membrane tension. Together, these studies underscore that shape anisotropy not only introduces local curvature variation but also governs the sequence of particle reorientation events necessary for complete membrane envelopment.^{22,25,26,30,32} In biological cells, particle–membrane interactions can involve more complex physics. While the wrapping process during phagocytosis can be adequately predicted using passive membrane elasticity and adhesion energetics,^{20,21,23,24,33,34} living cells frequently employ active processes to drive membrane deformation and facilitate the uptake of exogenous substances beyond what passive interactions alone can achieve. Protrusive forces from actin polymerization and contractile forces generated by myosin motors actively remodel the membrane to promote particle internalization.^4,35–40 Furthermore, membrane proteins, such as BAR-domain scaffolds and curvature-sensitive complexes, facilitate uptake dynamics by locally inducing curvature and recruiting cytoskeletal elements.^41–47 While recognizing the biological importance of these active membrane-deformation mechanisms, our study intentionally employs a simplified modeling framework that excludes active forces and cytoskeletal coupling. This passive model enables precise quantitative comparisons with controlled experiments involving synthetic vesicles or cells lacking active endocytic machinery, thus providing a foundational understanding that can guide future incorporation of active mechanisms.

The choice of numerical method for passive wrapping simulations depends on the specific objectives of the study and the particular aspects of membrane dynamics of interest (see Table S1). Energy-based simulations using Helfrich bending theory and numerical methods such as Monte Carlo predict wrapping phase diagrams and equilibrium vesicle shapes by minimizing membrane bending and adhesion energies.^4,24,29,45 However, these thermodynamic approaches fall short of revealing the kinetic pathway of wrapping. Specifically, they yield only equilibrium configurations and do not resolve the time-dependent dynamics of particle–vesicle interactions. This limitation is particularly significant for anisotropic and irregular shapes, whose orientation and dynamic responses to membrane interactions critically influence cellular uptake, as emphasized by experimental studies. Addressing these kinetic pathways thus requires a dynamic simulation framework capable of computing membrane-mediated forces and iteratively updating the translation and rotation of irregularly shaped particles.

We extend our previously developed force-based model⁴⁸ to simulate the interaction dynamics of vesicles with anisotropic particles, including those lacking analytical shape representations. Within a unified modeling framework, both deformable membranes and arbitrarily shaped rigid particles are represented using triangulated surface meshes. Two distinct numerical schemes are implemented to model adhesion between particle and membrane surfaces, and their accuracy in predicting wrapping energetics is systematically compared. Furthermore, we characterize the role of particle inertia in translation and rotational dynamics, explicitly capturing the force–torque coupling mediated by the membrane. Finally, we demonstrate the versatility of our approach by simulating interactions between non-spherical vesicles and particles with complex geometric features, including significant curvature variations and convex–concave transitions. This framework offers a versatile tool for exploring membrane–particle interactions beyond idealized geometries, supporting future research in both environmental biophysics and nanomedicine, particularly regarding shape-dependent uptake mechanisms at the cellular scale.

2. Methods

2.1. Triangulated surface model for vesicle

The process by which a vesicle membrane wraps around a particle is governed by the competition between bending energy E_bend and adhesion energy E_adh. The bending energy captures the cost of deforming the membrane, which is described by the Canham–Helfrich (CH) energy functional. This model is integral to predicting membrane deformation under various physical conditions and is described as follows^49,50Here, H, G, and H₀ represent the mean, Gaussian, and spontaneous curvatures of the membrane, respectively. Bending moduli κ_b and Gaussian curvature moduli κ_g characterize the strength of the membrane against different deformation modes. This study does not consider vesicle topology changes and thus the Gaussian curvature term is neglected according to the Gauss–Bonnet theorem.

In addition to the bending energy, a vesicle model must include area and volume constraints to maintain its mechanical stability and realistic membrane behavior. The area constraint energy ensures that the vesicle's surface area remains constant, reflecting the lateral incompressibility of the lipid bilayer (the per area lipid varies little without extreme tension).⁵¹ Similarly, the volume constraint energy penalizes variations in internal fluid volume, thereby maintaining the osmotic pressure difference between the vesicle and its surroundings.³⁴ These constraints work together to control the vesicle shape under physiologically relevant conditions, such as variations in membrane tension or osmotic concentrations. We employed simple quadratic forms for these two energies, given by

Here, κ_a is the area expansion modulus, while A_t and A denote the preferred and current surface areas of the vesicle membrane, respectively. In a similar fashion, V_t and V represent the target and current vesicle volumes, respectively. κ_v acts as a modulus that controls the strength of osmotic pressure imbalance when the vesicle is deformed from its natural shape by particle interactions.³⁴

The membrane shape corresponding to the lowest energy state can be theoretically derived by varying the CH energy functional. However, solving the resulting shape equation is often intractable due to complex geometric attributes. To overcome this, the membrane is modeled as a discretized two-dimensional surface using a triangulated mesh, which enables the application of discrete differential geometry to evaluate the membrane energy functional. Specifically, the membrane surface S is represented by a mesh network composed of N_v vertices, N_e edges, and N_f triangular faces. Each vertex v_i corresponds to a point on the surface embedded in a three-dimensional laboratory frame. All vertices collectively define the global shape of the membrane. A triplet of connected vertices defines each triangular face. The numerical derivations and formulations of geometric variables are detailed in the following references.^49,52,53 Details of vesicle mesh generation and energy calculations for the discretized membrane can be found in our previous work.⁴⁸

2.2. Triangulated surface representation of arbitrarily shaped particles

Analytical representations are often inadequate for particles with complex geometries, such as those exhibiting widely varying curvature, concave regions, and sharp edges or corners. Such irregular features typically cannot be captured by simple closed-form expressions.⁵⁴ Due to these limitations, triangulated surface meshes (similar to the vesicle surface) provide a flexible solution. By approximating a surface with numerous small facets, a mesh can closely fit virtually any shape and even accommodate complex topographies with intricate surface curvature variations.⁵⁵ Five particles with diverse shapes are simulated using triangulated meshes in this study. Spherical particles are generated by three-level icosahedron subdivision to achieve near-uniform triangulation. To enable direct comparisons with previous simulations and experiments, cube- and rod-like particles are generated using a MATLAB mesh generator⁵⁶ according to the equations of superellipsoids,^21,24x⁶ + y⁶ + z⁶ = R_p⁶ with R_p = 0.2 and [(x² + y²)/a²]^n/2 + (z/b)ⁿ = 1, with n = 6, a = 0.15, and b = 0.3, respectively. We also simulate cubical particles with rounded corners through a Minkowski sum of a cube and a sphere with a cube edge length of 0.6 and a rounded edge/corner radius of 0.2. The mesh for a bowl-shaped particle is created based on an initial oblate spheroid mesh with a major axis of 0.33 and a minor axis of 0.24. We carve out a smooth concave “bowl” region on its lower hemisphere by displacing the z positions of the vertices with z ≤ 0 following . Here, D_max is the dimple depth, which is set to 0.27, and R_bowl = 0.33. For a tetrahedral particle, we generate an ideal tetrahedral mesh using an edge length of 0.6 and then smooth the corners and edges using the libigl mesh library.⁵⁷

2.3. Mesh-based schemes for defining membrane–particle interaction

The membrane–particle interaction energy can be written as the surface integralHere, d(x) is the separation distance between the membrane at point x and the particle surface. V(d) is the interaction potential density acting between the two interacting surfaces. As the membrane surface is discretized into triangular faces, the total energy is computed by summing the contributions from the interacting vertices, weighted by their associated Voronoi areas A_i^33,58Here, d_i is the distance from the membrane vertex v_i to the particle surface. The adhesion between membrane and particle is modeled using Morse potential, taking the following form⁵⁹

V(d_i) = U[e^(−2d_i/ρ) − 2e^(−d_i/ρ)] (6)

In this expression, non-negative U is the adhesion strength, and positive ρ defines the interaction range. This potential has a minimum value of −U at d_i = 0, indicating attraction that lowers the total free energy (see Fig. 1). To reduce the computational cost, we introduce a numerical cutoff r_c (on the order of potential range) and include only those membrane vertices for which d_i < r_c in the calculation. We also define a dimensionless adhesion strength u = UR_p²/κ_b for spherical particles, where R_p is the radius. In the case of a non-spherical particle, we replace R_p with an effective radius () obtained from the total particle mesh surface area A_p. This R_eff is then used to rescale the adhesion strength, yielding u_mod = UR_eff²/κ_b for arbitrarily shaped particles.

Fig. 1 Nanoparticle–vesicle interaction potential. Morse potential profiles show how the interaction range is controlled by ρ. The inset schematically depicts a nanoparticle (yellow) approaching the vesicle membrane (gray) with the interaction distance d_i defined between the particle and membrane surfaces.

We determine d_i for an arbitrarily shaped particle represented by a triangulated mesh using two alternative geometric approaches. In the vertex-to-vertex scheme, each membrane vertex is paired with the nearest vertex on the particle mesh (a nearest-neighbor mapping) to estimate d_i. In the vertex-to-surface scheme, d_i is computed as the shortest Euclidean distance from the membrane vertex to any triangular face of the particle mesh. While the two approaches define the local membrane–particle separation differently for complex particle shapes, they allow us to impose adhesion between the particle and membrane, which drives wrapping and uptake. The implementation details of both methods are described as follows.

Vertex-to-vertex adhesion scheme. To compute the adhesion energy between the vesicle and particle surface, the Morse potential is applied based on the pairwise distance between the vertices of the two triangulated surfaces. This approach, referred to later as the vertex-to-vertex scheme, draws inspiration from prior studies on vesicle–vesicle aggregation in fluid flows.⁶⁰ This method establishes virtual bonds between nearest neighboring pairs of vesicle and membrane vertices (see Fig. 2), maintaining a one-to-one correspondence. The resulting bond distances are then used to calculate the adhesion energy and force. Notably, we ensure that bonds are formed only between mutually nearest vertex pairs, thereby avoiding redundant or overlapping interactions. The detailed algorithm for the vertex-to-vertex model is presented in the Appendix.

Fig. 2 Vertex-to-vertex interaction scheme. The schematic illustrates the vertex-to-vertex scheme for calculating the interaction between a triangulated rigid particle (orange mesh) and a biconcave vesicle (gray mesh). Each yellow arrow represents a one-to-one “bond” formed between a membrane vertex and its nearest neighbor on the particle mesh. For visual clarity, the vesicle and particle are not shown at their relative physical sizes.

Vertex-to-surface adhesion scheme. The vertex-to-surface scheme computes the shortest Euclidean distance from a given membrane vertex to the triangulated surface of the particle.^61,62 This approach measures the true minimal separation between two triangulated surfaces, independent of mesh resolution or vertex alignment. It ensures that even if the closest approach occurs between a vertex and a face (or edge) of the other mesh (rather than between two vertices), the distance is accurately captured.

Specifically, we describe how to find each vertex's closest point on the target surface to obtain the unsigned distance.⁵⁷ Consider a vesicle vertex v_i at the position P and a closed triangulated surface S¹ (see Fig. 3), P is orthogonally projected onto the plane of a triangle. If the projection point Q resides within the triangle (determined via barycentric coordinates or a cross-product test), then the normal distance |P–Q| is a candidate. If Q lies outside the triangle, we instead project P onto each of the three edges of the triangle, clamping the projection to the segment endpoints. Thus, the closest point is either this edge projection or the nearer vertex. The Euclidean distance from P to each of these boundary candidates for a triangle is computed and the smallest is retained. After evaluating all triangles, the global minimum of these candidate distances yields the closest point C on S¹ and the corresponding triangle, and the unsigned distance is |P–C|.

Fig. 3 Vertex-to-surface interaction scheme. For each membrane vertex P (yellow sphere), the closest point C (red sphere) on the particle mesh is identified by minimizing |P–C|. The interaction force is calculated based on the shortest distance |P–C|, acting along the line connecting P and C.

As the negative distance branch of the Morse potential is important for avoiding the particle unphysically penetrating the membrane, we define the sign of the distance between a vesicle vertex v_i and the particle surface S¹. This is achieved by taking the dot product of an angle-weighted pseudomonal⁶³n_α with the vector from the surface to the point, P−C. The sign of the dot product indicates the positional relationship: n_α(P−C) > 0 if P lies outside the surface S¹; n_α(P−C) < 0 if P lies inside S¹; n_α(P−C) = 0 if P lies exactly on S¹. This criterion works regardless of whether the closest point C is on a face, an edge, or coincides with a vertex of the mesh,⁶³ reliably determining the sign of the distance even at sharp features where a normal is not uniquely defined.

2.4. Time integration scheme

To simulate the motion of the vesicle membrane, the velocity-Verlet algorithm has been chosen as the time integration scheme, consistent with approaches used in molecular dynamics packages such as LAMMPS.⁶⁴ The motion of each vertex on the triangulated vesicle surface is governed by a deterministic form of the Langevin equation:^65,66where m is the vertex mass. The total deterministic force F^tot_i includes contributions from bending elasticity, area and volume constraints, and membrane-particle adhesion, with the detailed calculations described in ref. 48. The damping term with the drag coefficient γ serves two primary purposes: it ensures the numerical stability of the integrator and approximates viscous dissipation from the surrounding fluid. This approach is commonly adopted in simulations of particle–membrane interactions at the cellular scale.^67,68

It is well documented that hydrodynamic drag on particles increases with confinement, such as proximity to a membrane or solid boundary. While Stokes’ law describes constant drag in unbounded fluids, the presence of nearby surfaces results in higher and position-dependent drag.^69,70 However, such variation primarily affects the absolute time scale of particle motion, rather than the mechanistic sequence of orientational transitions or the equilibrium wrapping states.⁶⁷ Accounting for position-dependent drag induced by the presence of the membrane would require either explicit solvent modeling (e.g., using dissipative particle dynamics or lattice Boltzmann methods) or analytical treatment of hydrodynamic interactions.⁷¹ These approaches significantly increase computational cost and are beyond the scope of the present work. Moreover, many mesoscopic simulations of nanoparticle wrapping studies have similarly employed a constant damping coefficient.^65,66 Therefore, we believe that the assumption of a constant friction coefficient is justified for capturing the role of particle inertia and adhesion strength in wrapping dynamics. We set the vertex mass and drag coefficient (m = γ = 1 in the simulation units) to model an intermediate damping regime.

Notably, this fully deterministic Langevin formulation⁶⁷ omits stochastic thermal noise, solvent hydrodynamics (beyond simple drag), and active cellular factors such as cortical tension and receptor–cytoskeleton coupling. This simplification allows the adhesion–bending energy balance to be treated in a computational and analytically tractable manner. However, we acknowledge that this approach has limited quantitative predictive power at the nanoscale, where thermal fluctuations on the order of k_BT (with k_B being the Boltzmann constant), hydrodynamic dissipation, and active forces can facilitate energy barrier crossing and influence uptake kinetics.

To explore the coupled dynamics of the vesicle–particle system, the translational and rotational motion of nanoparticle is explicitly modeled using rigid body dynamics. A similar Langevin dynamics also governs the translation motion of the particle's center of mass, consistent with the vesicle vertex dynamics. The equation is given by

Here, m_p is the total mass of the rigid particle, r_com is the particle's center of mass position, and γ_p is the friction coefficient accounting for translational damping. F^p_j is the reaction force to membrane adhesion acting on the particle surface. For each membrane vertex v_i experiencing a force F^adh_i, we apply F^p_j = −F^adh_i at the corresponding point j on the particle, which is either the nearest mesh vertex (Fig. 2) in the vertex-to-vertex scheme or the closest surface point C (Fig. 3) in the vertex-to-surface scheme. Thus, gives the net force driving particle translation, while the drag term damps its motion. We fix m_p/γ_p = 1 to ensure a particle damping regime consistent with the membrane vertices.

The rotational dynamics is governed by the total torque generated from off-center reactions applied to the particle surface. The torque about the particle's center of mass is calculated as

In addition to the torque arising from the reaction force to membrane adhesion, we incorporate an effective rotational drag torque to model viscous dissipation acting on the particle. The rotational drag torque is proportional to the particle's angular velocity in the body frame and is scaled by its principal moments of inertia.⁶⁴ The drag torque is calculated in the body frame as

τ^b_d = −[small gamma, Greek, tilde]_rot[I_xxω^b_x, I_yyω^b_y, I_zzω^b_z] (10)

Here, [small gamma, Greek, tilde]_rot is the rotational damping coefficient. I_xx, I_yy, and I_zz are the principal moments of the inertia tensor I, while ω^b_x, ω^b_y, and ω^b_z are the components of the angular velocity vector in the body frame. This body-frame drag torque is then transformed into the laboratory frame using the rotation matrix R associated with the particle's current orientation

τ^s_d = Rτ^b_d (11)

This laboratory-frame drag torque is added to the total torque balance when updating the particle's rotational dynamics, providing essential rotational damping similar to linear damping applied to the particle's center of mass and vesicle vertices. Such damping prevents any artificial or ongoing spinning that may arise from numerical noise or integration drift and is crucial for achieving physically meaningful equilibrium states.

Subsequently, the angular momentum L and angular velocity ω are calculated from the total applied drag (τ_tot = τ − τ^s_d) on the particle over a time step Δt

L = τ_totΔt (12)

ω = I⁻¹L (13)

The orientation is tracked using a quaternion q = q₀ + q₁i + q₂j + q₃k, which is updated through

where ⊗ denotes the quaternion product. The updated quaternion is then normalized, followed by the reconstruction of the body frame to obtain new positions of all particle vertices. This quaternion-based formulation avoids singularity in rotation kinematics and ensures stable, smooth integration of rotational motion.^64,72,73 The detailed simulation parameters are given in the following table.

A mapping between the simulation units and physical units can be established by referencing comparable GUV experiments. The characteristic length scale is set by the vesicle radius l₀ = 10 μm, reflecting the typical GUV size. The characteristic time scale is τ = 0.05 s, based on experimentally observed uptake time.²⁶ The energy unit E₀ is defined by the bending modulus κ_b of typical phospholipid membranes, as bending energy represents the dominant energetic contribution during vesicle deformation. With this unit system, all other simulation parameters, including the area expansion modulus, the volume constraint modulus, membrane and particle masses, damping coefficients, and interaction potentials, can be consistently converted to physical units via dimensional analysis (Table 1). Notably, the resulting membrane vertex drag coefficient closely matches the estimate given by Stokes’ law for water.

Table 1 Standard simulation parameters

Parameters	Numerical values	Physical values
Bending modulus (κb)	0.01^52	20kBT^74,75
Area expansion modulus (κa)	1.0	8.3 × 10^−5 nN μm^−1
Volume constraint modulus (κv)	2.0 (for controlled-volume)^52	1.6 × 10^−5 nN μm^−2
	0.0 (for uncontrolled-volume)
Spontaneous curvature (H0)	0	0 μm^−1
Vesicle surface area (At)	4π	1.3 × 10^3 μm^2
Morse potential range (ρ)	0.01	100 μm
Simulation time step (Δt)	0.001	5 × 10^−5 s
Vesicle mass (mv)	10 242	2 × 10^−10 kg
Membrane vertex drag coefficient (γ)	1.0	4.1 × 10^−6 nN s μm^−1
Particle mass (mp)	1.0	2 × 10^−14 kg
Particle translational drag coefficient (γp)	1.0	4.1 × 10^−6 nN s μm^−1
Particle rotational damping coefficient ([small gamma, Greek, tilde]rot)	1.0	20 s^−1

---

3. Results and discussion

3.1. Comparison of interaction schemes for meshed particles

To evaluate how different numerical schemes affect the particle–membrane interaction, we characterize the energy profile of a spherical vesicle wrapping an external or internal nanosphere of radius R_p = 0.3 with reduced adhesion strength u = 2. Both the particle and vesicle meshes are generated by successive subdivisions of an icosahedron, yielding 2562 vesicle vertices (N_v) and 602 particle vertices (N_p). In these simulations, the particle is held stationary, and the wrapping energy is quantified when a steady-state configuration is reached. A harmonic umbrella potential is employed to sample the states with high wrapping fractions, which are not naturally attainable.^20,48 The wrapping fraction χ is calculated from the ratio of current adhesion energy to the theoretical value corresponding to full wrapping, χ = E_adh/UA_p. The results of the meshed sphere using the vertex-to-vertex and vertex-to-surface schemes are compared against the data in our previous study by modeling a parametric sphere of the same size.⁴⁸ Notably, no volume constraint is applied for consistency.

As shown in Fig. 4, the vertex-to-vertex scheme clearly overestimates both bending and total energies compared to the energetics predicted by the parametric sphere. This significant discrepancy is attributed to the calculation of surface distance. As the two triangulated surfaces approach each other, the vertex-to-vertex bonds shorten. However, as the membrane mesh locally conforms to the particle mesh, the computed bond lengths cannot reduce to zero even when membrane vertices are brought onto the triangular faces of the particle mesh. This residual distance results in the evaluation of the Morse potential and adhesion forces at artificially larger distances, leading to inaccurate energy profiles. This effect is more pronounced when there is a high local curvature generated at high χ. Regardless of whether particles are inside or outside the vesicle, the total energy sharply increases when the wrapping fraction exceeds 0.6. Neck formation at this high wrapping fraction introduces very large local curvature at the membrane–particle interface, causing the vertex-to-vertex distance to deviate increasingly from actual surface distance. The exact threshold at which this diverging behavior occurs depends on the resolution of the membrane and particle meshes. Finer meshes may extend this limit somewhat, but the qualitative constraint remains.

Fig. 4 Quantitative comparison between the two interaction schemes for triangulated rigid particles. Equilibrium (a) bending and (b) total energies of a spherical vesicle interacting with a spherical particle of dimensionless radius of 0.3 and interaction potential range ρ = 0.01, plotted against the wrapping fraction χ of the particle. Different branches correspond to particles positioned inside and outside the vesicle. Dashed lines represent results for an ideal sphere from ref. 48, while square and circular symbols denote simulations using triangulated spheres with vertex-to-vertex and vertex-to-surface interaction schemes, respectively.

In contrast, the vertex-to-surface scheme renders accurate adhesion forces, even at higher wrapping fractions where local membrane curvature also varies largely. The predicted energy curves closely align with those of the parametric sphere. Such an accurate evaluation of surface distance is essential for capturing contact mechanics and wrapping dynamics in our simulations. Due to this markedly improved agreement, the vertex-to-surface scheme is adopted for all subsequent simulations.

To further benchmark our numerical framework, we conducted simulations replicating the system in Yu et al.,²⁴ who reported normalized bending energy as a function of the wrapping fraction for a cube-like particle of size R_p = 0.2. In Fig. S1, we compare their results with those from our current model for two distinct cases: (i) a stationary particle with all translational and rotational motion suppressed, and (ii) a freely moving particle with both degrees of freedom enabled. Our model shows good agreement with Yu et al.'s data at low and high wrapping fractions. However, for intermediate wrapping states (χ ∼ 0.3–0.7), our simulations predict a lower equilibrium bending energy.

This discrepancy arises from fundamental differences in simulation methodology. In Yu et al.'s study, the membrane region adhering to the particle was held fixed, and thus the vesicle shape evolved with a fixed contact line. In contrast, our approach dynamically evolves the membrane and particle meshes independently, without constraining the adhesion interface/boundary. In other words, given a certain wrapping fraction, the membrane in our model can adaptively cover different regions of particles of the same area. This relative mobility between the meshes allows the system to explore energetically more efficient wrapping pathways, which contributes to the lower bending energies observed in our simulations at intermediate wrapping degrees.

3.2. Effect of particle inertia on interaction dynamics

In our Langevin-dynamics framework, the particle-to-vesicle mass ratio m_p/m_v (where m_v = mN_v is the total vesicle mass) governs the particle's inertia relative to the interacting membrane. This affects force–torque coupling, particularly how applied torques translate into particle rotation. Heavier particles tend to retain their initial orientation, while lighter ones respond more readily to vesicle deformation. To examine this effect, we simulate cubical particles with three mass ratios—0.0001, 0.0005, and 0.001—interacting with a biconcave vesicle with its reduced volume constrained at 0.65 (defined as the ratio of the vesicle volume to that of a sphere having the same surface area).^76,77 A rescaled adhesion strength u_mod = 4.0 is applied to enable a robust wrapping transition. The particle is initially oriented with one corner contacting the vesicle's north pole (“corner-attack” configuration).

As shown in Fig. 5, the two lighter particles (m_p/m_v = 0.0001 and 0.0005) undergo multiple reorientations during wrapping, while the heaviest particle (m_p/m_v = 0.001) retains its original orientation throughout the simulation. The lightest particle exhibits the most dynamic behavior, rapidly transitioning to a “face-attack” configuration (Fig. 5(a)). Because the top dimple of the biconcave vesicle is intrinsically shallow, this orientation smooths local curvature variations, reducing bending energy. However, as wrapping progresses, further advancement in the face-attack configuration requires bending the membrane around multiple corners, which becomes energetically unfavorable. Around t ∼ 500, the particle reverts to the corner-attack orientation, leveraging the dimple's natural concavity to reduce bending costs. Later, at t ∼ 800, it returns to the face-attack pose, enabling full wrapping of side faces, which is maintained thereafter (see SI Video 1). These rapid transitions reflect the particle's low inertia.

Fig. 5 Mass-dependent rotational dynamics of cubical particles interacting with a biconcave vesicle. The particle's instantaneous orientation is described by the rotation from the laboratory frame to the particle's body frame, defined by the principal axes of its moment of inertia tensor. This rotation is parameterized using the Euler vector θe, where the unit vector e indicates the direction of an axis of rotation and positive θ is the rotation angle following the right-hand rule. The direction of e is further expressed by its azimuthal and polar angles in the laboratory frame. Time evolutions of (a) rotation angle, (b) azimuthal angle, and (c) polar angle are shown for particles with particle-to-vesicle mass ratios m_p/m_v = 0.0001 (yellow), 0.0005 (green), and 0.001 (blue). Insets show corresponding vesicle–particle configurations at selected time points.

The intermediate-mass particle (m_p/m_v = 0.0005) initially remains in the corner-attack configuration and transitions only once to the face-attack orientation later in the simulation (Fig. 5(b)). This delayed reorientation suggests that inertia suppresses the early-stage transition observed for the lightest particle, but becomes less influential as wrapping progresses and adhesion forces dominate. In contrast, the heaviest cubical particle (m_p/m_v = 0.001) shows no significant reorientation and remains in the corner-attack configuration throughout (Fig. 5(c)). Its high inertia prevents the transition to a face-attack pose, limiting membrane wrapping. As shown in Fig. S2, the wrapping fraction for the heaviest particle plateaus around 0.55, whereas the lighter particles reach values near 0.8. This behavior aligns with recent experimental observations of poly(L-lactic-co-glycolic) acid nanoparticle uptake by macrophages, where heavier particles were internalized more slowly, with longer kinetic half-times and frequent arrest in partially wrapped states.⁷⁸ Overall, these results demonstrate how particle inertia influences rotational dynamics and wrapping efficiency through force–torque interactions with a deformable membrane. Given our focus on resolving rigid body dynamics, we use the lightest particle m_p/m_v = 0.0001 in subsequent simulations. The effects of particle and vesicle mesh densities are studied for this mass ratio. The wrapping kinetics shows convergent behavior when particle mesh density is higher than N_p = 702 (Fig. S3a). When vesicle density is changed, the overall wrapping pathway and final configuration remain the same, but the kinetics shows variations (see Fig. S3b). This difference is due to the bond density and formation sequence between the particle and membrane surfaces being altered for different mesh densities.

3.3. Cubical and rod-like particles interacting with a spherical vesicle

We examine the interaction between a cubical particle and a spherical vesicle to understand further how adhesion strength influences the final wrapping states of anisotropic particles. A series of simulations is performed with increasing adhesion strength, ranging from u_mod = 2.0 to 10.0, as shown in Fig. 6(a). At low adhesion strength (u_mod = 2.0), the particle retains its initial orientation with only weak membrane interaction and minimal wrapping. As the adhesion strength increases to intermediate values, the cubical particle reorients into the face-attack configuration, thereby increasing the effective contact area and lowering local adhesion energy. Interestingly, at higher adhesion strengths (u_mod = 8.0–10.0), the particle becomes fully wrapped by the vesicle while reverting to the corner-attack orientation. These transitions between corner-attack and face-attack configurations are consistent with the dynamics previously observed in Fig. 5.

Fig. 6 Effect of adhesion strength on anisotropic particles interacting with a spherical vesicle without the volume constraint. Final configurations with particles having different initial orientations at different rescaled adhesion strengths: (a) cubical particle with the corner contact, (b) rod-like particle with side-wise contact, and (c) rod-like particles with tip-wise contact. The side-wise rod achieves full wrapping at u_mod = 6.0, whereas both the cubical particle and tip-wise rod require u_mod = 8.0 for complete engulfment.

We extend our analysis to rod-like particles introduced in two distinct initial orientations: side-wise (lying horizontally along the vesicle surface) and tip-wise (aligned vertically). The interaction dynamics obtained from our simulations are compared with experimental observations reported by van der Ham et al.²⁶ (Fig. S4), showing excellent agreement. For side-wise particles, the membrane readily wraps upon initial contact, with the rods remaining horizontal in the partially wrapped state. As the membrane contact advances and begins to cover the edge of the rod tip, the particle gradually pivots upward along its long axis, eventually achieving full engulfment (Fig. S3a and b). This gradual reorientation is accompanied by steady increases in both membrane bending energy and wrapping fraction at u_mod ≥6.0 (Fig. S5a and c), indicating entry into the engulfment phase at u_mod = 6.0. The vesicle's reduced volume decreases concurrently (Fig. S4b), reflecting internal accommodation of the particle.

For tip-wise particles, the initial wrapping proceeds more slowly due to the small contact area and high energy cost associated with bending around the tip (Fig. S5c). Rapid engulfment begins once the particle breaks symmetry and transitions into a slanted orientation (Fig. S4c, d). Compared to side-wise particles, this slanted entry leads to sharper increases in bending energy, reduced volume, and wrapping fraction observed at u_mod ≥ 8.0.

Fig. 6(b) and (c) illustrate how adhesion strength influences the final configuration. Side-wise particles achieve full wrapping at lower adhesion strengths (u_mod = 6.0) compared to tip-wise particles, which require u_mod = 8.0. At lower adhesion strengths, particles remain partially wrapped while largely retaining their initial orientations. At high adhesion strengths (u_mod = 8.0 to 10.0), the vertical alignment of the particle and rotational symmetry of the vesicle are no longer maintained. This rapid change in curvature at the small catenoidal neck region is also reflected in the bending energy curves for high adhesion strengths. Together, these results demonstrate that the final orientation and wrapping state of anisotropic particles are governed by a combined effect of particle shape and adhesion strength.

3.4. Cubical particle interaction dynamics with a cigar-shaped vesicle

Capturing the wrapping dynamics is key to understanding particle uptake pathways in vesicles and cells. We examine the uptake of cubical particles by a cigar-shaped vesicle, focusing on morphological evolution across a range of rescaled adhesion strengths (u_mod = 2.0–10.0). The vesicle is modeled with a volume constraint corresponding to a reduced volume of 0.7,^76,77 resembling tubular endosomes and Weibel–Palade bodies—elongated organelles that support long-range cargo transport and targeted release.^79,80 Two initial contact locations are considered: the saddle-shaped waist with negative Gaussian curvature, and the highly convex north pole (Fig. 7). In both cases, the particle begins in a corner-attack orientation for consistency.

Fig. 7 Interaction dynamics of cubical particles with a cigar-shaped vesicle. Time evolution of the wrapping fraction for a cubical particle contacting the (a) saddle-shaped waist and (c) convex pole of the vesicle at different rescaled adhesion strengths. Time is normalized by the maximum duration required to reach the steady-state configuration, which depends on the rescaled adhesion strength. Filled and open star markers indicate the final states for contact at the waist and pole, respectively. Corresponding simulation snapshots are shown in panels (b) and (d).

Fig. 7 illustrates how adhesion strength and initial contact location affect uptake. Particles contacting the waist achieve full wrapping at lower adhesion strength (u_mod = 6.0) compared to pole-attached particles, which require u_mod = 10.0. This difference is attributed to the fact that the saddle-shaped waist promotes membrane bending more readily than the convex pole. For waist-attached particles, wrapping proceeds through three stages: an initial lag phase with a low wrapping fraction, a growth phase marked by reorientation to edge or face contact, and a final burst phase with rapid engulfment. These stages correspond to sharp increases in both bending energy and wrapping fraction at u_mod ≥ 6.0 (Fig. S4). The reduced volume initially decreases to accommodate the particle, then partially recovers due to the volume constraint, resulting in a capsule-like final shape (Video S4). In contrast, pole-attached particles remain partially wrapped at lower adhesion strengths, arrested in the lag phase for u_mod ≤ 4.0. At u_mod = 6.0 to 8.0, only one face attaches to the pole region (Fig. 7(d)), and full engulfment is observed only at u_mod = 10.0. During this process, the particle reorients from its initial corner-attack pose to a face-attack configuration, then returns to corner-attack after complete engulfment (Video S5). Substantial increases in bending energy are seen only at u_mod = 10.0. Correspondingly, reduced volume drops significantly without recovery, indicating a more energetically costly membrane deformation compared to waist interaction (Fig. S5). These results show that adhesion strength, initial contact geometry, and local membrane curvature together dictate the efficiency, dynamics, and morphology of particles wrapping by non-spherical vesicles.

3.5. Interaction between bowl-shaped and tetrahedral particles with a biconcave vesicle

To further test the robustness of our computational framework, we simulate the wrapping of two geometrically complex particles: a bowl-shaped particle with opposing concave and convex surfaces, and a tetrahedral particle with sharp facets and pronounced curvature variation. A biconcave vesicle is used in both cases to provide a membrane with regions of both positive and negative curvature, mimicking red blood cell geometry.

In the first case, the bowl-shaped particle is placed with its rim contacting the vesicle surface (Fig. 8(a)). Despite a high adhesion strength (u_mod = 10.0), wrapping remains shallow initially due to curvature mismatch and the associated bending energy cost. Over time, the particle shifts into the vesicle's dimple, where local curvature enables greater contact. Notably, the membrane opposite the adhering region bulges outward, forming a convex shape—consistent with previous observations of vesicle remodeling during spherical particle uptake.⁴⁸ When the particle is oriented with its concave face against the vesicle's convex waist (Fig. 8(b)), complementary curvature allows full coverage of the concave surface with minimal membrane deformation. However, further wrapping around the particle rim is hindered by the steep bending penalty required to curve around the rim, even at u_mod = 10.0.

Fig. 8 Biconcave vesicle interaction with bowl-shaped and tetrahedral particles. (a) Bowl-shaped particle (u_mod = 6.0) initially contacting the vesicle with its rim sinks into the concavity. (b) Bowl-shaped particle (u_mod = 10.0) minimally adheres to the vesicle waist. (c) Tetrahedral particle (u_mod = 10.0) rolls into the vesicle concavity. Images progress from left to right over simulation time.

The final simulation features a tetrahedral particle positioned near the flatter side of the vesicle (Fig. 8(c)). At u_mod = 10.0, full wrapping occurs spontaneously. The particle first approaches with a corner, then reorients to align a flat face for greater contact. As wrapping continues, it flips into an inverted pose that allows three faces to be enclosed, ultimately achieving full engulfment. Further rotation of the particle inside the vesicle follows. These simulations demonstrate the framework's capability to resolve complex wrapping behaviors driven by particle shape, adhesion strength, and membrane curvature.

4. Conclusion

In this work, we developed a force-based, dynamic simulation framework to resolve the translational and rotational dynamics of arbitrarily shaped particles interacting with fluidic vesicles. Both the deformable membrane and rigid particles are modeled using triangulated surface meshes. The method integrates the CH bending energy functional with a deterministic Langevin equation to describe the deformation of vesicles and the rigid-body motion of particles. Adhesive interactions between the triangulated meshes are implemented via two numerical schemes: vertex-to-vertex (nearest-neighbor bonds) and vertex-to-surface (closest-point projections). While the vertex-to-vertex scheme performs well for modeling vesicle aggregation⁶⁰ or weak particle adhesion,^81–83 we demonstrated that the more accurate distance calculations in the vertex-to-surface scheme are essential for simulating substantial membrane wrapping during particle uptake.

Dynamic simulations of a cubical particle interacting with a biconcave vesicle show that higher particle inertia impedes the reorientations necessary for full engulfment. In contrast, lighter particles undergo multiple reorientations to achieve complete uptake. Our simulations of rod–vesicle interactions with different initial configurations reveal sequential orientational transitions closely matching the experimental observations. Overall, our results demonstrate that membrane wrapping is governed by the interplay between particle orientation, local membrane curvature, and adhesion strength. Particles adapt their orientations during uptake to reduce energetic barriers and promote favorable membrane deformation. When curvature compatibility is high, this adaptive behavior leads to complete engulfment; otherwise, wrapping is arrested in partial or metastable states. These principles hold across a range of particle and vesicle shapes, underscoring general mechanisms of anisotropic particle uptake. To conclude, this geometry-agnostic framework not only captures steady-state membrane deformations and energy landscapes but also resolves the time-resolved trajectory of particle entry. This work provides a versatile tool for studying environmental colloidal particles and guiding the design of anisotropic nanocarriers in biomedical applications at cellular scales.

Author contributions

D. A. R., J. R., and X. Y. performed research; D. A. R. and J. R. analyzed data; and D. A. R. and X. Y. wrote the manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

The C++ code and data for this article are available at the following repository: https://doi.org/10.5281/zenodo.16906753.

Supplementary information is available. See DOI: https://doi.org/10.1039/d5sm00567a

Appendix

Algorithm: Vertex-to-vertex bond identification

Input: V_v (set of vesicle mesh vertices, size N), V_p (set of particle mesh vertices, size M)

Output: b̲o̲n̲d̲s̲ (list of mutual nearest-neighbor vertex index pairs)

Initialize nearestParticle and nearestVesicle to store nearest neighbor vertex indices

for each vesicle vertex v_i ∈ V_vdo:

 Compute the shortest distance for v_i, d_i = min_vj∈Vp|v_i − v_j|

 Store the index of the particle vertex corresponding to d_i, nearestParticle[i] = j

end for

for each particle vertex v_j∈V_pdo:

 Computer the shortest distance for v_j as d_j = min_vi∈Vv|v_i − v_j|

 Store the index of the particle vertex corresponding to d_j, nearestVesicle[j] = i

end for

Initialize empty list bonds

for each vesicle vertex index i = 1 to Ndo:

 if the pair is mutually nearest, nearestVesicle(nearestParticle[i]) = = i

 then Append (i,nearestParticle[i]) to b̲o̲n̲d̲s̲

Return b̲o̲n̲d̲s̲

Acknowledgements

X. Y. gratefully acknowledges funding from the U.S. National Science Foundation for supporting this work through award 2034855/2448213. Computing time was provided by the Center for Computational Research at the University at Buffalo and the Theory and Computation facility of the Center for Functional Nanomaterials (CFN), which is a U.S. Department of Energy Office of Science User Facility, at Brookhaven National Laboratory under contract no. DE-SC0012704. We would also like to acknowledge helpful discussions with Emad Pirhadi.

References

1. G. J. Cannon and J. A. Swanson, The macrophage capacity for phagocytosis, J. Cell Sci., 1992, 101, 907–913.
2. R. S. Flannagan, V. Jaumouillé and S. Grinstein, The Cell Biology of Phagocytosis, Annu. Rev. Pathol.: Mech. Dis., 2012, 7, 61–98.
3. S. D. Conner and S. L. Schmid, Regulated portals of entry into the cell, Nature, 2003, 422, 37–44.
4. R. K. Sadhu, S. R. Barger, S. Penič, A. Iglič, M. Krendel, N. C. Gauthier and N. S. Gov, A theoretical model of efficient phagocytosis driven by curved membrane proteins and active cytoskeleton forces, Soft Matter, 2023, 19, 31–43.
5. S. Han, E. Pirhadi, X. Yong and S. Shin, PFOA induces fission of phase-separated phospholipid vesicles, Chem. Commun., 2025 10.1039/D5CC03353E.
6. E. Pirhadi, J. M. Vanegas, M. Farin, J. W. Schertzer and X. Yong, Effect of Local Stress on Accurate Modeling of Bacterial Outer Membranes Using All-Atom Molecular Dynamics, J. Chem. Theory Comput., 2023, 19, 363–372.
7. S. Kefer, O. Miesbauer and H.-C. Langowski, Environmental Microplastic Particles vs. Engineered Plastic Microparticles—A Comparative Review, Polymers, 2021, 13, 2881.
8. X. Yong and K. Du, Effects of Shape on Interaction Dynamics of Tetrahedral Nanoplastics and the Cell Membrane, J. Phys. Chem. B, 2023, 127, 1652–1663.
9. E. J. Ewins, K. Han, B. Bharti, T. Robinson, O. D. Velev and R. Dimova, Controlled adhesion, membrane pinning and vesicle transport by Janus particles, Chem. Commun., 2022, 58, 3055–3058.
10. A. H. Bahrami and T. R. Weikl, Curvature-Mediated Assembly of Janus Nanoparticles on Membrane Vesicles, Nano Lett., 2018, 18, 1259–1263.
11. A. H. R. Koch, S. Morsbach, T. Bereau, G. Lévêque, H.-J. Butt, M. Deserno, K. Landfester and G. Fytas, Probing Nanoparticle/Membrane Interactions by Combining Amphiphilic Diblock Copolymer Assembly and Plasmonics, J. Phys. Chem. B, 2020, 124, 742–750.
12. X. Zhang, G. Ma and W. Wei, Simulation of nanoparticles interacting with a cell membrane: probing the structural basis and potential biomedical application, NPG Asia Mater., 2021, 13, 52.
13. J. Agudo-Canalejo, Engulfment of ellipsoidal nanoparticles by membranes: full description of orientational changes, J. Phys.: Condens. Matter, 2020, 32, 294001.
14. J. Agudo-Canalejo and R. Lipowsky, Critical Particle Sizes for the Engulfment of Nanoparticles by Membranes and Vesicles with Bilayer Asymmetry, ACS Nano, 2015, 9, 3704–3720.
15. J. Midya, T. Auth and G. Gompper, Membrane-Mediated Interactions Between Nonspherical Elastic Particles, ACS Nano, 2023, 17, 1935–1945.
16. N. Nambiar, Z. A. Loyd and S. M. Abel, Particle Deformability Enables Control of Interactions between Membrane-Anchored Nanoparticles, J. Chem. Theory Comput., 2024, 20, 1732–1739.
17. B. Li and S. M. Abel, Shaping membrane vesicles by adsorption of a semiflexible polymer, Soft Matter, 2018, 14, 185–193.
18. M. V. Baranov, M. Kumar, S. Sacanna, S. Thutupalli and G. Van Den Bogaart, Modulation of Immune Responses by Particle Size and Shape, Front. Immunol., 2021, 11, 607945.
19. J. A. Champion and S. Mitragotri, Role of target geometry in phagocytosis, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 4930–4934.
20. A. H. Bahrami, Orientational changes and impaired internalization of ellipsoidal nanoparticles by vesicle membranes, Soft Matter, 2013, 9, 8642.
21. S. Dasgupta, T. Auth and G. Gompper, Shape and Orientation Matter for the Cellular Uptake of Nonspherical Particles, Nano Lett., 2014, 14, 687–693.
22. H. Deng, P. Dutta and J. Liu, Entry modes of ellipsoidal nanoparticles on a membrane during clathrin-mediated endocytosis, Soft Matter, 2019, 15, 5128–5137.
23. S. Dasgupta, T. Auth and G. Gompper, Wrapping of ellipsoidal nano-particles by fluid membranes, Soft Matter, 2013, 9, 5473–5482.
24. Q. Yu, S. Othman, S. Dasgupta, T. Auth and G. Gompper, Nanoparticle wrapping at small non-spherical vesicles: curvatures at play, Nanoscale, 2018, 10, 6445–6458.
25. H. Tang, H. Zhang, H. Ye and Y. Zheng, Receptor-Mediated Endocytosis of Nanoparticles: Roles of Shapes, Orientations, and Rotations of Nanoparticles, J. Phys. Chem. B, 2018, 122, 171–180.
26. S. van der Ham, J. Agudo-Canalejo and H. R. Vutukuri, Role of Shape in Particle-Lipid Membrane Interactions: From Surfing to Full Engulfment, ACS Nano, 2024, 18, 10407–10416.
27. X. Liu, T. Auth, N. Hazra, M. F. Ebbesen, J. Brewer, G. Gompper, J. J. Crassous and E. Sparr, Wrapping anisotropic microgel particles in lipid membranes: Effects of particle shape and membrane rigidity, Proc. Natl. Acad. Sci. U. S. A., 2023, 120, e2217534120.
28. D. Paul, S. Achouri, Y.-Z. Yoon, J. Herre, C. E. Bryant and P. Cicuta, Phagocytosis Dynamics Depends on Target Shape, Biophys. J., 2013, 105, 1143–1150.
29. E. Alizadeh-Haghighi, A. Karaei Shiraz and A. H. Bahrami, Membrane-mediated interactions between disk-like inclusions adsorbed on vesicles, Front. Phys., 2022, 10, 1020619.
30. C. Huang, Y. Zhang, H. Yuan, H. Gao and S. Zhang, Role of Nanoparticle Geometry in Endocytosis: Laying Down to Stand Up, Nano Lett., 2013, 13, 4546–4550.
31. A. Azadbakht, B. Meadowcroft, T. Varkevisser, A. Šarić and D. J. Kraft, Wrapping Pathways of Anisotropic Dumbbell Particles by Giant Unilamellar Vesicles, Nano Lett., 2023, 23, 4267–4273.
32. X. Yi, X. Shi and H. Gao, A Universal Law for Cell Uptake of One-Dimensional Nanomaterials, Nano Lett., 2014, 14, 1049–1055.
33. A. H. Bahrami, R. Lipowsky and T. R. Weikl, Tubulation and Aggregation of Spherical Nanoparticles Adsorbed on Vesicles, Phys. Rev. Lett., 2012, 109, 188102.
34. Q. Yu, S. Dasgupta, T. Auth and G. Gompper, Osmotic Concentration-Controlled Particle Uptake and Wrapping-Induced Lysis of Cells and Vesicles, Nano Lett., 2020, 20, 1662–1668.
35. H. Ni and G. A. Papoian, Membrane-MEDYAN: Simulating Deformable Vesicles Containing Complex Cytoskeletal Networks, J. Phys. Chem. B, 2021, 125, 10710–10719.
36. K. Popov, J. Komianos and G. A. Papoian, MEDYAN: Mechanochemical Simulations of Contraction and Polarity Alignment in Actomyosin Networks, PLoS Comput. Biol., 2016, 12, e1004877.
37. A. Sciortino, H. A. Faizi, D. A. Fedosov, L. Frechette, P. M. Vlahovska, G. Gompper and A. R. Bausch, Active membrane deformations of a minimal synthetic cell, Nat. Phys., 2025, 21, 799–807.
38. S. Sadhukhan, S. Penič, A. Iglič and N. S. Gov, Modelling how curved active proteins and shear flow pattern cellular shape and motility, Front. Cell Dev. Biol., 2023, 11, 1193793.
39. R. K. Sadhu, M. Luciano, W. Xi, C. Martinez-Torres, M. Schröder, C. Blum, M. Tarantola, S. Villa, S. Penič, A. Iglič, C. Beta, O. Steinbock, E. Bodenschatz, B. Ladoux, S. Gabriele and N. S. Gov, A minimal physical model for curvotaxis driven by curved protein complexes at the cell’s leading edge, Proc. Natl. Acad. Sci. U. S. A., 2024, 121, e2306818121.
40. R. K. Sadhu, A. Iglič and N. S. Gov, A minimal cell model for lamellipodia-based cellular dynamics and migration, J. Cell Sci., 2023, 136, jcs260744.
41. M. Fošnarič, S. Penič, A. Iglič, V. Kralj-Iglič, M. Drab and N. S. Gov, Theoretical study of vesicle shapes driven by coupling curved proteins and active cytoskeletal forces, Soft Matter, 2019, 15, 5319–5330.
42. Y. Ravid, S. Penič, Y. Mimori-Kiyosue, S. Suetsugu, A. Iglič and N. S. Gov, Theoretical model of membrane protrusions driven by curved active proteins, Front. Mol. Biosci., 2023, 10, 1153420.
43. W. Pezeshkian and S. J. Marrink, Simulating realistic membrane shapes, Curr. Opin. Cell Biol., 2021, 71, 103–111.
44. A. Iglič, T. Slivnik and V. Kralj-Iglič, Elastic properties of biological membranes influenced by attached proteins, J. Biomech., 2007, 40, 2492–2500.
45. B. Li and S. M. Abel, Shaping membrane vesicles by adsorption of hinge-like nanoparticles, J. Chem. Phys., 2024, 160, 194901.
46. L. Mesarec, W. Góźdź, V. K. Iglič, S. Kralj and A. Iglič, Closed membrane shapes with attached BAR domains subject to external force of actin filaments, Colloids Surf., B, 2016, 141, 132–140.
47. A. Frost, V. M. Unger and P. De Camilli, The BAR Domain Superfamily: Membrane-Molding Macromolecules, Cell, 2009, 137, 191–196.
48. D. A. Redwan, K. Du and X. Yong, Probing wrapping dynamics of spherical nanoparticles by 3D vesicles using force-based simulations, Soft Matter, 2024, 20, 4548–4560.
49. C. Zhu, C. T. Lee and P. Rangamani, Mem3DG: Modeling membrane mechanochemical dynamics in 3D using discrete differential geometry, Biophys. Rep., 2022, 2, 100062.
50. W. Helfrich, Elastic Properties of Lipid Bilayers: Theory and Possible Experiments, Z. Naturforsch., C: J. Biosci., 1973, 28, 693–703.
51. J. Yang and J. Kim, Phase-field simulation of multiple fluid vesicles with a consistently energy-stable implicit–explicit method, Comput. Methods Appl. Mech. Eng., 2023, 417, 116403.
52. X. Bian, S. Litvinov and P. Koumoutsakos, Bending models of lipid bilayer membranes: Spontaneous curvature and area-difference elasticity, Comput. Methods Appl. Mech. Eng., 2020, 359, 112758.
53. A. Guckenberger, A. Kihm, T. John, C. Wagner and S. Gekle, Numerical–experimental observation of shape bistability of red blood cells flowing in a microchannel, Soft Matter, 2018, 14, 2032–2043.
54. J. Zhao, S. Zhao and S. Luding, The role of particle shape in computational modelling of granular matter, Nat. Rev. Phys., 2023, 5, 505–525.
55. S. Zhao and J. Zhao, Revolutionizing granular matter simulations by high-performance ray tracing discrete element method for arbitrarily-shaped particles, Comput. Methods Appl. Mech. Eng., 2023, 416, 116370.
56. P.-O. Persson and G. Strang, A Simple Mesh Generator in MATLAB, SIAM Rev., 2004, 46, 329–345.
57. A. Jacobson, D. Panozzo, C. Schüller, O. Diamanti, Q. Zhou and N. Pietroni, libigl: A simple C++ geometry processing library, Google Scholar, 2013.
58. A. Šarić and A. Cacciuto, Mechanism of Membrane Tube Formation Induced by Adhesive Nanocomponents, Phys. Rev. Lett., 2012, 109, 188101.
59. M. Raatz, R. Lipowsky and T. R. Weikl, Cooperative wrapping of nanoparticles by membrane tubes, Soft Matter, 2014, 10, 3570–3577.
60. D. A. Fedosov, W. Pan, B. Caswell, G. Gompper and G. E. Karniadakis, Predicting human blood viscosity in silico, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 11772–11777.
61. S. Runser, R. Vetter and D. Iber, SimuCell3D: three-dimensional simulation of tissue mechanics with cell polarization, Nat. Comput. Sci., 2024, 4, 299–309.
62. B. A. Payne and A. W. Toga, Distance field manipulation of surface models, IEEE Comput. Graph. Appl., 1992, 12, 65–71.
63. J. A. Baerentzen and H. Aanaes, Signed distance computation using the angle weighted pseudonormal, IEEE Trans. Vis. Comput. Graph., 2005, 11, 243–253.
64. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, LAMMPS – a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comput. Phys. Commun., 2022, 271, 108171.
65. P. Iyer, G. Gompper and D. A. Fedosov, Dynamic shapes of floppy vesicles enclosing active Brownian particles with membrane adhesion, Soft Matter, 2023, 19, 3436–3449.
66. P. Iyer, G. Gompper and D. A. Fedosov, Non-equilibrium shapes and dynamics of active vesicles, Soft Matter, 2022, 18, 6868–6881.
67. F. Frey, F. Ziebert and U. S. Schwarz, Dynamics of particle uptake at cell membranes, Phys. Rev. E, 2019, 100, 052403.
68. K. Xiao, R. Ma and C.-X. Wu, Wrapping dynamics and critical conditions for active nonspherical nanoparticle uptake, Phys. Rev. E, 2023, 107, 054401.
69. P. Czajka, J. M. Antosiewicz and M. Długosz, Effects of Hydrodynamic Interactions on the Near-Surface Diffusion of Spheroidal Molecules, ACS Omega, 2019, 4, 17016–17030.
70. R. Dimova and B. Pouligny, in Methods in Membrane Lipids, ed. A. M. Dopico, Humana Press, Totowa, NJ, 2007, pp. 227–236 DOI:10.1007/978-1-59745-519-0_15.
71. M. Sadeghi and F. Noe, Large-scale simulation of biomembranes incorporating realistic kinetics into coarse-grained models, Nat. Commun., 2020, 11, 2951.
72. F. Farivar, H. Zhang, Z. F. Tian and A. Gupte, CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles, Powder Technol., 2020, 360, 1017–1027.
73. D. Baraff, Physically based modeling: Rigid body simulation, SIGGRAPH Course Notes, ACM SIGGRAPH, 2001, vol. 2, p. 2-1.
74. W. Rawicz, K. C. Olbrich, T. McIntosh, D. Needham and E. Evans, Effect of Chain Length and Unsaturation on Elasticity of Lipid Bilayers, Biophys. J., 2000, 79, 328–339.
75. I. R. Cooke, K. Kremer and M. Deserno, Tunable generic model for fluid bilayer membranes, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 011506.
76. U. Seifert, K. Berndl and R. Lipowsky, Shape transformations of vesicles: Phase diagram for spontaneous- curvature and bilayer-coupling models, Phys. Rev. A: At., Mol., Opt. Phys., 1991, 44, 1182–1202.
77. S. Svetina and B. Žekš, Shape behavior of lipid vesicles as the basis of some cellular processes, Anat. Rec., 2002, 268, 215–225.
78. P. Kattel, S. Sulthana, J. Trousil, D. Shrestha, D. Pearson and S. Aryal, Effect of Nanoparticle Weight on the Cellular Uptake and Drug Delivery Potential of PLGA Nanoparticles, ACS Omega, 2023, 8, 27146–27155.
79. A. Sakane, T.-a Yano, T. Uchihashi, K. Horikawa, Y. Hara, I. Imoto, S. Kurisu, H. Yamada, K. Takei and T. Sasaki, JRAB/MICAL-L2 undergoes liquid–liquid phase separation to form tubular recycling endosomes, Commun. Biol., 2021, 4, 551.
80. F. Ferraro, M. L. da Silva, W. Grimes, H. K. Lee, R. Ketteler, J. Kriston-Vizi and D. F. Cutler, Weibel-Palade body size modulates the adhesive activity of its von Willebrand Factor cargo in cultured endothelial cells, Sci. Rep., 2016, 6, 32473.
81. S. Hillringhaus, A. K. Dasanna, G. Gompper and D. A. Fedosov, Importance of Erythrocyte Deformability for the Alignment of Malaria Parasite upon Invasion, Biophys. J., 2019, 117, 1202–1214.
82. S. Hillringhaus, A. K. Dasanna, G. Gompper and D. A. Fedosov, Stochastic bond dynamics facilitates alignment of malaria parasite at erythrocyte membrane upon invasion, eLife, 2020, 9, e56500.
83. A. K. Dasanna, S. Hillringhaus, G. Gompper and D. A. Fedosov, Effect of malaria parasite shape on its alignment at erythrocyte membrane, eLife, 2021, 10, e68818.

---