Non-equilibrium shapes and dynamics of active vesicles †

Active vesicles, constructed through the confinement of self-propelled particles (SPPs) inside a lipid membrane shell, exhibit a large variety of non-equilibrium shapes, ranging from the formation of local tethers and dendritic conformations, to prolate and bola-like structures. To better understand the behavior of active vesicles, we perform simulations of membranes modelled as dynamically triangulated surfaces enclosing active Brownian particles. A systematic analysis of membrane deformations and SPP clustering, as a function of SPP activity and volume fraction inside the vesicle is carried out. Distributions of membrane local curvature, and the clustering and mobility of SPPs obtained from simulations of active vesicles are analysed. There exists a feedback mechanism between the enhancement of membrane curvature, the formation of clusters of active particles, and local or global changes in vesicle shape. The emergence of active tension due to the activity of SPPs can well be captured by the Young– Laplace equation. Furthermore, a simple numerical method for tether detection is presented and used to determine correlations between the number of tethers, their length, and local curvature. We also provide several geometrical arguments to explain different tether characteristics for various conditions. These results contribute to the future development of steerable active vesicles or soft micro-robots whose behaviour can be controlled and used for potential applications.


Introduction
Many phenomena in biological cells, including cytoskeletal dynamics, 1,2 cytoplasmic streaming, 3 and cell migration, 4,5 are intrinsically out of equilibrium. 6Such processes and functions inevitably rely on force generation by internal active elements, whose primary example is a filament-based cytoskeleton driven by the dynamic growth/shrinkage of its filaments and the activity of molecular motors. 2,7,8Non-equilibrium force generation by the cytoskeleton can lead to dynamic changes in cell shape, such as the formation of filopodia and lamellipodia 9,10 and highly branched elongated structures in neurons and glia cells. 11,12Furthermore, non-equilibrium processes can contribute to different physical effects, such as active membrane fluctuations 13,14 in cells, and the emergence of collective behavior within a group of biological or artificial microswimmers. 15,16o tackle the complexity of biological cells in a more systematic and controlled way, the field of synthetic biology has recently emerged. 18,19In particular, the bottom-up approach, where a synthetic cell is built from scratch by combining different cell-mimicking components, 19,20 is appealing to physicists, as the composition and complexity of such systems can precisely be controlled.Examples include chemical factories that imitate different metabolic processes, 21 membrane compartments with cytoskeletal filaments, [22][23][24] and growth and spontaneous division of droplet-based or vesicle-based compartments. 25,26The development of cell mimics also includes artificial systems built from non-cellderived components.][34] A recent study combining experiments and simulations 17 has realized active vesicles driven from inside by diffusiophoretic self-propelled particles (SPPs).Determined by the balance between forces exerted by active particles and membrane bending rigidity and tension, several different regimes of vesicle shapes are observed, see Fig. 1.For weak particle activity, SPPs cannot strongly deform the membrane, resulting in a quasi-spherical shape with active fluctuations.As the particle activity is increased at small to moderate particle densities, tethers are formed by single SPPs or their clusters.For large particle densities and high propulsion strengths of SPPs, bola-like and/or prolate vesicle shapes can form with two or more satellite vesicles.A similar active system has been realized by enclosing swimming bacteria inside a lipid vesicle, though the density of bacteria was kept relatively low. 316][37] For instance, a non-zero spontaneous curvature induced by the adsorption of amphipathic peptides or BAR domain proteins can lead to string-of-pearls-like and tubular protrusions, 37,38 reminiscent of tether-like extensions of active vesicles.However, such structures in equilibrium are stationary and do not possess active force generation and dynamic shape changes, which occur in non-equilibrium systems driven by the underlying active processes.A recent example for a non-equilibrium vesicle system is the active shape oscillations of giant vesicles driven by dynamic changes in spontaneous curvature. 39For the case of vesicles with enclosed active particles, their highly dynamic shapes largely depend on the strength of active forces and their distribution, and the dynamics of the active components, so that fixed spontaneous curvature and reduced volume play at most a secondary role for such systems. 17everal physical mechanisms govern the behavior of active vesicles, including the curvature-induced clustering of active particles 40 near the membrane, and generation of active tension through the swim pressure. 41The accumulation of active particles near the membrane is similar to wall accumulation of SPPs, which occurs due to a slow rotational diffusion so that SPPs spend a significant time at the walls. 40,42Local membrane curvature further enhances accumulation of SPPs at the membrane, and results in the formation of particle clusters, which are able to collectively alter the shape of active vesicles. 17,27,29,30At large enough SPP densities, motilityinduced phase separation can also take place, [43][44][45] leading to large clusters at the membrane surface.In order to better understand how far active vesicles are driven out of equilibrium, it is necessary to study the interplay between vesicle shape changes, local membrane properties, and the dynamics of active particles for various SPP densities and propulsion strengths.
In this study, we perform a systematic simulation investigation and analysis of active vesicle shapes and the clustering characteristics of SPPs.We demonstrate how varying particle activity and density strongly affect membrane properties such as local curvature and tension.Furthermore, changes in vesicle shapes are directly correlated with alterations in the dynamics and clustering of SPPs, resulting in a feedback mechanism between the enhancement of membrane curvature and the formation of active-particle clusters.In the tethering regime, the tether properties such as length and local curvature correlate with the number of tethers.Most of the non-equilibrium aspects of active vesicles can well be rationalized by simple theoretical arguments, providing a good explanation for the non-equilibrium behavior of active vesicles.[48] 2 Methods and models An active vesicle is represented by a closed fluid membrane of spherical topology with radius R, enclosing a number of active Brownian particles.The activity of the particles is described by the dimensionless Peclet number Pe = sv p /D t , where s is the particle diameter, v p the propulsion velocity, and D t the translational diffusion coefficient.Note that Pe is a measure of the propulsion force f p of SPPs, with v p = f p /g p and D t = k B T/g p , where g p is the translational friction coefficient, so that Pe = f p s/k B T. Particle volume fraction within the vesicle is given by f = N p (s/2R) 3 , where N p is the number of enclosed SPPs.All simulation parameters are given in Table 1.

Model of self-propelled particles
SPPs are modeled as spherical active Brownian particles (ABPs) without hydrodynamic interactions.Each ABP experiences a propulsion force f p that acts along an orientation vector e i .The force results in a propulsion velocity hv p i = f p /g p , where g p is the friction coefficient.The orientation vector e i is subject to random rotation e ˙i = z i Â e i , where z i is a Gaussian random process with hz i (t)i = 0 and hz i (t)z j (t 0 )i = 2D r d ij d(t À t 0 ) with a rotational diffusional coefficient D r .D r is related to the ABP size s and translational diffusion coefficient D t as D t = D r s 2 /3.Furthermore, the ABPs repel each other and the membrane, which is implemented through the repulsive part of the 12-6 Lennard-Jones potential, with the potential minimum at r m = s for ABP-ABP interactions and r m = 0.5s for ABP-membrane interactions.The potential cut-off is set at r c = 2 1/6 r m , so that no attractive interactions are present.

Membrane model
The vesicle is modeled by a dynamically triangulated membrane of N v linked vertices. 49,50The interaction between linked vertices is controlled via a tethering potential 49,51 that is a combination of attractive and repulsive parts Here, k b is the bond stiffness, l min and l max are the minimum and maximum bond lengths, and l c 0 and l c 1 are the potential cutoff lengths.
The bending elasticity is modeled by the Helfrich curvature energy, 52 where k c is the bending rigidity and % c is the mean local curvature at the membrane surface element dA.In the discretized form, it becomes 53,54 where n i is the unit normal at the membrane vertex i, s i ¼ P jðiÞ s ij r ij is the area corresponding to vertex i (the area of the dual cell), j(i) corresponds to all vertices linked to vertex i, and s ij = r ij (cot y 1 + cot y 2 )/2 is the length of the bond in the dual lattice, where y 1 and y 2 are the angles at the two vertices opposite to the edge ij in the dihedral.In practice, since the dihedral terms corresponding to s ij are additive, the local curvature at each vertex can be calculated by summing over contributions from all triangles containing that vertex.The area conservation is imposed locally to each triangle by the potential where N t = 2(N v À 2) is the number of triangles, A l = A 0 /N t is the targeted local area (A 0 is the total membrane area), A i is the instantaneous local area, and k a is the local-area conservation coefficient.Also, there is a constraint applied on the total volume V of the vesicle via the potential where k v is the volume constraint coefficient and V 0 is the desired total volume.Membrane fluidity is modelled by a stochastic flipping of bonds following a Monte-Carlo scheme.The bond shared by each pair of adjacent triangles can be flipped to connect the two previously unconnected vertices. 49,54The flipping is performed with a frequency o and probability c.

Equation of motion
The system evolves in time according to the Langevin equation where m is the mass of membrane particle or SPP, € r i and r _ i represent the second and first time derivatives of particle positions, r i is the spatial derivative at particle i, and U tot is the sum of all interaction potentials described above.The effect of a viscous fluid is mimicked by the friction co-efficient g, whose value can be different for membrane particles and SPPs.Thermal fluctuations are modelled as a Gaussian random process x i with hx i (t)i = 0 and hx i (t)x j (t 0 )i = d ij d(t À t 0 ).Inertial effects are minimized by performing the simulations in the over-damped limit with m and g such that The positions and velocities of all particles are integrated using the velocity-Verlet algorithm. 56Results

Dynamic vesicle shapes
In the absence of volume constraint at Pe = 0 (i.e., passive particles), the vesicle has a nearly spherical shape with thermal membrane fluctuations for the entire range of volume fractions f.As Pe is increased, active fluctuations appear in addition to the thermal undulations, representing the fluctuating regime (see Fig. 1).However, in the fluctuating regime for Pe t 100, the average shape of the vesicle remains nearly spherical, independently of f, since the propulsion force of SPPs or their clusters is too weak to induce tethering or dramatic changes in the vesicle shape.For Pe \ 100 and f t 0.1, dynamic tether formation becomes possible by single SPPs or their clusters.At Pe \ 100 and large enough f \ 0.1, several large clusters of SPPs may form at the membrane, which can pull in different directions and induce a spontaneous symmetry breaking to prolate and bola-like shapes.Note that the shapes in the prolate/bola regime are also very dynamic.
Due to the rotational diffusion of the SPPs and the dynamic adjustment of membrane shape in reaction to the pushing forces generated by the active particles, vesicle shapes are highly dynamic (e.g., persistent tether formation and retraction).This is illustrated in Fig. 2 and corresponding Movie S1 (ESI †).Tether formation starts when the orientation of a SPP is essentially along the membrane normal.For times t t t r = 1/D r , the SPP moves ballistically and pulls a straight tether.The particle velocity slows down with increasing tether length, as the friction imposed by the tether increases linearly with tether length.At times t C t r , the pulling direction of the SPP starts to deviate from the initial tether orientation, resulting in tether curving as it elongates further.Now, the retraction force of the tether given by k/R b (R b is the tether radius) in the absence of a membrane tension shortens the tether length by pulling it toward a straight configuration.At the same time, the anchoring point of the tether at the vesicle may move in the force direction.Finally, on even larger time scale t 4 t r , there is a significant probability that the propulsion direction of the SPP reverts toward the vesicle surface, so that the SPP would move back through the tether accompanied by tether retraction.
When there are more than one SPP in the tether, the dynamics becomes even richer, as shown in Fig. 2(e) and (f) and Movie S1 (ESI †).In addition to several SPPs pulling together to extend the tether, some active particles may have changed their orientation at times t 4 t r , and either move back along the tether or on their way back to the vesicle pull in a direction roughly perpendicular to the tether orientation, thereby inducing a kink in the tether shape or a branching of the tether.
In order to characterize and understand this nonequilibrium behavior and the resulting dynamic shapes of active vesicles, we perform a detailed analysis of their various properties.Note that the analysis of the membrane properties and particle clustering focuses on simulations without a volume constraint, whereas for the analysis of tethers, simulations both with and without a volume constraint are considered.indicate considerable variations in the local curvature induced by thermal fluctuations of the membrane and changes in the vesicle shape due to SPPs.The distribution peak for the ''fluctuating'' phase is centered at hC eq i C 1/R, and shifts towards larger values from the fluctuating to the tethering regime.The tethering regime is associated with the development of regions with a high curvature.This leads to a non-zero skewness h(C À hCi) 3 i/h(C À hCi) 2 i 3/2 of the distribution, as shown in the inset of Fig. 3(a).In the 'bola' regime, the distribution of local curvature is narrower than for the other cases due to smaller fluctuations of the vesicle.In this case, membrane fluctuations become suppressed as a result of high density of SPPs and the development of significant active tension. 27,28t low f, an increase in Pe causes a transition from the fluctuating to the tethering regime.Tethers are regions of a  high curvature due to their small radii, causing an overall increase in the mean local curvature hCi.Consequently, the transition into the tethering regime at Pe B 100 is accompanied by a steep increase in hCi, as shown in Fig. 3(b).Note that the values of hCi at each Pe are averaged over various f values lying in the fluctuating and tethering regimes, i.e. f o 0.18.

Squared distance from the center of mass.
To differentiate between active vesicles in the tethering regime and in the fluctuating/bola regimes, a useful quantity is the squared distance R i 2 /hR i i 2 of the membrane vertices from the center-ofmass of the vesicle.This quantity characterizes well a deviation of the vesicle shape from a sphere, as tethers cause the mean value of R i 2 /hR i i 2 to increase.Fig. 3(c) shows that the distribution in the tethering regime is significantly wider than those in the bola and fluctuating regimes.Furthermore, the distribution for the tethering regime exhibits a pronounced tail due to the formation of multiple long tethers.The distribution for the bola regime also has a tail, although it decays much faster than that in the tethering case.The distribution for the fluctuating regime is nearly symmetric with no significant tails.

Membrane tension.
The membrane tension l of active vesicles can be divided into two contributions: active and passive (equilibrium) tension, i.e. l = l a + l eq .The passive tension is computed as the membrane tension (see Appendix A) for an equilibrium simulation of the vesicle at f = 0 (i.e., no SPPs), and is close to zero for a vesicle in equilibrium, i.e. l = l eq C 0. Even though l eq C 0, it is not exactly zero due to numerical errors and approximations used in the calculations.Therefore, only the deviations from the equilibrium are considered for the membrane tension, knowing that the tension of the vesicle in equilibrium is zero. 57,58In the presence of active particles, the membrane tension increases and is generated as a consequence of the swim pressure 41 due to the active Brownian particles.
For a fixed Pe, an increase in particle density f leads to a linear increase in the vesicle tension, as shown in Fig. 4(a).However, Peclet number affects the slope, so that it increases with increasing Pe, and saturates for Pe 4 100.A similar effect on membrane tension is observed for a fixed particle density, when Pe is increased, see Fig. 4(b).The dependence of l on Pe is also linear with a slope nearly independent of f.
To better understand these trends, we estimate the active tension originating from the swim pressure via the Young-Laplace equation as where w r 1 is the active tension weight which accounts for deviations of the particle force direction from orientation along the membrane normal.Its value changes with both Pe and f, i.e. w = w(Pe,f).With the proportionality Pe p f p and f p N p , the active tension can be described as where l = l/l 0 is the vesicle tension normalized by l 0 = R 2 k B T/ (ps 4 ).This theoretical argument supports the linear dependence of l in Fig. 4(a) and (b) obtained from simulations.
To investigate the dependence of w on the parameters involved, we determine the fraction r of SPPs located at the membrane for different Pe and f. r is calculated as the number of particles N p,mem , which are in direct or indirect contact with the membrane (within a distance of 1.75s, i.e. accounting for up to three layers of SPPs), divided by the total number of SPPs as r = N p,mem /N p .All values of r are averaged over 10 statistically independent timeframes of simulations.The values of w are obtained by fitting a linear function in eqn (9) to the data in Fig. 4 .This approximation allows us to qualitatively extract the effect of different Pe on w by measuring the slope w(Pe,f 0 ).The same method can be employed for the variation of l with Pe.Fig. 5(a) and (c) shows that both w and r follow similar trends as a function of Pe, suggesting that the larger is the This journal is © The Royal Society of Chemistry 2022 particle fraction at the membrane, the larger the contribution of SPPs to active tension.Note that for low particle densities f o 0.005, w can be larger than unity [see Fig. 5(a)], because the distribution of SPPs at the membrane is very inhomogeneous, which makes the analysis based on the Young-Laplace equation questionable.However, for f \ 0.005, w values remain in the range between 0.7 and 1.0, indicating that the majority of SPPs exert a force along a direction close to the normal direction of the membrane.This is due to the fact that the rotational diffusion of SPPs is small (or Pe is large), so that the SPPs primarily push in the direction of the membrane normal.Fig. 5(b) and (c) shows that nearly all SPPs are located at the membrane for large Pe, independently of f.However, at low Pe and large f, a considerable fraction of SPPs remain away from the membrane, as the SPPs are able to leave the membrane due to larger rotational diffusion (low Pe).Moreover, at low Pe, most SPPs can never reach the membrane surface due to the comparable effects of diffusion and propulsion.

Active particle clustering
The mobility of SPPs and the membrane-mediated interactions lead to the formation of particle clusters.We perform the analysis of cluster sizes and their dependence on Pe and f.Furthermore, the mobility of single SPPs is characterized by their fixed-time displacements, and the shape of clusters is described through cluster asphericity.
3.3.1 Cluster sizes.Clusters of SPPs are identified by using a clustering algorithm, where particles located within a cutoff distance r cut are considered as a part of the same cluster.In this analysis, r cut is selected to be 1.1s.Fig. 6(a) and (b) presents the distributions of cluster sizes for different Pe and f.At low f, only small clusters with sizes of N c C 10 are formed due to a limited number of particles.Nevertheless, an increase in Pe leads to a larger probability to form larger clusters as seen in Fig. 6(a), since at low f, particles primarily cluster at membrane locations with significant deformation (or large local curvature).SPPs simply spend more time at the locations of large membrane curvature due to longer escape times, so that large curvatures promote SPP accumulation and clustering.This process can also be described as a feedback-loop mechanism, where the arrival of additional particles to a cluster generates larger cluster forces which result in stronger membrane deformation, aiding SPP accumulation and cluster growth.Fig. 6(b) also shows a secondary peak in the cluster-size distribution at large enough f \ 0.1.At large f, SPPs cover a substantial portion of the membrane area, yielding a few large clusters in the prolate regime.Note that, even though small clusters are still present at large f, their probability reduces in comparison to the cases with low f, since the majority of SPPs is located within the large clusters.
3.3.2SPP mobility.In order to determine changes in particle mobility as a function of Pe and f, the fixed-time displacement d of SPPs is computed as a distance traveled by each active particle during a fixed time interval Dt = 0.1t r .For a free SPP in three dimensions, the crossover time t* from the diffusive to the ballistic regime can be estimated as t* C 18t r /Pe 2 .For Pe = 100, this implies Dt C 54t* which is well beyond the diffusive regime, with an expected fixed time displacement The results in Fig. 6(c) and (d) suggest that while a small fraction of SPPs exhibits ballistic motion (indicated by the distribution tail in Fig. 6(d)), the majority of SPPs moves much slower than a free active particle due to the interactions of SPPs with the membrane and each other.The slowing down of SPPs occurs due to their small rotational diffusion such that when the SPP is near the membrane, the orientation vector of the particle mostly points in the direction of the membrane normal.An approximation based on two-dimensional diffusion of a SPP at the membrane surface yields the estimate which is still significantly lower than the distribution peaks in  mobility along the membrane surface becomes hindered due to steric repulsion between SPPs within large clusters.Note that the distributions for f = 0.02 and f = 0.06 in Fig. 6(d) are nearly identical, as the cluster sizes for the both particle densities are small and lie within a similar range (N c A (1,40)).As the cluster size grows with increasing f [see Fig. 6(a)] for a fixed Pe, the lateral mobility of SPPs decreases.Therefore, the results in Fig. 6(c) and (d) suggest that the average cluster sizes determine the peak position of the mobility distribution, while the value of Pe affects the distribution tail.
3.3.3Cluster shapes.We characterize cluster shapes by their asphericity, as defined in Appendix B. Clusters whose asphericity is close to zero have a near-spherical shape, whereas asphericity close to unity implies a rod-like shape.Note that only clusters with more than four particles are taken into consideration in this analysis.Fig. 6(e) and (f) shows a strong dependence of the cluster-shape distribution on particle density for Pe = 100 (fluctuating regime).For low particle densities (f t 0.06), most clusters have an elongated shape; nearly all particles have a direct contact with the membrane surface, so that only single-layer flat clusters of sizes in the range N c A (5,30) are formed as seen in the insets of Fig. 6(e).At large SPP densities (f \ 0.18), large clusters develop, see Fig. 6(a), whose shape is closer to a filled sphere, as they generally consist of several layers of SPPs and are distributed over a considerable area of the vesicle, see the insets of Fig. 6(f).However, several small clusters with N c B 10 may still be present at large f, whose asphericity is close to unity as in the case for low particle densities.Therefore, asphericity distribution in Fig. 6(f) has two peaks at large f.For f = 0.36, almost all clusters have near zero asphericity, since hardly any small clusters remain.

Membrane tethers
Membrane tethers are detected by using an automated algorithm based on the positions of membrane vertices and SPPs, as described in Appendix C. Various tether properties, such as tether size, length, curvature and number, are calculated from simulations both with and without volume constraints.
3.4.1 Tether length distribution.Fig. 7 presents a distribution of tether lengths L teth , which are computed from simulation data as the distance between the vertices in the tether closest and farthest away from the vesicle center of mass.Since active vesicles are highly dynamic and the tethers constantly change their length, the data are collected every Dt = 0.1t r within a time frame 0.4t r o t o 5t r from simulations representing the tethering regime with Pe A (50,800).Note that the range of Pe starts from a lower value than the critical Pe for tether formation in the state diagram in Fig. 1, because simulation data for other reduced volumes, where tether formation may occur at a lower Pe, 17 are also included in the analysis.
To rationalize the computed distribution of L teth , a simple model of linear tether growth is considered as where a 1 accounts for multiple factors such as tether friction and effective pulling force of the SPPs due to clustering and varying propulsion direction, a 2 represents the existence of a minimum tether length required for the automated algorithm to detect it, and t* A [0,1] is the dimensionless time.
First, we randomly sample the lengths L 0 teth;th ¼ Pet Ã for various Pe A (50,800) and t* A [0,1].Then, the required coefficients a 1 and a 2 are calculated as a 1 ¼ stdðL teth Þ=stdðL 0 teth;th Þ and a 2 ¼ hL teth i À a 1 hL 0 teth;th i, where std(*) denotes standard deviation.Finally, the estimated tether lengths are given by L teth;th ¼ a 1 L 0 teth;th þ a 2 .Note that in this model, the range of t* can be taken arbitrarily as a 1 absorbs any scaling of the time, and a 2 accounts for a minimum tether length so that t* can start from zero.However, this is not true for the selection of Pe values, as changing the range of Pe causes a qualitative change in the distribution of L teth,th .Fig. 7 shows a very good agreement between the distributions computed from the simulation data and the model above, suggesting that the linear tether extension model provides a reasonable description for the tether growth and captures the principal mechanism.
3.4.2Geometrical estimates of tether curvature and length as a function of its size.To rationalize dependencies between various tether properties, two classes of tethers are assumed.For tethers formed by less than 10 SPPs, the front (or head) of a tether is relatively small and can be neglected in comparison with a long cylindrical stalk.Such tethers are characterized by the radius R b of the cylindrical stalk.For tethers formed by large clusters, we can assume that the head constitutes the major part of the tether, while the cylindrical stalk can be neglected.These tethers are characterized by the radius R h of the tether head.For a cylindrical tether, the radius in equilibrium is determined by the vesicle tension l as 59 , where k c is the membrane bending rigidity.The particle diameter s is another length scale that is relevant for tether formation; however, it should not affect the radius R b of the stalk for long enough tethers.The average equilibrium tether radius can now be estimated as R b ' ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi k c =2h li p ' 2R=32, where the calculated mean vesicle tension h li C 1.9l 0 is averaged over all simulations with tethers.Note that R b decreases with increasing Pe, as membrane tension increases.Also, it is important to mention that there is a minimal tether radius in our simulations, which is determined by the bond length l b of the discretized membrane; this corresponds to R b \ l b C R/45.
In addition, we define the tether size as the number of vertices N v,teth within the tether which is proportional to the area of the tether A teth as N v,teth = B 1 A teth .The constant B 1 can be calibrated through the relation between the total number of vertices in the membrane and the total membrane area A as B 1 = N v /A = 2387R À2 .Fig. 8(a) presents the relation between the average local curvature hC teth i of the tether and tether size, which is well described by a power law.The value of the powerlaw exponent can be estimated through dimensional analysis. .This rather crude approximation qualitatively captures the dependence of hC teth i on tether size in Fig. 8(a) with an exponent of 0.9 for tethers formed by less than 10 particles and an exponent of 0.4 for tethers formed by more than 10 particles.
For tethers with less than 10 SPPs, a positive linear correlation of tether length with size is presented in Fig. 8(b).We assume that the ''head'' of cylindrical tethers on average has a constant size, which implies a fixed offset.We then obtain the equation where L teth is the tether length.We assume that R b C 2R/32 and Here, the number of tethers and the curvature are averaged over every Dt = 0.1t r , as the active system is highly dynamic and a single snapshot is not sufficient for the analysis.A positive linear correlation of the mean curvature with the number of tethers is found, since the tethers are regions of enhanced local curvature.The slope of the linear fit in Fig. 8(c) can be estimated as follows.With the assumption that all tethers have some mean length hL teth i, a mean radius R h and R b for the tether head and body respectively, the average mean local curvature hCi of the vesicle can be expressed as where N teth is the number of tethers, A = 4pR 2 is the total vesicle area, A teth,b = 2pR b (hL teth i À 2R h ) is the area of the tether body, A teth,h = 4pR h 2 is the area of the tether head, A teth = A teth,b + A teth,h is the total tether area, and hC 0 i is the average mean local curvature of the reduced main vesicle body with area A 0 = 4pR 02 = A À A teth N teth .
In the limit of A teth N teth /A { 1, eqn (14) reduces to where hC eq i = 1/R is the equilibrium mean vesicle curvature.The values of hL teth i = 0.53R (measured from the distribution of lengths), R b C from the calculated mean vesicle tension h li C 1.9l 0 averaged over all simulations with tethers) and R h = 6R/32 (estimated from the calculated mean cluster size hN c i C 9), imply the slope and the value of hCi/hC eq i at N teth = 0 to be 0.2 and unity, respectively.A least squares fit to the data in Fig. 8(c) results in the values of 0.19 and 0.99, which agree very well with the estimate above.Note that despite the tether number, radius, and length being very dynamic quantities, this theoretical argument based on averages captures the data very well.

Summary and conclusions
Unlike static vesicle shapes obtained in equilibrium through energy minimization, active vesicles exhibit highly dynamic shape changes such as persistent tether formation and retraction.Therefore, active vesicles require a statistical physics approach, where stationary states of such systems and the variability in system properties can be characterised by the corresponding distributions and their moments.We have performed simulations and a systematic analysis of structural and dynamic properties of vesicles with enclosed active particles.Distributions of local curvature and the departure of vesicle geometry from a spherical shape are strongly affected by the non-equilibrium nature of this active system.As a result, dramatic and dynamic changes in vesicle shape have a substantial effect on its membrane properties.For instance, there exists a direct relation between active tension of the vesicle and the fraction of SPPs located at the membrane surface, which can be quantified via the Young-Laplace equation.Furthermore, enhanced local curvature of the vesicle favors the accumulation of SPPs, leading to a feedback loop between the evolution of local curvature and particle accumulation.This has been confirmed through the analysis of dynamic SPP clustering at the membrane, including variations in particle displacement, cluster size and shape as a function of Pe and f.At f t 0.1 and large enough Pe, small SPP clusters with N c B 10 are formed, and trigger the formation of various tether structures.At large particle densities, the cluster size distribution develops the second peak due to the formation of very large SPP clusters which cover a substantial area of the membrane.Note that when large clusters are formed, the number of small clusters generally becomes reduced.The formation of large clusters is similar to the motility-induced phase separation found in non-equilibrium systems of active particles. 45,60Distributions of particle mobility suggest a competition between SPP activity leading to an increase in the lateral mobility of particles along the membrane surface and steric repulsion between SPPs, which plays an increasing role in the formation of large clusters.Finally, we have employed an automated method for tether detection, which allows a subsequent analysis of tether number, size, shape, and local curvature for different simulated conditions.Most of the correlations between different tether properties can be rationalized by considering two classes of tether geometries: (i) tethers with a long cylindrical body whose head can be omitted, and (ii) tethers with a large head formed by a cluster of SPPs.The distribution of tether lengths can be rationalized well by considering a linear growth model of the tethers for a given range of Pe.These theoretical arguments provide qualitative understanding of tether structure and variability.
Despite the fact that most of the non-equilibrium aspects of active vesicles can well be rationalized, several observations still remain unclear.The analysis of SPP mobility at the membrane surface shows that the lateral motion of active particles is majorly diffusive.This raises an interesting question whether motility-induced phase separation can take place under certain conditions and how the local curvature of the membrane influences the possible coexistence of low-density and highdensity phases.Furthermore, it would be interesting to study the effects of particle shape and hydrodynamic interactions on the behavior of active vesicles, as such systems can be a precursor for active synthetic micro-robots capable of performing specific functions.For instance, the ability to change shape and induce membrane deformations is of importance for the motility of active systems, 33,46 the internalization of drugdelivery particles by cells, 61 and the development of biomimetic systems which can mimic phagocytic functions. 48verages must be performed for the equivalence between the virial stress and the mechanical stress.Since our system has both two-body and three-body interactions, the determination of the local stress tensor is not immediately obvious.However we are only interested in the trace of the stress tensor, which in two dimensions characterizes the tension of the vesicle.The sum over virial contributions from three body forces (area/volume constraints) is given by where indices i, j and k correspond to the vertices of the triangle and f av,a i is the force on vertex i in the direction a due to the area and volume constraints.This contribution is then divided equally between the three vertices.For bond forces f b,a i and f b,a j acting on vertices i and j, the virial contribution V b is divided equally among the two vertices and is calculated as The total virial contribution from particle forces at each vertex is therefore V = V av /3 + V b /2.The tension in the system is calculated as a spatial and temporal average of the local stresses where the factor two is due to the dimensionality of the system, V i (t) is the virial contribution at vertex i at time t due to interparticle forces and a i is the area of the dual cell, which is approximated by considering that each neighbouring triangle to the vertex contributes roughly one-third to the area.The contribution from momentum transfer is approximated by using the equipartition theorem.
where R peak and std(R) are the peak position and standard deviation of the distance distribution of membrane vertices from the vesicle center of mass, respectively.Vertices that are located beyond this cutoff may belong to a tether.The coefficient 1/4 in front of std(R) is chosen by a trial-and-error procedure such that a large fraction of tethers is retained, while the detection of non-tether-like membrane undulations is minimized.Note that a too large cutoff may exclude the base region of tethers.An adaptive threshold is necessary as a single cutoff value fails to properly detect tethers in different simulations due to a significant complexity of active vesicle shapes.Then, a density-based clustering algorithm is used to determine membrane regions which represent different tethers.
After the identification of possible tether regions, false tether detections have to be eliminated.First, existing SPP clusters are identified and only membrane regions that are associated with some SPP clusters are retained.This association is made if the distance between the membrane patch and a SPP cluster is within a cutoff distance of 2.2R/32 (compare to the particle radius of s/2 = 2R/32).In this step, no minimum cluster size is defined, so that even single SPPs are considered as clusters.This is a valid assumption since even single active particles are able to pull tethers at large enough Pe.The association of SPP clusters with possible tether regions does not guarantee that those regions are tethers, as enhanced membrane undulations are induced by active particles.Therefore, in the second check, the standard deviation s j = std(R i ) of a possible tether region j (R i are vertex distances from the vesicle's center of mass) is considered.For non-tether-like membrane undulations, a low standard deviation is expected, because such regions are typically flat with moderate deviations from the mean value hR i i.However, for a tether, a number of vertices will lie considerably above and below the mean hR i i, leading to a large standard deviation.As a rule of thumb, only regions with s j 4 2.5R/32 are retained as possible tethers, and this cutoff value of s j was selected manually to minimize false detection.
Finally, the last check is performed to eliminate 'global' shape changes of a vesicle misidentified as tethers.These are typically very large regions with a high tension due to a large number of SPPs near them.Such membrane regions are eliminated by considering the 'size' of such regions N v,teth using a selected threshold of N v,teth o 10 4 , i.e. we impose an order of magnitude restriction based on the total size of the vesicle N v = 3 Â 10 4 which ensures that the identified region is not a 'global' shape change.We also restrict the algorithm for f o 0.18 to avoid the bola/prolate regime.Fig. 9

illustrates this
This journal is © The Royal Society of Chemistry 2022 multistep procedure for of tethers, which allows us to automatically identify tethers in most simulations.Other properties of identified tethers such as mean tension, size of the associated cluster, tether length and tether size, can then easily be computed.

Fig. 1
Fig. 1 Phase diagram of various vesicle shapes, including fluctuating, tethering, and bola/prolate regimes, as a function of Peclet number Pe and the volume fraction f of active particles.The solid black lines are the theoretical phase boundaries of the tethering regime.The grey region indicates the crossover between the fluctuating and the bola/prolate regimes.In all simulations of the phase diagram, the membrane area is conserved, while volume conservation is not imposed, so that the dimensionless area-to-volume ratio can adjust freely.Adapted from ref.17.
An energetically favorable bond flip is accepted with a probability of p = 1.For an energetically unfavorable flip, the resulting change in energy due to an attempted bond flip DU = DU att + DU rep + DU A determines the probability of the flipping as c = exp[ÀDU/ k B T].The resulting membrane fluidity can be characterized by a 2D membrane viscosity for the selected frequency o and flipping probability c. 51,55

3. 2
Vesicle membrane characteristics 3.2.1 Local curvature.Fig. 3(a) shows distributions of the reduced local curvature % C = C/hC eq i for different activities and classes of vesicle shapes, where hC eq i is the average local mean curvature measured from an equilibrium simulation of the vesicle, i.e. with no enclosed active particles.The distributions This journal is © The Royal Society of Chemistry 2022

Fig. 3
Fig. 3 (a) Local curvature distributions for the tethering (T) (Pe = 200, f = 0.02), bola (B) (Pe = 200, f = 0.18) and fluctuating (F) (Pe = 25, f = 0.06) regimes.The inset shows the effect of global shape changes on the skewness h(C À hCi) 3 i/h(C À hCi) 2 i 3/2 of the distributions.(b) Mean local curvature hCi as a function of Pe.The values are averaged over different f with f o 0.18.(c) Distribution for the squared distances hR i 2 i/hR i i 2 of the membrane from the vesicle center of mass for the tethering, bola and fluctuating regimes.The tail is pronounced for the bola and tethering regimes with the latter having a significantly slow decay.

Fig. 2
Fig. 2 Several snapshots of a simulation illustrating tether pulling and retraction for parameters Pe = 400, f = 0.002 and a reduced volume n = 0.8.(a and b) For t o t r , SPPs move ballistically and pull a few tethers.(c) When t C t r , SPPs start changing their direction of propulsion which leads to tether curving.(d) For t 4 t r , some of the tethers are retracted since some SPPs may move back toward the vesicle.(e and f) For SPP clusters, some particles continuously leave and join different tethers.See also Movies S1 (ESI †).
(a) and (b).Note that in general, it is difficult to decouple the effects of Pe and f on w.As an approximation, we use w(Pe,f) C w(Pe,f 0 ) for the variation of l with f, where f 0 is an unknown parameter which minimizes the function f ðf 0 Þ ¼ P Pe;f ðwðPe; fÞ À wðPe; f 0 ÞÞ 2 = P Pe;f

Fig. 4
Fig. 4 Linear increase in the reduced vesicle tension l as a function of (a) f and (b) Pe.The values of the vesicle tension are normalized as l = l/l 0 , with l 0 = R 2 k B T/(ps 4 ), and have further been re-scaled by the corresponding value of f and Pe.The solid lines are linear fits to the data.

Fig. 5
Fig.5(a) Dependence of the active tension weight w on Pe and f.The values of w are obtained using a linear fitting of the function l = wPef.Variation in fraction r of SPPs located at the membrane surface as a function of (b) Pe and (c) f.The value of r is calculated as the number of particles N p,mem that are in a direct or indirect contact with the membrane divided by the total number of particles, i.e., r = N p,mem /N p .

Fig. 6 (
c) and (d).Thus, despite a considerable slowing down of the SPPs at the vesicle membrane, the fixed-time displacements are on average larger than d diff due to SPP activity.The value of Pe does not significantly affect the peak of d distributions [see Fig. 6(c)], but mainly alters the distribution tail.For a fixed Pe, an increase in f shifts the peak of mobility distribution to lower values of d [see Fig. 6(d)], as the lateral

Fig. 6
Fig. 6 Distributions of various SPP cluster properties.Cluster-size distributions (a) for different Pe at f = 0.06 and (b) for different f at Pe = 100.Fixedtime displacement d/R distributions, characterizing SPP mobility for various (c) Pe and (d) f for Dt = 0.1D r For tethers formed by clusters with N c o 10, N v,teth B R b and hC teth i B R b À1 , implying a power-law dependence hC teth i B N v,teth À1 .However, for tethers formed by large clusters, N v,teth B R h 2 , so that hC teth i B N v,teth À0.5 32 is determined by considering that on average clusters of size hN c i occupy the tether head, where we have used hN c i = 3.2 as determined from the simulation data.Then, S 2 = B 1 (4pR h (R h À R b )) C 170 can be estimated.The linear fit based on the determined slope S 1 and intercept S 2 is presented in Fig. 8(b) and follows well the simulation data.3.4.3Average membrane curvature with number of tethers.Fig. 8(c) shows the mean local curvature averaged over the whole vesicle as a function of the number of tethers.

Fig. 8
Fig. 8 (a) Mean local tether curvature hC teth i/C eq as a function of tether size N v,teth .The various colors represent tethers formed by SPP clusters of different sizes.The solid black line is a fit for tethers formed by clusters with N c o 10, while the dashed black line is for tethers formed by clusters with N c 4 10.The lines are least-squared fits to the data using the function hC teth i/C eq = AN v,teth ÀB + C, whose form is motivated by the geometrical

Table 1
Parameters used for 3D simulations of active vesicles both in model and physical units.N t = 2N v À 4 is the number of triangular faces in the vesicle discretization.Also, the table includes equilibrium value of local membrane curvature