## Abstract

Molecular dynamic simulations and experiments have recently demonstrated how cylindrical nanoparticles (CNPs) with large aspect ratios penetrate animal cells and inevitably deform cytoskeletons. Thus, a coupled elasticity–diffusion model was adopted to elucidate this interesting biological phenomenon by considering the effects of elastic deformations of cytoskeleton and membrane, ligand–receptor binding and receptor diffusion. The mechanism by which the binding energy drives the CNPs with different orientations to enter host cells was explored. This mechanism involved overcoming the resistance caused by cytoskeleton and membrane deformations and the change in configurational entropy of the ligand–receptor bonds and free receptors. Results showed that deformation of the cytoskeleton significantly influenced the engulfing process by effectively slowing down and even hindering the entry of the CNPs. Additionally, the engulfing depth was determined quantitatively. CNPs preferred or tended to vertically attack target cells until they were stuck in the cytoskeleton as implied by the speed of vertically oriented CNPs that showed much faster initial engulfing speeds than horizontally oriented CNPs. These results elucidated the most recent molecular dynamics simulations and experimental observations on the cellular uptake of carbon nanotubes and phagocytosis of filamentous *Escherichia coli* bacteria. The most efficient engulfment showed the stiffness-dependent optimal radius of the CNPs. Cytoskeleton stiffness exhibited more significant influence on the optimal sizes of the vertical uptake than the horizontal uptake.

## 1. Introduction

Nanoparticles (NPs) usually enter cells via endocytosis driven by the binding energy between diffusive receptors on the cell membranes and ligands on the surface of the NPs. This remarkable capability has resulted in the proposition of NPs as potential candidates for site-specific drug-delivery systems [1]. However, incomplete phagocytosis of long, rigid biopersistent asbestos fibres and carbon nanotubes (CNTs) with large aspect ratios has been hypothesized to cause the release of inflammatory mediators and even cell death [2,3]. Therefore, understanding how cells interact with cylindrical nanoparticles (CNPs) is greatly significant to advance our understanding of fundamental biological and cellular immunological recognition as well as improve various practical applications of drug delivery, pharmacology and nano-toxicology [4,5].

Numerous experiments and simulations have recently explored the mechanics and reasons for the penetration of one-dimensional nano-materials and micro-materials into living cells, including macrophages. Shi *et al*. [6] experimentally and theoretically demonstrated that CNPs, such as CNTs, initially enter the cells through the tip. This entry is assumed to be caused by tip recognition following rotation initiated by the asymmetric elastic strain at the tube–bilayer interface. To determine whether such shape- and orientation-dependent mechanisms are also manifested by the phagocytosis of bacterial filaments, Möller *et al*. [7] experimentally investigated the effects of shape and micro-environments on the phagocytosis of filamentous *Escherichia coli* bacteria by macrophages. They found that complete uptake occurs only if one of the terminal bacteria filament poles enters the cell first. Abdolahad *et al*. [8] used an array of vertically aligned multi-walled CNTs (radius = 32.5 nm) to distinguish healthy and cancerous cells by evaluating the entrapment of different cells. They found that the entrapment fraction of cancer cells with higher metastatic grades was significantly greater than those with lower metastatic grades, because the former are less stiff than the latter. These previous studies have shown that cellular uptake of one-dimensional nano-materials or micro-materials with high-aspect ratios can be dependent on orientation and cytoskeleton deformation.

The entry of long CNPs into cells inevitably deforms the cytoskeleton. The long CNPs can penetrate very deeply into cells, therefore, their interaction with the cytoskeleton is unavoidable. Obataya *et al*. [9] used an atomic force microscope with a long nanoneedle (radius = 100−150 nm) to indent living cells. They demonstrated that the loading force curve for an indentation depth of up to 2 µm is still consistent with the classic Hertz model in contact mechanics. Similarly, Beard *et al*. [10] 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 many viruses including the herpes simplex virus type 1 [11]. Likewise, they confirmed that the Hertz model can fit the loading force curve very well. The Hertz model is derived based on the contact problem of two elastic bodies. Thus, these results clearly prove that during indentation the living cell behaves as a deformable elastic solid instead of a membrane, implying that deformation of the cytoskeleton plays a dominant role in resisting CNP intrusion.

The cellular uptake of CNP processes via endocytosis for long, short or spherical NPs, regardless of size, accompany cytoskeleton remodelling, considering that many endocytic proteins can be potentially linked to the actin cytoskeleton, either directly or indirectly [12–15]. This phenomenon eventually leads to a close functional connection between the actin cytoskeleton and the internalization step of endocytosis. Therefore, many virologists and physiologists [16–18] have pointed out that the physical barrier imposed by the cortical actin meshwork on the endocytosis process should be overcome following the stimulation of actin cytoskeleton remodelling. This phenomenon is observed in HIV-1 viral particles such that the ERM protein family supplies functional linkage between integral membrane proteins and the cytoskeleton in mammalian cells to regulate membrane protein dynamics and cytoskeleton rearrangement [19,20]. After an HIV-1 Env protein binds to a host cell, a cytoskeleton-dependent clustering of infection receptors (CD4, CXCR4 and CCR5) occurs. These receptors are essential for efficient membrane fusion and subsequent entry of HIV-1 into the target cells [20]. A number of theoretical models have been proposed to understand the underlying biophysical mechanism for the cellular uptake of NPs based on these recognitions. Regarding diffusion kinetics, although neglecting the effect of cytoskeleton deformation, Gao *et al*. [21] and Shi *et al*. [22] showed that in the endocytosis of cylindrical and spherical NPs, optimal sizes exist because of the power balance between the kinetics of receptor diffusion and thermodynamic interaction of the NPs and the absorbing membrane. In addition, Sun & Wirtz [23] and Li *et al*. [24] developed an elastic model to investigate the resistance during the process of viral particle entry into host cells by neglecting receptor diffusion but considering a balance between energetic forces. They concluded that the resistance to engulfment is dominated by the membrane or cytoskeleton deformation depending on the size of the virus, engulfing stages and Young's modulus of the cell. Most recently, Yi *et al*. [25] investigated the adhesive wrapping of a soft elastic lipid vesicle by a lipid membrane analogous to the cellular uptake of elastic particles using variational methods and establishing free energy functionals. They identified five distinct wrapping phases. Zhang *et al*. [26] and Yuan *et al*. [27] analysed the equilibrium interaction between a group of NPs and the membrane. They also showed that although non-equilibrium processes of each individual NP are not considered, the predicted optimal particle size for maximal cellular uptake is close to that for the shortest wrapping time as provided by Gao *et al*. [21].

In spite of the aforementioned development, the understanding of cellular uptake of CNPs is still limited. Existing models have not investigated how cytoskeleton deformation may affect the equilibrium state and dynamic process of the cellular uptake behaviour. In this study, an elasticity–diffusion model is proposed to examine the engulfment of CNPs, considering the coupled effects on cytoskeleton and membrane deformations, receptor diffusion, ligand–receptor binding and CNP orientations. The majority of previous theories on the cellular uptake of NPs has focused on NP–membrane interaction. By contrast, the present model consists of a membrane-covered elastic solid with mobile molecular binders diffused along the membrane under the ligand binding action of NPs. The effects of the sizes and orientations of NPs, as well as stiffness and receptor densities of host cells, on the dynamic engulfing process of CNPs was investigated according to the proposed model.

## 2. Theoretical model

An elastic half space covered by a cell membrane embedded with diffusive mobile receptors was used to illustrate the proposed elasticity–diffusion model (figure 1*a*). The membrane engulfed horizontally and vertically oriented elastic CNPs coated with compatible ligands with a depth, *h*, below the initially undeformed position of the membrane. Each CNP with length, *L*, and radius, *R*, were capped by two hemispheres with radius, *R*, at its two ends. The ligands with a fixed constant density, *ξ*_{L}, were assumed to be immobile and uniformly distributed on the particle surface. By contrast, the receptors with an initial uniform density of *ξ*_{0} were mobile and underwent diffusive motion within the plane of the cell membrane. When the CNP with either orientation came into contact with the cell, the binding energy of the receptor–ligand complex drove the engulfment of the CNP to overcome the resistance caused by the change in the configuration entropy in the receptors by diffusion and the deformations of the membrane and cytoskeleton. Moreover, the free energy of interaction was lowered further to maintain the engulfing process. Thus, the receptors diffused to the contact site and bound with the ligands on the surface of the particle.

Figure 1*b* schematically shows the coordinate system at the centre of the CNP–cell contact zone. The half size of the contact area was *s* = *a*(*t*). At the contact zone, the receptor density was assumed to reach the ligand density, i.e. *ξ*(*s*, *t*) = *ξ*_{L}, *s* < *a*(*t*). During the engulfing process, the edge of the CNP–cell contact zone, *s* = *a*(*t*), and the distribution of membrane receptors outside the contact zone, *ξ*(*s*, *t*), continued to vary. The difference was determined by solving the dynamic problem influenced by ligand–receptor binding and elastic deformations of the membrane and the cytoskeleton.

### 2.1. Horizontally oriented cylindrical nanoparticles

The engulfment of the horizontally oriented CNP with high-aspect ratio, i.e. was considered first. In this case, receptor diffusion along the axial direction was neglected and mainly maintained along the direction perpendicular to the axial direction. The capped CNP with total surface area, 2*πR*(*L* + 2*R*), was also approximately replaced by a non-capped cylinder with the same diameter, surface area and length, 2*R* + *L*. In this situation, the CNP, which may or may not be capped, only had a negligible effect on receptor diffusion and deformation energy of contact because of its high-aspect ratio.

During the engulfing process, a global number of conservation conditions for the receptors could be expressed in terms of the receptor density, *ξ*(*s*, *t*), as [21,28]
2.1where *a*(*t*) is the half-width of the contact region. Outside the contact region, i.e. *s* ≥ *a*(*t*), the receptor diffusion flux is given by
2.2where *D* is the diffusion coefficient.

Within the contact region, *s* < *a*(*t*), *ξ*(*s*, *t*) = *ξ*_{L} and *j*(*s*, *t*) = 0. The continuity equation
2.3can be substituted into equation (2.1) to yield [21,28]
2.4where *ξ*_{+} ≡ *ξ*(*a*^{+}, *t*) and *j*_{+} ≡ *j*(*a*^{+}, *t*) denote receptor density and flux in front of the contact edge, respectively. Boundary conditions at *s* → ∞ were assumed to be *ξ*(*s*, *t*) → *ξ*_{0} and *j*(*s*, *t*) → 0.

Upon initial CNP–cell contact, the formation of ligand–receptor complexes drove the engulfing process to overcome the resistance caused by membrane and cytoskeleton deformation and the change in the configuration entropy of receptors. During this process, the binding energy of ligand–receptor complexes in the contact region could be expressed as , where *e*_{RL}*k*_{B}*T* is the binding energy of each single ligand–receptor bond. The deformation energy of the membrane can be written as , where *H*_{s} = 1/*R* and *H*_{c} = 1/2*R* are the mean curvatures of the spherical and cylindrical surfaces, respectively. *κk*_{B}*T* and *γk*_{B}*T* are the bending modulus and surface tension of the membrane, respectively. The free energy caused by the configurational change of all receptors on the membrane can be estimated as [21,28] . Here *k*_{B}*T* ln *ξ*_{L}/*ξ*_{0} and *k*_{B}*T* ln *ξ*/*ξ*_{0} are the energy per receptor associated with the loss of configuration entropy of the bonds and free receptors, respectively.

Cytoskeleton deformation also significantly contributes to the energy of the system [18]. Recently, Liu *et al*. [30] further addressed this issue by performing living cell indentation with a long nanoneedle (radius 90 nm). As shown in figure 2, when we compare their experimental loading force curve to the theoretical prediction based on linear elasticity, good agreement is found until cell membrane penetration. These experimental studies clearly prove that during deep indentation by nanoneedles living cells can somehow behave as linear elastic solids, implying that the Hertzian contact model based on linear elasticity might be an appropriate framework to describe the interaction between the CNPs and the cells, and under this situation, the complicated finite deformation formulation may not be necessary.

Therefore, by treating the cell as an elastic solid and within the limit where the cell is much larger than the NP, the energy of cytoskeleton deformation during the interaction between the CNP and the cell, based on the contact theory [29], can be estimated as
where 1/*E** = (1 − *μ*_{c}^{2}/*E*_{c}) + (1 − *μ*_{n}^{2}/*E*_{n}) is the combined elastic modulus. *μ*_{c} and *E*_{c} are Poisson's ratio and Young's modulus of the cell, respectively. The corresponding values for the NP are *μ*_{n} and *E*_{n}. Variable *p* is an effective contact force related to the engulfing depth by [29]
2.5where the contact depth, *h*, as shown in figure 1*a*, is geometrically related to the half-width of the contact region *a*(*t*) as
2.6Based on the above analysis, the free energy functional for the cellular uptake of the CNP can be integrated as
2.7Differentiation of equation (2.7) with respect to time yields
2.8A power balance equation for the horizontally oriented CNP can be obtained by equating the decrease rate of the free energy to the energy dissipated during receptor diffusion.
2.9which was obtained by inserting equations (2.5) and (2.6) into equation (2.8).

### 2.2. Vertically oriented cylindrical nanoparticles

In the engulfment of vertically oriented CNPs, the associated receptor diffusion was considered as an axisymmetric diffusion problem with the following governing equation outside the contact region
2.10The receptor conservation condition in this case was similar to that in equation (2.4). The time-dependent contact area can be expressed as *A*(*t*) = 2*πRh*(*t*) = *πa*^{2}(*t*). The binding energy caused by the formation of ligand–receptor bonds can be described as . The deformation energy of the membrane as a function of *a*(*t*) can be expressed in a sectional type as

2.11

The deformation energy of the cytoskeleton was estimated based on the contact theory [29,31] using . In addition, the energy contribution caused by the loss of configurational entropy of receptors can be given by .

Eventually, the total free energy functional was obtained as 2.12Differentiating equation (2.12) with respect to time leads to 2.13where the mean curvature is identified as 2.14Similarly, the power balance on the boundary of contact region can be given by 2.15

## 3. Results and discussion

### 3.1. Procedure for numerical solutions of governing equations

Obtaining the analytical solutions of equations (2.3) and (2.10) is a challenging task. In this section, these equations were numerically solved using the finite difference method similar to that described in previous studies [21,28]:

(1) At time step 0, i.e.

*t*= 0, the initial and boundary conditions on receptor distribution and diffusion flux were set to*ξ*(*s*, 0) =*ξ*_{0}and*ξ*(∞,*t*) =*ξ*_{0},*j*(∞, 0) = 0, respectively.(2) At time step 1, i.e.

*t*= Δ*t*, where Δ*t*was small enough so that the effect of cytoskeleton deformation was negligible. The power balance relations in equations (2.9) and (2.15) were respectively reduced to 3.1and 3.2In this situation, solutions of equations (2.3) and (2.10) at*t*= Δ*t*could be obtained analytically as [21,28] 3.3and 3.4where ,*A*_{h}and*A*_{v}are unknown constants of integration. Inserting equations (3.3 and (3.4), as well as [21,28], into equation (2.5) yields 3.5and 3.6where*α*is the speed factor [21,28]. Then, substituting equations (3.3)–(3.6) into the power balance relations of equations (3.1) and (3.2) yields 3.7and 3.8where and . The speed factor,*α*, was obtained by solving equations (3.7) and (3.8). Once*α*was determined, then the receptor density,*ξ*(*s*, Δ*t*), at this time step was obtained from equations (3.3)–(3.6). Furthermore, the diffusion flux,*j*(*s*,Δ*t*), and subsequently the engulfing speed, , was derived from equations (2.2) and (2.4), respectively.(3) At time step 2, i.e.

*t*= 2Δ*t*, the engulfing boundary was expressed as . Inserting*a*(2Δ*t*) into the power balance relations (equations (2.9) and (2.15)) provided another boundary condition of receptor density*ξ*_{+}(2Δ*t*) at*s*=*a*(2Δ*t*). Using this boundary condition and that specified in step 1, the receptor density,*ξ*(*s*, 2Δ*t*), was obtained by solving the diffusion equations (2.3) and (2.10) using the finite difference method. The diffusion flux,*j*(*s*,2 Δ*t*), and subsequently the engulfing speed, , were derived from equations (2.2) and (2.4).(4) Step 3 was repeated until the prescribed time step

*N*, i.e.*t*=*N*Δ*t*, was reached. Finally,*ξ*(*s*,*t*) was obtained, where*t*≤*N*Δ_{t}.

### 3.2. Faster engulfment of vertically oriented cylindrical nanoparticles

The bending modulus, *κ*, of a typical cell membrane range from 10 to 20 *k*_{B}*T* [32], and the surface tension is approximately 0.005 *k*_{B}*T*/nm^{2} [23]. Young's modulus of the cell is mostly on the order of 10 kPa or less [33]. In this study, Poisson's ratio of the cell was considered to be 0.5 [23]. The binding energy of a single closed ligand–receptor bond *e*_{RL} is approximately 15 *k*_{B}*T* at *T* = 300 K [34]. The diffusion constant of a receptor on the membrane is approximately 10^{4} nm^{2} s^{−1} [34,35]. A previous study [21] has shown that the ratio should range from 0.01 to 0.1. In general, Young's modulus of viral particles and most artificial NPs are much larger than that of a cell. Young's modulus of CNTs is on the order of 1 TPa [36]. Hence, the CNPs were assumed to behave like rigid bodies during cellular uptake. Ligands distributed on the surface of each CNP had typical density values of 5 × 10^{3} µm^{−2} [21].

The relevant parameters adopted in this study on cellular uptake of NPs are summarized in table 1.

Figure 3 plots the numerically determined normalized engulfment depth, *h*/*R,* as a function of engulfing time for CNPs with radius, *R* = 50 nm, and length, *L* = 500 nm, under horizontal (figure 3*a*) and vertical (figure 3*b*) orientations. Figure 3 shows that the vertically oriented CNP enters the cell with an initial uptake speed close to 0.0133 s^{−1}, which is much faster than that of the horizontally oriented CNP at 2.9 × 10^{−4} s^{−1}.

In the authors' previous study [24], based on a continuum model, for the enveloped virus to enter the host cells, the binding energy of the receptor–ligand complexes drove the engulfment of NPs to overcome the resistance alternatively dominated by membrane deformation and cytoskeleton deformation at different engulfing stages. In this study, the coupling effect between receptor diffusion and elastic deformation was considered. The resistance to CNP engulfment was still found to be initially dominated by membrane deformation, and only later by cytoskeleton deformation. Figure 3 clearly demonstrates that cell stiffness almost did not affect the initial uptake speed, which started to significantly influence the speed only after a certain engulfing depth. This phenomenon may imply that for the engulfment of long CNPs or large NPs, the cytoskeleton process can be crucial if the penetrations are sufficiently deep. By contrast, in tiny spherical NPs, the membrane process rather than cytoskeleton remodelling plays the dominant role.

Figure 3 also shows that a stiffer cell or a stiffer part of a cell can more easily resist the engulfing process after a certain engulfing depth. As has been shown in figure 3*b*, the entry of vertical CNPs can even be completely stopped after a certain engulfing time. This phenomenon mainly depends on the competition between different energy contributions. At higher cell stiffness, the elastic deformation energy of the cytoskeleton will balance more binding energy of ligand–receptor bonds. This behaviour decreases the effective energetic driving force, causing receptor diffusion for engulfment to slow down.

### 3.3. Faster cylindrical nanoparticle engulfment by softer cells or softer parts of a cell

To differentiate cell stiffness, the normalized receptor density on the boundary of the contact region, *ξ*_{+}/*ξ*_{0}, was numerically determined and plotted in figure 4 as a function of engulfment depth. This normalized receptor density reflects the maximum inward diffusion flux at time *t* (figure 5). Figure 4 also shows that for a small engulfing depth, the normalized receptor density is small and is almost not influenced by cell stiffness. This observation implies that the engulfing process at this stage possesses large inward receptor diffusion flux, fast engulfing speed and is not influenced by cytoskeleton deformation. However, for a relatively large engulfing depth, the normalized receptor density at the contact edge or the maximum inward diffusion flux, starts to be strongly influenced by cell stiffness (figure 4). At this stage of the engulfing process, the larger cell stiffness corresponds to the larger normalized receptor density at the contact edge, or the smaller maximum inward diffusion flux, implying a slower engulfing process. When Young's modulus of the cell becomes large enough, the diffusion flux eventually becomes zero, corresponding to *ξ*_{+}/*ξ*_{0} = 1 (figure 4). This condition also indicates that when a virus comes into contact with a stiff part of the cell, endocytosis may be stopped since receptors can no longer be recruited by diffusion. Hence, viruses may seek soft parts of a cell to attack, or viruses that attack soft parts of a cell could more successfully infect the cell.

For the cellular uptake of horizontally oriented CNPs, critical engulfment depth, *h*_{cr}, corresponding to the maximum receptor density on the engulfing boundary exists (figure 4). The deformation energy of the cytoskeleton increases along with the engulfing depth. When the engulfment depth is smaller than the critical value, the increased energy leads to a decreased energetic driving force. When the engulfment depth becomes larger than the critical value, the increasing rate of binding energy of ligand–receptor complex becomes larger than that of the deformation energy of cytoskeleton. This phenomenon results in an increased driving force to hasten the engulfing process as shown in the inset of figure 3*a*. Such critical *h*_{cr} can be numerically calculated by minimizing the seventh term, *f* = (*p*(*h*) sin *a*(*h*)/*R*)/(2*k*_{B}*T*(*L* + 2*R*)), in equation (2.10) with respect to *h*, where d*f*/d*h* = 0. The critical depth, *h*_{cr} = 76.8 nm, for the CNP with radius, *R* = 50 nm, only depends on the shape of the CNP.

For the vertical uptake of CNPs, abrupt changes were observed on the normalized receptor edge density when the engulfing depth reached *R*, because of the sharp change in the surface curvature of CNPs (figure 3*b*).

### 3.4. Effects on cellular uptake by initial receptor density

Figure 6 shows the engulfing depth as a function of time for various initial receptor densities. For both horizontal and vertical uptakes of CNPs, the receptor density significantly influenced the overall uptaking process. A larger receptor density provided a faster uptaking process. However, when the receptor density theoretically approaches infinity, an unrealistic ultrafast uptaking process occurs. This result reveals that the well-accepted receptor–diffusion-mediated mechanics model of endocytosis might not be applicable in the case when receptor density is high enough such that receptor recruiting through diffusion is energetically unfavourable during the engulfing process. Essentially, the process of cellular uptake of CNPs is rather complicated. This process cannot be fully understood based on theoretical models employing only the membrane process. For example, as shown in figure 7, at least three different regimes were hypothesized for the uptaking process of CNPs in terms of different initial receptor densities on the membrane.

#### 3.4.1. Regime 1

If the normalized initial receptor density, *ξ*_{0}/*ξ*_{L}, is much smaller than 1, the uptaking process of CNPs will be dominated by receptor diffusion. This juncture determines whether adhesion is sufficient to overcome the configurational entropy change of receptors and membrane bending at the initial stage.

#### 3.4.2. Regime 2

If *ξ*_{0}/*ξ*_{L} becomes larger but still smaller than 1, the uptaking process will be controlled by both receptor diffusion and the deformations of the membrane and cytoskeleton. In the case of intermediately large *ξ*_{0}/*ξ*_{L}, the initial adhesion can easily overcome the configurational entropy change of the receptors during diffusion to challenge membrane bending at the initial stage and later cytoskeleton deformation. This phenomenon also indicates that the uptaking process involves the full interplay between the receptor diffusion, membrane deformations and cytoskeleton deformations. The theoretical model in this regime should consider these three factors. Additionally, a pure membrane model or a coupled membrane–diffusion model may be valid only at the initial stage of the uptake or only for the uptake of small NPs with negligible cytoskeleton deformations. For the complete uptake of sufficiently long CNPs, a fully coupled model accounting for receptor diffusion is necessary for both membrane and cytoskeleton deformations.

#### 3.4.3. Regime 3

If *ξ*_{0}/*ξ*_{L} is larger than 1, receptor recruitment through diffusion is no longer energetically favourable. Adhesion will easily overcome membrane deformation and cytoskeleton deformation. In this case, the dynamic uptaking process will be controlled by the creeping of the cytoskeleton (remodelling). The membrane model is no longer valid in this regime. A detailed analysis of this regime has recently been provided by Wang *et al*. [37].

### 3.5. Effect of cytoskeleton stiffness on the optimal size and maximum engulfing depth of cylindrical nanoparticles

Similar to the previous studies [21], optimal sizes for both horizontal and vertical uptake of NPs exist. Small CNPs have large surface curvatures, so engulfing them implies the need to overcome high membrane bending energies. By contrast, engulfing large NPs implies the need to overcome large deformation energies of the cytoskeleton.

Figure 8 shows the inverse of the engulfing time at engulfing depth, *h* = *R*, as a function of Young's modulus of cytoskeleton and the radius of the CNPs. For the horizontal uptake of CNPs, optimal sizes corresponding to the fastest engulfing process exist (figure 8*a*). These optimal sizes are almost independent of Young's modulus. However, for the vertical uptake of NPs, Young's modulus of the cytoskeleton obviously influences the optimal sizes (figure 8*b*). A relatively fast increase in the optimal sizes along with Young's modulus can be seen in figure 8*b*.

Figure 9 shows the engulfing time, *h* = *R*, as a function of the radius of the CNPs under different Young's modulus of the cytoskeleton. These optimal sizes correspond to the smallest engulfing times and are influenced by cell stiffness. Interestingly, these optimal sizes are also approximately 20 and 30 nm, close to those previously predicted by the receptor-mediated membrane model [21] and experimental observations [38].

Figure 10 displays the maximum engulfing depth as a function of Young's modulus of cells for the vertical engulfment of CNPs with *R* = 50 nm. It can be seen from figure 10 that, for softer parts of cells, the maximum engulfing depth is larger. When Young's modulus is below 2 kPa, this final engulfing depth can even reach 400 nm, eight times *R*.

## 4. Conclusion

A coupled elasticity–diffusion model was established to study the cellular uptake of CNPs with horizontal and vertical orientations while considering the effects of ligand–receptor binding, receptor diffusion, membrane deformation and cytoskeleton deformations. The proposed model showed that the effect of cytoskeleton deformation should be considered to elucidate the uptaking processes of long CNPs or large spherical NPs. Membrane models could be considered a good substitute only when NPs are small enough. A CNP with vertical orientation exhibited a much faster initial uptaking speed than that with horizontal orientation at the initial stage of the uptaking process, where deformation of the membrane was more important than that of the cytoskeleton. Cytoskeleton stiffness significantly influenced the engulfing process only after a certain engulfing depth. Larger cytoskeleton stiffness corresponded to a slower engulfing process. Optimal sizes of CNPs engulfed at different orientations were observed. Cytoskeleton stiffness showed more significant influence on the optimal size for vertical uptake than for horizontal uptake.

Based on the proposed model, the phenomenon identified by Shi *et al*. [6] in which CNPs, such as CNTs, enter the cells through the tip first, was hypothesized to be true only at the initial stage of the uptaking process. In most cases, vertically oriented CNPs probably stuck at a certain length, which may have been engulfed by the soft part of the cytoskeleton. However, this attachment was subsequently stopped by the stiff part of the cytoskeleton.

## Funding statement

This research is supported by grants from the National Natural Science Foundation of China (11472119, 11032006, 11121202), National Key Project of Magneto-Constrained Fusion Energy Development Program (2013GB110002), the Fundamental Research Funds for the Central Universities (lzujbky-2013-1).

- Received September 12, 2014.
- Accepted October 27, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.