Nanodroplets at Membranes Create Tight-Lipped Membrane Necks via Negative Line Tension

The response of biomembranes to aqueousphase separation and to the resulting water-in-water droplets has been recently studied on the micrometer scale using optical microscopy and elasticity theory. When such a droplet adheres to the membrane, it forms a contact area that is bounded by a contact line. For a micrometer-sized droplet, the line tension associated with this contact line can usually be ignored compared with the surface tensions. However, for a small nanoscopic droplet, this line tension is expected to affect the membrane-droplet morphology. Here, we use molecular simulations to study nanodroplets at membranes and to gain insight into these line tension effects. The latter effects are shown to depend strongly on another key parameter, the mechanical tension experienced by the membrane. For a large membrane tension, a droplet adhering to the membrane is only partially engulfed by the membrane, and the membrane− droplet system exhibits an axisymmetric morphology. A reduction of the membrane tension leads to an increase in the contact area and a decrease in the interfacial area of the droplet, initially retaining its axisymmetric shape, which implies a circular contact line and a circular membrane neck. However, when the tension falls below a certain threshold value, the system undergoes a morphological transition toward a non-axisymmetric morphology with a non-circular membrane neck. This morphology persists until the nanodroplet is completely engulfed by the membrane and the membrane neck has closed into a tight-lipped shape. The latter morphology is caused by a negative line tension, which is shown to be a robust feature of membrane−droplet systems. A closed membrane neck with a tight-lipped shape suppresses both thermally activated and protein-induced scission of the neck, implying a reduction in the cellular uptake of nanodroplets by pinocytosis and fluid-phase endocytosis. Furthermore, based on our results, we can also draw important conclusions about the time-dependent processes corresponding to the surface nucleation and growth of nanodroplets at membranes.

A queous two-phase (or biphasic) systems, 1 which are intimately related to water-in-water emulsions, 2 have been frequently used in biochemical analysis and biotechnology 3,4 to separate biomolecules, organelles, and cell membranes. The classic example are aqueous solutions of PEG and dextran, which undergo phase separation when the weight fractions of the polymers exceed a few percent. 5 More recently, water-in-water droplets have also been observed in the form of membraneless organelles 6 and biomolecular condensates 7 in living cells. So far, the response of membranes and vesicles to such water-in-water droplets has been studied by optical microscopy and elasticity theory for fluid interfaces and membranes. 8−11 In these studies, the droplets and vesicles as well as the membrane necks had linear dimensions in the micrometer range, which implies that one could ignore the line tension associated with the droplets' contact lines on the membranes.
However, when the phase separation leads to individual, well-separated droplets, it proceeds via nucleation and growth, starting from droplets with a linear dimension in the nanometer range. When such a nanoscopic droplet is nucleated at a membrane surface, the line tension is expected to play an important role. Likewise, the engulfment of nanodroplets by cellular membranes represents an essential step of pinocytosis and fluid-phase endocytosis, which has been recently used to deliver imaging agents 12,13 or drugs 14−16 to biological cells. Such pinocytic processes are also involved in the uptake of hydrocarbons by microorganisms. 17,18 Here, we study such nanodroplets in contact with membranes by coarse-grained molecular simulations using dissipative particle dynamics. 19, 20 We start with a nanodroplet that adheres to a bilayer membrane as displayed in Figure 1.
The initial axisymmetric morphology is stabilized by the mechanical tension experienced by the membrane, a tension that we control by varying the lateral size of the simulation box for a certain, fixed number of lipid molecules. The droplet consists of the liquid phase α and is bounded by two surface segments: the α−β interface between the droplet and the liquid bulk phase β as well as the droplet's contact area with the membrane corresponding to the αγ membrane segment in Figure 1. When we reduce the mechanical tension, the droplet's contact area increases and the area of the α−β interface decreases. Initially, the droplet-bilayer morphology remains axisymmetric and the contact line retains its circular shape. However, when the tension falls below a certain threshold value, the system undergoes an unexpected transition to a non-axisymmetric morphology. The latter morphology persists until the nanodroplet is completely engulfed by the membrane and the membrane neck has closed into a tight-lipped shape. A detailed analysis of the force balance along the contact line reveals that these nonaxisymmetric morphologies arise from negative values of the line tension.
For liquid mixtures without membranes, the concept of line tension was introduced by Gibbs 21 who already pointed out that line tensions can be positive or negative (see, e.g., ref 22). In contrast, interfacial tensions must always be positive to ensure thermodynamic stability. Negative values of the line tension have been observed for sessile liquid droplets on solid surfaces, 23 for lense-shaped droplets between two bulk liquids, 24 and in simulations of Lennard−Jones fluids. 25 Negative line tensions have also been found for Plateau borders in foams. 26 None of these previous studies provided evidence for a morphological transition as described here for membrane-engulfed nanodroplets.
The formation of a tight-lipped membrane neck provides an example for a liquid nanostructure that undergoes a morphological transition with a spontaneously broken symmetry. In addition, the elongated, tight-lipped neck has important consequences for the pinocytosis and cellular uptake of nanodroplets. Indeed, such a neck shape implies a large increase in the free energy barrier for thermally activated scission of the neck. Likewise, a tight-lipped neck can hardly be cleaved by the known protein-based mechanisms for membrane scission. Indeed, our current models for membrane scission by proteins such as dynamin 27 and ESCRT 28,29 are all based on the assumption that the membrane neck has a circular shape.

RESULTS AND DISCUSSION
Bilayer Membranes Exposed to Three Liquid Phases. We studied bilayer membranes in contact with nanodroplets by molecular simulations using dissipative particle dynamics (DPD) with three different sets of force parameters f ij , as described in the Methods section. The liquid phases are built up from a binary mixture of two types of beads, liquid-A and liquid-B beads. This mixture undergoes phase separation into the A-rich phase α and the B-rich phase β (see Figure S1) and provides a relatively simple model for aqueous two-phase systems with small solutes, such as polyethylene glycol (PEG) and salt, 1 or for two different liquids such as water and oil that demix above a certain concentration of one molecular component. Oil-in-water nanoemulsions are stabilized by surfactants, which reduce the interfacial tension of the oil− water interface and shield the droplets against hydrophobic interactions. In our coarse-grained model, the surfactants are taken into account by a reduced force parameter f AB , which implies a reduced interfacial tension; see Figure S1.
The bilayer membranes are assembled from coarse-grained lipids as described previously; 30 see the Methods section. To obtain symmetric bilayers, the lipid head (H) beads and the lipid chain (C) beads were taken to experience the same interactions with the A and B liquid beads, corresponding to the equalities f AH = f BH and f AC = f BC between the force parameters f ij ; see Tables 1 and 2. To study αγ membrane segments with bilayer asymmetry, the interactions of the lipid H beads with the A beads were chosen to be different from those with the B beads, corresponding to f AH ≠ f BH ; see Table  3.
Each bilayer membrane was exposed to three liquid phases, denoted by α, β, and γ, as depicted in Figure 1. The α phase formed the nanodroplet adhering to the membrane, the β Figure 1. Nanodroplet of the α phase (dark blue) adhering to a bilayer membrane (green). The nanodroplet coexists with the bulk phase β (white with blue dots). The liquid phase γ (white) above the membrane corresponds to an inert spectator phase. The membrane-droplet morphology involves three surface segments: the α−β interface between the α droplet and the β phase as well as two membrane segments, αγ and βγ, which are in contact with the α droplet and the β phase, respectively. The simulation box is a cuboid with lateral size L ∥ . The image represents a cross-section parallel to the yz plane as indicated by the right-handed orthonormal trihedron (red−green−blue) in the lower left corner.

ACS Nano
Article phase coexisted with the α droplet, and the γ phase was an inert spectator phase, consisting only of B beads. As shown in Figure 1, the resulting membrane-droplet morphology involves three surface segments: the α−β interface between the α droplet and the β phase as well as two membrane segments αγ and βγ, which are in contact with the α droplet and the β phase, respectively. The bilayer membrane displayed in Figure 1 experiences a significant tension that prevents this membrane from increasing its contact area with the nanodroplet, thereby engulfing the droplet and decreasing the area of the α−β interface. Such an engulfment process was observed as soon as we reduced the membrane tension by decreasing the lateral size L ∥ of the simulation box; see Figure 2. This reduction of L ∥ was performed in such a way that the number of lipid molecules within the bilayer, the numbers of liquid-A and liquid-B beads, the bulk pressure of the liquid phases, and the volume L ∥ 2 L z of the simulation box remained constant. Because of the latter constraint, the reduction of L ∥ leads to an increase in the perpendicular box size L z ; see Figure 2.
Engulfment of Nanodroplets by Symmetric Bilayers. We first investigated the behavior of symmetric bilayers corresponding to the force parameter set DPD-1, as described by Table 1. We started from a membrane-droplet morphology as in Figure 1 and then reduced the lateral box size L ∥ and the concomitant membrane tension over a time period of 4 μs. The resulting evolution of the membrane-droplet morphology is shown in Figure 2 and in Movie S1.
Inspection of Figure 2a shows that the initial, axisymmetric droplet was only partially engulfed by the membrane and formed an extended α−β interface with the β phase. As we reduced the lateral box size L ∥ , the area A αβ of the α−β interface decreased and the contact area A αγ between the membrane and the droplet increased; see also Figure S2. During the initial reduction of the membrane tension, the droplet-bilayer morphology remained axisymmetric and the contact line retained its circular shape. However, when we reached a certain threshold value of the mechanical tension, the system underwent an unexpected transition to a nonaxisymmetric morphology. Close to this transition, the contact line exhibited strong shape fluctuations; see Figure 2b,c and Movie S2. As we further reduced the membrane tension, the non-axisymmetric morphology persisted until the nanodroplet was completely engulfed by the membrane and the membrane neck had been closed into an unusual, tight-lipped shape; see Figure 2d. The sum of the contact area A αγ and of the interfacial area A αβ is equal to the total surface area of the α droplet. As shown in Figure S2, the surface area A αγ + A αβ remained almost constant during the whole shape evolution in Figure 2, implying that the overall shape of the α droplet stayed close to a sphere even though, locally, the droplet shape has a ridge along the tight-lipped membrane neck. The details of the simulation protocol used to generate the morphological evolution in Figure 2 as well as Movies S1 and S2 are described in the Methods section.
Fluid-Elastic Parameters of Symmetric Bilayers. To understand the origin for the tight-lipped shape of the membrane neck as shown in Figure 2d, we now consider the different parameters that determine the free energy of the membrane-droplet system. The explicit form of this free energy is given by eq 4 in the Methods section. The two membrane segments αγ and βγ experience the mechanical tensions Σ αγ and Σ βγ , each of which includes both the overall lateral stress applied via the prescribed lateral area A ∥ ≡ L ∥ 2 and the adhesive energy between the membrane and the adjacent liquid phases. 11 For symmetric bilayers, both membrane segments have negligible spontaneous curvatures, and both segments have the same bending rigidity. The α−β interface has the interfacial tension Σ αβ and the interfacial area A αβ , which implies the interfacial free energy Σ αβ A αβ . In addition, the contact line with length L αβγ contributes the line free energy λ αβγ L αβγ , which is proportional to the line tension Table 2. Parameter Set DPD-2 for Symmetric Bilayers with Variable Force Parameters f AB and f AH , Which Were Chosen As Indicated by the 16 Red Crosses in Figure 6 The symbols have the same meaning as in Table 1, and the f ij values are again in units of k B T/d. a Force parameters f ij in units of k B T/d, used to study different affinity contrasts Δ aff as defined in eq 3. All parameters have the same values as those in Table 1 except for the parameter f AH , which was taken to be 22.5, 25, and 27.5 k B T/d, corresponding to the affinity contrasts Δ aff = 0.1,0, and −0.1. The engulfment process is driven by a reduction in the lateral box size L ∥ for a fixed volume of the simulation box, from L ∥ = 130 d at t = 0 to L ∥ = 120 d at t = 4 μs, where the bead diameter d is of the order of 1 nm. This decrease of L ∥ reduces the mechanical tension Σ βγ within the membrane segment βγ. The panels in the top row show bottom views of circular membrane segments (yellow and green) around the α−β interface (blue), separated by the contact line, which is circular at t = 0 but strongly noncircular after t = 3 μs, and has closed into a tight-lipped shape after t = 4 μs. The cross-sections in the middle and bottom row are taken along the red and green dashed lines in the top row. The same time evolution is also shown in Movie S1.

ACS Nano
Article λ αβγ . The latter tension has been ignored in previous studies of membrane wetting, which focused on giant vesicles. 11,31 In contrast to the interfacial tension Σ αβ , which is always positive as required by thermodynamic stability, the sign of the line tension λ αβγ can be positive or negative. 22 A negative line tension acts to elongate the contact line and provides the main driving mechanism for the formation of non-axisymmetric buds, as illustrated in Figure 2c,d and explained further below. In principle, the closed neck of the membrane could be further stabilized by effectively attractive forces between the two membrane segments that are in close contact along the elongated contact line. To examine this possibility, we performed additional simulations of two planar membranes that were in close contact but found no evidence for such attractive interactions.
Axisymmetric Membrane-Droplet Morphologies. An axisymmetric shape of the membrane-droplet system is uniquely defined by its one-dimensional shape contour, which can be obtained from any cross-section that contains the axis of rotational symmetry. Such a parametrization is not possible for non-axisymmetric shapes that are intrinsically two-dimensional, and it is then much more difficult to compute the different free energy contributions in eq 4. Therefore, we first studied axisymmetric shapes that we obtained by increasing the lateral box size L ∥ for fixed bead number and fixed volume of the simulation box, thereby increasing the mechanical tensions Σ βγ and Σ αγ within the two membrane segments. For the symmetric bilayers studied here, axisymmetric shapes were obtained for Σ βγ ≳ 0.6 Σ αβ . For these tensions, the nanodroplet was partially engulfed by the membrane with a circular contact line, at which the α−β interface and the membrane form the intrinsic contact angle θ α *; see Figure 3. 11,31 In the latter figure, we also define two additional quantities that characterize the contact line of an axisymmetric shape: the radius R co of the contact line and the angle ψ co between the membrane contour and the projected contact line, which is perpendicular to the axis of rotational symmetry.
The first variation of the free energy with respect to the position of the contact line leads to several boundary conditions, one of which describes the balance of the tangential force components at the contact line. For symmetric bilayers with negligible spontaneous curvatures, this boundary condition has the form of: which involves, apart from the geometric quantities θ α *, ψ co , and R co as defined in Figure 3; the mechanical tensions Σ αγ and Σ βγ of the two membrane segments; the interfacial tension Σ αβ of the liquid−liquid interface; and the line tension λ αβγ of the three-phase contact line.
To determine the three surface tensions Σ αγ , Σ βγ , and Σ αβ , we use the mechanical rather than the thermodynamic definition of these tensions, as discussed in more detail in the Methods section. Thus, all three tensions are calculated from the components P xx , P yy , and P zz of the local pressure tensor, which defines the stress profile 32,33 s(r, z) = P zz − 1 2 (P xx + P yy ), as displayed in Figure 4a. All geometric quantities can be directly obtained from the simulation data; see Figure 4b,c. We are then left with only one unknown parameter in the force balance eq 1, the line tension λ αβγ , which can be computed by inserting the measured values of Figure 3. Cross-section through an axisymmetric shape of a partially engulfed droplet (compare to Figure 1) with the axis of rotational symmetry provided by the z axis and the cylindrical coordinates r, z, and φ. The cross-section contains the contour of the membrane shape consisting of the αγ (red) and βγ (green) segments, the circular contour of the α−β interface between the α droplet (blue) and the β phase (white), and the projected contact line corresponding to the horizontal black line. Important geometric parameters that enter the force balance eq 1 are the intrinsic contact angle θ α * between the membrane contour and the interface contour, the angle ψ co between the membrane contour and the projected contact line, and the contact line radius R co . in units of k B T/d 3 . This profile varies from s = −0.5 (dark blue) to s = +1 (yellow). Away from the membrane, the three liquid phases are characterized by s = 0 (light blue). The tensions Σ αβ , Σ αγ , and Σ βγ are computed from the stress profiles in the red, black, and yellow boxes, respectively. (b) Density profile ρ H (r, z) of lipid head groups in units of 1/d 3 , which varies from ρ H = 0 (dark blue) away from the membrane to ρ H = 1.5 within the lower headgroup layer of the αγ membrane segment. Note that the headgroup density in the latter segment exceeds the one in the βγ segment, which implies that the αγ segment is compressed compared to the βγ segment. (c) Profile of the density difference ρ A (r, z) − ρ B (r, z), which varies from the negative bulk density −ρ B = −3/d 3 of the liquid-B beads (dark blue) to the positive bulk density ρ A = +3/d 3 of the liquid-A beads (yellow).

ACS Nano
Article the other parameters into eq 1. Next, we will describe in some detail how we obtained the different mechanical tensions and the intrinsic contact angle from the simulations.
Mechanical Tensions for Axisymmetric Shapes. The different mechanical tensions are obtained by integrating the stress profile s(r, z) in Figure 4a across the α−β interface and across the two membrane segments. To obtain well-defined normal and tangential stresses for the curved α−β interface as well as for the curved αγ and βγ membrane segments, we calculate the stress profile locally for those spatial regions in which the α−β interface, and the two membrane segments are almost planar and essentially normal to the z axis. Those regions are indicated by the red, yellow, and black boxes in Figure 4a. The regions are discretized using a cubic lattice of mesh points with lattice constant d. We first calculate the stress profile s(z,r) by integrating over the angular coordinate (or azimuth). The different mechanical tensions are then obtained by integrating the latter stress profile; see the Methods section, which leads to the tension values displayed in Figure 5a as a function of the base area A ∥ = L ∥ 2 ; for the precise numerical values of the different tensions, see Tables S1 and S2.
The interfacial tension Σ αβ measured for the membranedroplet system is very close to the tension calculated for a planar α−β interface; see the broken horizontal line in Figure  5a as well as Figure S1b. Therefore, the measured values of Σ αβ were not affected by the curvature of the α−β interface. Inspection of Figure 5a shows that the segment tensions Σ αγ and Σ βγ increase linearly with the projected area A ∥ = L ∥ 2 . The tension values obtained in this manner are essentially independent of the exact limits of integration for the z integrals as long as these limits are well-separated from the membrane and the α−β interface.

Intrinsic Contact Angles for Axisymmetric Shapes.
For axisymmetric membrane-droplet shapes, the intrinsic contact angle has a constant value along the contact line. To measure this angle, we determined the locations of the different surfaces from the density profiles depicted in Figure  4b,c. We used the density profile ρ H of the headgroup beads, as shown in Figure 4b, and image segmentation to define the locations of the interfaces separating the lower membrane leaflet from the α droplet and β phase, respectively; see the Methods section for more details. The location of the α−β interface was obtained from the crossing criterion ρ A = ρ B for the densities ρ A and ρ B of the A and B beads, see Figure 4c. To obtain the location of the contact line, the interfaces separating the lower membrane leaflet from the two liquid phases α and β were simultaneously fitted with basis-splines, while the α−β interface was fitted to a circular segment as in Figure 3. The intersection of these two fitting curves defined the contact line, and thus, the contact line radius R co and the tilt angle ψ co introduced in Figure 3. Finally, the intrinsic contact angle θ α * was obtained from the tangents of the two fitting curves at the contact line. The results for the intrinsic contact angle θ α * are displayed in Figure 5b as a function of the mechanical tension Σ βγ of the βγ membrane segment, both for the droplet volume V α,1 and for V α,2 = 2 V α,1 as defined in the Methods section, and for the precise numerical values of the contact angle θ α *; see Tables S3 and S4. The simulation data in Figure 5b show that, for a given mechanical tension Σ βγ , the intrinsic contact angle has the same value, within the accuracy of our simulations, for both volumes V α,1 and V α,2 .
Negative Values of the Line Tension. Using the different mechanical tensions and the geometric parameters of the axisymmetric membrane-droplet morphologies as obtained from the simulation data for mechanical tensions Σ βγ ≳ 0.6 Σ αβ of the βγ membrane segment, we can now compute the line tension λ αβγ of the contact line from the force balance eq 1. All line tension values obtained in this way are negative as shown in Figure 5c, where the line tension λ αβγ is plotted as a function of the mechanical segment tension Σ βγ . The numerical values of λ αβγ are given in Tables S1 and S2. Inspection of Figure 5c shows that a larger membrane tension Σ βγ , imposed by a larger box size L ∥ , leads to a more negative value of the line tension. Furthermore, a linear extrapolation of these data to vanishing segment tension Σ βγ = 0 leads to the line tension λ αβγ = −6.87 k B T/d, indicating that the line tension remains negative even for tensionless membranes.
For the range 0.85 Figure 5c, the line tension varies within the interval −12 Using the thermal energy k B T = 4 × 10 −21 J at room temperature and the bead diameter d = 1 nm, we obtain the interval −4.9 × 10 −11 N ≲ λ αβγ ≲ −2.9 × 10 −11 N, which is comparable with the three-phase contact line tensions that have been theoretically estimated in the absence of a membrane. 34,35 For such membraneless droplets, the experimentally deduced values of λ αβγ vary over a much wider range but several experimental studies have also found negative line tensions with a similar order of magnitude. 23,24 Negative Line Tension Creates Tight-Lipped Membrane Necks. The negative value of the line tension explains the elongation of the contact line and the formation of the tight-lipped membrane neck as shown in Figure 2d and Movie S1. Indeed, the contact line with length L αβγ reduces the free energy of the system by L αβγ λ αβγ < 0, a reduction that Figure 5. (a) Mechanical tensions Σ αβ , Σ αγ , and Σ βγ as a function of the base area A ∥ = L ∥ 2 of the cuboid-shaped simulation box for droplet volume V α,2 . The interfacial tension Σ αβ for a large planar α−β interface is also displayed as the horizontal dashed line. (b) Intrinsic contact angle θ α *. (c) Line tension λ αβγ as well as (d) interfacial free energy Σ αβ A αβ > 0 and line free energy λ αβγ L αβγ < 0 as functions of the mechanical segment tension Σ βγ . In panels b− d, the black and red data points were obtained for small droplet volume V α,1 and large droplet volume V α,2 = 2 V α,1 , respectively. The black dashed line in panel c corresponds to the functional dependence λ αβγ = −0.711 d Σ βγ −6.87 k B T/d and represents a linear fit to the data for droplet volume V α,1 . The error bars indicate the standard deviations around the average values.

ACS Nano
Article becomes more significant for larger L αβγ values and thus favors the elongation of the contact line.
The symmetry breaking also increases the bending energy of the membrane but this increase is overcompensated by the negative line free energy; see Figure S3 and Movie S2. The increase in bending energy arises primarily from the highly curved membrane segment along the contact line, because the overall shape of the completely engulfed droplet remains close to a sphere as follows from Figure S2. For the tight-lipped membrane neck in Figure 2d, these highly curved membrane segments resemble hemicylinders with curvature radius π = ⊥ R /2 me , which is of the order of the membrane thickness me . As shown in the Methods section, the associated bending energy increase ΔE be is proportional to the bending rigidity κ of the membrane and leads to the effective line tension: corresponding to the superposition of the "bare" line tension λ αβγ and the bending energy contribution ΔE be /L αβγ . The tight-lipped shape of the closed membrane neck will be energetically favorable compared to the closed axisymmetric shape as long as λ eff < 0. For the symmetric bilayers studied here, the bending rigidity has the value κ ≃ 12.6 k B T, which was calculated from the area compressibility modulus as in refs 36 and 37. This κvalue falls within the range of experimental values measured for a typical phospholipid such as POPC at room temperature. 38 Using the bending rigidity κ ≃ 12.6 k B T together with the bilayer thickness ≃ d 5 me and the length L αβγ ≃ 70 d of the contact line in Figure 2d, we obtain the estimate ΔE be / L αβγ ≃ 3.1 k B T/d for the positive line tension contribution from the highly curved membrane segments. This estimate is somewhat larger than the numerically calculated value ΔE be / L αβγ ≃ 1.74 k B T/d obtained for the tight-lipped membrane neck in Figure 2d, see section S5 in the Supporting Information. However, the "bare" line tension λ αβγ as plotted in Figure 5c is always smaller than −6.87 k B T/d. It then follows from eq 2 that the effective line tension λ eff is negative for all of these λ αβγ values and leads to a tight-lipped membrane neck as observed in our simulations. Furthermore, because the bending energy increase ΔE be is proportional to κ, we also conclude from eq 2 that the effective line tension remains negative for significantly larger κ-values up to about 28 k B T.
Negative Line Tensions for an Enlarged Parameter Set. The negative values of the line tension as depicted in Figure 5c apply to symmetric bilayers with vanishing spontaneous curvature as obtained for the parameter set DPD-1 in Table 1. To find out whether the negative sign of the line tension is a robust property of membrane-droplet systems, we next studied symmetric bilayers for the enlarged parameter set DPD-2 in Table 2. In the latter set, we varied the two force parameters f AB and f AH in a systematic manner, retaining the symmetry condition f BH = f AH . The two parameters f AB and f AH are likely to have the largest effect on the force balance eq 1 because they determine the properties of the three interfaces that meet at the contact line.  Figure 6a, the line tension was obtained by two-dimensional interpolation.
In Figure 6b, we replot the λ αβγ values as a function of the interfacial tension Σ αβ , which is determined by the force parameter f AB ; see Figure S1b. For each of the four values of f AH used in Figure 6a, we obtain an essentially linear dependence of the line tension on the interfacial tension. Indeed, each set of data is well fitted by a linear expression of the form c 1 Σ αβ +c 0 with c 0 < 0 and c 1 < 0; see inset of Figure  6b. Thus, these linear expressions remain negative for all positive values of the interfacial tension Σ αβ , i.e., for the whole physically meaningful range of this tension. Therefore, the linear fits of our data predict that the line tension is negative for all possible values of the interfacial Σ αβ .
Interfacial Tensions in Nanoemulsions and Aqueous Two-Phase Systems. In Figure 6b, the interfacial tension Σ αβ is given in units of k B T/d 2 , which is about 4 mN/m at room temperature, comparable with the interfacial tensions of oil-in-water droplets stabilized by surfactants in nanoemulsions. 39,40 Droplets in aqueous two-phase systems, however, exhibit interfacial tensions Σ αβ that can vary over a fairly wide range. The highest interfacial tensions measured in aqueous PEG−dextran solutions, for example, were of the order of 1 mN/m, whereas the lowest tensions were only about 0.2 × 10 −3 mN/m, reflecting the vicinity of a critical demixing point. 5,11 The linear extrapolation of our data as displayed in Figure 6b should provide reliable estimates for Σ αβ ≃ 1 mN/m but may become unreliable for Σ αβ ≃ 10 −3 mN/m. Therefore, we predict that the line tension λ αβγ is also negative for aqueous two-phase systems, provided that one considers two-phase coexistence sufficiently far from the critical demixing point. For the displayed range of parameters, the line tension λ αβγ was always found to be negative and to have a value between −5 k B T/d (yellow) and −14 k B T/d (dark blue); see the vertical bar with color code. All simulations were performed for the smaller droplet volume V α,1 , lateral box size L ∥ = 140 d, and force parameter f BH equal to f AH as in Table 2. (b) Line tension λ αβγ as a function of interfacial tension Σ αβ , which is directly determined by the force parameter f AB (see Figure S1b) for different values of the force parameter f AH . The four straight lines, corresponding to the linear expressions in the inset, provide good fits to the data for all four values of f AH . Furthermore, linear extrapolation of the data to small values of Σ αβ suggests that the line tension remains negative even in the limit of small interfacial tensions. The basic tension scale k B T/d 2 is about 4 mN/m.

Article
To obtain reliable predictions for interfacial tensions that are much smaller than k B T/d 2 ≃ 4 mN/m, we would have to study much larger simulation boxes. Indeed, if we decreased the interfacial tension Σ αβ by decreasing the force parameter f AB , the width of the α−β interface would increase, and the interface would become more and more fuzzy. More precisely, hyperscaling 41 implies that the interfacial width grows as Σ αβ k T/ B as we approach a critical demixing point. Now, to determine the membrane-droplet geometry in a reliable manner, we need to consider droplet sizes that are large compared to this interfacial width. To study interfacial tensions Σ αβ that are of the order of 0.1 k B T/d 2 , for example, we would have to simulate droplets with linear dimensions that are increased by a factor 10 1/2 ≃ 3.3 or, equivalently, with volumes V α that are increased by a factor 10 3/2 ≃ 32. In principle, such droplet sizes are accessible to our simulation approach but they would be computationally very expensive.
Engulfment of Nanodroplets by Asymmetric Bilayers. So far, we discussed the engulfment of nanodroplets by symmetric bilayers for which the lipid heads and chains experience the same interactions with the liquid-A and liquid-B beads. As a consequence, the α and the β phase have the same adhesion free energy per unit area. If the membrane were planar and the droplet sufficiently large so that we could ignore the line tension λ αβγ , the symmetry condition f BH = f AH would lead to the intrinsic contact angle θ α * = 90 • . As shown in Figure 5b, the nanodroplets studied here exhibit contact angles in the range 80 • ≲ θ α * ≲ 100 • , with deviations from 90 • that arise from the line tension.
We will now describe the engulfment of nanodroplets by asymmetric bilayers that are obtained for force parameters f BH ≠ f AH , corresponding to the parameter set DPD-3 in Table 3, and show that these systems lead to negative line tensions as well. The different interactions of the liquid-A and liquid-B beads to the lipid head beads H will be characterized by the affinity contrast: For symmetric bilayers with f BH = f AH as discussed above, the affinity contrast Δ aff vanishes. Positive values of Δ aff are obtained for f BH > f AH , which describe lipid head groups that prefer liquid-A beads over liquid-B beads. In contrast, negative values of Δ aff correspond to f AH >f BH and thus to a stronger affinity between the lipid head beads and the liquid-B beads. For planar membranes and sufficiently large droplets, affinity contrasts Δ aff > 0 and Δ aff < 0 would lead to contact angles θ α * < 90 • and θ α * > 90 • , respectively.
It is important to note that a nonzero affinity contrast Δ aff ≠ 0 has two important consequences. First, this contrast affects the overall adhesion energy of the membrane-droplet system. Indeed, for positive and negative affinity contrast Δ aff , the system tries to maximize the contact area A αγ and the noncontact area A βγ , respectively, for a fixed total number of lipid molecules and fixed droplet volume V α . Second, a nonzero affinity contrast also implies that the αγ membrane segment acquires a spontaneous curvature.
During pinocytosis and fluid-phase endocytosis of nanodroplets, the droplets originate from the exterior solution and adhere to the outer bilayer leaflets. To compare our results with pinocytic and endocytic processes, it will be convenient to use the convention that the spontaneous curvature m αγ is negative if the αγ bilayer segment prefers to bulge toward the spectator-phase γ, which then represents the interior solution. Likewise, the mean curvature of a membrane patch is also taken to be negative if this patch bulges toward the γ phase. This convention implies that the αγ membrane segment in Figure 1 has a negative mean curvature. As described in section S6 of the Supporting Information and illustrated by Figures S4 and S5, the spontaneous curvature m αγ of the αγ membrane segment can be determined by studying planar bilayers, with one leaflet exposed to the α phase and the other leaflet exposed to the γ phase. The latter method leads to positive and negative values of m αγ for positive and negative affinity contrasts, respectively, because the membrane segment prefers to enlarge the area of the bilayer leaflet in contact with the A-rich phase α for Δ aff > 0 and the area of the other leaflet in contact with the A-poor phase γ for Δ aff < 0.
The affinity contrasts Δ aff = +0.1 and Δ aff = −0.1, for example, generate the spontaneous curvatures m αγ = +0.047/d and −0.044/d, which are quite large. Indeed, for bead diameter d ≃ 1 nm, we obtain |m αγ | ≃ 1/(20 nm) in both cases. For such strongly asymmetric bilayers, the force balance condition along the contact line, which is given by eq 1 for symmetric bilayers, becomes more complicated and involves additional terms that depend on the local membrane curvatures close to the contact line. 11 In fact, minimization of the total free energy as described by eq 4 leads to a discontinuity of the mean curvature across the contact line, which is, however, difficult to determine by simulations. Thus, instead of deducing the line tension λ αβγ from the force balance equation for asymmetric bilayers, we directly studied the droplet engulfment by such bilayers. Using the same simulation protocol as in Figure 2, we observed elongated contact lines and tight-lipped membrane necks for all spontaneous curvatures within the range −1/(20 nm) ≤ m αγ ≤ 1/(20 nm) as illustrated in Figure 7. Therefore, we conclude that the line tension is negative for all of these m αγ values.
Increased Stability of Tight-Lipped Necks against Membrane Scission. To cleave a membrane neck by membrane scission, one has to create two hydrophobic edges across the membrane segment that forms the neck. The corresponding free energy barrier is proportional to the combined length L ed of the two edges. For a closed circular  Figure  2d. (c) Positive affinity contrast Δ aff = +0.1, corresponding to lipid heads that prefer the A liquid beads and to an asymmetric bilayer with positive spontaneous curvature m αγ = +1/(22.7 nm). In all cases, the membrane neck closes into a tight-lipped shape.

ACS Nano
Article neck, the neck has a diameter that is comparable to the bilayer thickness ≃ d 5 me , which implies the neck perimeter π me ≃ 15.7d and the combined edge length L ed = π 2 me ≃ 31d. However, the tight-lipped neck displayed in Figure 2d has an increased perimeter of about 70 d and cleavage of the latter neck leads to the combined edge length L ed ≃ 140 d. Therefore, the free energy barrier for the tight-lipped neck is about 140/31 = 4.5 times larger than the barrier for a circular neck.
In general, a neck can be cleaved by thermal fluctuations or by some active protein machinery. For thermally activated scission, the scission rate depends exponentially on the free energy barrier. As a consequence, the increased barrier for cleavage of the tight-lipped neck leads to a strong reduction of thermally activated scission. However, during protein-induced scission, proteins such as dynamin 27 form a constrictive ring around the neck, whereas other proteins such as ESCRT 28,29,42 form an adhesive cone within the neck. In all models that have been used to describe these proteinmediated scission processes, the neck is taken to have a circular shape. In fact, according to these standard models for protein-induced scission, tight-lipped necks as observed in our simulations (see Figures 2d and 7) can hardly be cleaved by any of these proteins. Therefore, a tight-lipped neck shape suppresses both thermally activated and protein-induced scission, which implies that such a neck may arrest the pinocytic process in the completely engulfed state.
Membrane Necks Arising from Wetting by Micrometer-Sized Droplets. Finally, we discuss the implications of our results for the formation of closed membrane necks between micrometer-sized water-in-water droplets in contact with GUVs as observed in ref 43 and theoretically studied in ref 11. When the aqueous solution within the vesicle undergoes phase separation, the coarsening of small droplets eventually leads to two coexisting droplets. The two droplets will be completely engulfed by the vesicle membrane provided that the combined volume of these droplets is sufficiently small, the membrane area is sufficiently large, or both. 11 Such a wetting morphology, corresponding to exocytic engulfment, implies a closed membrane neck that replaces the α−β interface between the two droplets. Likewise, when the exterior aqueous solution undergoes phase separation, droplets of the minority phase can adhere to the outer leaflet of a vesicle membrane. Such an adhering droplet can also become completely engulfed by the membrane, provided that the droplet volume is sufficiently small, the membrane area is sufficiently large, or both. The latter wetting morphology provides an example for endocytic engulfment and again involves a closed membrane neck that now replaces the α−β interface between the droplet and the aqueous bulk phase.
If these closed membrane necks are assumed to be circular, they are governed by a stability condition that is completely analogous to the condition for two-domain vesicles, with the line tension of the domain boundary replaced by the contact line tension λ αβγ . 11 However, in view of the simulation results presented here, the assumption of a circular neck does not apply in general. Indeed, it follows from Figure 6b that the line tension λ αβγ is negative for interfacial tensions Σ αβ of the order of millinewtons per meter. The latter tension values apply, e.g., to oil-in-water emulsions stabilized by surfactants 39,40 as well as to aqueous two-phase systems sufficiently far from the critical demixing point. 5 Indeed, for the oil-inwater emulsion droplets studied in ref 39., the surfactant monolayers at the oil−water interface reduced the interfacial tension to about 2 mN/m. Therefore, when such an oil-inwater or a water-in-water droplet adheres to a lipid membrane, the corresponding contact line should have a negative line tension and should attain an elongated, tight-lipped shape when the membrane tension falls below a certain threshold value.
As we approach the critical demixing point of an aqueous two-phase system, the interfacial tension Σ αβ of water-in-water droplets becomes ultralow and eventually vanishes. 5 The linear extrapolation of our simulation data to small Σ αβ values is consistent with a negative line tension for the whole range of physically meaningful tensions Σ αβ > 0 (see the straight lines in Figure 6b), but we cannot exclude the possibility that the line tension depends on the interfacial tension in a nonlinear manner and becomes positive for sufficiently small Σ αβ . For a positive line tension, the membrane neck will be axisymmetric and governed by the stability condition derived in ref 11.

CONCLUSIONS
In summary, we have shown that the line tension of a nanodroplet adhering to a bilayer membrane is negative for a wide range of parameters. As the mechanical membrane tension is reduced, the negative line tension leads to a morphological transition from axisymmetric to non-axisymmetric membrane-droplet morphologies (Figure 2 and Movie S1). The latter type of morphology persists until the nanodroplet is completely engulfed by the membrane and the membrane neck attains an elongated, tight-lipped shape, as displayed in Figures 2d and 7 for symmetric and asymmetric bilayers, respectively. The formation of such a membrane neck provides an example for a liquid nanostructure that undergoes a morphological transition with a spontaneously broken symmetry.
We also argued that this unusual neck shape impedes pinocytosis and fluid-phase endocytosis of nanodroplets. Examples are provided by oil-in-water nanoemulsions 39,40 as well as aqueous two-phase systems sufficiently far from the critical demixing point, 5 which are both characterized by interfacial tensions Σ αβ of the order of 1 mN/m, in the same range as investigated here (Figure 6b). When such oil-in-water or water-in-water droplets interact with GUVs, they may be completely engulfed by the vesicle membrane but the resulting membrane neck will be stabilized against scission by its tightlipped shape. When these droplets interact with cells, however, the cellular uptake often involves membrane-bound proteins such as clathrin. Indeed, both clathrin-dependent 12,14,16 and caveolae-dependent 12 pathways have been discussed for pinocytosis and cellular uptake of nanoemulsion droplets.
One interesting question that requires further studies is how these protein-dependent pathways of pinocytosis are affected by a negative line tension of the nanodroplet. If the negative line tension leads to a tight-lipped membrane neck as observed here in the absence of proteins, scission and cellular uptake of the completely engulfed droplet will be strongly suppressed. However, it is also conceivable that the protein coat can control the engulfment process in such a way that the closed neck becomes circular even for a negative line tension. Another open issue is the sign of the line tension for ultralow interfacial tensions Σ αβ ≃ 1 μN, as found for aqueous twophase systems close to their critical demixing point 5 as well as

ACS Nano
Article for membrane-less organelles and biomolecular condensates. 6,44 Linear extrapolation of our simulation data predicts a negative line tension for the whole physically meaningful range Σ αβ > 0 (Figure 6b), but we cannot exclude the possibility that the line tension becomes positive for sufficiently small values of Σ αβ .
For rigid nanoparticles, clathrin-dependent endocytosis and cellular uptake depend non-monotonically on the particle size, with an optimal diameter of about 50 nm. 45 This size dependence can be understood in terms of the spontaneous curvature m pro generated by the protein coat, which implies a preferred radius of 1/|m pro | ≃ 40 nm for the clathrin-encaged endocytic vesicle. 46 Because of this preferred size of the endocytic vesicle, the uptake of nanodroplets should also be characterized by a non-monotonic dependence on the droplet size. The latter prediction is accessible to experimental studies in close analogy to those performed in ref 45. In the present paper, we studied nanodroplets with a constant, time-independent volume, but, based on our simulation results, we can also draw important conclusions for the closely related time-dependent processes corresponding to surface nucleation and growth of nanodroplets at membranes. A single droplet that has been nucleated at a membrane will initially attain a partially engulfed morphology as in Figure 1 and will keep this morphology during its diffusion-limited growth provided the membrane tension is sufficiently large. For small membrane tension, however, we will observe a competition between the engulfment and the growth of the droplet. If the growth is slow, the droplet can only grow up to a certain size before it is completely engulfed by the membrane. This competition becomes particularly interesting if we consider multisite nucleation of several droplets because the tensions of the membrane segments along the contact line of one droplet will now depend on the engulfment states of all droplets. Thus, the first few droplets that are nucleated at a membrane may become completely engulfed but the nucleation of additional droplets will effectively increase the membrane tension, which tends to open the closed membrane necks up again and to allow the droplets to continue their growth. We are currently extending our computational approach to address these nucleation and growth phenomena in a systematic manner as will be described in a subsequent paper.

METHODS
Dissipative Particle Dynamics. In our molecular simulations, DPD method, 19,20 which is based on discrete particles or beads that represent groups of liquid molecules as well as groups of atoms within the head groups and hydrophobic tails of the lipid molecules. All beads have the same diameter, d, which provides the basic length scale of the system. The basic energy scale is taken to be the thermal energy k B T. A pair of beads of types i and j repel each other by shortranged conservative forces of magnitude f ij . The same beads also interact by a random and a dissipative force connected via the dissipation-fluctuation theorem that act as a thermostat. DPD has been previously used to elucidate the properties of bilayer membranes at molecular and nanoscopic scales. 20,30,47−52 Our simulations were based on the molecular model of ref 30, in which all lipid molecules consist of three hydrophilic H beads and two hydrophobic chains, each consisting of six hydrophobic C beads. In addition, two types of liquid beads, liquid-A and liquid-B beads, were used to model the coexisting liquid phases in contact with the membrane: the α phase is enriched in liquid-A beads, while the β phase consists primarily of liquid-B beads.
Choice of DPD Force Parameters. In general, the DPD force parameters are chosen in such a way that they reproduce certain nanoscopic and mesoscopic properties of the system. Thus, the forceparameter values f AA = f BB = 25 k B T/d ensure that the liquid density is equal to 3/d 3 as appropriate for water. The DPD force parameter f AB determines the interfacial tension Σ αβ of the α−β interface. We measured this tension for the slab geometry shown in Figure S1a. For 40 k B T/d ≤ f AB ≤ 60 k B T/d, the interfacial tension lies within the interval 0.8 k B T/d 2 ≤ Σ αβ ≤ 2.6 k B T/d 2 ; see Figure S1b. The force parameters f HH , f HC , and f CC were chosen as in ref 30 to obtain a linear relationship between the mechanical membrane tension and the molecular area per lipid, as observed experimentally, and to ensure that flip-flops of lipids between the two leaflets are suppressed over microsecond time scales.
The complete parameter set DPD-1 that we used to obtain the simulation results in Figures 2, 4, and 5 is given in Table 1. The data in Figure 6 correspond to the extended parameter set DPD-2 given in Table 2. Both parameter sets DPD-1 and DPD-2 lead to symmetric bilayers with negligible spontaneous curvature and intrinsic contact angles close to 90 • . We also studied asymmetric bilayers, using the parameter set DPD-3 with f BH ≠ f AH ; see Table 3.
Geometry of Simulation Box and Boundary Conditions. The simulation box was taken to be a cuboid with its edges parallel to the Cartesian coordinates x, y, and z; the base of the cuboid with area A ∥ = L ∥ 2 was taken to be parallel to the xy plane. Initially, a planar bilayer of lipid molecules was assembled parallel to the xy plane, together with a hemispherical nanodroplet of liquid-A beads that formed a circular contact area with the bilayer. The remaining volume of the box was filled with liquid-B beads. The bilayer consisted of 13225 lipids per leaflet. A pair of droplet volumes consisting of 21 707 and 43 414 A beads were simulated, both of which were much smaller than the combined number of about 3 445 000 A and B beads. The corresponding droplet volumes V α,1 = 6837.6 d 3 and V α,2 = 2 V α,1 = 13675.3 d 3 can be estimated from the volume of a spherical α droplet immersed in the β phase. All simulations were performed with the LAMMPS package; 53 density profiles were analyzed using the Matlab image processing toolbox; 54 and VMD 55 was used for visualization.
During the simulations, the liquid-A beads were forced to stay on one side of the membrane using a planar and semipermeable wall that repelled the A beads by a purely repulsive Lennard−Jones potential with cutoff r c = 1.122 d. The semipermeable wall could adjust its z coordinate to ensure mechanical equilibrium between the β and γ phases with P β = P γ = 23.7 k B T/d 3 . In this way, we obtained two distinct bulk phases β and γ on the two sides of the membranes: the β phase contained some liquid-A beads, whereas the γ phase contained only liquid-B beads. Periodic boundary conditions were used in the all three Cartesian directions. In addition, the lateral diffusion of the nanodroplet was compensated by lateral displacements of the system so that the center-of-mass of the droplet remained at x = y = 0. All beads were taken to have the same mass m b . A time step of 0.01 τ was used, where the basic time scale corresponds to about 1 ns and the bead diameter d to about 1 nm as estimated in refs 47 and 30 from the lateral diffusion constant of the lipid molecules and from the thickness of the bilayer membrane, respectively.
On the time scales of our simulations, the lipid bilayers are essentially impermeable to the liquid-A and liquid-B beads. This impermeability can be inferred from the density profiles of the liquid-A beads as shown in Figure S5. This figure displays the bead densities of the α and γ phases, which are now separated by two planar bilayers. The α phase consists primarily of liquid-A beads, whereas the γ phase contains no such beads at the beginning of the simulations. If the bilayers were permeable to the A beads, these beads should eventually show up in the γ phase. However, we do not find any such beads, even after 25 μs.
Mechanical Definition of Surface Tensions. The tension of an interface between two liquid phases can be defined in two ways. First, it can be defined using the thermodynamic limit and an expansion of ACS Nano Article the system's free energy in powers of the system size L. The interfacial tension is then obtained from the term proportional to L 2 (for three-dimensional systems). Second, it can be defined mechanically by the work expended to increase the area of the interface. One then obtains the interfacial tension from the integral over the local stress profile s(r, z), which is derived from the pressure tensor. Even though these two definitions look quite different, they are in fact equivalent (see, e.g., ref 22). The situation for membranes is somewhat different. On the one hand, a sufficiently large membrane segment that is under mechanical tension can always lower its free energy by rupture or poration. Therefore, we cannot perform the thermodynamic limit of a membrane under tension. On the other hand, the membrane tension can still be defined by the local stress profile 32 and can then be calculated for relatively small membrane segments as considered here.
Control Parameters. For a given lateral size L y = L x ≡ L ∥ of the simulation box, we first used an NPT ensemble with the Berendsen barostat 56 and adjusted the perpendicular extension L z of the box in such a way that the pressure P zz attained the standard DPD value P zz = 23.7 k B T/d 3 corresponding to the bulk liquid density 3/d 3 . We then kept the perpendicular extension L z at the adjusted value and continued our simulations in the corresponding NVT ensemble. For the simulations of tense membranes, the lateral box size L ∥ was increased in a stepwise manner using L ∥ = 130, 135, 140, 145, and 150. For each L ∥ value, we collected data from two separate simulations with 3500 τ. The three-dimensional pressure and density profiles were calculated on a cubic grid with mesh size equal to the bead diameter d using Irving−Kirkwood contours and the algorithm described in refs 33 and 57 for the local pressure.
The sequence of snapshots as shown in Figure 2 and in Movie S1 was based on a more-elaborate simulation protocol with a timedependent lateral box size L ∥ = L ∥ (t). The system was first equilibrated for L ∥ = 130 d as described in the previous paragraph. We then switched from the NVT to the NPT ensemble to keep the bulk pressure P zz at its standard DPD value P zz = 23.7 k B T/d 3 as we decreased the lateral size L ∥ and, thus, the overall lateral stress acting on the membrane using a small and constant rate dL ∥ /dt. The constant bulk pressure and the constant bead number imply a constant volume of the simulation box. Therefore, the perpendicular extension L z of the simulation box increased with decreasing L ∥ as can be seen in Figure 2.
For the snapshots displayed in Figure 2, the box compression rate was taken to be dL ∥ /dt = −0.0025 d/ns, which reduced the lateral size L ∥ from 130 d to 120 d within 4 μs. Additional simulations with different values of the rate dL ∥ /dt confirmed that the sequence of observed morphologies does not depend on the precise value of dL ∥ / dt as long as |dL ∥ /dt| ≲ 0.01 d/ns.
Free Energy of Membrane−Droplet System. The two membrane segments αγ and βγ with areas A αγ and A βγ experience the mechanical tensions Σ αγ and Σ βγ , each of which includes both the overall lateral stress applied via the prescribed lateral area A ∥ ≡ L ∥ 2 and the adhesion free energy between the membrane and the adjacent liquid phase. 11 The αβ interface with area A αβ contributes the interfacial free energy Σ αβ A αβ and the volume term P αβ V α with the pressure difference P αβ ≡ P α −P β > 0 between the pressures P α and P β within the α droplet and the β phase. In addition, the contact line with length L αβγ contributes the line free energy λ αβγ L αβγ , which is proportional to the line tension λ αβγ .
The local membrane shape is described by the mean curvature M, which tries to adapt to the spontaneous curvatures m αγ and m βγ of the two membrane segments. The total free energy E of the membrane-droplet system is then given by: with the bending rigidities κ αγ and κ βγ of the two membrane segments. The basic energy scale is taken to be the thermal energy k B T, which is equal to 4 × 10 −21 J at room temperature. For negligible spontaneous curvatures m αγ and m βγ , the first variation of E with respect to the position of the contact line leads to the force balance calculations (eq 1).
Bending Energy of Tight-Lipped Membrane Necks. For a uniform membrane with negligible spontaneous curvature and uniform bending rigidity κ, the contribution from the bending energy in eq 4 reduces to 2κ∫ dA M 2 , which represents the surface integral over the square of the local mean curvature M. We now apply this expression to closed membrane necks. If the neck closed in an axisymmetric manner, the diameter of the circular contact line would be comparable to the membrane thickness me and its perimeter to π me . When the neck closes into a tight-lipped shape as in Figure 2d, the contact line resembles two straight and parallel line segments that are connected by two highly curved line segments, each of which has the curvature radius π /2 me . Therefore, the two straight segments of the contact line have the combined length π Δ = − αβγ αβγ L L me . Along these straight segments, the membrane has a roughly hemicylindrical shape, see Figure 2d. The radius R ⊥ of these hemicylinders is comparable to the membrane thickness me , which implies that their mean curvature ≃ M 1/(2 ) me and their area π ≃ Δ αβγ A L me . The bending energy of these hemicylinders is then given by: per unit length of the contact line. Mechanical Tensions from Stress Profiles. For the axisymmetric membrane-droplet morphologies, we calculated the mechanical tension from the two-dimensional stress profiles s(z,r), using trapezoidal numerical integration. The integration domains are chosen in such a way that the membrane segments and the α−β interface are essentially flat; see the black, red, and yellow boxes in Figure 4a. The mechanical segment tensions Σ αγ and Σ βγ and the interfacial tension Σ αβ of the α−β interface are then obtained by integrating the stress profile using a cubic lattice of discrete sites with a lattice constant d. The profile s(z, r) was obtained by integrating the local stress s(x, y, z) over the angular coordinate (or azimuth). This angular integration was performed by summing up all values of the local stress within an annulus with inner radius r − d/2 and outer radius r + d/2.
Leaflet−Liquid Interfaces from Density Profiles. The threedimensional density profile ρ(x, y, z) is first calculated as a function of the Cartesian coordinates x, y, and z using a cubic lattice. For axisymmetric morphologies, one can then calculate the density as: which depends on the cylinder coordinates = + r x y 2 2 and z but is independent of the azimuth φ. The sum runs over the number N r of lattice points within an annulus of radius r ± d/2. Twodimensional density profiles are then generated by two-dimensional interpolation within the (r, z) plane, as illustrated in Figure 4b for L || = 130 d.
The location of the leaflet−liquid interfaces is based on iso-density stripes, which are defined by a certain range of H bead density ρ H as given by 0.15/d 3 ≤ ρ H ≤ 0.3/d 3 . To extract the spatial location of these stripes, image files are generated from the two-dimensional density profiles ρ(r, z) as obtained from eq 6, and the Matlab image processing toolbox 54 is used for image segmentation of these densityprofile images.

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acsnano.8b06634.

ACS Nano
Article Additional details, figures, and tables on phase separation, contact area and interfacial area of nanodroplets, numerical results for symmetric bilayers, contributions to total free energy, numerical computation of bending energies, as well as affinity contrast and spontaneous curvature, movie captions (PDF) A video showing time-dependent engulfment of nanodroplets by lipid bilayers (AVI) A video showing shape and energy fluctuations of the nanodroplet partially engulfed by the bilayer membrane (AVI)