Role of Ligand Distribution in the Cytoskeleton-Associated Endocytosis of Ellipsoidal Nanoparticles

Nanoparticle (NP)–cell interaction mediated by receptor–ligand bonds is a crucial phenomenon in pathology, cellular immunity, and drug delivery systems, and relies strongly on the shape of NPs and the stiffness of the cell. Given this significance, a fundamental question is raised on how the ligand distribution may affect the membrane wrapping of non-spherical NPs under the influence of cytoskeleton deformation. To address this issue, in this work we use a coupled elasticity–diffusion model to systematically investigate the role of ligand distribution in the cytoskeleton-associated endocytosis of ellipsoidal NPs for different NP shapes, sizes, cytoskeleton stiffness, and the initial receptor densities. In this model, we have taken into account the effects of receptor diffusion, receptor–ligand binding, cytoskeleton and membrane deformations, and changes in the configuration entropy of receptors. By solving this model, we find that the uptake process can be significantly influenced by the ligand distribution. Additionally, there exists an optimal state of such a distribution, which corresponds to the fastest uptake efficiency and depends on the NP aspect ratio and cytoskeleton stiffness. We also find that the optimal distribution usually needs local ligand density to be sufficiently high at the large curvature region. Furthermore, the optimal state of NP entry into cells can tolerate slight changes to the corresponding optimal distribution of the ligands. The tolerance to such a change is enhanced as the average receptor density and NP size increase. These results may provide guidelines to control NP–cell interactions and improve the efficiency of target drug delivery systems.

Recently, Schubertová et al. [33] focused on the specific case of receptors being overexpressed in cancer cells and demonstrated that spherical NPs with uniform ligand distribution yield the fastest wrapping speed by the cell membrane through performing a molecular dynamic simulation, in which the receptors can directly react with the ligands on the surface of NPs without diffusion. While the receptors have relatively low density, NP endocytosis is generally limited by receptor diffusion [4,8,34]. In the case of receptordiffusion-mediated endocytosis of NPs, we have previously examined the effect of ligand distribution on wrapping cylindrical NPs by cell membrane and determined an almost uniform ligand distribution as the optimal distribution associated with the highest cellular uptake efficiency [35]. The above studies have focused on describing the effect of the ligand distribution on the internalizations of NPs with spherical and cylindrical shapes, both of which have surfaces of constant Gaussian curvature. In addition, cylindrical NPs have also been used to investigate how aspect ratios may affect the cell-NP interactions [7,36,37]. However, for NPs with complex shapes, it is still unclear how non-constant surface curvature may couple other factors like the distribution patterns of ligands and cytoskeleton remodeling to influence their cellular uptake.
Endocytosis includes a number of biophysical processes by which cells internalize NPs. It has been recognized that the cytoskeleton plays an essential role in several of these processes [38]. The molecular structure, on the basis of the association of the cytoskeleton to endocytosis, is made up of many endocytic proteins potentially linked to the actin cytoskeleton so that the NP entry via endocytosis essentially accompanies cytoskeleton remodeling [38][39][40][41][42][43][44]. Such a point of view has been confirmed through the observations of HIV viral particle entry into host cells, where a large number of protein-protein interactions are indicative of the links between the endocytic machinery and the actin cytoskeleton [45,46]. Most recently, while taking advantage of today's nanotechnology, an experiment has been performed on measuring the dynamic process of internalization of a single human enterovirus 71 with size 26 nm via endocytosis by a force tracing technique on atomic force microscopy (AFM) [47]. It is found that if cytochalasin B, a cell-permeable mycotoxin that strongly inhibits network formation by actin filaments, is injected into the AFM liquid chamber during force tracing measurement, most of the force signal disappeared, which indicates that the wrapping of the NP cannot occur without the association of the actin filament network [47]. From a mechanics point of view, the involvement of cytoskeleton remodeling in the endocytic process embodies mechanical resistance through elastic deformation [21][22][23]. This deformation of the cytoskeleton during the NP-cell interaction, at the length scale of about ten nanometers, usually follows the linear theory of continuum mechanics, as confirmed by the experiments on nanoneedle (radius from 20 to 150 nm) indentation to living cells [48][49][50].
Although tremendous progress in understanding cell-NP interactions has been achieved, the role of ligand distribution in the cytoskeleton-associated endocytosis of non-spherical NPs remains unclear. In this work, we propose a coupled elasticity-diffusion statistical dynamic model to investigate the influence of ligand distribution on the cellular uptake of ellipsoidal NPs, where we assume that the uptake process is driven by the binding energy between diffusive receptors on the cell membrane and ligands on the NP surface to overcome resistances from membrane deformation, cytoskeleton deformation, and changes in the configuration entropy of receptors. Figure 1 shows the schematic for the interaction between an ellipsoidal NP and a cell via cytoskeleton-associated endocytosis. The cell is modeled as an elastic half-space covered by a cell membrane embedded with diffusive mobile receptors. The ellipsoidal NP is considered as rigid and coated with immobile ligands. For simplicity, we consider the NP as an ellipsoid of revolution with equatorial and polar semi-axes, R a and R b , as shown in Figure 1a. We assume that the mobile receptors are distributed evenly with constant density ξ 0 prior to the engulfment. The ligand density ξ L (h), as a function of the central axis h, can be non-uniform. And the depth of the NP engulfment is denoted as h d . As long as the receptors diffuse to binding sites and bind with the ligands on the NP surface, the receptor density within the contact area becomes identical to the ligand density. The distribution of the receptors on the membrane surface then become non-uniform, which is denoted as ξ(s, t) at time t. During the cytoskeleton-associated endocytosis of NPs, engulfment is driven by the formation of ligand-receptor complexes to overcome the resistance from the membrane and cytoskeleton deformations, and the changes in the configuration entropy of receptors. The energies that determine the internalization process include: (1) , favorable binding energy of ligand-receptor; (2) , unfavorable energy of membrane deformation; (3) , unfavorable energy attributed to changes in the configuration entropy of receptors; and (4) , unfavorable energy of cytoskeleton deformation. The total energy for cellular uptake of NPs can then be expressed by,

Materials and Methods
The binding energy during the wrapping process can be obtained as, where is the contact area between the NP and cell, and is the chemical energy released by each receptor-ligand binding event. The axially symmetric system causes the infinitesimal area element to be expressed in the following form, The fluid membrane theory of Canham-Helfrich [51] states that during NP-cell interaction, the energy of membrane deformation is, where and are the bending modulus and surface tension of the membrane, respectively. Once the membrane is in close contact with the NP, its mean curvature, , becomes equal to that of the corresponding position of the NP, and can be expressed as (see Appendix A for details), During the cytoskeleton-associated endocytosis of NPs, engulfment is driven by the formation of ligand-receptor complexes to overcome the resistance from the membrane and cytoskeleton deformations, and the changes in the configuration entropy of receptors. The energies that determine the internalization process include: (1) F 1 , favorable binding energy of ligand-receptor; (2) F 2 , unfavorable energy of membrane deformation; (3) F 3 , unfavorable energy attributed to changes in the configuration entropy of receptors; and (4) F 4 , unfavorable energy of cytoskeleton deformation. The total energy for cellular uptake of NPs can then be expressed by, The binding energy during the wrapping process can be obtained as, where A w is the contact area between the NP and cell, and e RL is the chemical energy released by each receptor-ligand binding event. The axially symmetric system causes the infinitesimal area element to be expressed in the following form, The fluid membrane theory of Canham-Helfrich [51] states that during NP-cell interaction, the energy of membrane deformation is, where κ and γ are the bending modulus and surface tension of the membrane, respectively. Once the membrane is in close contact with the NP, its mean curvature, H, becomes equal to that of the corresponding position of the NP, and can be expressed as (see Appendix A for details), To estimate the entropic contribution to the energy associated with the receptor distribution, we adopted the model based on the theory of ideal gas [4,52]. The unfavorable energy of changes in the configuration entropy of the bonds and free receptors is given as, where T is the absolute temperature in the order of 300 K.
In the cytoskeleton-associated endocytosis, the entry of NPs into the cells inevitably deforms the cytoskeleton. Obataya et al. [48] used an atomic force microscope with a long nanoneedle (radius = 100 nm to 150 nm) to perform indentations on living cells. They demonstrated that the loading force curve for an indentation depth of up to 2 µm remains consistent with the classic Hertz model in contact mechanics. Similarly, Beard et al. [49] used a nanoneedle probe with a radius of only 20 nm and recorded the loading force curve during the indentation of a corneocyte. The corneocyte is usually a target of numerous viruses, including the Herpes simplex virus Type 1 [53]. This study [49] confirmed that the Hertz model can well fit the loading force curve. As the Hertz model is derived based on the contact problem of two elastic bodies, these results therefore clearly indicate that living cells during indentation behave as a deformable elastic solid, implying that cytoskeleton deformation plays an important role in resisting NP intrusion. Based on the Hertz contact theory for rotating solids [54], the deformation energy of cytoskeleton can be derived as, where R 2 a /R b is the curvature radius at the initial contact point, and the combined elastic modulus is defined as 1/E * = 1 − µ 2 c /E c + 1 − µ 2 n /E n , in which µ c and E c represent the Poisson ratio and Young modulus of the cell, and the corresponding values for NPs are µ n and E n . In the limit wherein NPs are stiffer than cells, the combined elastic modulus can be simplified to E * = E c / 1 − µ 2 c . The differentiation of free energy in Equation (1) with respect to time yields, (8) where µ = 1 + ln ξ/ξ 0 is the local chemical potential of a receptor, ξ + L ≡ ξ L (h d ) and ξ + ≡ ξ(a + , t) denote the receptor-ligand bond density and receptor density, respectively.
dh is the contact area, and the function f (h) associated with NP shape can be given as (see Appendix B for details), Receptors on the membrane are driven by the density gradient to the adhesion region, and this driving force is, Further, the diffusion flux can be given as, with the diffusion coefficient D.
Finally, owing to membrane fluidity, the free receptors outside the contact region are mobile within the membrane, and the receptor density can be governed by the diffusion equation as follows [4,52], The conservation condition of receptors on the membrane can be given as [4,52], Substitution of the continuity equation ∂ξ ∂s into Equation (13) under the boundary conditions ξ(s, t) → ξ 0 , j(s, t) → 0 as s → ∞ yields [4,52], where j + ≡ j(a + , t) denotes the flux in front of the contact edge. By incorporating the conservation relation and integrating by parts, the integral term in Equation (8) takes the form [4,52], The first term on the right side represents the energy transport across the front at s = a(t) due to receptor diffusion. The second term on the right in Equation (15) is precisely the rate of energy dissipation as a result of receptor diffusion along the cell membrane.
Similar to the Griffith crack growth criterion, the decreasing rate of free energy should equal to the energy dissipated during receptor diffusion for a power-balanced process [52]. Hence, the energy balance equation at the contact edge features the form, Given the ligand distribution on the NP's surface, the diffusion Equation (12) can be numerically calculated by adopting the finite difference method subjected to the boundary conditions ξ(∞, t) → ξ 0 and ξ + ≡ ξ(a + , t). Under the contact edge condition, wherein the contact area reaches the total ellipsoid surface area as A t = πa 2 = 4π R 2 a + 2R a R b /3, the NP is totally internalized, and the wrapping time is determined.

Results
To simulate the complex situation of non-uniform ligand distribution, we denote a distribution pattern of the ligand density with an approximate harmonic distribution, where the dimensionless constant c represents the degree of non-uniformity of ligand distribution. For example, the case of c = 0 means that the ligand is uniformly distributed with a density of ξ L0 . Under such a situation, the total number of ligands is independent of the constant c and is denoted by ξ L0 A t . In Figure 2, we illustrate how such a degree of non-uniformity of the ligand distribution c, determines the ligand distribution on an NP with different shapes.

Influence of Ligand Distribution on the Cellular Uptake of NPs with Different Shapes
We exclude the effects of the cytoskeletal deformation and set the initial receptor density as 0.01 to systematically investigate the effect of NP shape on uptake. Figure  3 shows the wrapping time as a function of the degree of non-uniformity of ligand distribution, c, for three NPs with different shapes. The wrapping time initially decreases and then increases as the constant increases; this trend indicates the existence of an optimal ligand distribution corresponding to the minimum wrapping time. This optimal distribution is influenced by the shape of the NP. For spherical NPs, the optimal ligand distribution approaches a uniform distribution, which is consistent with previous predictions by Li et al. [35] and Schubertová et al. [33]. However, when the NP shape changes from spherical to oblate, the optimal ligand distribution transitions from uniform to polarized at the edge of the oblate NP (that is, ℎ = ). For the cellular uptake of elongated NPs, the optimal ligand distribution becomes densely distributed at both ends of the elongated NP. For each case of the oblate and elongated NPs, a large amount of binding energy from the receptor-ligand binding events is required to consume the high bending energy of the local cell membrane with large deformation. Figure 3 also demonstrates that wrapping times are only slightly different for a large range of the constant c. When the degree of non-uniformity of ligand distribution, , is higher or lower than a threshold value, the wrapping process becomes difficult to complete. Then, we consider the typical values of the binding energy of a single bond e RL = 15 k B T [55], bending modulus of the cell membrane κ = 20 k B T [56], surface tension of the cell membrane γ = 0.005 k B T/nm 2 [22], diffusion coefficient of the receptors on the membrane D = 10 4 nm 2 /s [4,55], and ξ L0 = 5000/µm 2 [4].

Influence of Ligand Distribution on the Cellular Uptake of NPs with Different Shapes
We exclude the effects of the cytoskeletal deformation and set the initial receptor density as 0.01ξ L0 to systematically investigate the effect of NP shape on uptake. Figure  3 shows the wrapping time as a function of the degree of non-uniformity of ligand distribution, c, for three NPs with different shapes. The wrapping time initially decreases and then increases as the constant c increases; this trend indicates the existence of an optimal ligand distribution corresponding to the minimum wrapping time. This optimal distribution is influenced by the shape of the NP. For spherical NPs, the optimal ligand distribution approaches a uniform distribution, which is consistent with previous predictions by Li et al. [35] and Schubertová et al. [33]. However, when the NP shape changes from spherical to oblate, the optimal ligand distribution transitions from uniform to polarized at the edge of the oblate NP (that is, h = R b ). For the cellular uptake of elongated NPs, the optimal ligand distribution becomes densely distributed at both ends of the elongated NP. For each case of the oblate and elongated NPs, a large amount of binding energy from the receptor-ligand binding events is required to consume the high bending energy of the local cell membrane with large deformation. Figure 3 also demonstrates that wrapping times are only slightly different for a large range of the constant c. When the degree of non-uniformity of ligand distribution, c, is higher or lower than a threshold value, the wrapping process becomes difficult to complete.

Effect of Ligand Distribution on the Cytoskeleton-Associated Endocytosis
We consider the internalization of spherical NPs under the influence of cytoskeleton deformation and determine how the cytoskeleton stiffness may affect the optimal ligand distribution associated with the smallest wrapping time. Figure 4 displays the relation between the wrapping time and the degree of nonuniformity of ligand distribution, , for spherical NPs with a radius of 50 nm for entry into cells. It can be seen that the cytoskeleton deformation can significantly affect the uptake. In the membrane-mediated endocytosis, the optimal distribution of the ligands coated on the spherical NP surface tends to be uniform. With an increase of Young's modulus of the cytoskeleton, the optimal ligand distribution becomes dense at the NP ends with a large curvature, in which the large deformation energy of the cytoskeleton plays a key role in resistance to the cellular uptake. When NPs are uptaken by stiff cells, since a large amount of energy for elastic deformation needs to be overcome, long wrapping times will be needed. Similar results for the oblate ellipsoid and elongated NPs can be seen in Figure 5. Regardless, both Figures 4 and 5 show that the range of the ligand distribution constant, , for effective NP uptake decreases as the Young modulus of the cytoskeleton increases. These findings imply that the internalization of spherical NPs into a soft cytoskeleton can not only shorten the wrapping time but also broaden the range of ligand distribution for effective NP uptake.

Effect of Ligand Distribution on the Cytoskeleton-Associated Endocytosis
We consider the internalization of spherical NPs under the influence of cytoskeleton deformation and determine how the cytoskeleton stiffness may affect the optimal ligand distribution associated with the smallest wrapping time. Figure 4 displays the relation between the wrapping time and the degree of nonuniformity of ligand distribution, c, for spherical NPs with a radius of 50 nm for entry into cells. It can be seen that the cytoskeleton deformation can significantly affect the uptake. In the membrane-mediated endocytosis, the optimal distribution of the ligands coated on the spherical NP surface tends to be uniform. With an increase of Young's modulus of the cytoskeleton, the optimal ligand distribution becomes dense at the NP ends with a large curvature, in which the large deformation energy of the cytoskeleton plays a key role in resistance to the cellular uptake. When NPs are uptaken by stiff cells, since a large amount of energy for elastic deformation needs to be overcome, long wrapping times will be needed. Similar results for the oblate ellipsoid and elongated NPs can be seen in Figure 5. Regardless, both Figures 4 and 5 show that the range of the ligand distribution constant, c, for effective NP uptake decreases as the Young modulus of the cytoskeleton increases. These findings imply that the internalization of spherical NPs into a soft cytoskeleton can not only shorten the wrapping time but also broaden the range of ligand distribution for effective NP uptake.

Effect of Ligand Distribution on the Cytoskeleton-Associated Endocytosis
We consider the internalization of spherical NPs under the influence of cytoskeleton deformation and determine how the cytoskeleton stiffness may affect the optimal ligand distribution associated with the smallest wrapping time. Figure 4 displays the relation between the wrapping time and the degree of nonuniformity of ligand distribution, , for spherical NPs with a radius of 50 nm for entry into cells. It can be seen that the cytoskeleton deformation can significantly affect the uptake. In the membrane-mediated endocytosis, the optimal distribution of the ligands coated on the spherical NP surface tends to be uniform. With an increase of Young's modulus of the cytoskeleton, the optimal ligand distribution becomes dense at the NP ends with a large curvature, in which the large deformation energy of the cytoskeleton plays a key role in resistance to the cellular uptake. When NPs are uptaken by stiff cells, since a large amount of energy for elastic deformation needs to be overcome, long wrapping times will be needed. Similar results for the oblate ellipsoid and elongated NPs can be seen in Figure 5. Regardless, both Figures 4 and 5 show that the range of the ligand distribution constant, , for effective NP uptake decreases as the Young modulus of the cytoskeleton increases. These findings imply that the internalization of spherical NPs into a soft cytoskeleton can not only shorten the wrapping time but also broaden the range of ligand distribution for effective NP uptake.

Effect of Ligand Distribution on the Cellular Uptake of NPs under Different Initial Recepto Densities
The relation between the wrapping time of NPs and the degree of non-uniformity ligand distribution, , for different initial receptor densities is displayed in Figure  which further identifies how the ligand distribution may influence the wrapping time u der the different initial receptor densities. The semi-axis and engulfing pattern are set = 50 nm and membrane-mediated endocytosis, respectively. It can be seen that t wrapping time decreases as the initial receptor density increases. Therefore, the role receptor diffusion becomes negligible for large initial receptor densities. Additional there exists a broad range of ligand distribution patterns determined by the coefficient corresponding to an almost identical wrapping time. The ranges of coefficient c are almo independent of the initial density of the receptors. However, the corresponding wrappi time decreases with the increase of the initial density of the receptors.

Effect of Ligand Distribution on the Cellular Uptake of NPs under Different Initial Receptor Densities
The relation between the wrapping time of NPs and the degree of non-uniformity of ligand distribution, c, for different initial receptor densities is displayed in Figure 6, which further identifies how the ligand distribution may influence the wrapping time under the different initial receptor densities. The semi-axis and engulfing pattern are set as R a = 50 nm and membrane-mediated endocytosis, respectively. It can be seen that the wrapping time decreases as the initial receptor density increases. Therefore, the role of receptor diffusion becomes negligible for large initial receptor densities. Additionally, there exists a broad range of ligand distribution patterns determined by the coefficient c, corresponding to an almost identical wrapping time. The ranges of coefficient c are almost independent of the initial density of the receptors. However, the corresponding wrapping time decreases with the increase of the initial density of the receptors.

Effect of Ligand Distribution on the Cellular Uptake of NPs under Different Initial Receptor Densities
The relation between the wrapping time of NPs and the degree of non-uniformity of ligand distribution, , for different initial receptor densities is displayed in Figure 6, which further identifies how the ligand distribution may influence the wrapping time under the different initial receptor densities. The semi-axis and engulfing pattern are set as = 50 nm and membrane-mediated endocytosis, respectively. It can be seen that the wrapping time decreases as the initial receptor density increases. Therefore, the role of receptor diffusion becomes negligible for large initial receptor densities. Additionally, there exists a broad range of ligand distribution patterns determined by the coefficient c, corresponding to an almost identical wrapping time. The ranges of coefficient c are almost independent of the initial density of the receptors. However, the corresponding wrapping time decreases with the increase of the initial density of the receptors.

Effect of Ligand Distribution on the Cellular Uptake of NP with Different Sizes
Size-dependent NP endocytosis has been demonstrated in previous studies [4][5][6]. However, how the ligand distribution and NP size may couple together to influence the cellular uptake of NPs remains unclear. We consider membrane-mediated endocytosis to address this issue. Figures 7 and 8 plot the wrapping time of NPs with different aspect ratios, such as λ = 1, λ = 5/3, and λ = 5/9, as a function of the degree of non-uniformity of ligand distribution c, for NPs with different sizes. The optimal ligand distributions for different NP shapes, such as c = 0 (uniform) for a spherical NP, c < 0 for an oblate NP, and c > 0 for an elongated NP, are determined. The results also show that the shortest wrapping time associated with the optimal ligand distribution increases with an increase in the NP's size. For large NPs, a long wrapping time is required to recruit numerous mobile receptors to bind with the ligands on the NP's surface.

Effect of Ligand Distribution on the Cellular Uptake of NP with Different Sizes
Size-dependent NP endocytosis has been demonstrated in previous studies [4][5][6]. However, how the ligand distribution and NP size may couple together to influence the cellular uptake of NPs remains unclear. We consider membrane-mediated endocytosis to address this issue. Figures 7 and 8 plot the wrapping time of NPs with different aspect ratios, such as = 1, = 5/3, and = 5/9, as a function of the degree of non-uniformity of ligand distribution , for NPs with different sizes. The optimal ligand distributions for different NP shapes, such as = 0 (uniform) for a spherical NP, < 0 for an oblate NP, and > 0 for an elongated NP, are determined. The results also show that the shortest wrapping time associated with the optimal ligand distribution increases with an increase in the NP's size. For large NPs, a long wrapping time is required to recruit numerous mobile receptors to bind with the ligands on the NP's surface.  For the wrapping of large NPs, there exists a very large range of c corresponding to a wrapping time almost identical to the optimal one. Under this situation, resistance to the cellular uptake from membrane deformation becomes negligible when compared to the energy associated with the receptor distribution, and a long wrapping time mainly results from the large area of NPs that should be wrapped.

Effect of Ligand Distribution on the Cellular Uptake of NP with Different Sizes
Size-dependent NP endocytosis has been demonstrated in previous studies [4 However, how the ligand distribution and NP size may couple together to influence cellular uptake of NPs remains unclear. We consider membrane-mediated endocytosi address this issue. Figures 7 and 8 plot the wrapping time of NPs with different asp ratios, such as = 1, = 5/3, and = 5/9, as a function of the degree of non-uniform of ligand distribution , for NPs with different sizes. The optimal ligand distributions different NP shapes, such as = 0 (uniform) for a spherical NP, < 0 for an oblate and > 0 for an elongated NP, are determined. The results also show that the shor wrapping time associated with the optimal ligand distribution increases with an incre in the NP's size. For large NPs, a long wrapping time is required to recruit numer mobile receptors to bind with the ligands on the NP's surface.  For the wrapping of large NPs, there exists a very large range of c corresponding a wrapping time almost identical to the optimal one. Under this situation, resistanc the cellular uptake from membrane deformation becomes negligible when compared the energy associated with the receptor distribution, and a long wrapping time mai results from the large area of NPs that should be wrapped. For the wrapping of large NPs, there exists a very large range of c corresponding to a wrapping time almost identical to the optimal one. Under this situation, resistance to the cellular uptake from membrane deformation becomes negligible when compared to the energy associated with the receptor distribution, and a long wrapping time mainly results from the large area of NPs that should be wrapped.

Discussion
It has been revealed that the optimal ligand distributions of NPs depend on the NP's shape and cytoskeleton stiffness in cytoskeleton-associated endocytosis. Such optimal ligand distributions become non-uniform when cytoskeleton deformation is taken into account, which is different from that of membrane-mediated endocytosis on the spherical or cylindrical NPs [35]. The non-uniform optimal ligand distribution results from the competition between the thermodynamic driving force and the kinetics of receptor diffusion. For ligand densities lower than their optimal values at local wrapping edges, insufficient binding energy prolongs the wrapping time. For ligand densities higher than their optimal values at the local wrapping edges, receptor recruitment can also increase the wrapping time, as the high energy dissipation is due to the large change in the configuration entropy of the receptors.
In the cellular uptake of spherical or cylindrical NPs via membrane-mediated endocytosis, the uniform distribution of ligands becomes the optimal one corresponding to the high uptake efficiency. By contrast, during the cellular uptake of ellipsoidal NPs via cytoskeleton-associated endocytosis, the mean curvature of the membrane at each wrapping edge is no longer a constant, but is determined by NP shape. In the power balance in Equation (18), the fourth term at the left hand is the energy contribution of cytoskeleton deformation, which is dependent on the engulfment depth. Therefore, the optimal ligand distribution turns to non-uniform and becomes strongly dependent on NP shape and cytoskeleton stiffness.
Bio-inspired methods from viruses are suitable for designing drug delivery systems. Thus, a biophysical understanding of NP-cell interactions is urgently needed. For spherical NPs, fast wrapping occurs in a large range of different ligand distribution patterns around a uniform distribution; this finding provides a physical insight into robust virus infection rather than the point of view of gene expression [57,58]. The optimal size (tens of nanometers) [4,6] and shape (sphere) [8] have been revealed from a physical optimization standpoint. In this study, we confirm that ligand distribution is another significant factor in determining the receptor-diffusion-mediated NP uptake into cells. The almost uniform ligand distribution of spherical viruses is possibly controlled by physical evolution and guarantees viral infectivity via receptor-mediated endocytosis.
We also examined the critical state as follows: By solving Equation (16), where ProductLog(.) is the Lambert-W function. In the gradient-driven diffusion process of mobile receptors along the cell membrane, effective wrapping requires ξ + < ξ 0 , When the receptor density in front of the contact edge ξ + is equal to the initial receptor density ξ 0 , the wrapping process begins to be terminated. Accordingly, the solution to Equation (18) in such a critical case, ξ + = ξ 0 , yields the critical ligand density in front of the contact edge, If the ligand density at the wrapping edge is larger than the critical value, then the wrapping process cannot occur due to the large unfavorable entropic contribution of the energy from the receptor configuration change.
Our study presents certain limitations. Although we provide a detailed description on the dependence of the optimal ligand distribution on the NP shape and cytoskeleton stiffness, we have treated the NPs as a rigid body. We did not consider in detail whether NP deformation significantly influences the optimal ligand distribution. We also disregard membrane tension [18], the kinetic reaction between receptor and ligand molecules [32,[59][60][61], as well as the NP concentration [62]. We apply these assumptions in this study so that we can easily focus on the effect of ligand distribution on the cellular uptake of ellipsoidal NPs via cytoskeleton-associated endocytosis.

Conclusions
We show the systematically distinct effects of ligand distribution on the dynamics of cytoskeleton-associated endocytosis of ellipsoidal NPs with different sizes under different initial receptor densities. NPs with the same number of ligands but different distributions can have a wide range of different dynamic processes of cellular uptake. We find that there exist optimal ligand distributions corresponding to the minimum wrapping times for NPs with different shapes. Unlike the findings in previous studies on the cell entry of spherical and cylindrical NPs with different ligand distributions via membrane-mediated endocytosis, such optimal ligand distributions become non-uniform and dependent on the NP shape and cytoskeleton stiffness. For example, the optimal distribution favors that the ligands are densely distributed in the region with large local curvature or large cytoskeleton deformation. These findings supply an insight into the physical understanding of NP-cell interactions and may provide guidelines for targeted drug delivery.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
In order to better describe the mean curvature at any point p corresponding to engulfment depth h, a Cartesian coordinate system is established at the center of an ellipsoid. As shown in Figure A1, the mean curvature at point p x p , y p of this ellipsoid can be expressed as, with two Gaussian principal curvatures c 1 and c 2 . Furthermore, c 1 is the Gaussian principal curvature of the blue intersecting line at point p shown in Figure A1, and c 2 is the Gaussian principal curvature of the red intersecting line which is orthogonal to the blue one. One of the principal curvatures along the symmetry plane can be determined as, c 1 = x p (θ)y p (θ) − x p (θ)y p (θ) x p (θ) 2 + y p (θ) 2 3/2 , Besides, the parameter equation of the symmetry plane and geometrical relation can be written as, where θ is the phase angle. Substitution of Equation (A3) into Equation (A2) yields, in which, = ⁄ . The other principal curvature at point is the reciprocal of the principal curvature radius ( Figure A1). As well, the line pn in Figure A1 can be given as, in which, λ = R a /R b . The other principal curvature c 2 at point p is the reciprocal of the principal curvature radius l pn ( Figure A1). As well, the line pn in Figure A1 can be given as, where angle α depends on the position of point p, and it can be obtained by, dy dx x=x p , y=y p = tan π 2 + α = − cot α.
According to the geometrical relation, The combination of Equations (A6) and (A7) and sin 2 α + cos 2 α = 1 gives, sin 2 α = 1 Hence, the other principal curvature c 2 can be deduced from Equations (A5) and (A8), which is, The substitution of Equations (A9) and (A4) into Equation (A1) yields the mean curvature as a function of engulfment depth h, About the curvature of an ellipsoid, there are some existing results [63] based on different analysis.

Appendix B
As shown in Figure A1, the infinitesimal wrapping area of the ellipsoidal NP by cell membrane can be expressed as, with infinitesimal arc length dx 2 + dy 2 . Differentiation of Equation (A3) with respect to engulfment depth h yields, The infinitesimal wrapping area can be rewritten as, Insertion of Equations (A12) and (A3) into Equation (A13) and subsequent integration yields,