Uniform and Janus-like nanoparticles in contact with vesicles : energy landscapes and curvature-induced forces

Biological membranes and lipid vesicles often display complex shapes with non-uniform membrane curvature. When adhesive nanoparticles with chemically uniform surfaces come into contact with such membranes, they exhibit four different engulfment regimes as recently shown by a systematic stability analysis. Depending on the local curvature of the membrane, the particles either remain free, become partially or completely engulfed by the membrane, or display bistability between free and completely engulfed states. Here, we go beyond stability analysis and develop an analytical theory to leading order in the ratio of particle-to-vesicle size. This theory allows us to determine the local and global energy landscapes of uniform nanoparticles that are attracted towards membranes and vesicles. While the local energy landscape depends only on the local curvature of the vesicle membrane and not on the overall membrane shape, the global energy landscape describes the variation of the equilibrium state of the particle as it probes different points along the membrane surface. In particular, we find that the binding energy of a partially engulfed particle depends on the ‘unperturbed’ local curvature of the membrane in the absence of the particle. This curvature dependence leads to local forces that pull the partially engulfed particles towards membrane segments with lower and higher mean curvature if the particles originate from the exterior and interior solution, respectively, corresponding to endoand exocytosis. Thus, for partial engulfment, endocytic particles undergo biased diffusion towards the membrane segments with the lowest membrane curvature, whereas exocytic particles move towards segments with the highest curvature. The curvature-induced forces are also effective for Janus particles with one adhesive and one non-adhesive surface domain. In fact, Janus particles with a strongly adhesive surface domain are always partially engulfed which implies that they provide convenient probes for experimental studies of the curvature-induced forces that arise for complex-shaped membranes.


Introduction
Understanding the interaction between nanoparticles and biomembranes is of key importance in areas such as drug delivery, nanotoxicity, viral infection or biomedical imaging. [1][2][3][4][5] The uptake of particles by cells requires the engulfment of the particles by the cell membrane, a process governed by the interplay between particle-membrane adhesion and membrane bending. 6,7 Experimentally, the fundamental features of this engulfment process can be studied using model systems such as lipid [8][9][10][11][12] vesicles or polymersomes. 13 Several computational studies have previously explored the engulfment process of both spherical [14][15][16][17][18][19][20] and nonspherical particles. [21][22][23][24][25] These studies have been impeded by the high dimensionality of the parameter space that describes particle-membrane interactions. Indeed, one may consider (i) different particle properties such as size, adhesiveness, shape and softness; (ii) particles with chemically uniform, Janus-like, or patchy surfaces; (iii) variations in the bending rigidity and spontaneous curvature of the membranes as well as in their phase behavior; or (iv) different geometrical and topological constraints on the membrane, such as closed vesicles with and without volume constraint, which lead to spherical, oblate, prolate, discocyte or stomatocyte-shaped vesicles, as well as planar membranes connected to area reservoirs with a certain membrane tension. Computational studies can typically explore only a small subspace or 'slice' of this high-dimensional parameter space, and are therefore limited in their ability to guide future experiments.
Recently, we overcame this limitation for the particular case of rigid spherical particles with a uniform adhesive surface. Such particles can in principle remain free, that is, unbound from the membrane; they can be partially engulfed, when the particles are bound to but not fully covered by the membrane; or they can be completely engulfed, when the entire particle surface is membrane-bound. In ref. 7, we developed a local stability analysis that allowed us to obtain exact analytical conditions for the energetic stability of free and completely engulfed particles. The two stability conditions depend on the particle size and adhesiveness, the material properties of the membrane, and the local curvature of the membrane, which is affected by the geometric constraints acting on the membrane. Combining the two stability conditions with numerical calculations, we found that small particles can exhibit only four distinct behaviors when in contact with a membrane: a particle can either (i) remain free, (ii) become completely engulfed, (iii) be partially engulfed, or (iv) display a bistability between free and completely engulfed states, which are then separated by an energy barrier. The boundaries between these four regimes directly follow from the stability limits associated with the two stability conditions. Subsequently, in ref. 26, we used the local stability analysis to show that complex-shaped vesicles with non-uniform curvature display different engulfment patterns when exposed to adhesive particles. One example of such a pattern is displayed in Fig. 1(a) for the case of a prolate vesicle. † In this particular example, particles are completely engulfed at the strongly curved poles (green regions) but partially engulfed in the weakly curved equatorial region (blue). However, the approach based on local stability analysis does not give us access to the binding energy and degree of engulfment in the case of partially engulfed particles, or to the height of the energy barrier between free and completely engulfed states in the bistable engulfment regime.
In the present paper, we develop an analytical expression for the energy of the particle-vesicle system to leading order in the ratio of particle-to-vesicle size. This expression allows us to determine the binding energies and degree of engulfment for partially engulfed particles as well as the barrier heights for particles in the bistable regime. In particular, we find that the binding energy of a partially engulfed particle depends on the local, 'unperturbed' curvature of the vesicle in the absence of the particle. As a consequence, partially engulfed particles experience curvature-induced forces towards regions of lower membrane curvature if the particles get in contact with the outer leaflet of the vesicle membrane as in Fig. 1(a), or towards regions of higher membrane curvature if the particles get in contact with the inner leaflet of this membrane. In contrast, completely engulfed particles do not experience such curvatureinduced forces.
The membrane curvature-sensing capability of nanoparticles is certainly of theoretical interest, and may be potentially useful in biomedical imaging or drug-delivery applications. However, nanoparticles with a chemically uniform surface are only partially engulfed if their size and adhesiveness are finely tuned with respect to the material properties and local curvature of the membrane. 7,26 In order to overcome this limitation, we will also consider Janus particles with one strongly adhesive and one non-adhesive surface domain, see Fig. 1(b), which provides a much more robust method to explore curvatureinduced forces. For such Janus particles, the membrane spreads only up to the domain boundaries between the two domains on the particle surface. As a consequence, these particles are always partially engulfed and, thus, subject to curvatureinduced forces.
The paper is organized as follows. First, in Section 2 we describe the combined model for the membrane's curvature elasticity and the membrane-particle adhesion used throughout the paper, and briefly summarize the stability limits for free and completely engulfed particles as obtained in ref. 7. These stability limits are then used, in Section 3, to obtain conditions on the energy of the system close to free and completely engulfed states, in the general case of particles of any size. The small particle approximation is finally developed and verified against numerical calculations in Section 4. From this small particle approximation, in Section 5 we obtain the energy landscapes, binding energies, barrier heights, engulfment patterns and curvature-induced forces experienced by uniformly adhesive particles in contact with vesicles. Finally, in Section 6 we describe the curvature-induced forces experienced by Januslike particles in contact with vesicles of complex shape.

Interplay of curvature elasticity and adhesion
The fate of a particle that comes into contact with a membrane is determined by the interplay between particle adhesion and membrane bending. 6,7 According to the spontaneous curvature model, 27,28 the bending energy of a membrane is given by Fig. 1 (a) A prolate vesicle (green-blue) in contact with uniform adhesive particles (grey) originating from the exterior solution. In this particular example, particles are completely engulfed at the strongly curved poles (green segments) and partially engulfed in the weakly curved equatorial region (blue). Partially engulfed particles experience a curvature-induced force towards regions of lower membrane curvature, whereas completely engulfed particles experience no such force. (b) For the same prolate vesicle, Janus particles with one strongly adhesive (grey) and one nonadhesive (red) surface domain are partially engulfed everywhere on the membrane, and therefore always experience curvature-induced forces towards regions of lower membrane curvature as indicated by the arrows.
This journal is © The Royal Society of Chemistry 2017 Soft Matter, 2017, 13, 2155--2173 | 2157 E be ¼ Ð dA2kðM À mÞ 2 , where the integration runs over the whole surface of the membrane or vesicle. The mean curvature of the membrane is defined as M = (C 1 + C 2 )/2, where C 1 and C 2 are the principal curvatures, and is taken to be positive or negative if the membrane bulges towards the exterior or interior compartment of the vesicle, respectively. Because we do not consider transformations that change the topology of the membrane, the contribution of the Gaussian curvature to the bending energy is constant and can be omitted. The elastic properties of the membrane are then described by its bending rigidity k and its spontaneous curvature m. The attractive interaction with the particle is included via the adhesive energy E ad = À|W|A bo , where |W| is the absolute value of the adhesive energy per unit area and A bo the area of the membrane that is bound to the particle. 29 The total energy of the system is then given by E = E be + E ad .
We will focus on the engulfment of rigid spherical particles with radius R pa by a vesicle with fixed total membrane area A and enclosed volume V. Unless otherwise stated, we will consider the case of endocytic engulfment, for which the particles originate from the exterior solution outside the vesicle; the exocytic case of particles originating from the interior of the vesicle is treated in Appendix B. With the exception of Section 6, we will focus on chemically uniform particle surfaces that are characterized by a spatially independent adhesive strength |W|. In Section 6, we will consider Janus particles with two distinct surface domains, only one of which is adhesive.

Relevant parameters for the vesicle-particle system
The vesicle-particle system has a typical geometry as shown in Fig. 2. In order to minimize the free energy of the system, one can decompose the membrane into a bound and an unbound segment. The bound membrane segment (green) is in contact with the particle and extends up to the contact line, which may in principle have any shape. Irrespective of the latter shape, we can define the bound area fraction q A bo /4pR pa 2 (1) which is a measure for the degree of engulfment, varying from q = 0 for a free particle to q = 1 for a completely engulfed particle, and can be regarded as the reaction coordinate for the engulfment process. Intermediate q-values with 0 o q o 1 correspond to partially engulfed states. The total membrane area A is equal to the sum of the area A bo of the bound membrane segment and the area A un of the unbound segment.
The total membrane area A defines the vesicle size which will be used as the basic length scale. The equilibrium shape of a free vesicle with area A and enclosed volume V is then completely determined in terms of the rescaled spontaneous curvature % m mR ve (3) and the rescaled volume which varies from v = 1 for a spherical vesicle to v = 0 for a completely deflated vesicle. For given values of these two parameters, the equilibrium shape of an axisymmetric free vesicle can be directly obtained by the numerical shooting method developed in ref. 28. When the membrane is in contact with a chemically uniform nanoparticle of adhesive strength |W|, the deviation of the unbound membrane segment from the nanoparticle is governed by the contact mean curvature which increases both with increasing particle size R pa and increasing adhesive strength |W|. The quantity M co corresponds to the equilibrium mean curvature of the unbound membrane segment along the contact line, 7 both for axisymmetric 29 and non-axisymmetric 30 contact lines. This curvature, which can be positive or negative, represents a material parameter that depends on the adhesive strength |W| as determined by the membrane-particle interactions and, thus, by the surface chemistry of the particle, on the bending rigidity k which varies with the membrane composition, and on the particle size R pa . The expression (5) applies to the endocytic case in which the particles originate from the exterior solution. The contact mean curvature for the exocytic case is given by eqn (109) in Appendix B. Note that relation (5) for the contact mean curvature can be rewritten in the form with the rescaled contact mean curvature Fig. 2 Non-axisymmetric geometry of a spherical particle (grey) with radius R pa adhering to a prolate vesicle. The total energy of the system can be separated into contributions E bo and E un from the bound (green) and unbound (black) membrane segments. The total membrane area is A = A bo + A un . As a reaction coordinate for the engulfment process, we choose the area fraction q = A bo /4pR pa 2 of the membrane-covered particle surface. which will be useful in order to express the adhesion energy |W|R pa 2 in terms of the other parameters.

Energies of bound and unbound membrane segments
The total energy E of the vesicle-particle system can also be decomposed into the contributions of the bound and unbound membrane segments such that E = E bo + E un . The bound membrane segment follows the shape of the particle which implies that the energy E bo of the bound segment is given by the expression which is linear in the bound area fraction q. When we replace the term |W|R pa 2 in (8) via relation (6), the energy of the bound membrane segment becomes with the rescaled particle radius On the other hand, the unbound membrane segment contributes only the bending energy where the integral extends over the area A un = A À A bo of the unbound segment. In order to obtain the shape of the unbound segment that minimizes E un for a given value of the bound area fraction q, one typically has to perform a functional minimization of (11), which leads to the well-known Euler-Lagrange or shape equations for a free membrane. 27 These have to be solved numerically, with appropriate boundary conditions ensuring a smooth matching with the bound membrane segment, as well as the conservation of the vesicle area and volume, as was done in ref. 7 for the special case of axisymmetric geometries. For the general case of non-axisymmetric geometries, the shape equations cannot be easily solved and one then has to perform time-consuming, direct energy minimization of triangulated membranes. 25,31 However, we will see in the following sections that important information about E un (q) at the boundary values q = 0 and q = 1 can be obtained analytically for the latter case as well.

Stability of free and completely engulfed particles
In recent work, 7 we studied the engulfment of rigid, spherical particles by membranes with non-zero spontaneous curvature. Examination of our numerical results coupled to theoretical considerations led us to discover exact analytical conditions for the stability of free and completely engulfed particle states, corresponding to q = 0 and q = 1, respectively. These conditions involve the particle size, the adhesive strength, the membrane's bending rigidity and spontaneous curvature as well as two local curvatures of the membrane in contact with the particle. The two stability conditions can be expressed in terms of the derivatives dE/dq of the total energy E(q) at the boundary values q = 0 and q = 1 of the bound area fraction q. The stability limit L fr of the free particle state is given by with % M MR ve where M is the mean curvature of the free membrane segment at the point of initial contact with the particle. The free particle state is stable for % M co r % M and unstable for % M co 4 % M. Likewise, the stability limit L ce of the completely engulfed particle state has the form with % M 0 M 0 R ve where M 0 is the mean curvature of the mother vesicle at the position of the narrow neck that connects the mother vesicle to the completely engulfed particle. The completely engulfed state is stable for % M co + % M 0 Z 2 % m and unstable for % M co + % M 0 o 2 % m. It is important to note that the curvatures that enter in (12) and (13) have a different character: the contact mean curvature % M co represents a material parameter, as explained after (5), whereas the mean curvatures % M and % M 0 represent geometric quantities related to the shape of the vesicle membrane. Because the contact mean curvature % M co increases both with particle size and adhesive strength, the stability limits (12) and (13) reflect the fact that small or weakly adhesive particles will tend to remain free, whereas larger or more adhesive particles tend to be completely engulfed.
If the vesicle-particle system is axially symmetric, the contact line is a circle and one may use the wrapping angle f as the reaction coordinate which varies from f = 0 to f = p. The bound area fraction q is then given by q ¼ 1 2 ð1 À cos fÞ and the stability limits (12) and (13) are then equivalent to d 2 E/df 2 | f=0 = 0 and d 2 E/df 2 | f=p = 0 as used in ref. 7. When a membrane segment comes into contact with a spherical particle of radius R pa , its mean curvature M must satisfy the inequality M 4 À1/R pa in order to avoid intersections between the membrane segment and the particle surface. Likewise, when we consider a completely engulfed particle, the mean curvature M 0 of the mother vesicle at the neck must fulfill the inequality M 0 o 1/R pa in order to avoid intersections between the membrane of the mother membrane and the particle surface.

Energy of the unbound membrane segment
We will now show how the stability limits (12) and (13) lead to expressions for the derivatives of the energy E un (q) at the boundaries of the energy landscape as defined by the bound area fractions q = 0 and q = 1. Using expression (9) for the energy E bo of the bound membrane segment, the derivative of the total energy E = E bo + E un becomes According to stability relations (12) and (13), the first derivative dE(q)/dq vanishes at the boundary values q = 0 and q = 1 which allows us to express the corresponding derivatives dE un (q)/dq in terms of the other parameters. Furthermore, using the alternative form of stability relations (12) and (13), we can express the contact mean curvature % M co in terms of the other curvatures which leads to and These two expressions for the first derivatives of E un (q) at q = 0 and q = 1 are direct consequences of the stability relations and involve no additional approximations. Furthermore, the actual values of the energy of the unbound segment at the boundaries q = 0 and q = 1 can be directly related to the bending energy of free vesicles. Indeed, the case q = 0 trivially corresponds to the free vesicle, and therefore E un (q = 0) = E fr (18) where E fr is simply the energy of the free vesicle with total membrane area A, enclosed volume V, and spontaneous curvature m. As mentioned before, in our continuum approach the completely engulfed state represents a limit shape with an ideal neck which costs no energy. 28,32,33 The energy of the unbound membrane segment of the mother vesicle after complete engulfment can therefore be calculated as where E fr 0 is now the energy of an 'effective' free vesicle representing the mother vesicle, with total membrane area A À 4pR pa 2 , enclosed volume V + 4pR pa 3 /3, and spontaneous curvature m.
4 Particle size small compared to the vesicle size We will now focus on the case of particles that are much smaller than the vesicle and use the rescaled particle size % R pa R pa /R ve as an expansion parameter. In the following sections, we will expand the energies of the membrane-particle system up to first order in % R pa . We will also restrict ourselves to the case of small spontaneous curvatures, which are of the order of the inverse vesicle size, i.e., m B O(1/R ve ) or % m B O(1). In the absence of large spontaneous curvatures, the mean curvature at any point of a free vesicle will also be of the order of the inverse vesicle size, and we can therefore always take . In addition, we will be interested in the behavior of the system in the proximity of the stability limits (12) and (13), which implies that the contact mean curvature is also of the order of the inverse vesicle size, In this limit, the energy (9) of the bound membrane segment can be rewritten as As described in Appendix A, it is not difficult to show that, under the conditions described in the previous paragraph, the difference between the energy E fr of the free vesicle before engulfment and the energy E fr 0 of the mother vesicle after complete engulfment is of second order in the particle size, i.e., Likewise, the difference between the mean curvature of the membrane segment in contact with the particle before and after complete engulfment at any given point on the vesicle surface is also of second order in % R pa , with We can now expand eqn (15)- (19) up to first order in the particle size, which leads to and where the asymptotic equality % M 0 E % M has been used in the last equation. Therefore, to first order in the rescaled particle size % R pa = R pa /R ve , the properties of the energy function E un (q) at the boundaries depend only on those of the free vesicle before engulfment, namely on the quantities E fr and M. In addition, E un (q) attains a particularly symmetric form in the limit of (1) and It is important to note that the derivation of relations (15)- (19) and (23)-(26) did not involve any assumptions concerning the symmetry of the particle-vesicle geometry. In particular, the use of the bound area fraction q as a reaction coordinate makes no assumption about the particular shape of the contact line between the particle and the membrane. The simplest functional form of the energy E un (q) that satisfies all four asymptotic equalities (23)-(26) is provided by a parabola with an extremum at q = 1/2, as given by We have compared this analytical interpolation formula with full numerical calculations for a variety of axisymmetric vesicle shapes that differed in the rescaled volume v, spontaneous curvature % m, and local mean curvature % M at the poles, see Fig. 3. For each of these shapes, we considered a particle at one of the poles and varied the location of the contact line and, thus, the bound area fraction q, imposing the condition that the unbound and bound membrane segments join smoothly along the contact line. 33 In all cases studied, the analytical expression (27) was confirmed by the numerical results to first order in the particle-to-vesicle size, % R pa R pa /R ve , see Fig. 3. In Fig. 3(a), we show how the full numerical results (colored solid lines) approach the first-order term in (27) as we reduce the relative particle size and, thus, higher order correction terms B % R pa 2 .
These higher order corrections depend on the global geometry of the vesicle, as given by its rescaled volume v and the particular branch of solutions (prolate, oblate, etc.) to which the vesicle shape belongs, whereas the first order term depends only on the local curvature % M as demonstrated in Fig. 3(b). The evolution of the energy of the unbound segment as a function of the bound area fraction is strongly dependent on the sign of the difference % M À % m between the mean curvature of the membrane at the contact point and the spontaneous curvature. Indeed, it follows from (27) that the energy of the unbound segment at q = 1/2 has the value corresponding to an energy maximum for % M 4 % m but to an energy minimum for % M o % m. Whenever % M = % m, the energy of the unbound membrane segment is independent of q up to first order in the particle size. The particular case % M = % m = 0 is consistent with the known results for particles interacting with tensionless planar membranes with no spontaneous curvature, 15 in which case the unbound segment forms a catenoidal shape with E un (q) 0.
In the small particle limit with M 0 E M, the two nonintersection constraints mentioned above can be combined into Because the small particle limit considered here implies that (29) are automatically fulfilled in this limit and the resulting constraints on % M do not play any role in the following calculations. 5 Results for uniform particle surfaces

Local energy landscapes
In the following sections, we will focus on the first order terms proportional to % R pa , corresponding to small particle-to-vesicle size ratios, and we will ignore higher-order terms O( % R pa 2 ).
Combining expressions (20) and (27) for the energies of the bound and unbound membrane segments, we obtain an analytical expression for the local energy landscape of engulfment as defined by that has the explicit form or Fig. 3 (a) Numerically calculated energies of the unbound segment as a function of bound area fraction q, for the endocytic engulfment of a particle at the pole of a discocyte vesicle with rescaled volume v = 0.62 and no spontaneous curvature, for three different values of the rescaled particle size % R pa R pa /R ve , see the inset. The rescaling of the energy follows the analytical expression (27) with % M= À1.063. As the relative particle size becomes smaller, the numerical curves approach the black dashed line q(1 À q) as described by (27); and (b) numerically calculated energies of the unbound segment at q = 1/2 as a function of the rescaled particle size % R pa for the discocyte in (a) and two other vesicle shapes with the rescaled volume v, and spontaneous curvature % m as given in the inset. At their poles, the oblate and the prolate shapes have mean curvatures % M= 0.263 and 1.530, respectively. Negative values of the rescaled particle size correspond to exocytic engulfment; see Appendix B. At % R pa = 0 all three curves have a common tangent (black dashed line) corresponding to the first-order term in eqn (28). The latter expressions are indeed purely local in the sense that they depend only on the mean curvature M of the mother vesicle at the attachment point, and not on the overall shape of the vesicle. In addition to M, the local energy landscape also depends on the contact mean curvature M co and on the spontaneous curvature m. In the following sections, we will consider particles with a certain radius R pa exposed to complex-shaped vesicles. Both the contact mean curvature and the spontaneous curvature are then material parameters with fixed values, whereas the mean curvature M varies along the vesicle membrane. Even though the local energy landscape DE(q) as given by eqn (32) is a simple quadratic function of q, this landscape describes a large variety of different cases depending on the relative magnitude of the three curvatures M co , m, and M. In order to discuss these different cases in a transparent manner, it will be useful to first consider the extrema of DE(q) within the physically meaningful interval 0 r q r 1 and to introduce some notation for specific q-values.
The quadratic q-dependence of expression (32) implies that the energy can have at most one extremum in the physically meaningful q-interval with 0 r q r 1. Indeed, the condition dDEðqÞ dq q¼q Ã ¼ 0 leads to only one solution as given by which may or may not lie within the physically meaningful interval 0 r q r 1. The extremum is located at the lower boundary value corresponding to the free particle state, and at the upper boundary value corresponding to the completely engulfed state. Furthermore, the extremum at q = q* represents the minimum of a noninverted parabola for M o m and the maximum of an inverted parabola for m o M. In order to distinguish these two cases, we will use the notation and For the boundary case m = M, local energy landscape (32) becomes linear in q and does not exhibit any extremum. It is useful to note that the M-dependence of q* as given by (33) stems from the particle-induced energy change of the unbound membrane segment as described by (27). The equilibrium state of the particle is provided by the state of minimal energy with bound area fraction q = q eq . Because the physically meaningful q-values must lie within the interval 0 r q r 1, the equilibrium state is provided by a boundary minimum at q = 0 or q = 1 whenever q*o0 or q* 4 1. Thus, the equilibrium value of the bound area fraction is given by q eq = q fr = 0 (38) corresponding to the free particle state for two different cases: for a non-inverted parabola with M o m and q min * o 0 as well as for an inverted parabola with m o M and q max * 4 1. Likewise, we have q eq = q ce = 1 (39) corresponding to the completely engulfed state for a non-inverted parabola with M o m and q min * 4 1 as well as for an inverted parabola with m o M and q max * o 0. Furthermore, the equilibrium state corresponds to the bound area fraction corresponding to a partially engulfed state. Finally, if the local energy landscape exhibits a maximum in the physically meaningful q-interval, we have two boundary minima at q = q fr = 0 and q = q ce = 1 which are separated by an energy barrier at In the latter case, the two boundary minima have the energies DE(q = 0) = 0 and as follows from eqn (32)  The different cases are illustrated in Fig. 4 where we plot the rescaled energy landscapes which have the dimensionless form

Stability and engulfment regimes
Up to first order in the particle-to-vesicle size ratio % R pa , the stability limits L fr and L ce of the free and completely engulfed state as given by (12) and (13) can be written as M = M co and M = 2m À M co , respectively, because % . These stability limits can be directly recovered from expression (33) by imposing q* = q fr = 0 and q* = q ce = 1, see relations (34) and (35). In general, expression (31) for the local energy landscape DE(q) leads to four qualitatively distinct landscapes corresponding to four distinct stability or engulfment regimes which are separated by the stability limits L fr and L ce : This journal is © The Royal Society of Chemistry 2017 (i) Free or non-engulfment regime F st in which the free state at q = q fr = 0 is stable and the completely engulfed state at q = q ce = 1 is unstable. The function DE(q) as given by eqn (32) must then increase monotonically with increasing q within the physically meaningful interval 0 o q o 1. This monotonic q-dependence applies both to a parabola with its minimum at q = q min * r 0, which is obtained for M co r M o m, and for an inverted parabola with its maximum at q = q max * Z 1 as obtained for m o M r 2m À M co . These two curvature inequalities can be combined into the condition for the mean curvature M, including the case M = m for which the local energy landscape DE(q) is linear in q. It is useful to rearrange the two inequalities for the local mean curvature M as given by eqn (47) in order to obtain the upper bounds for the contact mean curvature M co . Examples of the energy landscapes of F-segments are displayed in Fig. 4 (ii) Complete engulfment regime C st in which the completely engulfed state at q = q ce = 1 is stable, whereas the free state at q = q fr = 0 is unstable. In this case, the function DE(q) in (32) must decrease monotonically with increasing q within the whole interval 0 o q o 1. The latter q-dependence applies both to an inverted parabola with its maximum at q = q max * r 0, which is obtained for m o M r M co , and to a parabola with its minimum at q = q min * Z 1 as obtained for 2m À M co r M o m. The two curvature inequalities can be combined into the condition 2m À M co r M r M co (C-segments) (iii) Bistable engulfment regime B st in which both the free state at q = q fr = 0 and the completely engulfed state at q = q ce = 1 are locally stable and separated by an intermediate energy barrier at q = q ba = q max *. Now, the function DE(q) in eqn (32) represents an inverted parabola with its maximum located within the interval 0 o q max * o 1. The parabola is inverted for M 4 m which implies that the inequalities q max * 4 0 and for the contact mean curvature M co ; and (iv) Partial engulfment regime P st in which both the free state at q = q fr = 0 and the completely engulfed state at q = q ce = 1 are unstable and the parabolic function DE(q) in eqn (32) for the contact mean curvature M co . When we specify the above curvature inequalities to the case of vanishing spontaneous curvature, m = 0, we conclude that concave membrane segments with M o À|M co | belong to the partial engulfment regime P st , whereas convex segments with M 4 +|M co | represent B-segments and exhibit bistable engulfment behavior.
The four engulfment regimes just described are visualized in Fig. 5. The upper panel displays these regimes within the (M co ,m)-plane for a fixed value of the local mean curvature M. The lower panel of Fig. 5 depicts the four regimes within the (M,m)-plane for a fixed value of the contact mean curvature M co . As previously mentioned, the local mean curvature M varies along the surface of a complex-shaped vesicle, whereas the spontaneous curvature m and the contact mean curvature M co represent material parameters (for fixed particle size R pa ).
The four stability and engulfment regimes F st , C st , B st , and P st have been previously derived in ref. 7 and 26, but this derivation was based on the numerical observation that, for particles much smaller than the vesicle, partially engulfed states are only stable if both free and completely engulfed states are unstable. Now, we see that this latter condition is a direct consequence of the energy landscape (32) as obtained to first order in the rescaled particle size % R pa = R pa /R ve . On the other hand, as we increase the magnitude of % R pa , deviations from the simple parabolic energy landscape become more and more apparent. In particular, numerical results show that new stability regimes appear, in which stable free and completely engulfed states can coexist with stable 'satellite' minima, corresponding to partially engulfed states with very low or very high bound area fraction, q \ 0 or q t 1. It is now clear that these deviations are a consequence of higher order terms in the rescaled particle size % R pa , which are expected to depend on the overall shape of the vesicle as well.

Binding energies and barrier heights
For partial engulfment, the energy landscape has a minimum at q* = q min * with 0 o q min * o 1 as shown in Fig. 4(a) and (f), corresponding to the curvature inequalities M o M co o 2m À M co and M o 2m À M co o M co , respectively. The binding energy of the partially engulfed particles is obtained by inserting expression (33) for the location q* = q min * of the energy minimum back into the energy landscape (31). As a result, we find the binding energy within the partial engulfment regime. For the special case M co = m, the free and the completely engulfed states have the same energy, DE(q = 1) = DE(q = 0) = 0 and q min Ã ¼ 1 2 as follows from eqn (32) and (33). In the latter situation, the binding energy attains the particularly simple form It is useful to note that the M-dependence of the binding energy DE(q min *) in (55) and (56) arises from the M-dependence of the bound area fraction q min * as given by (33) and, thus, from the particle-induced energy change of the unbound membrane segment as in (27). For bistable engulfment, the energy landscape has a maximum at q* = q max * with 0 o q max * o 1, as shown in Fig. 4(e) and (j), corresponding to the curvature inequalities M 4 2m À M co 4 M co and M 4 M co 4 2m À M co , respectively. By inserting eqn (33) with q* = q max * into eqn (31), we obtain the energy barrier height It follows from relationships (56) and (58) that the magnitude of the typical binding energy or barrier height is proportional to the particle size R pa and, thus, decreases with decreasing particle size. This dependence is important because the binding energies and barriers heights obtained here are only relevant as long as their magnitude is large compared to the thermal energy. For lipid bilayers with a negligible spontaneous curvature and a typical bending rigidity k = 20k B T, we have |DE(q* = 1/2)| 4 k B T as long as |M|R pa 4 k B T/4pk C 1/250 with the segment curvature |M| of the order of the inverse vesicle size 1/R ve . As a consequence, if the particle size R pa is smaller than about R ve /250, the binding energies and energy barriers will be 'blurred out' by thermal fluctuations.

Engulfment patterns
In ref. 26, we showed how the dependence of the engulfment behavior on the local curvature of the membrane leads to the formation of engulfment patterns when a complex-shaped vesicle is exposed to many particles. If a vesicle has nonuniform curvature, that is, if the value of M varies along the surface of the vesicle, the vesicle surface may be divided into F-, C-, Bor P-segments according to the local membrane curvature, depending on whether a particle coming into contact with the membrane at that point belongs to the free or nonengulfment regime F st , the complete engulfment regime C st , the bistable regime B st or the partial engulfment regime P st . One example of an engulfment pattern is shown in Fig. 6, corresponding to a discocyte vesicle which is both axially symmetric and up-down symmetric. Inspection of this figure shows that the corresponding engulfment pattern consists of two P-segments around the north and south poles, a B-segment or -belt around the equator, and two C-segments or -rings which connect the P-segments with the B-belt.
The shape of any axisymmetric vesicle can be parametrized by the arc length s of the shape contour. 28 For the discocyte vesicle shape in Fig. 6, the arc length varies from s = 0 at the north pole to s = L = 2.90R ve at the south pole. In order to eliminate one parameter, we now define the rescaled arc length The vesicle in Fig. 6 has the rescaled volume v = 0.62 and zero spontaneous curvature, m = 0, and the membrane-particle adhesion is characterized by the contact mean curvature % M co = 1/2. Because % M co 4 % m = 0, the general classification described in Section 5.2 implies that such a vesicle does not exhibit any F-segments and that the rescaled mean curvature % M = MR ve satisfies % M o À % M co = À1/2 for the two P-segments, (60) and Thus, the boundaries between the Pand the C-segments, which are located at s 1 and s 4 , are characterized by the mean curvature % M = À1/2. Likewise, the boundaries between the C-segments and the B-belt have mean curvature % M = +1/2, see Fig. 7(a). Furthermore, for the discocyte shape in Fig. 6, the local energy landscape as given by (32) has the explicit form which depends on the mean curvature % M = MR ve of the membrane segment. As we move along the membrane surface of the vesicle, the mean curvature % M changes and so does the Fig. 6 Engulfment pattern on a discocyte vesicle with rescaled volume v = 0.62 and zero spontaneous curvature, m = 0, when exposed to small particles with contact mean curvature M co = 1/(2R ve ). The axisymmetric vesicle displays two partially engulfed P-segments (blue) around the north and south poles, two completely engulfed C-segments (green), and one bistable B-segment (yellow) which forms a belt around the equator of the vesicle. local energy landscape DE(q) which attains its maximal or minimal value for the bound area fraction as follows from (55) and (33). The functional dependence of the rescaled segment curvature % M on the rescaled arc length s is displayed in Fig. 7(a). Both the north and the south pole are characterized by the negative mean curvature % M = À1.063 which implies that the local energy landscape (63) has the form of a non-inverted parabola with its minimum located at the bound area fraction q min * = 0.735 which is equal to the equilibrium area fractions q eq (s = 0) and q eq (s = 1) of the partially engulfed states at the two poles, see Fig. 7(b). As we move away from the two poles, the mean curvature % M increases towards zero and the minimum of the local energy landscape is shifted towards larger values of q min * = q eq (s) until we reach q min * = q eq (s 1 ) = q eq (s 4 ) = 1 corresponding to the mean curvature % M = À1/2. The membrane loci with the latter curvature define the boundaries between the Pand C-segments. Moving further towards the equator of the discocyte, the location q min * of the energy minimum satisfies q min * 4 1 and, thus, moves out of the physically meaningful interval 0 o q o 1. The particle then attains a completely engulfed state corresponding to the boundary minimum at q eq (s) = 1.
As we move across one of the two C-segments towards the equator, the mean curvature % M vanishes at certain arc length values, s 1 + z and s 4 À z with z = 0.06. For these membrane loci with % M = 0, the energy landscape DE(q) as given by (63) becomes a linear and monotonically decreasing function of q. Therefore, as we approach % M = 0 from negative % M-values, the location q min * of the energy minimum moves out to infinity. For % M 4 0, on the other hand, the local energy landscape has the form of an inverted parabola, and the location q max * of its maximum moves out to minus infinity as we approach % M = 0 from positive % M-values. With increasing % M 4 0, the maximum of the energy landscape shifts towards less negative values of q max * until we reach the boundary value q max * = 0 for membrane curvature M = 1/2. The membrane loci with the latter curvature define the boundaries between the two C-segments and the B-segment around the equator. For the latter segment, we have two (meta)stable particle states with area fractions q eq = 0 and q eq = 1 as well as an energy barrier at q ba = q max *, as shown in Fig. 7(b). The location q ba of the energy barrier reaches its largest value q ba = 0.375 at the equator where the membrane curvature has its largest value % M = 2.006, as follows from (65).

Global energy landscapes
Because the local energy landscape DE(q) as given by (32) depends explicitly on the mean curvature M, the variation of M along the membrane surface implies a corresponding variation in the local energy landscape and, thus, in the equilibrium value q eq of the bound area fraction. The latter variation is illustrated in Fig. 7(b) for the case of the discocyte shape in Fig. 6. Because of its axisymmetry, this shape can be parametrized by the rescaled arc length s and the engulfment pattern can then be expressed Fig. 7 Different quantities that determine the engulfment pattern of the discocyte shape in Fig. 6 as a function of the rescaled arc length s: (a) rescaled mean curvature % M of the vesicle membrane, (b) bound area fraction q eq , and (c) global energy landscape DE eq = DE(q eq ). The latter two quantities correspond to the equilibrium state of a particle in contact with the vesicle membrane at s. For the bistable engulfment regime with s 2 = 0.26 o s o s 3 = 0.74, the free particle states with q = q fr = 0 and DE fr DE(q fr ) are metastable and separated from the equilibrium states by energy barriers at q = q ba with barrier heights DE ba DE(q ba ). The color of the lines corresponds to the different membrane segments defined in Fig. 6. in terms of the function q = q eq (s) as in Fig. 7(b). We will now use the latter function to define the global energy landscape DE eq (s) DE(q eq (s)) (66) with the local energy landscape DE(q) as given by eqn (32). For the discocyte shape, the global energy landscape DE eq has the form displayed in Fig. 7(c) which implies that this vesicleparticle system attains its lowest energy state (or ground state) when the nanoparticle is partially engulfed at the north or south pole of the vesicle. Thus, for % m = 0 and % M co = 1/2, the energy of the completely engulfed states exceeds the energies of the partially engulfed states which implies that the energy increase of the unbound membrane segments overcompensates the energy decrease of the bound membrane segment as we move from a Pto a C-segment. Expression (67) also applies to the completely engulfed states of the B-segment around the equator. For the latter segment, the free states with q eq = 0 are metastable as well but their energy DE(q eq = 0) = 0 is always larger than the energy DE(q eq = 1) of the completely engulfed states. The two (meta)stable states are separated by an energy barrier of height DE(q ba ) as given by (64) which attains its largest value DE(q ba ) = 14.212k % R pa at the equator. When we defined the global energy landscape by eqn (66), we implicitly assumed that the bound area fraction q of the particle quickly adapts to its local equilibrium value q eq as the particle moves along the surface of the vesicle. For axisymmetric vesicles, we then obtain global energy landscapes that depend only on the rescaled arc length s, as illustrated in Fig. 7(c) for a discocyte shape. These one-dimensional energy landscapes can be generalized in two ways. First, we may consider an arbitrary, non-axisymmetric shape of the free vesicle which can be parametrized by a three-dimensional vector, -X = -X(s 1 ,s 2 ), that depends on two properly chosen surface coordinates, s 1 and s 2 . The latter coordinates determine both the local mean curvature, M = M(s 1 ,s 2 ), of the vesicle membrane and the position of the particle on this membrane. Second, we may vary the bound area fraction q of the particle independently of its position. As a consequence, the most general energy landscape depends on the three parameters s 1 , s 2 , and q and this three-dimensional landscape DE = DE(s 1 ,s 2 ,q) can be directly obtained from the local landscape as given by eqn (32) when we insert the local mean curvature M = M(s 1 ,s 2 ) into this latter equation.

Curvature-induced forces
The binding energy of a partially engulfed particle depends on the 'unperturbed' curvature M of the adjacent membrane segment, as described by eqn (55) and illustrated by the blue segments of the global energy landscape in Fig. 7(c). Because the binding energy (55) decreases with decreasing curvature M, the particle experiences a curvature-induced force that pulls it towards points of lower mean curvature of the vesicle membrane. For an axisymmetric shape, the curvature-induced force is given bỹ where t s is the unit tangent vector of the shape contour at arc length s = Ls. The global energy landscape DE eq (s) depends on the rescaled arc length s only via the bound area fraction q eq which varies within a P-segment but is constant within all other segments, as illustrated in Fig. 7(b). The latter property is true in general because F-, C-, or B-segments are characterized by boundary minima with q eq = 0 or q eq = 1 and the corresponding energies DE(0) = 0 or DE(1) = 16pkR pa (m À M co ) do not depend on s. On the other hand, within a P-segment, the s-dependence of q eq arises only from the s-dependence of the local mean curvature M = M(s) as described by eqn (55). Therefore, the gradient of the global energy landscape has the value @DE eq ðsÞ @s ¼ 0 for F-; C-; or B-segments (69) and the s-dependent form @DE eq ðsÞ @s ¼ @M @s @DE eq @M ¼ @M @s CðMÞ for P-segments (70) with where the latter expression depends on s via M = M(s), see the blue segments in Fig. 7(a). Note that C(M) is always positive because P-segments are characterized by the curvature inequalities M o M co o 2m À M according to eqn (54). A combination of relations (68) and (70) then implies that a partially engulfed particle experiences the curvature-induced forcẽ on an axisymmetric vesicle. For an arbitrarily shaped vesicle, the curvature-induced force has the form with the surface gradient operator r S . When we parametrize the vesicle shape by the vector-valued function -X(s 1 ,s 2 ), we obtain the two tangent vectors -X i q -X/qs i with i = 1,2 that span the plane tangential to the vesicle surface, the metric tensor g ij = -X i Á -X j , and the inverse metric tensor (g ij ) (g ij ) À1 . The action of the gradient operator r S onto any scalar function U(s 1 ,s 2 ) is then given by 34 with an implicit summation over repeated indices. In response to the curvature-induced force f M , as given by eqn (72) or (73), partially engulfed particles undergo biased diffusion towards membrane segments with lower curvature, whereas completely engulfed particles do not experience any curvature-induced force and diffuse freely on Band C-segments. However, if a completely engulfed particle diffuses from a Bor C-segment onto an For P-segment, it will unbind or partially unwrap from the membrane. Once on the P-segment, the particle will be pulled towards a membrane segment of minimal mean curvature. At temperature T, the probability to find the particle at a certain position, -X (s 1 ,s 2 ), on the membrane is then proportional to the Boltzmann factor exp[ÀDE eq ( -X)/k B T].

Results for Janus particles
For particles with a chemically uniform surface, partial engulfment is only found if the particle size and adhesiveness, which determine the contact mean curvature M co , are finely tuned so that M o M co o 2m À M as in (54). To soften this constraint and to enlarge the P st -regime, we will now consider Janus particles with surfaces that are divided into one adhesive and one non-adhesive domain. The adhesive surface domain is taken to form a spherical cap with area A 1 which implies a circular domain boundary of the adhesive surface domain and leads to the adhesive surface area fraction Because the area fraction q of the bound membrane cannot exceed the area fraction of the adhesive surface domain, the accessible q-range is now further restricted to 0 o q r Q 1 (Janus particles).
The adhesive domain is taken to have the adhesive strength |W 1 | and the contact mean curvature which grows as M co;1 $ ffiffiffiffiffiffiffiffiffi ffi W 1 j j p for large |W 1 |. A membrane segment with local mean curvature M will start to spread over such an adhesive surface domain provided the free state is unstable which implies dDE/dq| q=0 o 0 or M co,1 4 M ( Janus, onset of adhesion) as follows from the local energy landscape DE(q) in eqn (32) with M co replaced by M co,1 . Inspection of Fig. 4 shows that the inequality dDE/dq| q=0 o 0 is satisfied within the partial and complete engulfment regimes of the uniform particle with |W| = |W 1 |, as illustrated in Fig. 4(a) and (f)-(h). The same conclusion follows from Fig. 5(a). Because the Janus particle is characterized by the restricted q-range in eqn (76), the equilibrium value q = q JP eq of the membrane's bound area fraction on the Janus particle is given by and by q JP eq = q eq for q eq o Q 1 where q eq denotes the equilibrium value of the bound area fraction of the uniform particle as before. We have tacitly assumed here that the contact line adapts to the domain boundary of the adhesive surface domain and, thus, attains a circular shape as well. In principle, the equilibrium shape of the contact line could be non-circular as for partially engulfed particles with a uniform surface. Therefore, we assume here that the energy cost of forcing the contact line to be circular is of higher order in R pa /R ve . Furthermore, irrespective of the detailed geometry of the contact line, we conclude that a Janus particle with Q 1 o 1 is always partially engulfed by a membrane segment if its mean curvature M o M co,1 as in eqn (78) and the membrane starts to spread on the adhesive surface domain of the particle.
If the uniform particle with |W| |W 1 | is completely engulfed as in Fig. 4(g) and (h), the bound area fraction is q eq = 1, corresponding to a boundary minimum at q = 1. Therefore, in this case, the Janus particle is characterized by the bound area fraction q JP eq = Q 1 . On the other hand, if the uniform particle with |W| |W 1 | is partially engulfed as in Fig. 4(a) and (f), the bound area fraction q eq = q min * o 1. In the latter case, the membrane will still cover the whole adhesive surface domain provided the area fraction Q 1 of this domain satisfies the condition Therefore, a membrane segment with local mean curvature M starts to spread onto the adhesive surface domain for M co,1 4 M as in eqn (78) and spreads over the whole domain for M co,1 Z U(M,Q 1 ) as in eqn (82). A combination of these two relationships implies that the adhesive surface domain is fully covered by the bound membrane segment with bound area fraction whereas this surface domain is only partially covered with The latter case is only possible for M o m corresponding to partial engulfment of the uniform particle with |W| = |W 1 |, compare Fig. 5(a). We now focus on full coverage of the adhesive surface domain as described by eqn (83) and derive global conditions that ensure such a full coverage by all membrane segments of the vesicle. The function U(M,Q 1 ), defined in eqn (82), monotonically increases and decreases with increasing M for and where M max and M min are the maximal and minimal mean curvatures of the complex-shaped vesicle under consideration.
In the following sections, we will consider Janus particles with a sufficiently large contact mean curvature M co;1 $ ffiffiffiffiffiffiffiffiffi ffi W 1 j j p so that one of the conditions (85)-(87), corresponding to the chosen value of Q 1 , is satisfied. The Janus particles are then always partially engulfed with the bound membrane segment covering the whole area A 1 of the adhesive surface domain, and the equilibrium value of the bound area fraction is given by q JP eq = Q 1 irrespective of the location of the Janus particle on the vesicle membrane.
Using the explicit expression (32) for the local energy landscape DE(q), the binding energy of such a Janus particle is then given by which becomes DE JP eq = 4pkR pa (M + m À 2M co,1 ) for Q 1 = 1/2 (90) corresponding to an adhesive hemisphere.
Using the same line of reasoning that led to eqn (72) for uniform particles, the curvature-induced force acting on a partially engulfed Janus particle now has the form for an axisymmetric vesicle with the M-independent factor For an arbitrarily shaped vesicle, we now obtain as follows from eqn (73) and (92). Thus, the Janus particle is also pulled towards regions of lower membrane curvature. It is interesting to note that, in contrast to the case of uniform particles as described by eqn (72), the curvature-induced force acting on the Janus particle with q = Q 1 does not depend on the contact mean curvature M co and depends only implicitly on the spontaneous curvature m via the mean curvature M. The curvature-dependent binding energy (89) and the curvatureinduced force as given by eqn (91) or (93) have far-reaching consequences. Janus particles in contact with membranes of non-spherical shape will always 'feel' the local mean curvature gradient and tend to diffuse towards the points of minimal mean curvature. Biological objects involving complex-shaped membranes include intracellular organelles such as the endoplasmic reticulum, 35 dividing cells, red blood cells, and nascent autophagosomes. The shapes of the latter three examples are analogous to prolate, discocyte and stomatocyte vesicles, as shown in Fig. 8(a)-(c). The corresponding energy landscapes experienced by Janus particles are displayed in Fig. 8(d). In the latter figure, we plot the binding energy difference defined by in terms of rescaled variables and is, thus, directly proportional to the mean curvature difference % M À % M min . The rescaled arc length s is again defined by eqn (59), where L denotes the total contour length of the respective vesicle shape. Thus, the parameter s varies from s = 0 at the north pole to s = 1 at the south pole for all vesicle shapes displayed in Fig. 8(a)-(c).
For prolates and discocytes, which are up-down symmetric, particles tend to localize at the equator and at the poles, respectively. For a prolate vesicle, particles can freely diffuse along the equatorial line because all points on this line have the same curvature. For a discocyte vesicle, particles will localize either at the north or at the south pole. In the latter case, the diffusion from one pole to the other is impeded by a large energy barrier located along the equatorial line. These results are consistent with the collective behavior of particles in contact with prolate and discocyte vesicles, as described in ref. 9.
The behavior of Janus particles on stomatocyte-shaped vesicles, which have no up-down symmetry, is more complex. As shown in Fig. 8(c), a stomatocyte shape consists of an outer and an inner membrane segment that are connected by a membrane neck. The inner segment corresponds to a membrane invagination or in-bud. The associated energy landscape DDE JP eq now exhibits two distinct minima: (i) a local, metastable minimum located at the convexly curved north pole; and (ii) a globally stable minimum located at the concavely curved south pole which also represents the pole of the invagination. The energy barrier between these two minima is located at s = s sto C 0.55, right outside of the membrane neck, as indicated by the black diamond in Fig. 8(c). In order to reach the inner membrane segment, a particle attached to the outer segment will have to overcome this energy barrier, whose height is approximately 2k % R pa for a stomatocyte with rescaled volume v = 0.4, as displayed in Fig. 8(c). This implies that particles can spontaneously overcome this barrier via thermal diffusion only if they are sufficiently small. A crossover value % R † pa for the particle radius can be estimated from the requirement that 2k % R † pa C k B T which implies % R † pa C 1/40, assuming a typical bending rigidity for lipid membranes k = 20k B T. Small particles with radii % R pa t % R † pa will not encounter any significant barrier, whereas large particles with % R pa c % R † pa will be trapped in the outer membrane segment for a long time. As soon as a Janus particle has overcome the energy barrier close to the membrane neck, it will be strongly pulled through the neck region towards the central part of the invagination as follows from the strong decay of the stomatocyte's energy landscape DDE JP eq for s 4 s sto , as depicted in Fig. 8(d).

Summary and conclusions
In this paper, we derived an analytical theory for the engulfment of rigid spherical nanoparticles by vesicle membranes that allowed us to calculate local energy landscapes for such particles at any point of the vesicle membrane. More precisely, we used the bound area fraction q of the particle surface, see eqn (1), as our basic reaction coordinate and considered the limit in which the particle size R pa is small compared to the vesicle size R ve . Up to first order in the size ratio R pa /R ve , we found the explicit form (32) for the local energy landscape DE(q) of a uniform particle in contact with a membrane segment of mean curvature M. The most important feature of this energy landscape is that it does not depend on the overall shape of the vesicle.
Another remarkable feature of the local energy landscape DE(q) as given by eqn (32) is that it has a simple quadratic dependence on the bound area fraction q, as confirmed by extensive numerical calculations, see Fig. 3. Furthermore, this simple functional form provides a complete classification for the complex engulfment behavior that arises from the interplay of mean curvature M, spontaneous curvature m, and contact mean curvature M co , as displayed in Fig. 4 and discussed in Section 5.2.
The local energy landscape as given by eqn (32) is intimately related to our previous results, described in ref. 7 and 26. In ref. 7, we made the empirical (numerical) observation that small particles interacting with vesicles exhibit four distinct engulfment regimes: a free regime F st in which only the free particle state is stable, a completely engulfed regime C st in which only the completely engulfed state is stable, a partially engulfed regime P st in which only the partially engulfed state is stable, and a bistable regime B st in which both the free and the completely engulfed states are stable and separated by an energy barrier. We then used this empirical observation to describe the interaction of nanoparticles with complex-shaped vesicles in ref. 26. Because of the local energy landscape derived in the present paper, see eqn (32) and Fig. 4, we now understand that these previous results directly follow from the energetics of the system up to first order in R pa /R ve . Deviations from this behavior will always be a signature of higher order terms in R pa /R ve , which are expected to depend on the global shape of the vesicle.
The physically meaningful part of the local energy landscape DE(q) is provided by all q-values within the interval 0 r q r 1. When we restrict the energy landscape to this q-range, we always find a state of lowest energy which determines the particle's equilibrium state with bound area fraction q = q eq . For partially engulfed states on P-segments, the equilibrium state corresponds to a smooth minimum with q eq = q min * and 0 o q min *o1. For all other engulfment regimes, the equilibrium state is provided by a boundary minimum with q eq = q fr = 0 or q eq = q ce = 1. As we move along the membrane of a complexshaped vesicle, the membrane's mean curvature M changes in a continuous manner and so does the local energy landscape as well as the associated equilibrium state q eq . We can then use the energy DE(q eq ) of the equilibrium particle states to define a global energy landscape DE eq for all points on the vesicle membrane. If the vesicle is axisymmetric, its shape can be parametrized by the rescaled arc length s of the shape contour, see eqn (59), and the global energy landscape is obtained via DE eq (s) = DE(q eq (s)) as in eqn (66). One example is provided by the discocyte shape in Fig. 6 with its global energy landscape depicted in Fig. 7(c). For F-, Cor B-segments of the vesicle membrane, the equilibrium states of the particles correspond to the boundary minima of the local energy landscapes with q eq = q fr = 0 or q eq = q fr = 0. The corresponding energies DE(0) = 0 and DE(1) = 16pk(m À M co )R pa as in eqn (42) do not depend on the local mean curvature M of the segments. As a consequence, the global energy landscape DE eq = const for these three types of membrane segments, see Fig. 7(c), and, thus, does not lead to any curvature-induced forces on the particles.
For partially engulfed particles on P-segments, on the other hand, the equilibrium states are provided by the smooth minima of the local energy landscapes with q eq = q min * and 0 o q min * o 1. The latter minima depend on the local mean curvature M as described by eqn (55) and, thus, change as the particle moves along a P-segment with variable mean curvature M. When the particles originate from the exterior solution (endocytosis), the binding energy DE eq = DE* of the partially engulfed particles is reduced when they move towards membrane segments with lower mean curvature M. In contrast, partially engulfed particles that originate from the interior solution (exocytosis), reduce their binding energy when they move towards membrane segments with higher mean curvature M, see Appendix B. The corresponding curvature-induced forces acting onto particles with a chemically uniform surface have the explicit form as given by eqn (72) for axisymmetric vesicles and by eqn (73) for arbitrarily shaped membranes.
These curvature-induced forces arise from the M-dependence of the binding energy in eqn (55), which can be traced back to the M-dependence of the particle-induced energy change of the unbound membrane segment, as described by eqn (27). Furthermore, a decomposition of the binding energy (55) into separate contributions from the bound and unbound segments shows that the curvature-induced forces are dominated by the particle-induced energy changes of the unbound segments. Therefore, an endocytic (exocytic) particle moves towards a membrane segment of lower (higher) mean curvature in order to lower the bending energy of the unbound membrane segment in the vicinity of the contact line.
For uniform particles, partial engulfment is only possible if the contact mean curvature M co , which depends on the particle size and adhesiveness, is finely tuned with the local mean curvature M and the spontaneous curvature m of the membrane, see eqn (54). In contrast, Janus particles with one strongly adhesive and one non-adhesive surface domain are predicted to be always partially engulfed if the contact mean curvature M co,1 of their strongly adhesive domain satisfies the relevant condition in (85)-(87). The corresponding curvature-induced forces as given by eqn (91) and (93) again pull the Janus particles towards points of minimal membrane curvature in the endocytic case, as illustrated by the energy landscapes in Fig. 8, or towards points of maximal membrane curvature in the exocytic case, as described in Appendix B.
Strictly speaking, the results described here were obtained for a single particle interacting with a vesicle membrane. However, these results should also apply to a finite concentration of such particles provided one can ignore the interactions between the membrane-bound particles. The latter situations is provided by dilute particle concentrations on F-, Cor B-segments of the vesicle membrane. On the other hand, if several partially engulfed particles are present on the vesicle surface, membranemediated particle-particle interactions are expected to be important when the particles come close to each other, 36,37 and particles may cluster into tubes 31 or buds. 38 The latter process provides an example of domain-induced budding in which the domains acquire a spontaneous curvature arising from the asymmetry of the membrane-bound particles. 39,40 The curvature-induced forces derived here for partially engulfed particles are conceptually similar to the capillary forces that have been recently described for particles at fluid interfaces [41][42][43][44] and tense membranes 45 of arbitrary shape. The underlying mechanisms are very different, however: in the latter systems, the forces are a consequence of interfacial or membrane tension and of the changes in interfacial or membrane area caused by the presence of the particle. In contrast, for deflated vesicles as considered here, the total membrane area is constant and tension plays no role, which implies that the curvature-induced forces described here are a consequence of pure bending.
Nonspherical vesicles with rescaled volume v o 1 undergo thermally-excited shape fluctuations that lead to deviations from the equilibrium shapes depicted in Fig. 6 and 8. One may then wonder whether these fluctuations can affect the engulfment behavior. The engulfment of a particle takes up excess membrane area or, equivalently, increases the rescaled volume v of the vesicle, as described by eqn (99) in Appendix A. Particle engulfment therefore acts to suppress membrane fluctuations, and this fluctuation effect implies a corresponding entropic contribution to the free energy of the system. However, entropic contributions to the free energy are proportional to the thermal energy k B T, whereas the energy contributions discussed here are proportional to the membrane's bending rigidity k, which has a typical value of 20k B T. In particular, the entropic contributions to the local energy landscapes should have the form Ck B T % R pa up to first order in the rescaled particle size % R pa , with C being a dimensionless coefficient of order one. In contrast, the local energy landscape as given by eqn (32) is proportional to 16pk % R pa C 1000k B T % R pa . Therefore, entropic contributions to this local energy landscape can be safely ignored in typical cases.
An important feature of the analytical theory presented here is that it can be used to describe particle engulfment at nonaxisymmetric locations of vesicles. Usually, the study of nonaxisymmetric membrane geometries requires numerical energy minimization of dynamically triangulated membranes, which is computationally rather costly. 25,31 It will certainly be useful to