Introduction

Clathrin-mediated endocytosis (CME) is a fundamental cellular process to transport lipids, membrane proteins and extracellular cargo molecules into the cell (McMahon and Boucrot, 2011; Sorkin and Puthenveedu, 2013; Lu et al., 2016; Kaksonen and Roux, 2018; Lacy et al., 2018; Mettlen et al., 2018). In mammalian cells, a small patch of flat plasma membrane is shaped into a spherical vesicle when CME occurs (Avinoam et al., 2015). Clathrin molecules are essential for the membrane remodelling process. They are made of three subunits that form a triskelion, which further assemble into a cage-like structure in vitro (Musacchio et al., 1999; Shraiman, 1997). The minimum cages contain 16 polygons (Fotin et al., 2004) and the most commonly observed ones are semi-regular icosahedral cages (Cheng et al., 2007; Dannhauser and Ungewickell, 2012; Heuser, 1980). Two main hypotheses are under debate regarding how the clathrin coat scaffolds the flat membrane into a spherical vesicle in vivo (Avinoam et al., 2015; Chen and Schmid, 2020; Frey and Schwarz, 2020; Kaksonen and Roux, 2018; Scott et al., 2018). The constant area model asserts that the clathrin molecules initially polymerize into a flat lattice with a regularly arranged hexagonal structure, and later reorganization of the bonds between adjacent clathrins results in formation of pentagons in the hexagonal lattice, which in turn leads to curvature generation of the clahtrin coat (Fotin et al., 2004; Heuser and Anderson, 1989) (Figure 1a). Adhesion of clathrin molecules with the substrate has been suggested to contribute a flattening force that prevents curvature generation. Release of the flattening force therefore could induce curvature of the clathrin coat with preloaded pentagons (Sochacki et al., 2021). The alternative constant curvature model asserts that the intrinsic curvature of the elements of the lattice is kept constant during the assembly and expansion of the lattice, therefore curvature generation occurs from the very beginning of clathrin assembly (Figure 1b).

Figure 1.

In order to distinguish between the two models, the dynamics of clathrin assembly and the geometry of membrane shapes are needed. Fluorescence microscopy, including light sheet and MINFLUX, has revealed the assembly dynamics of the clathrin coat (Aguet et al., 2016; Ferguson et al., 2016), as well as other proteins that participate in endocytosis (Sirotkin et al., 2010; Taylor et al., 2011; Kukulski et al., 2016; Kaksonen et al., 2006; Balzarotti et al., 2017; Gwosch et al., 2020; Schmidt et al., 2021), while electron tomography has been able to resolve membrane shapes during endocytosis (Avinoam et al., 2015; Kukulski et al., 2012). However, neither of the methods can capture both spatial and temporal information at the same time. Under conventional fluorescence microscopy, the clathrin-coated pits appear as diffraction-limited spots due to their small size, which is typically ∼ 30 − 150nm in mammals (McMahon and Boucrot, 2011) and yeast (Kukulski et al., 2012), and shape information of the membrane is completely lost (Cocucci et al., 2012; Loerke et al., 2009; Picco et al., 2015; Sirotkin et al., 2010; Kural and Kirchhausen, 2012). On the other hand, super-resolution fluorescence microscopy has been able to reveal the protein organization at the endocytic pit (Mund et al., 2018; Cocucci et al., 2012; Arasada et al., 2018) and to reconstruct the shape of the clathrin coat (Sochacki et al., 2017; Scott et al., 2018; Mund et al., 2023) from averaging over ensembles of endocytic sites. Correlative light and electron microscopy (CLEM) method has exploited the fluorescence of fiducial markers to locate endocytic sites while resolving membrane shapes using electron tomography. However, both super-resolution and CLEM requires sample fixation, therefore, one can identify multiple endocytic sites at the same time and perform the average, yet unable to trace a single endocytic site over time. The temporal information is nevertheless lost.

As a result of the incomplete information obtained by existing experimental methods, both hypotheses have experimental support. Experiments that apply electron microscopy to resolve the membrane shapes of endocytic pits favor the constant area model (Avinoam et al., 2015; Bucher et al., 2018; Sochacki and Taraska, 2019; Sochacki et al., 2021). However, super-resolution imaging combined with analysis of the fluorescence intensity of the clathrin coat is inclined towards the constant curvature model (Willy et al., 2021). In addition, it was argued that the energetic cost of bond reorganization in a regular hexagonal lattice in the constant area model may be too large to be fulfilled (Frey et al., 2020; Frey and Schwarz, 2020; Kirchhausen et al., 2014).

Extensive theoretical efforts have been dedicated to model membrane morphology during endocytosis (Fu and Johnson, 2023), of which molecular dynamics simulations (Varga et al., 2020) and continuum mechanics (Hassinger et al., 2017; Ma and Berro, 2021; Walani et al., 2015; Agrawal and Steigmann, 2008; Rangamani et al., 2013) are two common approaches. Hybrid models were also broadly applied to gain higher resolution than continuum mechanics and lower computing expense than molecular dynamics (Fu et al., 2019, 2021). However, most theoretical investigations have focused on how mechanical properties, such as membrane tension and bending rigidity of the clathrin coat, influence the membrane morphology. The process of curvature generation is either neglected or taken for granted. Only few of them have addressed the difference between the constant area model and the constant curvature model.

In fact, the constant curvature and constant area models are only two extreme models for clathrin assembly during endocytosis and any change in area or curvature are possible at any time point during endocytosis. In this paper, we extend the classic Helfrich theory for membrane deformation to incorporate energy terms associated with clathrin assembly and curvature generation, and compare geometric features calculated by theory with that extracted from experimental data. The negative gradient of the total free energy defines a pathway that neither fits the constant area model nor the constant curvature model. We find that a pathway that is close to the constant area model fits electron tomograms of the endocytic pits the best. Our study also offers experimental suggestions to distinguish between the two main hypotheses.

Models and methods

We model the membrane patch of the CCP (clathrin-coated pit) as a surface which is rotationally symmetric with respect to the z-axis. The shape of the membrane is parameterized with the meridional curve {r(s), z(s)}, where s denotes the arc length along the curve. The bending energy of the membrane (together with the clathrin coat) assumes the Helfrich model (Helfrich, 1973)

where k denotes the bending rigidity of the CCP, C₁ and C₂ denote the principal curvatures of the surface, C₀ denotes the intrinsic curvature of the membrane induced by the clathrin coat. To model a finite area of the clathrin coat, we assume the intrinsic curvature C₀ spatially varies as

where a denotes positions on the membrane. Here we choose a to be the surface area calculated from the tip of the membrane, and C₀ equals c₀ for area a < a₀, and rapidly drops to zero when a > a₀. The parameter α controls the sharpness of the drop. In the constant area model, we vary the intrinsic curvature c₀ but keep the coating area a₀ constant, while in the constant curvature model, we vary the coating area a₀ but keep the intrinsic curvature c₀ constant. As a result of clathrin coat, the bending rigidity k also varies as

where k_coat and k_bare denotes the bending rigidity of the clathrin-coated membrane and bare membrane, respectively. Besides the bending energy, the membrane tension contributes to the free energy in the form of

where σ_e denotes the membrane tension at the base and A denotes the surface area of the membrane patch within a fixed radius of R_b. The membrane tension σ_e and bending rigidity k defines a characteristic length (Derényi et al., 2002). The total free energy E_tot = E_b + E_t is a functional of the membrane shape. Given a coating area a₀ and an intrinsic curvature c₀, we numerically solve the variational equations of the energy functional to obtain membrane shapes that minimizes E_tot. More detailed descriptions of the model can be found in Appendix 1.

Results

Difference between the constant area model and the constant curvature model in terms of membrane morphology

Vesiculation requires assembly of a clathrin coat on the membrane, as well as curvature generation from the clathrin coat. A vesiculation process defines a pathway in the phase space (a₀, c₀) of the clathrin coat area a₀ and the intrinsic curvature c₀ of the coat. The constant area model and the constant curvature model are pathways that are made of a vertical line and a horizontal line. Besides these two extreme cases, there is a continuum spectrum of pathways with simultaneously increasing coating area a₀ and intrinsic curvature c₀ that could lead to vesiculation. Along the pathway, the membrane evolves from a flat shape to a dimple shape, and finally to an Ω-shape, as shown in Figure 2. Hereafter we use the maximal tangential angle Ψ_max of the membrane as an indicator of the progression of vesiculation - when the membrane is flat, Ψ_max = 0°, and when the membrane becomes spherical, Ψ_max = 180°. In our simulation, the neck becomes extremely narrow before Ψ_max = 180°. Therefore, we consider vesiculation occurs when Ψ_max reaches 150° (Figure 2a).

Figure 2.

First, we analyze the difference between the two models in terms of membrane morphology evolution along their pathways. We fit a circle around the membrane tip and use the radius R_t of the circle to characterize the curvature of the membrane at the tip. When R_t is plotted against the maximal angle Ψ_max, both models show that R_t decreases with increasing Ψ_max (Figure 2b and c). We stress that even though the intrinsic curvature c₀ is fixed in the constant curvature model, it does not imply the tip radius along the vesiculation pathway is a constant. When the coating area is small, the geometric curvature at the membrane tip remains small and differs from the intrinsic curvature. Tip radius in the constant curvature model decays more steeply with Ψ_max than in the constant area model. This difference becomes obvious when one plots the ratio of the tip radius at Ψ_max = 150° and Ψ_max = 30° (Figure 2b and c insets). For the constant area model, the ratio approximately equals to a constant 0.268 regardless of the coating area a₀, which agrees well with the analytical result (See Appendix 5), while for the constant curvature model, the ratio stays above and approaches 1 at a certain c₀ value.

Another difference between the two models is the evolution of the projected area of the clathrin coat on the substrate. We use the maximal radius R_coat of the membrane within the clathrin-coated area as the indicator of the projected area (Figure 2a). In the constant area model, R_coat decreases with increasing Ψ_max, while in the constant curvature model, R_coat increases with Ψ_max and reaches a plateau (Figure 2d and e). The ratio R_coat (150°)/R_coat (90°) is about 0.732 in the constant area model and around 1 in the constant curvature model (Figure 2d and e, insets). The analytical curves of R_t and R_coat against Ψ_max also fit perfectly with numerical solutions in the constant area model (Figure 2b and d, compare dotted and solid curves). Our calculations therefore demonstrate that the two models exhibit clear differences in the evolution of R_t and R_coat at the beginning of endocytosis (i.e. when Ψ_max is small) which can be determined from shapes of endocytic pits obtained experimentally.

To demonstrate how the coating area a₀ and the intrinsic curvature c₀ of the clathrin coat influence the membrane morphology, for each pair of (a₀, c₀), we calculate the corresponding membrane shapes and plot the contour lines for Ψ_max which indicate the stage of endocytosis (Figure 2f). The contour line with Ψ_max = 150° represents the critical line where vesiculation occurs. The line can be well fitted by the analytical expression (Figure 2f, thick black line, see Appendix 4). It implies that a small intrinsic curvature of the clathrin coat is able to induce vesiculation of the membrane with a large clathrin coat. We find that the distances between contour lines for higher values of Ψ_max is smaller than those for lower values of Ψ_max. It means that at the late stage of endocytosis, a small change in a₀ and c₀ could result in a more dramatic change in the membrane shape than that at the early stage. This trend is demonstrated clearly from the orange dots in Figure 2f, which correspond to shapes in Figure 2a.

Note that in order to produce the phase diagram (Figure 2f), it requires that the bending rigidity of the clathrin coated area k_coat is significantly larger than k_bare in the uncoated area. If k_coat is comparable with k_bare, there exists a region in the phase diagram in which a single (a₀, c₀) corresponds to three possible membrane shapes (See Appendix 2—Figure 1a-d). Physically, it implies a discontinuous transition in the membrane shape along a path that passes through this region, and a gap in the maximal angle Ψ_max would appear. Because in experiments, a wide spectrum of Ψ_max are observed and no gap in Ψ_max is found (Avinoam et al., 2015), we keep k_coat much greater than k_bare for the rest of the paper. In this regime, the membrane shapes evolve continuously along any pathway that connects the origin (0, 0) with a point on the critical vesiculation curve.

Vesiculation needs free energy sources to drive clathrin assembly and curvature generation

In the previous section, we take curvature generation in the constant area model and clathrin assembly in the constant curvature model for granted, so that the coating area a₀ and the intrinsic curvature c₀ are imposed, such that the physical forces behind clathrin assembly and curvature gen-eration are ignored. However, the bending energy E_b and the tension energy E_t typically increase with a₀ and c₀ along a vesiculation pathway. Therefore, vesiculation will be energetically unfavorable if no additional free energy sources are provided. In this section, we extend the model to include free energy terms for the assembly of the clathrin coat and its reorganization for curvature generation.

To describe the assembly of the clathrin coat in the constant curvature model, we introduce

where μ denotes the effective surface binding energy density of clathrin molecules with the membrane. This term reduces the free energy with increasing a₀, therefore, driving the assembly of clathrin. We can identify three types of free energy curves for different assembly strength μ: (1) When μ is small, the total free energy E_tot = E_b + E_t + E_a as a function of the coating area a₀ has two local minima, with the lowest one at a small a₀ and the other one at the maximum a₀ where vesiculation occurs (Figure 3a, red curve). The minima are separated by an energy barrier that is significantly higher than the thermal energy k_BT and the clathrin coat would assemble to a small area and halt. (2) With increasing μ, the lowest free energy minimum is shifted to the vesicula-tion point, but the energy barrier still exists and the clathrin coat remains small (Figure 3a, orange curve). (3) Vesiculation could happen for large enough μ such that the energy barrier vanishes and the free energy E_tot monotonically decreases with a₀ (Figure 3a, green). Based on the above analysis of the energy landscape we construct the phase diagram of the constant curvature model with clathrin assembly in the phase space of (c₀, μ) and classify the points into four types. Besides the three types mentioned above, when the intrinsic curvature c₀ is small, increasing a₀ to its maximum value (10⁵nm²) cannot produce vesiculation. (Figure 3b, gray region). The critical assembly energy density μ at which the energy barrier vanishes is found to increase with the intrinsic curvature c₀ (Figure 3b, interface between the green region and the orange region), which implies that a larger assembly strength of clathrin coat μ is needed to complete vesiculation if the clathrin coat has a higher intrinsic curvature c₀. When comparing the contour lines of the energy barrier ΔE_tot = 1k_BT and ΔE_tot = 10k_BT (Figure 3b, dotted curve and dash-dotted curve), the gap between them increases with c₀, which means that the energy efficiency is reduced with c₀ in the sense that, for larger c₀, a larger increase in μ is needed to reduce the same amount of free energy.

Figure 3.

We next consider curvature generation in the constant area model. As the molecular mechanisms of curvature generation of the clathrin coat remains debated, we introduce a phenomenological model in which the free energy has the general form,

where v denotes the strength of curvature generation, m and n are two positive numbers that are associated with the molecular mechanisms of curvature generation. The free energy E_c in Equation 6 decreases with increasing c₀, therefore, driving curvature generation. We set m = 1 such that E_c is proportional to the coating area. Note that m cannot be zero, otherwise, E_c only depends on the intrinsic curvature c₀ and can be nonzero even when the coating area is zero. As for the power n of the intrinsic curvature c₀, we set n = 1 or 2 (called Model(1,1) and Model(1,2), respectively). Physically, Model(1,2) implies cooperativity in the curvature generation such that the reduction of free energy per increase of unit curvature is proportional to the current curvature, i.e., ΔE_c ∝ −c₀Δc₀, while in Model(1,1), the reduction of free energy per increase of unit curvature is independent of current curvature. For Model(1,1), when v is small, the total free energy E_tot = E_b + E_t + E_c as a function of the intrinsic curvature c₀ has two minima, the lowest one at a small positive c₀ and the other one at the maximum c₀ where vesiculation occurs (Figure 3c, red curve). Further curvature generation is strictly limited by the high energy barrier (sometimes more than 100k_BT) between the two minima. With increasing v, the lowest minimum shifts to the vesiculation point, but the energy barrier still prevents curvature generation (Figure 3c, orange line). For a large enough v, the free energy monotonically decreases with c₀ and the curvature generation proceeds until vesiculation occurs (Figure 3c, green curve). When the coating area a₀ is very small, vesiculation fails to occur even when the intrinsic curvature is increased to its maximum value (0.125nm⁻¹) (Figure 3c, gray region). In the phase space of (a₀, v), the critical value of v where the energy barrier vanishes increases with the coating area a₀ (Figure 3d, interface between the orange region and the green region), which implies that a larger clathrin coat needs a stronger strength of curvature generation to complete vesiculation.

Model(1,2) has similar free energy landscape as Model(1,1) (Compare Figure 3c and e, d and f). However, in Model(1,2), for very small v, the lowest free energy minimum is strictly pinned at c₀ = 0, which implies no spontaneous curvature generation. In contrast, the minimum is at a small positive c₀ in Model(1,1), which indicates slight curvature generation.

Determination of the vesiculation pathway from the energy landscape

In this section, we combine the assembly energy Equation 5 and the curvature generation energy Equation 6 together and calculate the total free energy E_tot (a₀, c₀) = E_b +E_t +E_a +E_c as a function of both the coating area a₀ and the intrinsic curvature c₀. A pathway from the origin can be constructed by the descent along the negative gradient of the free energy landscape −∇E_tot. In Figure 4a we show the free energy landscape for a fixed assembly strength and varied reorganization strength for model(1,2). When v is small, the energy contour lines near the origin are kinked, which represents an energy barrier that prevents the path from going up, i.e. from generating curvature. The path extends horizontally and terminates on the a₀ − axis (Figure 4a, first column, red curve). With increasing v, the kinked contour lines shift towards larger a₀ and the path can be lifted up to in the middle and drops to the a₀ − axis in the end (Figure 4a, second column, orange curve). Beyond a critical v, the energy barrier vanishes and the path bends up and terminates on the vesiculation curve (Figure 4a, third column, green curve). Further increasing μ leads to the path lifting up at a smaller coating area a₀ (Figure 4a, fourth column, green curve). The path goes horizontally first and is later lifted up, which resembles the path of the constant area model. The three types of pathways are classified into three colored regions in the phase diagram of (μ, v), which represent complete vesiculation (green), partial vesiculation (orange) and no vesiculation (red), respectively (Figure 4c).

Figure 4.

The free energy landscapes of Model(1,1) dramatically differs from Model(1,2). When v is small, the energy gradient is strongly biased towards the horizontal direction, and the path extends horizontally with little or no curvature generation (Figure 4b, first column). For an intermediate v, the path first goes towards the top right direction until and then slowly bends down and extends towards large coating area along a valley formed in the energy landscape, which corresponds to membrane shapes with a small dimple (Figure 4b, second column). For large enough v, the path shoots nearly straightly towards the top right direction before it reaches the vesiculation line (Figure 4b, third column). Further increasing v makes the path more straight and terminates at a smaller coating area (Figure 4b, fourth column). The (μ, v) phase diagram shows the parameter regions that lead to complete, partial or no vesiculation for Model(1,1) (Figure 4d).

Comparison between different models with the experimental data

The constant area model and the constant curvature model represent two extreme pathways of membrane vesiculation. We have found constant-area-like pathways in Model(1,2) and straight-line-like pathways in Model(1,1). In order to understand which model is the most plausible, we compare membrane shapes predicted by the models with the membrane profiles obtained by electron microscopy in (Avinoam et al., 2015). The fitting error ϵ of a vesiculation path reflects the relative difference between the model-predicted geometric features along the path and the rolling median of the corresponding experimental data (Figure 5 and Appendix 3). The fitting geometric features include neck width, tip radius and invagination depth (Figure 5c). We draw the corresponding optimum energy paths that minimize the fitting error (Figure 5b), and compare the best model-predicted shapes with the experimental ones (Figure 5d).

Figure 5.

For Model(1,2) and Model(1,1), the fitting parameters include the polymerization strength and the reorganization strength which together determine the vesiculation pathway, as well as the characteristic length L₀ which scales the size of the membrane. For Model(1,2), the best fits are obtained for and L₀ = 30nm (Figure 5a). The resulting path moves horizontally at first and then bents up vertically (Figure 5b), which resembles the behavior of the constant area model. For Model(1,1), the best fitting parameters are and L₀ = 40nm (Figure 5a). The fitting path is almost a straight line towards the vesiculation line (Figure 5b). The optimum fitting error of Model(1,2) (ϵ = 0.14) is slightly better than that of Model(1,1) (ϵ = 0.17).

We also perform the fitting procedure to the constant area model and find the optimum parameters are a₀ = 1.69 × 10⁴nm² and L₀ = 30nm. For the constant curvature model, the best fitting parameters are c₀ = 0.043nm⁻¹ and L₀ = 50nm. The minimum fitting error of the constant curvature model (ϵ = 0.28) is exactly twice as large as that of the constant area model (ϵ = 0.14) (Figure 5a and b). So considering the fitting error and the pathway in (a₀, c₀) phase diagram, we raise the conclusion that the experimental vesiculation process probably favors constant-area-like pathways.

When comparing the model-predicted geometric features with the rolling median of the experimental data, we find that the four models fit almost equally well the experimental data for the neck width(Figure 5c left). Model(1,2) and the constant area model predict very similar results such that the curves almost overlap with each other (Figure 5c, red curve and orange curve). The predictions of these two models match the rolling median of the experimental data best. The constant curvature model strongly deviates from the rolling median of the experimental tip radius, particularly in the early stage of vesiculation when Ψ_max < 90°.

To compare the axisymmetric membrane shapes predicted by the models with the non-axisymmetric membrane profiles obtained with electron microscopy, we symmetrize the experimental data with a procedure (Appendix 3—Figure 1). Then, we average the symmetrized profile within an interval of Ψ_max ∈ [Ψ₀ −10°, Ψ₀ +10°] and overlay the averaged profile with model-predicted shapes for Ψ_max = Ψ₀ (Figure 5d). At the early stage when the membrane exhibits a dimple shape (0° < Ψ_max < 60°), the membrane morphology predicted by the constant curvature model is distinguishable from the other three models, particularly when looking at the tip radius. At the late stage when the membrane exhibits an Ω-shape, i.e., Ψ_max > 90°, the difference in shape between models is mainly manifested in the invagination depth. The constant curvature model and Model(1,1) mainly predict a deeper invagination depth than the symmetrized experimental profile, while the constant area model and Model(1,2) usually give much better fitting.

Discussion

Three types of clathrin coats

In this paper, we have constructed a physical model to describe how curvature generation and clathrin assembly are interrelated during the vesiculation process in CME. Previous experiments have reported three groups of clathrin coated pits, which are plaques, abortive pits and pits that lead to vesiculation, according to their structural and dynamic properties (Maupin and Pollard, 1983; Ehrlich et al., 2004; Loerke et al., 2009; Saffarian and Kirchhausen, 2009; Kirch-hausen, 2009; Lampe et al., 2016).

In Figure 4, we show that depending on the clathrin assembly strength μ and its reorganization strength v, the pathway might end up with three possible final shapes: (i) a flat membrane with no curvature generation, (ii) a nearly flat membrane with small curvature generation, (iii) a spherical cap that leads to vesiculation. They essentially correspond to the three types of clathrincoated pits found in experiments. Based on the phase diagram of Model(1,2) (Figure 4c), the difference between the three groups comes from the difference in the assembly and reorganization strengths of clathrin molecules. Furthermore, Model(1, 2) predicts that at the boundary between the type (iii) region and the type (ii) region, the reorganization strength v increases with the assembly strength. This result has important implications to explain why large plaques are commonly observed in experiments. The large area of the plaques are due to the strong assembly strength μ. However, for these plaques to go to vesiculation, a strong reorganization energy v is also needed. Therefore, the combination of a strong μ and weak v leads to the formation of large plaques. Model(1,2) predicts that a plaque or an abortive pit can be transformed into a vesicle by either increasing the reorganization strength or reducing the assembly strength (Figure 6). The former ends up with a large vesicle and the latter ends up with a small vesicle. This can be used as a test of our theory with experiments to modify the binding affinity of clathrin molecules with adaptor proteins on the membrane. Weakening the affinity might increase the portion of vesicles and reduce the portion of plaques, though the vesicles would become smaller.

Figure 6.

Cooperativity in the curvature generation process

In Figure 5, we show the best fitting results for all the four models and find that Model(1,2) produces better fits than Model(1,1), which suggests the existence of cooperativity in the curvature generation process. In particular, if curvature generation is driven by breaking bonds in the hexagonal lattice, cooperativity implies that the number of newly broken bonds is proportional to the number of already broken bonds. Because of this cooperativity, at the early stage of endocytosis bonds are broken slowly and clathrin assembly dominates over curvature generation. At the late stage of endocytosis, an increasing number of bonds are broken and curvature generation could happen rapidly and dominate over clathrin assembly. Altogether, this cooperativity leads to a constant-area-like behavior. Similar effect have been reported in (Mund et al., 2023).

The difference in membrane morphology between the different models is most salient at the early stage of endocytosis

When we compare the model predictions, we find that the difference in membrane morphology between models is not as big as expected, which might explain why it has been difficult to distinguish between the constant area and the constant curvature models for so long. For instance, the neck width vs. Ψ_max and the invagination depth vs. Ψ_max are similar for all the four models (Figure 5c, left and right). The best fitting error of the four models calculated from the geometric features are relatively close, except for the constant curvature model, which gives the worst fit (Figure 5a). The models are mainly distinguishable from the tip radius vs. Ψ_max plot at the early stage of endocytosis when the membrane is nearly flat (Figure 5c, middle), i.e. for shapes with small Ψ_max. However, published experimental shape at early stages of endocytosis are sparse. Our result hints that in order to distinguish between the models, collecting membrane shapes at the early stage is necessary and the relation of tip radius vs. Ψ_max is the key geometric feature to tell the models apart.

The projected area of clathrin coat in the plane of the plasma membrane could distinguish between the two models

In Figure 2d and e we have shown that the coat radius R_coat, which represents the projected area of the clathrin coat in the plane of the plasma membrane, as a function of Ψ_max exhibit opposite trends for the constant curvature model and the constant area model. The results suggest that in experiments the projected area for the constant area model would first increase and then decrease over time, finally reaching a plateau. However, for the constant curvature model, the projected area would increase over time and finally reach a plateau without a decreasing phase. This result suggests another method to distinguish between the two models via the projected area measurement. The idea has been used in a study where the clathrin-coated pit was imaged with platinum replica and cryoelectron microscopy and tomography (Sochacki et al., 2021). The results support a constant-area-like model, consistent with the prediction of our Model(1,2), in which the dome structures have a slightly larger coating area than flat structures. On the other hand, another study has used the super-resolved live cell fluorescence imaging with TIRF to measure the growth of the clathrin coat over time (Willy et al., 2021). The authors found a smooth drop in the projected area of clathrin coat over time. However, based on a computer simulation of clathrin assembly, they concluded that the smooth drop of the projected area was the result of a constant-curvature-like model because a constant-area-like model would produce a sharp drop. We attribute the difference between their model and our model to the fact that they model the clathrin coat as a discrete lattice while we use a continuum mechanics method. More importantly, in their model, the moment at which curvature generation occurs was arbitrarily imposed at 80% of clathrin triskelions. If the transition were chosen to occur with fewer triskelions, e.g. 40%, the sharp drop in the project area might not happen in the constant-area-like model. Furthermore, the authors used the number of triskelions to monitor the progress of endocytosis which terminates when the triskelions reach the maximum number. This choice might bias towards the constant-curvature-like model because the vesiculation may not happen at all when the clathrin numbers reaches its maximum.

The bending rigidity of the coated area should be much larger than the uncoated area

Comparison of our model to experimental data demonstrates that the relative bending rigidities of the coat and the membrane are constrained. Indeed, if k_coat /k_bare < 50, the model predicts an abrupt change (or a gap) in Ψ_max at the end of vesiculation (See Appendix 2—Figure 1), i.e., a snap-through transition (Hassinger et al., 2017). If a gap in Ψ_max existed, we would expect that the distribution of experimental shapes to be discontinuous, with no or very few data corresponding to a certain range of Ψ. However, in the experiments (Avinoam et al., 2015), the endocytic pits shapes are continuously distributed across the Ψ_max spectrum, indicating the ergodicity of Ψ_max during the endocytic process. Our calculation suggests that the clathrin coat is about 50 times stiffer than the membrane.

Acknowledgements

We thank Prof. Ori Avinoam and Marko Kaksonen for generously sharing their data with us. R.M. acknowledges financial support from Fundamental Research Funds for Central Universities of China under Grant No. 20720240144. Part of this work was funded by NIH R01 grant GM115636 awarded to J.B.

Appendix 1

Detailed formula derivation

We assume the membrane shape is axisymmetric and is parameterized with its meridional curve {r(u), z(u)}, with u = 0 corresponding to the membrane tip and u = 1 corresponding to the membrane edge. Our aim is to derive the shape equations via minimizing the bending energy of the membrane under certain geometric constraints (Jülicher and Seifert, 1994; Seifert et al., 1991; Zhong-Can and Helfrich, 1987).

Hereafter we use f ‘ ≡ df/du to denote the derivative of an arbitrary function f with respect to the rescaled arclength u. It should be noticed that if {r(u), z(u)} describes a membrane shape, {r[g(u)], z[g(u)]} describes the same membrane shape if the functiong maps the interval [0, 1] to [0, 1], e.g., g(u) = u². In order to fix the issue, we introduce and impose that h is a constant. By this definition, we essentially let u = s/S, where s is the arclength calculated from the membrane tip and S is the total arclength. The constant h = S is an unknown parameter to be determined via solving the shape equations. In order to simplify the form of the bending energy, we introduce the tangent angle Ψ(u) which satisfies the geometric relation:

and

The two principal curvatures can be expressed as

The bending energy of the membrane then reads

Note that the bending rigidity k[a(u)] and the intrinsic curvature C₀[a(u)] are functions of the area a(u), which fullfils the equation

This relationship means that we study an inhomogeneous membrane which is locally incompressible in its area. In order to impose the geometric relation Eqs. (7),(8) and (11), we introduce three Lagrangian multipliers γ(u), η(u) and σ(u) to the integral

The total free energy to be minimized under the geometric constraints reads

in which ε reads

The variation of the functional E_tot in Equation 13 reads

which contains both the bulk terms (first line) and the boundary terms (second line). The Euler-Lagrange equations can be obtained by the vanishing of the former 7 bulk terms, which are reduced to

and

as well as Eqs. (7), (8) and (11). Note that the vanishing bulk term of δh gives us

which is a boundary condition (not one of the differential equations).

Next, we need 9 boundary conditions (BCs) to numerically solve 6 first-order differential equations, 1 second-order differential equation together with 1 unknown parameter. At the membrane tip u = 0, we can easily obtain r(0) = 0, Ψ(0) = 0, a(0) = 0 by geometric relations. Then, ∂ε/∂h = 0 given by Equation 20 is satisfied at u = 0. At the membrane base u = 1, we can also easily get r(1) = R_b, z(1) = 0 by geometric relations, where R_b is the distance from the boundary u = 1 to the axisymmetric axis. All boundary terms in Equation 15 vanish simultaneously because geometric relations lead to δf = 0, except δz at u = 0 and δΨ at u = 1. This forces ∂ε/∂z’ = 0, which is, η(0) = 0, and ∂ε/∂Ψ’ = 0, and leads to

Surface tension is fixed to σ = σ_e at the base. The definition of the characteristic length L₀ is , depicting the radius of a tube of membrane elongated by a point force.

In summary, all the BCs are listed below

Appendix 2

Gap of maximum angle

All of our results in the main text were computed under the assumption of a very large ratio k_coat /k_bare = 50 to prevent a discontinuous jump of Ψ_max. When k_coat /k_bare = 5, the gap of Ψ_max is observable in some (a₀, c₀) sets (Appendix 2—Figure 1). For the constant area model, the gap is observed around Ψ_max = 125° when a₀ is above a limit value and increases slightly with increasing a₀ (Appendix 2—Figure 1a). The solid and dotted lines deviate further with an increasing a₀. The dotted free energy curve with larger a₀ generates a Gibbs triangle (Appendix 2—Figure 1c). The bottom point of the Gibbs triangle is the phase change point and corresponds to the c₀ value where the gap of Ψ_max is situated on. When a₀ is less than the limit value, the free energy curve is smooth and has no obvious phase change point. Result in the constant curvature model is qualitatively similar to the constant area model (Appendix 2—Figure 1b and d). If Ψ_max gap exists, the dotted line and solid line intersect at three multiple-solution points with the same c₀ and a₀.

The (a₀, c₀) phase diagram shows how both arguments effect on Ψ_max gap (Appendix 2—Figure 1e).

The orange lines shows the deviation and variation trend between the dotted lines and the solid lines in Appendix 2—Figure 1a and b. Note that the deviation of any path that starts from the origin and terminates on the vesiculation boundary can be described by orange lines, rather than just constant curvature paths or constant area paths. A physical vesiculation pathway just goes across three orange lines directly and terminates at the vesiculation boundary while a numerical one turns back at the upper line, then turns back at the lower line, and finally completes vesiculation. Note that some section of the upper dotted curve overlaps with the numerical vesiculation boundary, because the system reaches vesiculation boundary before passing through the phase change point. Finally, the introduction of the assembly energy and the reorganization energy have no effect on Ψ_max gap, and just change the shape of the Gibbs triangle.

Boundary radius value

We set throughout our study. Values of tip radius and neck radius appear almost identical for different R_b in both models (Appendix 2—Figure 2). However, membrane heights have distinct differences for different R_b when Ψ_max is in the middle range. This difference is likely due to fact the uncoated area contains a smaller region able to generate curvature and does not contribute in lifting the membrane center at this stage. From to , the curves are virtually identical, even for membrane height. Since a larger R_b requires a larger number of mesh points in the simulations, we chose to balance computation time and numerical accuracy.

Appendix 2—figure 1.

Appendix 2—figure 2.

Appendix 3

Symmetrization algorithm

In our investigation to determine the model parameters from experimental data, we need to symmetrize the shapes extracted from electron tomograms using the paired coordinates (r_i, z_i)(Appendix 3—Figure 1a is an example). Firstly, we define the vertical line that crosses the maximum z-value of the shape as its axisymmetric axis, and normalize the left- and right-poritons of shapes using the same rescaled mesh points u = s_i/S_i. Secondly, we average r and z with the same u on both sides to achieve the symmetrized shape (Appendix 3—Figure 1b).

Rolling median calculation

We note (x₁, y₁), (x₂, y₂), …, (xi, y_i), …, (x_N, y_N) the coordinates of all N data points in a given figure (Fig. 5c), sorted by independent variable xi in an ascending order. We calculate the median points of ten consecutive neighbouring points (xi, y_i), …, (xi₊₉, y_i+9) by x and y, respectively. We connect consecutive median points to plot the rolling median line. The same procedure is performed for each of the three parameter figures.

The symmetrized experimental shapes are grouped by their Ψ_max value in eight intervals of equal width ranging from Ψ_max = 0° to Ψ_max = 160° (Figure 5d). A single range of Ψ_max contains n shapes {r, z}₁, …, {r, z}_n. The {u}_n meshes being exactly the same for each n (See previous section), we calculate the average values in each Ψ_max interval as

to obtain eight average shapes . Then, we draw the shapes of four theoretical models with the midpoint Ψ_max values for each interval to compare our theoretical shapes with the averaged experimental shapes.

Error calculation

We compare four models with the rolling median lines, and draw the relationship between total relative fitting error and parameters (Figure 5a and c). In the Ψ_max interval where the-oretical and experimental data points are well-defined, we use an interpolation method to map their Ψ_max into the same mesh points. The total relative fitting error is calculated by

where N_i is the number of well-defined experimental data points in the i-th parameter figure, and is the k-th theoretical value in the i-th para figure, and ^{is the k-th experi-} mental value in the i-th para figure. We choose the parameter sets with minimum ϵ in four theory models. Then, we compare the curves and shapes obtained using the best-fit parameters with the experimental data. We only show the error-parameter figure with the best-matched L₀ value in Fig. 5a.

Normalization

We impose Ā as the dimensionless variable of A and A_unit as the normalizing units, such that A/A_unit = Ā. Units of model variables are listed in Table 1. Other variables are derived variables and are not listed in the table.

Table 1.

Appendix 3—figure 1.

Appendix 3—table 1.

Appendix 4

Critical vesiculation curve

We approximate the system as a combination of a spherical cap and a plane (Appendix 4—Figure 1, cross-section view). The undeformed shape of the coating area is a spherical cap with curvature c₀ and the uncoated area is flat. For simplicity, We postulate the spherical cap deforms uniformly and the plane has no deformation. From the geometric relationship, the area can be written as

where A_s is the area of spherical cap and A_p is the area of plane. Correspondingly, their surface bending energy density expressions are

The uncoated area has no deformation so it has no bending energy. The total free energy includes the surface tension, so

where λ is the Lagrange multiplier to give geometric restriction on A_s = a₀, and σ = σ_e is a constant surface tension of uncoated area. Physically, the system prefers the shape that minimizes E_tot, i.e.

Eliminating λ and R and ignoring solutions where θ < 0, which are not physiclaly possible, and θ = π, which can only be achieved after passing an energy barrier, we obtain

Then, we postulate that the rigidity of the solid shell is very large, giving 8πk_coat ≫ a₀σ. We can then use the approximation:

Then, we consider a closed sphere as the vesiculation state, where θ = π. Finally, we substitute , to obtain the dimensionless relation

which is the analytical vesiculation boundary of this simplified model.

Appendix 5

Model fitting

Using the approximated model (Appendix 4—Figure 1), we provide an analytical result to distinguish between the constant area model and the constant curvature model (Figure 2b-e). This result is confirmed by numerical analyses and can be used to differentiate the two models from experimental data. The area of a spherical cap is

In the constant area model, R and θ vary during the vesiculation process but a₀ remains contant. Equation 32 gives

Therefore, we calculate the ratios R_t (150°)/R_t (30°) ≈ 0.268 and R_coat (150°)/R_coat (90°) ≈ 0.732, which both are close to the results from the numerical simulations. In the constant curvature model, θ and a₀ vary during the vesiculation process but R remains constant, so

which holds true for any θ. So, in the analytic model, R_t (150°)/R_t (30°) = 1 and R_coat (150°)/R_coat (90°) =

The latter ratio matches the numerical results well but the former ratio doesn’t because of the very large tip deformation that significantly differ from a spherical cap in some shapes determined numerically. Using different pairs of Ψ_max values, R_t and R_coat ratios are different in both models, and these ratios can be used as indicative variables to distinguish between the two models from the experimental shapes.

In addition, the simplified model postulates R = R_t, θ = Ψ_max and Equation 32 leads to

For R_coat the analytical expressions are

Equation 35 and Equation 36 both fits the numerical results well (Figure 2b-e), proving the validity of the spherical cap approximation.

Appendix 5—figure 1.