Mutual remodeling of interacting nanodroplets and vesicles

Liquid-liquid phase separation within the cytosol leads to the formation of protein-enriched droplets inside cells. These droplets known as biomolecular condensates have ultra-low interfacial tensions and fulfill a vast range of functions inside cells. Biomolecular condensation can take place at the plasma membrane and generate mechanical forces on membranes as a result of membrane wetting. But little is known about the wetting of membranes by biomolecular condensates. In this study, we utilize energy minimization to explore a wide range of parameters and determine the dependence of membrane wetting phenomena on interfacial tension, bending rigidity, line tension, and spontaneous curvature. We observe that interacting nanodroplets and vesicles mutually remodel one another. In addition, we determine the parameter regimes for which the droplet-membrane systems exhibit axisymmetric and non-axisymmetric contact lines. Our results provide insights into understanding intracellular processes and physical mechanisms based on the mutual remodeling of droplets and membranes. Biomolecular condensates are membrane-less organelles performing various functions inside cells which behaviour can be understood in terms of liquid-liquid phase separation and wetting. Here, the authors characterize the low interfacial tension regime of nanodroplets during endocytic and exocytic engulfment within an elastic membrane, study the role of the contact line symmetry, and show that nanodroplets and vesicles mutually remodel one another

W etting is a ubiquitous phenomenon in nature with implications ranging from daily-life examples to industrial applications 1 . Moreover, wetting is essential to understand the physicochemical properties of interfaces 1,2 . The typical example is the wetting of a rigid substrate, which can be described via Young's equation as a force balance along the contact line 1,2 . The interaction of liquid droplets with flexible substrates such as lipid membranes 3,4 , is a generalized example of wetting, which represents a relatively unexplored research field. One example is provided by the liquid-liquid phase separation (LLPS) of polymeric solutions within lipid vesicles 5 . Such aqueous two-phase systems enclosed in lipid vesicles lead to a wide range of phenomena, such as complete wetting to partial wetting transition, vesicle budding, and membrane tubulation 3,6 .
Analogous wetting phenomena are also relevant in the context of cell biology. LLPS of protein solutions within the cytosol leads to the formation of biomolecular condensates (also known as membrane-less organelles), which behave like liquid droplets 7 . Such biomolecular condensates have a very low interfacial tensions in the range of 7-9 1 − 100 μN m −1 , similar to the interfacial tension of polymeric condensates 10 . LLPS in vivo is expected to have many functional consequences 11 , which are not yet fully understood. Biomolecular condensates can interact with distinct biomolecular entities inside cells. For instance, biomolecular condensates inside the cell nucleus have been shown to change the structure of the genome by forming a capillary bridge among chromatin fibers 12 . In addition, protein condensation is observed to generate mechanical forces on the membrane and cause membrane remodeling and endocytosis 13 .
Furthermore, a recent study suggests an active involvement of protein phase separation in the formation of cell-cell tight junctions 14 . Likewise, due to the interaction of biomolecular condensates with membranes, a cluster of synaptic vesicles has been observed to be relocated inside phase-separated synapsin protein 15 . Such synaptic vesicles eventually are drawn for exocytosis in synapses to deliver neurotransmitters 15 . All these examples can be better understood in the context of membranewetting phenomena.
A theory of membrane wetting has been developed based on the membrane elasticity and variational calculus 3,16,17 . The obtained force balance is a generalized form of Young's equation which couples the elasticity of the membrane to the capillary force arising from the liquid-liquid interface. The membranedroplet force balance links the mechanical segment tensions to the intrinsic contact angle of the droplet, the interfacial tension, the three-phase line tension, as well as the spontaneous curvature of membrane segments. On the micron scale, the effect of line tension can be neglected 3,6,16,17 . However, at the nanometer scale, the line tension makes a significant contribution to the total energy of the system 4 , as demonstrated by molecular dynamics (MD) simulations. For a large range of parameters, the line tension is found to be negative in sign and leads to spontaneous symmetry breaking of the contact line shapes, which eventually creates a tight-lipped membrane neck 4 . The system studied by MD simulations is characterized by high interfacial tension, with magnitudes that correspond to those of oil-water interfaces, which are of the order of mN m −1 . These tension values are three orders of magnitude larger than those observed for biomolecular condensates. MD simulations become, however, unfeasible when the interfacial tension becomes too low. Indeed, as we decrease the interfacial tension, the liquid-liquid interface becomes rather diffuse, and we would need to simulate larger droplets to obtain reliable results for which MD simulations are computationally very expensive to perform.
In this work, we explore the low interfacial tension regime through an alternative approach based on energy minimization, which is computationally less expensive. We examine the endocytic and exocytic engulfment of liquid nanodroplets with low interfacial tension interacting with elastic vesicle membranes. We explore the effect of nanodroplet size, interfacial tension, spontaneous curvature, and line tension on membrane-wetting phenomena and nanodroplet-vesicle mutual remodeling. In addition, we study how the fluid-elastic parameters affect the shape transformation of the vesicle and the wetting energy. Finally, we determine the boundary between axisymmetric and non-axisymmetric contact line geometries within the three-dimensional parameter space obtained for vanishing spontaneous curvature.

Results and discussion
Mutual remodeling of vesicles and nanodroplets. We studied the wetting of a vesicle with radius R ve = 1μm and volume-toarea ratio of vesicle v = 0.7 (as defined in Supplementary Note 1 and Supplementary Equation (2)) interacting with nanodroplets of variable size. The vesicles can interact with nanodroplets formed in the exterior solution which is equivalent to the scenario of droplet endocytosis (see Fig. 1a) and with nanodroplets formed inside the vesicle, which is equivalent to the scenario of droplet exocytosis (see Fig. 1b). For α droplet endocytosis, the vesicle interior liquid phase is called γ and the exterior phase is called β. For α droplet exocytosis, the vesicle interior liquid phase is called β and the exterior phase is called γ (Supplementary Figure 1). In a non-wetting state, the liquid nanodroplet attained a spherical geometry, as expected for a droplet in equilibrium. The vesicle morphology, however, transformed to a prolate shape for the volume-to-area ratio of v = 0.7, (leftmost shape in Fig. 1), for which the prolate morphology is the minimum energy shape. Such a shape transformation is known to occur when we vary the volume v at zero spontaneous curvature 18 (see Supplementary Note 2). In the wetting state, we initialized the minimization process where the nanodroplet partially wetted the vesicle and then we minimized the total energy of the nanodroplet-vesicle systems (see Supplementary Note 3). All nanodroplet-vesicle morphologies shown in all figures represent minimum energy shapes. Even though different initialization can lead to several coexisting stable and meta-stable solutions (see Supplementary Notes 4 and 5).
In response to the capillary force arising from the αβ interface of nanodroplets, the vesicle deforms and consequently reduces the interfacial energy of the nanodroplet, until the bending energy counteracts the capillary forces and prevents further vesicle deformation. This interplay between interfacial and bending energies leads to mutual remodeling of both nanodroplets and vesicles. For droplet endocytosis (see Fig. 1a), and relatively small nanodroplet sizes, r α = R α /R ve = 0.1, nanodroplets wet the vesicle at the equator which leads to a gain in interfacial energy overcompensating the increased cost of bending energy. If the size of the nanodroplet is increased by a factor of two, to r α = 0.2, wetting induces a large morphological transformation by changing the vesicle morphology from prolate to oblate. By further increasing the nanodroplet size, to r α = 0.3, the vesicle engulfs the nanodroplet and back-transforms to the prolate shape with the droplet being engulfed by a spherical in-bud at the equator. If the nanodroplet size is further increased, beyond r α > 0.4, the vesicle cannot fully engulf the nanodroplet and it takes the form of a stomatocyte. For nanodroplet exocytosis (see Fig. 1b), nanodroplets wet the vesicle at the poles which enables the system to reduce its interfacial energy to a large extent without causing much of a penalty from the bending energy. In the latter case, the nanodroplet forms an out-bud, the size of which is determined by the size of the droplet, and the daughter vesicle keeps its original prolate shape. For a nanodroplet size of r α = 0.7, the vesicle takes the shape of a dumbbell. However, for large nanodroplet sizes, the wetting-induced out-bud takes the form of a prolate which connects to the daughter vesicle via its pole. Our result demonstrates that a local perturbation of a vesicle by a nanodroplet can propagate over the whole membrane geometry and globally change the shapes of interacting nanodroplets and vesicles. Similar shape transformations of vesicles have been reported for the system of interacting nanoparticles and vesicles 19 as well as the adhesion of the vesicle to flat surfaces 20 . The wetting-induced shape transformation is the result of global energy minimization, where the cost of vesicle deformation is overcompensated by the gain of interfacial energy. The endocytic and exocytic processes of growing nanodroplets as in Fig. 1 can be regarded as the heterogeneous nucleation and growth of nanodroplets inside and outside of the model biomembranes. Such processes can be readily performed in an in-vitro experimental setup.
Morphology diagram after complete engulfment. The morphology of the vesicle after complete engulfment of the nanodroplet is determined by the volume-to-area ratio of the two daughter vesicles, where the first subcompartment engulfs the nanodroplet and the second subcompartment encloses the second aqueous phase. For droplet endocytosis, the vesicle volume that contains the γ phase effectively increases by the nanodroplet volume while its membrane area decreases by the surface area of the nanodroplet, see Supplementary Equation (3). The endocytosis volume-to-area ratio of enclosing γ phase v γ is described by the Supplementary Equation (3) after complete engulfment of the nanodroplet is shown in Fig. 2a, and the volume-to-area ratio of nanodroplet v α is shown in Fig. 2b. Energy of the complete engulfment is calculated from the morphology diagram of the vesicle at zero spontaneous curvature 21 , as the sum of bending energy of two daughter vesicles i.e., (E αγ + E βγ )/8πκ (Fig. 2c). The corresponding nanodroplet-vesicle morphology examples are illustrated in Fig. 2g. For nanodroplet exocytosis, the vesicle volume, which contains β and α phases, effectively decreases by the nanodroplet volume, and its membrane area decreases by the surface area of the nanodroplet, see Supplementary Equation (4). The exocytosis volume-to-area ratio of enclosing β phase v β is described by the Supplementary Equation (4) after complete engulfment of the nanodroplet is shown in Fig. 2d, and the volume-to-area ratio of nanodroplet v α is shown in Fig. 2e. Energy of the complete engulfment is calculated from the morphology diagram of the vesicle at zero spontaneous curvature 21 , as the sum of bending energy of two daughter vesicles i.e., (E αγ + E βγ )/8πκ (Fig. 2f). The corresponding nanodroplet-vesicle morphology examples are shown in Fig. 2h. In the dark gray regions, only partial engulfment is possible because the area of the vesicle membrane is not sufficiently large, however, in the light gray region engulfment is not possible because the vesicle volume is too small (Fig. 2). The yellow, blue, green, and red lines in the morphology diagram (Fig. 2a, d, e), correspond to the volume-toarea ratio of v i = 0.5915, 0.6515, 0.98 and 1.0 (i = γ for nanodroplet endocytosis and i = β for nanodroplet exocytosis). These volume-to-area ratio isolines set the boundary between minimum energy shapes 21 (Supplementary Figure 2), where a prolate is the global minimum for the volume range 0.6515 < v i < 1.0, oblate in the range of 0.5915 < v i < 0.6515, and stomatocytes for v i < 0.5915, respectively.
The liquid nature of the nanodroplets does not set any area constraint on the shape of the daughter vesicles after complete engulfment of the nanodroplet, see Eq. (1). This provides the freedom that membrane area segments A αγ and A βγ are partitioned in a way so that the total bending energy of the two daughter vesicles reaches a minimum, under the constraint of constant vesicle area A ve = A αγ + A βγ and constant nanodroplet volume and constant vesicle volume Eq. (1). We used this feature of the nanodroplets to calculate the area A αγ which corresponds to the minimum energy morphology for the two daughter vesicles, see Supplementary Equation (5). For nanodroplet endocytosis, there is a small region in the morphology diagram (denoted by nested prolates in Fig. 2b), where the nanodroplet deviates from the spherical shape and attains a prolate shape for the αγ segment of the membrane. Likewise, the βγ segment of the membrane attains a prolate shape and the final morphology is two nested prolates, see the diamond symbol in Fig. 2g. For nanodroplet exocytosis, there is a larger region in the morphology diagram, where the nanodroplet deviates from the spherical shape and attains prolate and oblate shapes for the αγ segment of the membrane. In the region denoted by double prolates in Fig. 2e, where both αγ and βγ segments of the membrane attain prolate shapes and the final morphology is two connected prolates that we call double prolates, see the diamond symbol in Fig. 2h. The energy change associated with the membrane neck fission, for both nanodroplet endocytosis and nanodroplet exocytosis, is in the same order of topological transformation which is estimated to be 4π from the integral of the Gaussian curvature on a sphere 22 (see Supplementary Note 6). For nanoparticles, nonetheless, the morphology diagram is different from nanodroplets 19 , because the area and shape are fixed before and after engulfment (see Supplementary Note 7).  (a) and (b) show the cutaway views of the nanodroplet-vesicle system for different sizes of endocytic and exocytic nanodroplets, respectively. The wetting energy values ΔE are listed below for each shape. The vesicle surface is shown in gray, the nanodroplet in light blue, and the αγ segment of the membrane, where the membrane is in contact with both α and γ phases, is highlighted with red color. The nanodroplet size increases from left to right. The vesicle radius is R ve = 1 μm, the vesicle volume-to-area ratio is v = 0.7, and the relative size of nanodroplet to vesicle is denoted by r α = R α /R ve , where R α is the nanodroplet radius and R ve is the vesicle size. The bending rigidity of the vesicle is κ = 20 k B T, and the interfacial tension of the nanodroplet is Σ αβ = 5 μNm −1 .
Shape transformation by droplet endocytosis. Here we study wetting-induced shape transformations of vesicles for different nanodroplet sizes and interfacial tensions. For the volume-to-area ratio v γ < 1 of the γ phase, for which complete engulfment is possible, we observed three different regions with different wetting and morphological features. First, for small nanodroplets and very low interfacial tensions, the nanodroplet partially wetted the vesicle and did not change the prolate morphology of the vesicle (the green region in Fig. 3a). For a range of parameters shown in the pink region in Fig. 3, wetting of the vesicle induced shape transformation from prolate to oblates or stomatocytes, and vesicle partially engulfed the nanodroplet. Two shapes have the same free energy on the transition line (dash-dotted brown lines in Fig. 3a). However, for interfacial tension larger than the critical value, where interfacial energy is equal to bending energy of the sphere Σ αβ = (8πκ)/A αβ , a complete engulfment is observed (blue region in Fig. 3a). Critical interfacial tension locates the limit shapes on the morphology diagram, where the membrane neck closes (black line in Fig. 3a). In the complete engulfment region, the nanodroplet is engulfed from the equator, and the daughter vesicle takes the prolate shape. For v γ > 1, we only observed partial wetting and partial engulfment. The latter is due to the lack of membrane area for complete engulfment (the light gray region in Fig. 3a). For very low interfacial tension, the nanodroplet wetted the vesicle from the pole and did not affect the morphology due to ultra-low interfacial energy. However, for large interfacial tension, the interfacial energy dominated and changed the morphology from prolate to stomatocytes to reduce The droplet α phase volume-to-area ratio v γ for droplet endocytosis, and (c) the total elastic energy ðE αγ þ E βγ Þ=8πκ for droplet endocytosis after complete engulfment. Above v γ > 1 only partial engulfment is possible, shown by the dark gray region. (d) The phase β (aqueous phase enclosed inside vesicle) volume-to-area ratio v β for nanodroplet exocytosis. (e) the droplet α phase volume-to-area ratio v α for nanodroplet exocytosis, and (f) The total elastic energy for nanodroplet exocytosis after complete engulfment. Above v β > 1 only partial engulfment is possible, shown by dark gray region and no engulfment is possible below v β < 0 shown by light gray region. (g) The typical morphology obtained for droplet endocytosis, the symbols corresponds to the parameter space in (a). (h) The typical morphology obtained for nanodroplet exocytosis, the symbols corresponds to the parameter space in (d). The vesicle radius is R ve = 1 μm, the bending rigidity of the vesicle is κ = 20 k B T, the spontaneous curvature is zero, and the interfacial tension of the nanodroplet is Σ αβ = 17.5 μNm −1 , respectively. the interfacial energy at the cost of an increase in bending energy. The nanodroplet-vesicle morphology for different regions in the morphology diagram is shown in Fig. 3b.
Once the vesicle is in contact with the nanodroplet, in response to the capillary force arising from the αβ interface, the nanodroplet deformed the membrane and partially wetted the vesicle surface. The Lagrange multiplier Σ ve is coupled to vesicle's area A ve and acts as mechanical tension, see Eq. (1) and 23 . The bending energy counteracted the capillary force in presence of the vesicle area and volume constraint. Our results demonstrate that bending rigidity κ and interfacial tension Σ αβ both contribute to the morphology that minimizes the total energy.
Geometry of three-phase contact line and membrane neck. So far we have studied two main energetic contributions to the total energy of the nanodroplet-vesicle system, specifically bending energy and interfacial energy. However, the three-phase contact line introduces an extra energetic contribution known as line energy, see Eq.(1). A negative sign of the line tension can lead to spontaneous symmetry breaking of the three-phase contact line 4 . However, even in presence of negative line tension, breaking the contact line symmetry involves a considerable bending and interfacial energy penalty. Thus all three energetic contributions determine whether the three-phase contact line αβγ is axisymmetric or non-axisymmetric. Here we study the parameter boundary between axisymmetric and non-axisymmetric solutions for the three-phase contact line. This parameter boundary is represented as a two-dimensional surface in three-dimensional parameter space, defined by line tension, interfacial tension, and bending rigidity. We constructed this dividing surface by exploring the boundary between axisymmetric and nonaxisymmetric contact line geometries (Fig. 4a). We first scan curves of constant bending rigidity where line tension and interfacial tension determine the three phase contact line symmetry. Then, a similar isoline for constant line tensions is obtained in the two-dimensional parameter space of interfacial tension and bending rigidity. By combining discrete data points (shown by red stars in Fig. 4a) in the three-dimensional parameter space and fitting two-variable polynomial, see Supplementary Note 8, the dividing surface is constructed (Fig. 4a). The typical non-axisymmetric and axisymmetric contact line geometries are shown in Fig. 4b and c, respectively. The interpretation of the parameter boundary is rather straightforward. If we are below the parameter boundary we obtain a nonaxisymmetric contact line geometry (Fig. 4b).
However, if we are above this surface, we obtain axisymmetric contact line geometries (Fig. 4c). Several examples can intuitively clarify the tendencies and the shape of this parameter boundary. For instance, for very high interfacial tension and intermediate values of the line tension and bending rigidity, we expect to obtain axisymmetric solutions for contact line geometry. This is because breaking the contact line symmetry would cost very high interfacial energy. Likewise, we expect axisymmetric shapes for very rigid membranes, even for droplets with low interfacial tension and for intermediate values of the line tension. This is because breaking the contact line symmetry would cost a substantial bending energy. On the other hand, for very negative line tension, one would anticipate a nonaxisymmetric contact line even for rigid membranes and droplets with relatively high interfacial tensions because the breaking of the contact line symmetry reduces the line energy, similar to the nanodropletmembrane system at the molecular scale 4 . The symmetry of the three-phase contact line is an important characteristic of the nanodroplet engulfment. On molecular scales, closure of the membrane neck tends to form a tight-lipped shape for contact lines with non-axisymmetric geometries 4 . Such tight-lipped shapes are expected to suppress neck scission and reduce the cellular uptake of nanodroplets 4 .

Conclusions
In summary, we studied the interaction of nanodroplets with vesicles in the context of membrane-wetting theory. The nanodroplets studied here have interfacial tensions comparable to those of biomolecular condensates. We showed that the interaction of the nanodroplets with vesicles can mutually transform the morphology of vesicles and nanodroplets by wetting the membrane based on the competition between interfacial energy and bending energy. In addition, we found a prominent effect of the Fig. 3 Morphology diagram for nanodroplet endocytosis. The volume-to-area ratio of the vesicle is constant v = 0.7 which implies that in the absence of the nanodroplet, the vesicle attains a prolate shape as the minimum energy shape 21 . The dash-dotted brown line represents the transition line where two shapes with different morphologies have the same free energy. Complete engulfment is only possible when γ phase volume-to-area ratio of the γ phase satisfies v γ < 1 provided that the interfacial tension is large enough. The green region corresponds to the partial engulfment region where the nanodroplet does not change the prolate morphology of the vesicle. In the pink region, shape transformations from prolates to oblates and stomatocytes are observed where the nanodroplet is partially engulfed. The solid black line indicates the critical interfacial tension when the interfacial energy equals the bending rigidity of spherical out-bud, Σ αβ = 8πκ/A αβ , where A αβ is the area of αβ interface. The complete engulfment and the limit shapes are only possible in the light blue region where the interfacial tension is larger than the critical value. The vertical dashed gray line corresponds to v γ = 1 and for larger values of r α only partial engulfment is possible, (gray region). For the morphology diagram in (a) the background color and the symbols correspond to the typical vesicle-droplet morphology in (b). The vesicle radius is R ve = 1 μm, the bending rigidity of the vesicle is κ = 20 k B T and the spontaneous curvature is zero.
spontaneous curvature on the membrane-droplet morphology (see Supplementary Note 9), capable of changing the membrane curvature in the segment which is in contact with the nanodroplet. We determined the boundary in the three-dimensional parameter space, which separates axisymmetric from nonaxisymmetric geometries. We explored all material parameters involved in the shape transformation of nanodroplet-vesicle system, including bending rigidity, spontaneous curvature, interfacial tension, and line tension. The available experimental methods for the measurement of line tension are still very challenging. But an observation of contact-line symmetry seems to be feasible using super-resolution microscopy at the nanometre scale 24 , similar to the observation of endocytosis by proteins condensates with apparently non-axisymmetric contact lines 13 . Such observations, combined with the theoretical approach presented here, suggest the possibility to estimate the range of physically meaningful line tension values.

Methods
Continuum model for membrane wetting. The total energy of the nanodropletvesicle system in a partially wetted state is described by including three main energetic contributions, namely, the bending energy of the vesicle, interfacial energy of the nanodroplet, and the line energy of the three-phase contact line, as described before in the literature 3,4 . Thus, the total energy of the partially wetted state, E W , of the nanodroplet-vesicle reads the subscripts α, β and γ refer to the three distinct liquid phases. The nanodroplet is named α phase with an interface with β phase. Thus, for droplet endocytosis, the enclosed liquid phase inside the vesicle is defined by γ phase and the external liquid phase is specified by β phase. While for nanodroplet exocytosis, the enclosed liquid phase inside the vesicle is defined by β phase and the external liquid phase is specified by γ phase, see Supplementary Figure 1. The first term on the right-hand side of the Eq. (1) shows the total elastic energies of the two membrane segments with areas A αγ and A βγ in contact with the α and β phase, respectively, where the total vesicle area is constant, i.e., A ve = A αγ + A βγ . These membrane segments are in principle characterized by the two bending rigidities κ αγ and κ βγ and the two spontaneous curvatures m αγ and m βγ . The local mean curvature of the membrane is denoted by M. The second and third terms in Eq. (1) are the interfacial energy of the nanodroplet and the line energy of the nanodroplet-vesicle three-phase contact line, respectively. Here Σ αβ and λ αβγ are the interfacial tension of αβ interface and the three-phase αβγ contact line tension, respectively. Finally, the last three terms in Eq. (1), employ three different Lagrange multipliers, Σ ve , P ve and P α , to impose constant vesicle area and volume as well as constant nanodroplet volume. The area of the nanodroplet is not conserved. To reduce the number of parameters, the bending rigidity is taken to be the same in both segments κ = κ αγ = κ βγ , and we assume that the spontaneous curvature in the membrane βγ segment is negligible m βγ = 0. Furthermore, in this study, we exclude the adhesion energy of nanodroplets, by considering the special case of zero adhesion strength.
Energy minimization. A discretized form of the continuum model Eq. (1), using a triangulated mesh is used to minimize the energy of the nanodroplet-vesicle system. Here, we used the Surface Evolver (SE) package version 2.70y for triangulation and energy minimization 25 of the nanodroplet-vesicle system. The four-fold symmetry is employed for calculations of systems with zero line tension. However, for the case of negative line tension, the whole surface is considered for energy minimization. The calculation of the mean curvature and the corresponding bending energy is done based on the SE method called "star_perp_sq_mean_curvature" 25 . In this method, the mean curvature at each vertex v is given by: where ∇ A v and ∇ V v are the area gradient and the volume gradient at vertex v. Then the bending energy is the sum over the vertices: where A v is the area of adjacent facets, m v is the spontaneous curvature (intrinsic curvature in SE terminology), and κ is the bending rigidity of the membrane (modulus in SE terminology).
To keep the membrane area constant, the "facet_area" method is used along with a prescribed Lagrange multiplier. Similarly, the "edge_length" method is used with zero or negative line tension (modulus in SE terminology) to take the line energy into account. The tension value is only prescribed in liquid facets where the area is not kept constant. The combination of gradient descent, conjugate gradient, and Hessian methods are used for minimization of the total energy of the system by taking all constraint and boundary conditions into account. The vertices along the three-phase contact line are subject to a specified contact line constraint. The contact line is parametrized on an elliptic curve ðx=R a Þ 2 þ ð y=R b Þ 2 ¼ 1, where the R a and R b are the major and minor axes. These two axes are treated as two independent optimization parameters in energy minimization. For zero line tension, the two axes are set to be equal R a = R b and treated as one single optimization parameter. However, for the case of negative line tension, initially, the two axes are set to be equal to avoid any initial Fig. 4 Parameter boundary and geometry of three-phase contact line. (a) Boundary surface in the three-dimensional parameter space defined by bending rigidity, interfacial tension and line tension: above this surface, the three-phase contact lines of the nanodroplet-vesicle are axisymmetric but non-axisymmetric below the surface. Red stars show the scanned data points. The surface is a polynomial (degree of three with two variables) fit Σ fit αβ ð λ αβγ ; κÞ over the red data points, see Supplementary Note 8 for more details. The color of the surface represents the critical value of interfacial tension Σ αβ . (b, c) The nanodroplet-vesicle shapes illustrate the morphology of nanodroplet endocytosis. The morphology of the nanodroplet-vesicle system with non-axisymmetric contact line geometry is observed below the boundary surface (b), while axisymmetric contact line geometries are observed above this surface (c). The panels in the top row of (b) and (c) show the top views of the elliptic and circular contact lines, respectively. The second and third rows show two cutaway views. The vesicle radius is R ve = 1 μm, and volume-to-area ratio is v = 0.75, the spontaneous curvature is zero and relative size of nanodroplet to vesicle is r α = 0.22, respectively. condition bias. In the latter case, two axes can independently change and the contact line can then either remain in a circular geometry or can deviate from it to attain an elliptic shape as a result of the energy minimization. Such contact line parametrization can be considered as the first perturbation mode caused by negative line tension. This can be extended to higher-order perturbations by different parametrization of the contact line.

Data availability
All relevant data are available from the authors upon reasonable request.

Code availability
All relevant codes are available from the authors upon reasonable request.