# Influence of the mechanical and geometrical parameters on the cellular uptake of nanoparticles: A stochastic approach

Sarah Iaquinta, Shahram Khazaie, Éléna Ishow, Christophe Blanquart, Sylvain Fréour, Frédéric Jacquemin

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Sarah Iaquinta, Shahram Khazaie, Éléna Ishow, Christophe Blanquart, Sylvain Fréour, et al.. Influence of the mechanical and geometrical parameters on the cellular uptake of nanoparticles: A stochastic approach. International Journal for Numerical Methods in Biomedical Engineering, 2022, 38 (6), 10.1002/cnm.3598. hal-03993822

HAL Id: hal-03993822

https://hal.science/hal-03993822

Submitted on 17 Feb 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Influence of the mechanical and geometrical parameters on the cellular uptake of nanoparticles: a stochastic approach

Sarah Iaquinta<fnref n="1" /> | Shahram Khazaie<fnref n="1" /> | Éléna Ishow<fnref n="2" /> | Christophe Blanquart<fnref n="3" /> | Sylvain

Fréour<fnref n="1" /> | Frédéric Jacquemin<fnref n="1" />

Correspondence

*Shahram Khazaie, Laboratoire GEM, IUT
de Saint Nazaire, 58 rue Michel Ange,
F-44600 Saint-Nazaire, France. Email:
shahram.khazaie@univ-nantes.fr

Funding information

i-Site NexT, Région Pays de la Loire and
CNRS (French National Centre for Scientific
Research) through the ENAMEL project,
Grant/Award Number: ANR-16-IDEX-0007

Summary

Nanoparticles (NPs) are used for drug delivery with enhanced selectivity and reduced side-effect toxicity in cancer treatments. Based on the literature, the influence of the NPs mechanical and geometrical properties on their cellular uptake has been studied through experimental investigations. However, due to the difficulty to vary the parameters independently in such a complex system, it remains hard to efficiently conclude on the influence of each one of them on the cellular internalization of a NP. In this context, different mechanical / mathematical models for the cellular uptake of NPs have been developed. In this paper, we numerically investigate the influence of the NP's aspect ratio, the membrane tension and the cell-NP adhesion on the uptake of the NP using the model introduced in<fnref n="1" /> coupled with a numerical stochastic scheme to measure the weight of each one of the aforementioned parameters. The results reveal that the aspect ratio of the particle is the most influential parameter on the wrapping of the particle by the cell membrane. Then the adhesion contributes twice as much as the membrane tension. Our numerical results match the previous experimental observations.

KEYWORDS:

cellular uptake, mechanical properties, shape, sensitivity analysis, surrogate modeling

1 | INTRODUCTION

The development of nanoparticles (NPs) for the delivery of anti-cancer agents is a rapidly expanding field<fnref n="2" />⁻⁸. The effectiveness of this strategy relies on the internalization of NPs by malignant cells to obtain an intracellular release of the active compounds. A better understanding of the mechanisms and the interactions involved is necessary to develop adapted NPs. Investigations have initially and mainly been conducted experimentally⁹⁻¹⁴, while numerical studies have been emerging more recently¹⁵⁻²³, many of them focusing on the influence of the NP's physico-chemical properties, illustrated in Figure 1. The cellular internalization pathway usually preferred by NPs is endocytosis²⁴,²⁵, among which one can distinguish specific from non-specific ones. In the latter, NPs can be internalized by the cells using non-specific interactions, also known as non-receptor-mediated uptake, to be distinguished from internalization involving interactions between ligands and receptors present on the cell surface. Thus, in the case of a non specific uptake, the cellular NP's internalization process initiates spontaneously after that the NP's adheres to the cell membrane⁴,²⁶,²⁷. In the literature, the cellular uptake of NPs is often simulated with methods based on molecular

<footnote n="1">¹Nantes Université, Ecole Centrale Nantes,
CNRS, GeM, UMR 6183, F-44600
Saint-Nazaire, France</footnote>

<footnote n="2">²Nantes Université, CNRS, CEISAM, UMR
6230, F-44000 Nantes, France</footnote>

<footnote n="3">³Nantes Université, Univ Angers, CHU
Nantes, INSERM, CNRS, CRCI2NA,
F-44000 Nantes, France</footnote>
dynamics<sup>4,15,20,28,29</sup> that allow to model molecules and atoms. However, models such as the ones described in<sup>1,30–34</sup> represent the phenomenon at the scale of the NP and model both the cell membrane and the NP with thin homogeneous lines.

In order to take into consideration the real-world manufacturing challenges relative to the control and repeatability of the NP's size, aspect ratio or surface adhesion<sup>13,35,36</sup>, it is necessary to identify the sensitivity of its uptake to those parameters. Indeed, this identification would allow the producers to prioritize some properties over others, taking into account the manufacturing and economic constraints among others.

In this study, we present a simplified model for the cellular uptake of an elliptic NP to investigate the influence of the mechanical properties of the cell membrane and of the aspect ratio of the particle, using a stochastic approach. In order to perform a purely mechanical investigation, we chose to model non specific NP's cellular internalization, focusing at the scale of the NP. The novelty of this work is (i) to clarify the models that already exist in the literature<sup>1,4,37</sup> and to present complementary results from the ones commonly exposed, (ii) to provide a biophysical explanation of them and (iii) to perform a stochastic modeling of the input parameters, to build a metamodel for the output random wrapping degree at equilibrium and subsequently to estimate the sensitivity of the results to each of them. The outline of this paper is the following. Section 2 presents the model along with its corresponding hypotheses, the adopted investigation strategy with a broad description of the stochastic approach and the subsequent sensitivity analysis (Section 3). Section 4 then provides an overview of the results, highlighting separately the influence of the cell membrane's mechanical properties and aspect ratio of the particle and the sensitivity of the endocytosis to these parameters. Section 5 discusses the results and hypotheses. Finally, Section 6 summarizes the results and discussions with the concluding remarks and perspectives for further work.

<figure><img src="image_3.png" /><figcaption>**FIGURE 1** The cellular uptake of a NP is influenced by its mechanical and geometrical properties.</figcaption></figure>

## 2 | METHODS

### 2.1 | Investigation strategy and hypotheses

In order to identify the influence of the mechanical properties of the cell membrane and of the aspect ratio of the particle, the following strategy is adopted. First, the influence of the cell membrane's mechanical properties on the uptake of a rigid circular NP is studied. An elliptic NP is then considered and the influence of its aspect ratio is investigated. A coupled study involving the aspect ratio of the NP and the cell's mechanical properties is afterwards conducted. Finally, a metamodel is built in order to perform a sensitivity analysis of different parameters on the cellular uptake of the NPs.

In this paper, the focus is put on rigid particles for two reasons: (i) reasonable computational cost and complexity allowing us to subsequently carry out a stochastic modeling and (ii) a moderate number of influencing parameters. Nonetheless, given that the rigidity of the NP is typically at least 2 orders of magnitude greater than that of the cell membrane<sup>20,36,38</sup>, the assumption of a rigid NP is acceptable in most cases. The model focuses on the non-receptor-mediated endocytosis of a rigid elliptic NP and does not apply to phagocytes, as their physiology is specific to engulf objects from the extra-cellular middle<sup>39</sup>. A 2D model of the system is used for the sake of simplification and to decrease the computation costs. However, the former could have been straightforwardly extended to 3D<sup>30</sup>.

The diameter of rigid NPs used for drug delivery is usually around 100 nm<sup>20</sup> while that of the cell is in the order of 10 µm<sup>40</sup>. Moreover, the cell membrane's thickness is about 4 nm<sup>41</sup>. These scales allow us to model the system at the scale of the NP as we can neglect the mechanical interactions due to the membrane's microstructure. Moreover, it also allows us to neglect the cell's curvature compared to the NP's one.

The model consists in evaluating the variation of total energy during the endocytosis and then identifying the wrapping degree at equilibrium.
## 2.2 | Nanoparticle - membrane system

The system is illustrated in Figure 2. It is divided into four regions: region 1 is the free part of the particle, region 2r and 2l are the free parts of the membrane, respectively on the right and left sides of the particle and region 3 is the contact region between the particle and the membrane.

<figure><img src="image_4.png" /><figcaption>FIGURE 2 Illustration of the nanoparticle-membrane interaction and the coordinates.</figcaption></figure>

The spatial coordinates are denoted with the tuple $(r, z)$ wherein the origin is located at the intersection between the regions 1, 2l and 3. Each region $i \in \{1, 2l, 2r, 3\}$ is parametrized by its arclength $s_i \in [0; l_i]$ and by the angle $\psi_i$ used to compute the curvature and subsequently the bending energy. The system is symmetric, therefore the length of the two free parts of the membrane are equal and will henceforth be denoted as $l_2$ ($l_{2l} = l_{2r} = l_2$). Animated figures of the system can be found in the Supplemental Material.

The lengths $l_1$ and $l_3$ are defined in terms of the wrapping degree $f$, which corresponds to the fraction of the particle wrapped by the membrane. Hence, $l_3 = pf$, $l_1 = p(1 - f)$ and thus $f = l_3(l_1 + l_3)^{-1}$, where $p$ is the circumference of the particle. The membrane is supposed to be flat far from the nanoparticle, as the effects of the latter on the membrane weaken. A convergence study was then conducted on $l_2$ to verify that setting $l_2 = 20a$, where $a$ is the relative radius of the ellipse, defined as the ratio between the circumference of the particle and $200\pi$ nm (which is the circumference of a circular particle with a radius of 100 nm), satisfies the hypothesis of membrane flatness at both extremities. See¹⁸ for more details.

## 2.3 | Evaluation of the variation of energy

The variation of energy is due to the bending of the membrane $\Delta E_b$, the adhesion between the particle and the membrane $\Delta E_{\bar{\gamma}}$ and the membrane tension $\Delta E_{\bar{\sigma}}$¹,³⁷. Hence, the former reads:

$$ \Delta E = \Delta E_b + \Delta E_{\bar{\gamma}} + \Delta E_{\bar{\sigma}} \quad (1) $$

The bending energy of the membrane is decomposed as:

$$ \Delta E_b = \frac{1}{2} \kappa_{2r} \int_{0}^{l_2} (\psi_{2r} - \tilde{c}_{2r})^2 ds_{2r} + \frac{1}{2} \kappa_{2l} \int_{0}^{l_2} (\psi_{2l} - \tilde{c}_{2l})^2 ds_{2l} + \frac{1}{2} \kappa_3 \int_{0}^{l_3} (\psi_3 - \tilde{c}_3)^2 ds_3, \quad (2) $$

where the overdot denotes the derivative with respect to $s$. $\kappa_i$ and $\tilde{c}_i$ denote the bending rigidity and the initial curvature of the region $i \in \{1, 2l, 2r, 3\}$, respectively. Since the membrane is considered initially flat, its initial curvature is zero, leading to $\tilde{c}_{2r} = \tilde{c}_{2l} = \tilde{c}_3 = 0$. Since we assume that the bending rigidity of each region $\kappa_i$ is independent of $s$ (homogeneous membrane), the former is put outside of the corresponding integrals. Hence, we have $\kappa_{2r} = \kappa_{2l} = \kappa_2$. The particle being considered rigid, it does not deform. Consequently, the bending of the region 3 is only due to the deformations of the membrane. Thus, the bending rigidity of the zone 3 is equal to that of the zone 2, i.e. $\kappa_3 = \kappa_2$. Furthermore, as defined in Figure 2, $\psi_{2l} = 2\pi - \psi_{2r}$. As such, equation 2 becomes:
$$ \Delta E_b = \kappa_2 \underbrace{\int_0^{l_2} \dot{\psi}_{2r}^2 ds_3}_{\Delta E_{b2}} + \underbrace{\frac{1}{2} \kappa_2 \int_0^{l_3} \dot{\psi}_3^2 ds_3}_{\Delta E_{b3}} \quad (3) $$

The bending energy of the free part of the membrane is computed based on the evolution of $\psi_{2r}(s_2)$ which could be analytically obtained using the following formula$^{42}$:

$$ \psi_{2r}(s_2) = 4 \arctan \left[ \tan \left( \frac{\psi_{2r}(0)}{4} \right) \exp \left( -\frac{s_2}{a} \sqrt{\frac{\bar{\sigma}}{2}} \right) \right] \quad (4) $$

The spatial coordinates $(r, z)$ in regions $2r$ and $2l$ are then evaluated as follows$^{18}$:

$$ r_{2r}(s_2) = r_{2r}(0) + s_2 - \sqrt{\frac{\bar{\sigma}}{2}} \frac{1 - \cos(\psi_{2r}(0))}{\coth\left(\sqrt{\frac{\bar{\sigma}}{2}} s_2\right) + \cos\left(\frac{\psi_{2r}(0)}{2}\right)} \quad (5a) $$

$$ r_{2l}(s_2) = r_{2r}(0) - r_{2r}(s_2) \quad (5b) $$

$$ z_{2r}(s_2) = z_{2r}(0) + \sqrt{\frac{8}{\bar{\sigma}}} \sin\left(\frac{\psi_{2r}(0)}{2}\right) \left[ 1 - \frac{\coth\left(\sqrt{\frac{\bar{\sigma}}{2}} s_2\right) + \cos\left(\frac{\psi_{2r}(0)}{2}\right)}{\coth\left(\sqrt{\frac{\bar{\sigma}}{2}} s_2\right) + \cos\left(\frac{\psi_{2r}(0)}{2}\right)} \right] \quad (6a) $$

$$ z_{2l}(s_2) = z_{2r}(s_2) \quad (6b) $$

In which $\psi_2(0)$ is the value of $\psi_2$ at the intersection between the regions 3 and $2r$ (represented with a square in Figure 2). The points ($s_{2r} = 0$) and ($s_3 = l_3$) are coincident. Consequently, $\psi_{2r}(s_2 = 0)$ can be calculated from $\psi_3(s_3 = l_3)$ which is well known from the definition of the ellipse (as well as $\psi_1, r_1, z_1, r_3$ and $z_3$). Then, $\Delta E_{b2}$ can be readily computed.

The adhesion energy is defined as $\Delta E_{\gamma} = -\gamma l_3$, in which $\gamma$ is the lineic adhesion force between the membrane and the particle. The tension energy is defined as $\Delta E_{\sigma} = \sigma \Delta l$, wherein $\sigma$ is the membrane tension and $\Delta l$ is the length of stretched membrane.

Finally, in order to ease further comparison, the adimensional energy variation $\overline{\Delta E}$ is introduced as $\overline{\Delta E} = \Delta E \frac{2a}{\kappa_2}$ and Equation 1 is reformulated as Equation 7, in which adimensional variables $\bar{\gamma}$ and $\bar{\sigma}$ are introduced as $\bar{\gamma} = \gamma \frac{2a^2}{\kappa_2}$ and $\bar{\sigma} = \sigma \frac{2a^2}{\kappa_2}$. A Python script to compute the evolution of $\overline{\Delta E}$ in terms of $f$ for a given set of input parameters can be found in the Supplemental Material.

$$ \overline{\Delta E} = \frac{a}{4} \underbrace{\int_{0}^{l_2} \dot{\psi}_3^2 ds_3}_{\overline{\Delta E_{b2}}} + \frac{a}{4} \int_{0}^{l_2} \dot{\psi}_{2r}^2 ds_2 - \frac{1}{4a} \bar{\gamma} l_3 + \underbrace{\frac{1}{4a} \bar{\sigma} (2l_2 + l_3 - r_{2r}(l_2) + r_{2l}(l_2))}_{\overline{\Delta E_{\bar{\gamma}}}} \quad (7) $$

In this paper, all the investigated parameters are supposed constant during endocytosis and their values are set based on experimental results from the literature: $\kappa_2 \approx 10^{-18} \text{N.m}$, $\gamma \approx 10^{-3} \text{N.m}^{-1}$ and $\sigma \approx 10^{-5} \text{N.m}^{-1}$.$^{4,37,43,44}$. The units commonly used for these variables are k$_B$T or erg for $\kappa_2$, k$_B$T.nm$^{-2}$ or erg.cm$^{-2}$ for $\gamma$ and dyne.cm$^{-1}$ for $\sigma$.$^{4,45,46}$ but we chose to convert them to the SI units in order to allow a better understanding by an unfamiliar audience.

An example of the contributions of $\overline{\Delta E}_b$, $\overline{\Delta E}_{\bar{\gamma}}$ and $\overline{\Delta E}_{\bar{\sigma}}$ on the total energy variation $\overline{\Delta E}$ are presented in Figure 3. With the particular set of parameters used to generate this figure, it appears that tension and adhesion contribute more to the total energy than bending.

## 2.4 | Phases

The next step of the study consists in evaluating the influence of the aspect ratio $\bar{r}$ (ratio between the semi-major and semi-minor axes) of the particle on the wrapping degree at equilibrium, denoted by $\tilde{f}$. Contrary to circular NPs, $\overline{\Delta E}(f)$ may have
<figure><img src="image_5.png" /><figcaption>FIGURE 3 Contribution of bending, adhesion and tension on the total energy variation, $\bar{\gamma} = 10$, $\bar{\sigma} = 2$.</figcaption></figure>

multiple local minima. It is then necessary to provide a clear definition of the equilibrium. According to the literature$^{31,32}$, there is no mechanism that could contribute to overcome the energy barriers following a local minimum, except the energy due to thermal fluctuations. However, the latter are too weak to overcome the energy barriers we met in this study. Consequently, the equilibrium is defined as the first local minimum.

Three phases are schematically defined in Figure 4. The first phase (no wrapping) gathers the configurations which lead to $\tilde{f} < 0.2$. For this range of wrapping degrees, one can consider that the endocytosis aborts soon after the particle touched the cell (no wrapping). The phase 3 (full wrapping) includes the configurations which lead to an equilibrium where the two sides of the free membrane cross, i.e. $\max r_{2l}(s_{2l}) \ge \min r_{2r}(s_{2r})$. The former is called full wrapping, since the particle is already engulfed by the cell. The phase 2 is composed of all the intermediate configurations, which lead to partial wrapping. A Python script to determine the wrapping phase at equilibrium for a set of input parameters is available in the Supplemental Material.

<figure><img src="image_6.png" /><figcaption>FIGURE 4 Illustration of the three wrapping phases: no wrapping (left), partial wrapping (middle) and full wrapping (right).</figcaption></figure>

To investigate the influence of the membrane tension $\bar{\sigma}$ and adhesion $\bar{\gamma}$, $\overline{\Delta E}(f)$ is computed for $f \in [0, 1]$ for a given aspect ratio $\bar{r}$, say 1 (a circular NP), and a range of values for $\bar{\sigma}$ and $\bar{\gamma}$. It is necessary to set ranges of values because of the variability of these parameters (considering the variety of cell types and the diversity of their membrane composition$^{47,48}$). Based on the literature$^{1,18,30,43}$ and considering a radius of 100 nm, the ranges were considered as $\bar{\sigma} \in [0.5, 5.5]$ and $\bar{\gamma} \in [1, 10]$. For each one of these testcases, the phases are identified, as illustrated in Figure 6.

From the phase diagram, the proportion of each phase is evaluated. Then, to investigate the influence of another parameter, here, the aspect ratio $\bar{r}$ of the elliptic particle, one generates the phase diagram and then computes the proportion of the three phases, denoted as $\psi_1$, $\psi i_2$ and $\psi i_3$, for each value of the investigated parameter and compare them in a phase comparison graph, as illustrated further in Section 4.

To investigate the influence of $\bar{r}$, the circumference of the NP was fixed to $2\pi$ and the values of the semi-minor and semi-major axes were defined based on the Ramanujan's approximation of the circumference of an ellipse$^{49}$, for a given aspect ratio. The
<table><caption>TABLE 1 Geometric parameters of the investigated ellipses, the circumference is set to 200π nm.</caption><thead><tr><th rowspan="2"></th><th colspan="4">vertical ellipse</th><th colspan="5">horizontal ellipse</th></tr><tr><th>1/6</th><th>1/4</th><th>1/3</th><th>1/2</th><th>1</th><th>2</th><th>3</th><th>4</th><th>6</th></tr></thead><tbody><tr><th>aspect ratio</th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>semi-major axis (x100) [nm]</th><td>0.25</td><td>0.37</td><td>0.47</td><td>0.65</td><td>1</td><td>1.30</td><td>1.41</td><td>1.46</td><td>1.5</td></tr><tr><th>semi-minor axis (x100) [nm]</th><td>1.5</td><td>1.46</td><td>1.41</td><td>1.30</td><td>1</td><td>0.65</td><td>0.47</td><td>0.37</td><td>0.25</td></tr></tbody></table>

geometric parameters used to perform the investigation are detailed in the Table 1. The results, necessary to build the response surface used to carry out sensitivity analyses detailed in the next section, are presented in Section 4.

## 3 | SURROGATE MODELING AND SENSITIVITY ANALYSIS

This section briefly introduces the uncertainty quantification (UQ) methods that will be used in this paper. Let $\tilde{F}$ denote the positive, second-order stochastic wrapping degree at equilibrium. The former is a function of a vector of three input random variables $X = \{\tilde{\Gamma}, \tilde{\Sigma}, \tilde{R}\}^\top$, i.e. $\tilde{F} = \Xi(X)$. To investigate the probabilistic content of the random variable $\tilde{F}$, classical sampling techniques such as Monte Carlo (MC) simulation can be used. However, the convergence is only guaranteed for a large number of datasets implying high computational costs. As such, metamodeling approach will be used as an appropriate alternative allowing us to build fast-to-evaluate surrogate models requiring significantly a lower number of calls to the numerical solver. For this purpose, we will use two popular surrogate models: Kriging$^{50-53}$ and Polynomial Chaos Expansion (PCE)$^{54-57}$ which are briefly discussed in Sections 3.1 and 3.2. Finally in Section 3.3, the Sobol indices are introduced to quantify the global sensitivities of $\tilde{F}$ with respect to each input random variable.

### 3.1 | Kriging

Kriging, also known as Gaussian process regression, is based on the decomposition of the random variable $\tilde{F}$ into a deterministic mean (trend) and a residual Gaussian process as $\Xi^K(\mathbf{x}) = \mathbf{a}^\top g(\mathbf{x}) + \sigma_Z Z(\mathbf{x}; \mathcal{L})^52$. In this equation $\mathbf{x}$ is a realization of the random vector $\mathbf{X}$, $E[\Xi^K] = \mathbf{a}^\top g(\mathbf{x})$ is the mean of $\tilde{F}$, represented by a vector $\mathbf{a}$ and a set of predefined shape functions $g_i, i \in \{0, ..., P-1\}$, and $Z(\mathbf{x}; \mathcal{L})$ is a standard Gaussian random field defined by its auto-correlation function (ACF) $R(\mathbf{x}, \mathbf{x}'; \mathcal{L}) = E[Z(\mathbf{x}; \mathcal{L})Z(\mathbf{x}'; \mathcal{L})]$. The vector $\mathcal{L} = \{\ell_1, \ell_2, \ell_3\}^\top$ is called the correlation distance that should be calculated. Note also that $\sigma_Z$ represents the standard deviation of the Kriging's residual.

The objective is to evaluate the unknown parameters using a sample of size $N$ (experimental design), $\mathbf{X} = \{\mathbf{x}_1, ..., \mathbf{x}_N\}^\top$ and the corresponding wrapping degrees at equilibrium $\tilde{F} = \{\tilde{f}_1, ..., \tilde{f}_N\}^\top$ in which $\tilde{f}_i = \Xi(\mathbf{x}_i)$, $i \in \{1, ..., N\}$ is the value of $\tilde{f}$ for $i^{th}$ realization of input parameters. The least-squares estimate of the $P \times 1$ coefficients vector $\mathbf{a}$ is $\hat{\mathbf{a}} = (G^\top R^{-1} G)^{-1} G^\top R^{-1} \tilde{F}$ in which the $N \times P$ matrix $G$ is defined as $G_{ij} = g_j(\mathbf{x}_i)$, $i \in \{1, ..., N\}$, $j \in \{1, ..., P\}$ and $R$ is the $N \times N$ correlation matrix whose components are $R_{ij} = R(|\mathbf{x}_i - \mathbf{x}_j|; \mathcal{L})$, $i, j \in \{1, ..., N\}$. The main steps of the Kriging method are summarized as follows. The first step consists in choosing the trend function. Then a model for the ACF is chosen$^{58}$ with a correlation length vector $\mathcal{L}$ to be determined. Then, using the experimental design $\mathbf{X}$ and the values of the wrapping degree $\tilde{f}_i$, a maximum likelihood optimization problem is solved to calculate $\mathcal{L}^{59}$. The latter is then used to obtain the coefficients of the trend $\hat{\mathbf{a}}$, the mean and the variance of the Kriging predictor (see$^{60}$ for more details).

Finally, we use a leave-one-out (LOO) cross-validation error estimator $\epsilon_{LOO}$ in order to quantify the quality of the metamodels. The former is calculated by comparing the predictions obtained using $N$ surrogate models $\mathbf{X}^{(-j)} = \{\mathbf{x}_i, i = 1, ..., N, i \neq j\}$ for each of which one point is removed from the initial experimental design, with the real value of $\tilde{F}$ at the excluded point$^{61}$.

### 3.2 | Polynomial Chaos Expansion (PCE)

PCE is a functional representation of the random quantity of interest wherein the latter is written as an infinite linear combination of some orthogonal polynomials with respect to the probability density functions (PDFs) of input parameters$^{62-64}$. The exact
infinite expansion is then truncated up to some degree p as:

$$ \Xi_p^{PCE}(\mathbf{x}) = \sum_{i=0}^{P-1} a_i \Psi_i(\zeta) = \mathbf{a}^\top \Psi(\zeta), \qquad (8) $$

in which $P = (p+3)$ polynomial basis functions $\Psi$ and deterministic coefficients $\mathbf{a}$ are used. An isoprobabilistic transformation $\mathcal{T}$ links the standard random vector $\zeta = (\zeta_1, \zeta_2, \zeta_3)$ to the random input vector $\mathbf{X}$, i.e. $\mathbf{X} = \mathcal{T}(\zeta)$. The type of polynomials is determined based on the PDFs of the input random parameters. For instance, Legendre and Jacobi polynomials are used for uniform and beta random variables. The basis functions are then constructed by multiplying the one dimensional bases. The expansion coefficients are then calculated using a least squares approach via $\hat{\mathbf{a}} = (\mathbf{G}^\top \mathbf{G})^{-1} \mathbf{G}^\top \hat{\mathbf{F}}$ wherein the $N \times P$ matrix $\mathbf{G}$ is defined as $G_{ij} = \Psi_j(\zeta_i)$. The mean and variance of the random wrapping degree at equilibrium are subsequently calculated in terms of the expansion coefficients via $\mathbb{E}[\Xi_p^{PCE}] = a_0$ and $\text{Var}[\Xi_p^{PCE}] = \sum_{i=1}^{P-1} a_i^2$, respectively. Similar to the previous case, a LOO cross-validation error estimation will be calculated to assess the quality of the model. Finally, the Sobol sensitivity indices are calculated for both types of metamodeling techniques as will be discussed in the next Section.

## 3.3 | Sensitivity analysis

To calculate the sensitivity indices⁶⁵, we use Quasi Monte Carlo (QMC) estimators for Kriging and direct estimations based on the expansion coefficients for PCE (see⁶⁶ for more details). A QMC-based estimator for the Sobol index of the random variable $X_i$ ($i \in \{1, 2, 3\}$) which describes the sensitivity of $\hat{\mathbf{F}}$ to $X_i$ writes:

$$ S_i = \frac{\sigma_i^2}{\sigma_{\hat{\mathbf{F}}}^2} \approx \frac{1}{N \sigma_{\hat{\mathbf{F}}}^2} \sum_{j=1}^{N} \hat{F}_j(S_2) (\hat{F}_j(S_1') - \hat{F}_j(S_1)), \qquad (9) $$

in which $S_1$ and $S_2$ are two independent QMC samples of size $N \times 3$ of input random parameters and $S_1^i$ is the matrix $S_1$ whose $i$-th column is replaced by the $i$-th column of $S_2$. Note that in (9) $\sigma_{\hat{\mathbf{F}}}^2 = \text{Var}[\hat{\mathbf{F}}]$ is the estimate of the total variance of the random wrapping degree at equilibrium $\hat{\mathbf{F}}$ and $\sigma_i^2 = \text{Var}[\mathbb{E}[\hat{\mathbf{F}}|X_i]]$ is the part of the variance of $\hat{\mathbf{F}}$ resulting from $X_i$. The sensitivity of the wrapping degree to the interactions between input random variables $X_i$ and $X_j$ is described by the second order indices $S_{ij}$ ($1 \le i \ne j \le 3$) excluding the first order effects:

$$ S_{ij} = \frac{\sigma_{ij}^2}{\sigma_{\hat{\mathbf{F}}}^2} = \frac{\text{Var}[\mathbb{E}[\hat{\mathbf{F}}|X_i X_j]]}{\text{Var}[\hat{\mathbf{F}}]} - S_i - S_j, \qquad (10) $$

which could be similarly generalized to higher order interactions. Hence, there are $2^3 - 1 = 7$ different Sobol indices for different level of interactions. However, instead of calculating all these indices, total Sobol indices are used to estimate the total contribution of each random variable $X_i$ ($1 \le i \le 3$) to the variance of the wrapping degree at equilibrium⁶⁷:

$$ S_{T_i} = S_i + \sum_{i \ne j} S_{ij} + \sum_{(j<k) \ne i} S_{ijk} + \dots \qquad (11) $$

A QMC-based estimator for the total Sobol index then reads:

$$ S_{T_i} \approx \frac{1}{2N\sigma_{\hat{\mathbf{F}}}^2} \sum_{j=1}^{N} (\hat{F}_j(\Xi_1^i) - \hat{F}_j(\Xi_1))^2. \qquad (12) $$

As mentioned earlier, for the PCE response surface, the sensitivity indices could be analytically estimated in terms of the coefficients of the expansion. The readers are referred to Sudret⁵⁷ for more details. The numerical applications and the corresponding results are discussed in Section 4.

# 4 | RESULTS

## 4.1 | Mechanical parameters $\bar{\gamma}$ and $\bar{\sigma}$

The mechanical parameters $\bar{\gamma}$ and $\bar{\sigma}$ tend to alter the evolution of $\Delta \bar{E}$ with respect to $f$, as illustrated in Figure 5. Indeed, a larger $\bar{\gamma}$ as well as a smaller $\bar{\sigma}$, tend to decrease the values of $\Delta \bar{E}$. The membrane tension does not influence $\Delta \bar{E}$ for $f < 0.2$, while
the effect of $\bar{\gamma}$ is significant for all values of $f$. This matches with the decomposition of $\overline{\Delta E}$ presented in Figure 3 in which one can see that $\overline{\Delta E}_b$ and $\overline{\Delta E}_\sigma$ (the energy variations in which $\bar{\sigma}$ contributes) do not change for low values of $f$, contrary to $\overline{\Delta E}_{\bar{\gamma}}$ which decreases significantly.

<figure><img src="image_7.png" /><figcaption>FIGURE 5 Effect of adhesion (left) $\bar{\gamma}$ with $\bar{\sigma} = 2$ and (right) membrane tension $\sigma$ with $\bar{\gamma} = 10$ on $\overline{\Delta E}(f)$.</figcaption></figure>

The conclusions from Figure 5 apply to the range of values of $\bar{\gamma}$ and $\bar{\sigma}$ tested and can be generalized to larger ranges of $\bar{\gamma}$ and $\bar{\sigma}$, as illustrated Figure 6, which reveals that larger values of $\bar{\gamma}$ and lower values of $\bar{\sigma}$ favor the occurrence of the phase 3. That means that a NP is more likely to be uptaken by a cell whose membrane tension is low and which adheres well to the NP.

These results match some of the previous studies conducted by Yi and Gao$^{1,15,18}$, Zhang$^{4,68}$ and the ones presented in the review from Ding$^{21}$ leading to similar results as ours, being that the uptake of NP is enhanced by low $\bar{\sigma}$ and large $\bar{\gamma}$. These conclusions were also drawn experimentally in$^{69}$ for instance, which showed that the adhesion between colloidal drug delivery systems and macrophages tend to increase their uptake. Similar conclusions are also presented in$^{70}$.

<figure><img src="image_8.png" /><figcaption>FIGURE 6 Phase diagram for a circular NP. Colors: values of $\tilde{f}$, from 0 (lightest) to 1 (darkest).</figcaption></figure>

An illustration of the effect of membrane tension and adhesion on the wrapping is given in Figure 7. In this figure, even though the adhesion does not affect the shape of the membrane during the endocytosis, it is clear that for low membrane tension, the two sides of the membrane tend to get in contact more easily, leading to full wrapping.
<figure><img src="image_9.png" /><figcaption>**FIGURE 7** Effect of adhesion (top) $\bar{\gamma}$ with $\bar{\sigma} = 2$ and membrane tension (bottom) $\bar{\sigma}$ with $\bar{\gamma} = 10$ on the wrapping for $f = 0.8$. The $r$ axis was truncated to $[-10; 10]$. For the same value of $f$, the NP is fully wrapped by the cell only for $\bar{\sigma} = 0.5$.</figcaption></figure>

## 4.2 | Aspect ratio of the elliptic particle

The aspect ratio $\bar{r}$ of the elliptic particle affects the evolution of $\overline{\Delta E}$ with respect to $f$, as depicted in Figure 8.

<figure><img src="image_10.png" /><figcaption>**FIGURE 8** $\overline{\Delta E}(f)$ for the aspect ratios introduced in the table 1. The mechanical parameters are $(\bar{\gamma}, \bar{\sigma}) = (10, 2)$.</figcaption></figure>

Except for the cases of circular and slightly vertical particle ($\bar{r} \in [1/2, 1]$), an energy barrier occurs at low values of $f$ for vertical ellipses and at $f \approx 0.5$ for horizontal ones. It is easy for the horizontal ellipse to reach $f = 0.5$, since it is basically the contact between the particle and the membrane, without even bending the latter. The vertical particle tends to stay at low values of $f$ because a sharper bent on the cell membrane is necessary to increase $f$ with the same amount.

However, the energy barrier to overcome is lower for vertical than for horizontal particles. As illustrated in Figure 9, this barrier for vertical ellipses is reduced for lower $\bar{\sigma}$ and higher $\bar{\gamma}$, while the energy barrier for horizontal ellipses is only slightly reduced by lowering $\bar{\sigma}$.

For a better understanding of the influence of $\bar{r}$ and to generalize these observations to a large range of values of $\bar{\gamma}$ and $\bar{\sigma}$, phase diagrams are generated for all aspect ratios (see 6 for the phase diagram generated for a circular NP). From the latter, the proportion of each phase is computed and finally compared in Figure 10, which confirms the previous assumptions. Indeed, the horizontal particles do not pass from phase 2 while the vertical ones remain in the phase 1. However, vertical particles that overcome the phase 1 tend to reach the phase 3, contrary to horizontal particles. For highly elongated particles neither the vertical nor the horizontal particles reach the phase 3.
<figure><img src="image_11.png" /><figcaption>**FIGURE 9** Influence of (left) $\bar{\gamma}$ with $\bar{\sigma} = 2$ and (right) $\bar{\sigma}$ with $\bar{\gamma} = 10$ on $\overline{\Delta E}(f)$ for elliptic particles. Dark lines: vertical particle ($\bar{r} = 1/4$); Light lines: horizontal particle ($\bar{r} = 4$).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**FIGURE 10** Proportion of phase 1, 2, 3 (respectively denoted as $\psi_1$, $\psi_2$, and $\psi_3$) in terms of $\bar{r}$.</figcaption></figure>

These results match the ones obtained experimentally by Champion *et al.*$^{13,14}$ who manufactured NPs with various aspect ratios and shapes and observed that the circular and vertical NPs are engulfed faster by macrophages than the ones with other shapes. However, their results show that the highly elongated vertical NPs could not be engulfed.

## 4.3 | Probabilistic modeling of the random wrapping degree at equilibrium

### 4.3.1 | Random input parameters

In this section the mechanical parameters of the NP $\bar{\gamma}$ and $\bar{\sigma}$, along with its aspect ratio $\bar{r}$ are modeled as random variables (denoted by ($\bar{\Gamma}$, $\bar{\Sigma}$, $\bar{R}$)). The objective is to first model these parameters via some appropriate probability distributions and then to investigate the probabilistic content of the random wrapping degree at equilibrium $\tilde{F} = \Xi(\bar{\Gamma}, \bar{\Sigma}, \bar{R})$. Having recourse to the maximum entropy principle$^{71}$, since the only available information are the lower and upper bounds, the best candidate for the Probability Density Function (PDF) (maximizing the entropy measure) is the uniform distribution. However, a uniform distribution for the ratio $\bar{R}$ will result in realizations that are biased towards more horizontal NPs since they are six times more likely to be produced. To overcome this issue, we used a transformed beta distribution to model $\bar{R}$, i.e. $\bar{R} \sim \beta(0.5, 1.7; 1/6, 6)$. The particularity of this PDF is that its median is 1 so that $P(\bar{R} < 1) = P(\bar{R} > 1) = 0.5$ almost surely. As such, the input parameters are modeled as three independent random variables $\bar{\Gamma} \sim U(1, 8)$, $\bar{\Sigma} \sim U(0.5, 5.5)$ and $\bar{R} \sim \beta(0.5, 1.7; 1/6, 6)$, the
interval of definition of these variables being defined in Section 2.4. The coefficients of variation (CoV) of $(\bar{\Gamma}, \bar{\Sigma}, \bar{R})$ being about (45%, 48%, 92%), implies significant fluctuation levels of the input parameters.

### 4.3.2 | Sampling of the random wrapping degree at equilibrium

Once the uncertainties related to each input parameter are characterized, 1000 realizations based on QMC sampling method are generated. The corresponding values of the wrapping degree at equilibrium, i.e. $\tilde{f}$, for each realization are then calculated. Figure 11 depicts the convergence of the average and the CoV of the random wrapping degree at equilibrium in terms of the number of samples along with the corresponding 95% two-sided confidence intervals (shaded zones in gray). The confidence bounds for the CoV are calculated by using bootstrap resampling technique. The right plot of Figure 11 shows a considerable fluctuation level of the wrapping degree at equilibrium (around 94%) which results from large variabilities of the input parameters. It seems that one should use more than 1000 realizations to guarantee the convergence of the first and second order statistics of $\tilde{F}$. This highlights the importance of using surrogate models allowing us to approximate the probabilistic content of $\tilde{F}$ using small sample sizes.

<figure><img src="image_13.png" /><figcaption>**FIGURE 11** Convergence of the mean (left) and coefficient of variation (right) of the random wrapping degree at equilibrium in terms of the sample size (thick solid black curves) along with 95% two-sided confidence intervals denoted by the gray regions.</figcaption></figure>

### 4.3.3 | Surrogate modeling

In this section, the PDF of the random wrapping degree at equilibrium $\tilde{F}$ is modeled using two different surrogate modeling approaches.

The Kriging model is constructed based on the following parameters. For the trend part, a quadratic function as $a^\top g(x) = a_0 + \sum_{i=1}^M a_i x_i + \sum_{i=1}^M \sum_{j=1}^M a_{ij} x_i x_j$ is used with coefficients to be determined. A von Kármán model with a shape parameter of $H = 1.5$ is used for the ACF of the underlying Gaussian process. A hybrid genetic algorithm$^{72}$ with upper and lower bounds of [0.001, 10] for each $\ell_i$ ($1 \le i \le 3$) is used to solve the optimization problem for the hyperparameter $\ell$.

On the other hand, to build the PCE response surface, the polynomial degree $p$ is chosen such that the corresponding LOO error is minimized. For this purpose, the latter is calculated for different values of $p$ from 1 to 30. Figure 12 shows that $p = 23$ minimizes the error ($\epsilon_{LOO} = 0.045$) and thus will be used to construct the PCE response surface. These metamodels could now be evaluated with a near-zero computational cost. Using $10^5$ realizations, the estimated PDFs of $\tilde{F}$ are compared to the histogram of the simulated data in Figure 13 which shows a bimodal distribution of the wrapping degree. We can then estimate the probability of occurrence of each phase. For instance, the probability of no wrapping $P(\tilde{F} < 0.2)$ (phase 1) is 0.033, 0.032 and 0.036 based on the simulation data, Kriging and PCE, respectively, which shows the accuracy of the metamodels (with relative errors of 3% and 9% for Kriging and PCE, respectively).

The Sobol sensitivity indices are finally estimated using both approaches. Kriging yields the following first-order indices 0.12, 0.06 and 0.41. On the other hand, based on the PCE method, we get the values of 0.15, 0.07 and 0.46 (for $\bar{\Gamma}$, $\bar{\Sigma}$ and $\bar{R}$, respectively) which reveals that if we consider the effect of the variability of each single input random variables on the variance of the random wrapping degree at equilibrium, the aspect ratio of the nanoparticle $\bar{r}$ is the most influential. Then the adhesion $\bar{\gamma}$
<figure><img src="image_14.png" /><figcaption>FIGURE 12 Variation of the LOO error estimate in terms of the degree *p* of the PCE.</figcaption></figure>

contributes twice as large as the membrane tension $\bar{\sigma}$. The total sensitivity indices obtained with PCE method for $\bar{\Gamma}$, $\bar{\Sigma}$ and $\bar{R}$ are 0.4, 0.24 and 0.74, respectively. The difference $S_{T_i} - S_i$ shows the influence of the interactions of different orders on the variance of $\bar{F}$. This difference is in the same order as for the first-order indices, i.e. the interactions of $\bar{R}$ with other inputs contribute the most on the variability of the wrapping degree at equilibrium. The readers can use open-source software like OpenTURNS (openturns.org) and UQLab (uqlab.com) to reproduce the results.

<figure><img src="image_15.png" /><figcaption>FIGURE 13 Estimations of the PDF of the random wrapping degree at equilibrium obtained by Kriging (dashed line) and PCE (solid curve) metamodels. The histogram corresponds to the PDF of the output data.</figcaption></figure>

# 5 | DISCUSSION

## 5.1 | Cell membrane reorganization during endocytosis

This paper only considers constant values of $\bar{\sigma}$ and $\bar{\gamma}$ during endocytosis, while they are actually variable with respect to $f$, as explained in$^{27,38,73-80}$. Indeed, the cell tends to adapt to the phenomenon by simultaneously increasing the amount of interactions at the contact region (and consequently increasing the adhesion) and reorganizing its actin network, helping the membrane deformation (and consequently reducing the tension). It may explain why the *in vivo* observations$^{13,14}$ concluded that vertical elliptic NPs are more likely to enter the cell than horizontal ones in almost all cases. Indeed, the variation of the mechanical properties leads to vanish the energy barrier in most cases.
## 5.2 | Heterogeneity and composition of the cell membrane

In this study, the NP was considered about 100 times smaller than the cell and the system was set at the NP's scale (around 100 nm). However, membrane's thickness is around 4 nm and is, among other things, composed of several trans-membrane proteins whose lengths can reach 16 nm⁴¹. Consequently, they could play a role in this process and directly interact with the particle. Moreover, the membrane is highly heterogeneous and is composed of various components (mainly carbohydrates, proteins and lipids)⁴⁷,⁸¹. The outer part of the membrane is also covered by a sugar layer called glycocalyx or fuzzy coat⁸²,⁸³, which alters the interaction with the extra-cellular middle by reducing the adhesion. The presence of glycocalyx is usually more important in cancer cells than in healthy cells²⁶, meaning that it is a parameter to take into account to model the adhesion between the cell and the NP. In addition, the composition of the cell membrane may change alongside its circumference. For instance, the presence of lipid rafts⁴⁸,⁸⁴ will locally increase the bending rigidity of the membrane. These observations about the structure of the cell membrane challenge the modeling of the region 3 as homogeneous alongside the cell. However, since the study was performed for a range of $\tilde{\gamma}$ and $\tilde{\sigma}$, our results remain valid as they account for the variability of the membrane's properties.

## 5.3 | Rigidity of the NP

Our model was built considering a rigid NP. As explained in the hypotheses section, this choice was made for simplification purposes and it appeared to fit with many of the existing NPs. However, NPs with low bending rigidity (in the same order of magnitude as the membrane) exist as well²⁰. Moreover, previous studies demonstrated that the rigidity of the NP is a prominent parameter of endocytosis¹,¹⁸. Consequently, additional investigations should be done accounting for the rigidity of the NP to provide results adapted to a larger range of NPs.

## 5.4 | Rotation of the elliptic NP during endocytosis

It was also found experimentally⁹ and numerically¹⁸,⁸⁵,⁸⁶ that if an elliptic particle touches the membrane by its horizontal side, it will tend to rotate to be engulfed vertically. However, this phenomenon appears to depend on the membrane tension, as demonstrated in¹⁵: the NP is wrapped in its vertical configuration if the tension is lower than a threshold value of $\tilde{\sigma}$ and in its horizontal configuration otherwise.

## 5.5 | 2D model

Zhou et al.³⁰ suppose that a 2D model suits for the study of the adhesion of a vesicle to a substrate. However, Yi¹ performed both 2D and 3D axisymmetric models of the cellular uptake of a NP and found the same results for circular NPs but different behavior of $\Delta E(f)$ for other shapes. Indeed, for non-circular particles, studying only one section of the NP is not enough to have a full picture of the wrapping since it would not appropriately account for the whole bending of the membrane.

## 5.6 | Steps of drug delivery

NPs are used for many medical purposes. Indeed, they can be used as markers for radiotherapy or even for imaging to diagnose cancers⁸⁷,⁸⁸. They may also be used to deliver drugs as in chemotherapy for instance⁸⁹. For the NPs to achieve their goal, they first need to reach the target cell. However, depending on their physico-chemical properties, these NPs may be cleared by defense mechanisms, as immunity cells, or by the liver or the kidneys for example⁹⁰. Supposing that the NP avoids these clearance processes, it still has to be internalized by the cell, as it has been detailed in this paper. However, the cellular wrapping of the NP is actually the very first step of endocytosis⁶. Indeed, it is not sure that the NP will not be rejected (via exocytosis) before reaching its target⁴⁸. Finally, the NP being a foreign object in the body, it remains essential to be aware of the potential hazards and toxicity it could cause⁹¹.

# 6 | CONCLUSIONS AND PERSPECTIVES

This work presented a 2D model of the cellular uptake of a rigid elliptic NP at the scale of the NP. Investigations have been conducted to evaluate the influence of the mechanical properties of the cell membrane coupled with the aspect ratio of the particle
on the cellular uptake of NP. Afterwards, a sensitivity analysis highlighted that the wrapping degree is twice more sensitive to the NP's aspect ratio than to the adhesion, which itself influences twice more than the membrane tension. Thus, slightly vertical NPs with large adhesion are the ones that are the most likely to be fully wrapped. Whenever these two features are difficult to be implemented in the manufacturing process of the NPs, we recommend manufacturers to focus mainly on the NP's aspect ratio since it is, according to our study, the parameter that influences the most the predictions of endocytosis.

As suggestions for future research, the influence of the deformability of the particle and the variation of the mechanical properties of the cell during wrapping should be investigated. One could also model the bending rigidity of the membrane as a random field in order to take into account the membrane's heterogeneity.

## References

1. Yi X, Shi X, Gao H. Cellular uptake of elastic nanoparticles. *Physical review letters* 2011; 107(9): 098101.

2. Kreuter J. Nanoparticles—a historical perspective. *International journal of pharmaceutics* 2007; 331(1): 1–10.

3. Davis ME, Chen Z, Shin DM. Nanoparticle therapeutics: an emerging treatment modality for cancer. *Nanoscience and technology: A collection of reviews from nature journals* 2010: 239–250.

4. Zhang S, Gao H, Bao G. Physical principles of nanoparticle cellular endocytosis. *ACS nano* 2015; 9(9): 8655–8671.

5. Briolay T, Petithomme T, Fouet M, Nguyen-Pham N, Blanquart C, Boisgerault N. Delivery of cancer therapies by synthetic and bio-inspired nanovectors. *Molecular Cancer* 2021; 20(1): 1–24.

6. Behzadi S, Serpooshan V, Tao W, et al. Cellular uptake of nanoparticles: journey inside the cell. *Chemical Society Reviews* 2017; 46(14): 4218–4244.

7. Sangtani A, Nag O, Field L, Breger J, Delehanty J. Multifunctional nanoparticle composites: progress in the use of soft and hard nanoparticles for drug delivery and imaging. *Wiley Interdisciplinary Reviews: Nanomedicine and Nanobiotechnology* 2017; 9(6): e1466.

8. Boucard J, Linot C, Blondy T, et al. Small Molecule-Based Fluorescent Organic Nanoassemblies with Strong Hydrogen Bonding Networks for Fine Tuning and Monitoring Drug Delivery in Cancer Cells. *Small* 2018; 14(38): 1802307.

9. Shi X, Dem Bussche vA, Hurt RH, Kane AB, Gao H. Cell entry of one-dimensional nanomaterials occurs by tip recognition and rotation. *Nature nanotechnology* 2011; 6(11): 714–719.

10. Chithrani B, Ghazani A, Chan W. Determining the size and shape dependence of gold nanoparticle uptake into mammalian cells. *Nano letters* 2006; 6(4): 662–668.

11. Florez L, Herrmann C, Cramer JM, et al. How shape influences uptake: interactions of anisotropic polymer nanoparticles and human mesenchymal stem cells. *Small* 2012; 8(14): 2222–2230.

12. Yoo JW, Mitragotri S. Polymer particles that switch shape in response to a stimulus. *Proceedings of the National Academy of Sciences* 2010; 107(25): 11205–11210.

13. Champion JA, Katare YK, Mitragotri S. Making polymeric micro-and nanoparticles of complex shapes. *Proceedings of the National Academy of Sciences* 2007; 104(29): 11901–11904.

14. Champion JA, Mitragotri S. Role of target geometry in phagocytosis. *Proceedings of the National Academy of Sciences* 2006; 103(13): 4930–4934.

15. Yi X, Shi X, Gao H. A universal law for cell uptake of one-dimensional nanomaterials. *Nano letters* 2014; 14(2): 1049–1055.

16. Das S, Du Q. Adhesion of vesicles to curved substrates. *Physical Review E* 2008; 77(1): 011907.

17. Bahrami AH. Orientational changes and impaired internalization of ellipsoidal nanoparticles by vesicle membranes. *Soft Matter* 2013; 9(36): 8642–8646.
18. Yi X, Gao H. Phase diagrams and morphological evolution in wrapping of rod-shaped elastic nanoparticles by cell membrane: a two-dimensional study. *Physical Review E* 2014; 89(6): 062712.

19. Vácha R, Martinez-Veracoechea F, Frenkel D. Receptor-mediated endocytosis of nanoparticles of various shapes. *Nano letters* 2011; 11(12): 5391-5395.

20. Huang C, Zhang Y, Yuan H, Gao H, Zhang S. Role of nanoparticle geometry in endocytosis: laying down to stand up. *Nano letters* 2013; 13(9): 4546-4550.

21. Ding H, Ma Y. Theoretical and computational investigations of nanoparticle-biomembrane interactions in cellular delivery. *Small* 2015; 11(9-10): 1055-1071.

22. Yuan X, Wang Y. Adhesion of carbon nanotubes on elastic substrates with finite thickness. *Journal of Applied Physics* 2018; 124(15): 155306.

23. Gupta R, Rai B. Effect of size and surface charge of gold nanoparticles on their skin permeability: a molecular dynamics study. *Scientific Reports* 2017; 7(1): 1-13.

24. Sahay G, Alakhova Y, Kabanov A. Endocytosis of nanomedicines. *Journal of controlled release* 2010; 145(3): 182-195.

25. Sahay G, Kim J, Kabanov A, Bronich T. The exploitation of differential endocytic pathways in normal and tumor cells in the selective targeting of nanoparticulate chemotherapeutic agents. *Biomaterials* 2010; 31(5): 923-933.

26. Abercrombie M, Ambrose E. The surface properties of cancer cells: a review. *Cancer research* 1962; 22(5 Part 1): 525-548.

27. Freund L, Lin Y. The role of binder mobility in spontaneous adhesive contact and implications for cell adhesion. *Journal of the Mechanics and Physics of Solids* 2004; 52(11): 2455-2472.

28. Shen Z, Ye H, Yi X, Li Y. Membrane wrapping efficiency of elastic nanoparticles during endocytosis: Size and shape matter. *ACS nano* 2018; 13(1): 215-228.

29. Spangler EJ, Upreti S, Laradji M. Partial wrapping and spontaneous endocytosis of spherical nanoparticles by tensionless lipid membranes. *The Journal of chemical physics* 2016; 144(4): 044901.

30. Zhou X, Liu J, Zhang S. Adhesion of a vesicle on an elastic substrate: 2D analysis. *Colloids and Surfaces B: Biointerfaces* 2013; 110: 372-378.

31. Deserno M, Bickel T. Wrapping of a spherical colloid by a fluid membrane. *EPL (Europhysics Letters)* 2003; 62(5): 767.

32. Deserno M. Elastic deformation of a fluid membrane upon colloid binding. *Physical Review E* 2004; 69(3): 031903.

33. Góźdź W. Deformations of lipid vesicles induced by attached spherical particles. *Langmuir* 2007; 23(10): 5665-5669.

34. Döbereiner H, Evans E, Kraus M, Seifert U, Wortis M. Mapping vesicle shapes into the phase diagram: A comparison of experiment and theory. *Physical Review E* 1997; 55(4): 4458.

35. Rebouillat S, Pla F. State of the art manufacturing and engineering of nanocellulose: a review of available data and industrial applications. 2013.

36. Vauthier C, Bouchemal K. Methods for the preparation and manufacture of polymeric nanoparticles. *Pharmaceutical research* 2009; 26(5): 1025-1058.

37. Helfrich W. Elastic properties of lipid bilayers: theory and possible experiments. *Zeitschrift für Naturforschung c* 1973; 28(11-12): 693-703.

38. Bao G, Suresh S. Cell and molecular mechanics of biological materials. *Nature materials* 2003; 2(11): 715-725.

39. Hillaireau H, Couvreur P. Nanocarriers' entry into the cell: relevance to drug delivery. *Cellular and molecular life sciences* 2009; 66(17): 2873-2896.
40. Cooper S. Control and maintenance of mammalian cell size. *BMC cell biology* 2004; 5(1): 1–21.

41. Zhao W, Tian Y, Cai M, et al. Studying the nucleated mammalian cell membrane by single molecule approaches. *PLoS One* 2014; 9(5): e91595.

42. Müller M, Deserno M, Guven J. Balancing torques in membrane-mediated interactions: exact results and numerical illustrations. *Physical Review E* 2007; 76(1): 011921.

43. Dietrich C, Angelova M, Pouligny B. Adhesion of latex spheres to giant phospholipid vesicles: statics and dynamics. *Journal de Physique II* 1997; 7(11): 1651–1682.

44. Yi X, Gao H. Cell membrane wrapping of a spherical thin elastic shell. *Soft Matter* 2015; 11(6): 1107–1115.

45. Dietrich C, Angelova M, Pouligny B. Adhesion of latex spheres to giant phospholipid vesicles: statics and dynamics. *Journal de Physique II* 1997; 7(11): 1651–1682.

46. Gao H, Shi W, Freund LB. Mechanics of receptor-mediated endocytosis. *Proceedings of the National Academy of Sciences* 2005; 102(27): 9469–9474.

47. Singer SJ, Nicolson GL. The fluid mosaic model of the structure of cell membranes. *Science* 1972; 175(4023): 720–731.

48. Doherty GJ, McMahon HT. Mechanisms of endocytosis. *Annual review of biochemistry* 2009; 78: 857–902.

49. Villarino MB. Ramanujan's Perimeter of an Ellipse. *arXiv preprint math/0506384* 2005.

50. Sacks J, Welch WJ, Mitchell TJ, Wynn HP. Design and analysis of computer experiments. *Statistical science* 1989: 409–423.

51. Stein ML. *Interpolation of spatial data: some theory for kriging*. Springer Science & Business Media. 2012.

52. Santner TJ, Williams BJ, Notz WI. *The design and analysis of computer experiments*. Springer Science & Business Media. 2013.

53. Rasmussen CE. Gaussian processes in machine learning. In: Springer. 2004 (pp. 63–71).

54. Ghanem R, Spanos P. *Stochastic finite elements: a spectral approach*. Dover publications. 2003.

55. Choi SK, Grandhi RV, Canfield RA, Pettit CL. Polynomial chaos expansion with Llatin hypercube sampling for estimating response variability. *American Institute of Aeronautics and Astronautics* 2004; 42(6): 1191–1198.

56. Berveiller M, Sudret B, Lemaire M. Stochastic finite element: a non intrusive approach by regression. *European Journal of Computational Mechanics/Revue Européenne de Mécanique Numérique* 2006; 15(1-3): 81–92.

57. Sudret B. Global sensitivity analysis using polynomial chaos expansions. *Reliability Engineering & System Safety* 2008; 93(7): 964–979.

58. Khazaie S, Cottereau R, Clouteau D. Influence of the spatial correlation structure of an elastic random medium on its scattering properties. *Journal of Sound and Vibration* 2016; 370: 132–148.

59. Lataniotis C, Marelli S, Sudret B. Gaussian process modelling using UQLab. *arXiv preprint arXiv:1709.09382* 2017.

60. Khazaie S, Wang X, Komatitsch D, Sagaut P. Uncertainty quantification for acoustic wave propagation in a shallow water environment. *Wave Motion* 2019; 91: 102390.

61. Blatman G, Sudret B. An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis. *Probabilistic Engineering Mechanics* 2010; 25(2): 183–197.

62. Ghanem RG, Spanos PD. Stochastic Finite Element Method: Response Statistics. In: Springer. 1991 (pp. 101–119).

63. Le Maître O, Knio OM. *Spectral methods for uncertainty quantification: with applications to computational fluid dynamics*. Springer Science & Business Media. 2010.
64. Savin E, Faverjon B. Computation of higher-order moments of generalized polynomial chaos expansions. *International Journal for Numerical Methods in Engineering* 2017; 111(12): 1192–1200.

65. Sobol I. Sensitivity estimates for non-linear mathematical models. *Mathematical Modeling and Computational Experiment* 1993; 1(4): 407-414.

66. Saltelli A, Annoni P, Azzini I, Campolongo F, Ratto M, Tarantola S. Variance-based sensitivity analysis of model output. *Design and estimator for the total sensitivity index*. *Computer Physics Communications* 2010; 181(2): 259-270.

67. Homma T, Saltelli A. Importance measures in global sensitivity analysis of nonlinear models. *Reliability Engineering & System Safety* 1996; 52(1): 1–17.

68. Yuan H, Zhang S. Effects of particle size and ligand density on the kinetics of receptor-mediated endocytosis of nanoparticles. *Applied Physics Letters* 2010; 96(3): 033704.

69. Vonarbourg A, Passirani C, Saulnier P, Benoit JP. Parameters influencing the stealthiness of colloidal drug delivery systems. *Biomaterials* 2006; 27(24): 4356–4373.

70. Beningo KA, Wang YL. Fc-receptor-mediated phagocytosis is regulated by mechanical properties of the target. *Journal of cell science* 2002; 115(4): 849–856.

71. Jaynes ET. Information theory and statistical mechanics. *Physical review* 1957; 106(4): 620.

72. El-Mihoub TA, Hopgood AA, Nolle L, Battersby A. Hybrid Genetic Algorithms: A Review. *Engineering Letters* 2006; 13(2): 124–137.

73. Staykova M, Arroyo M, Rahimi M, Stone HA. Confined bilayers passively regulate shape and stress. *Physical review letters* 2013; 110(2): 028101.

74. Yi X, Gao H. Kinetics of receptor-mediated endocytosis of elastic nanoparticles. *Nanoscale* 2017; 9(1): 454–463.

75. Yuan H, Li J, Bao G, Zhang S. Variable nanoparticle-cell adhesion strength regulates cellular uptake. *Physical review letters* 2010; 105(13): 138101.

76. Noguchi H, Takasu M. Adhesion of nanoparticles to vesicles: a Brownian dynamics simulation. *Biophysical journal* 2002; 83(1): 299–308.

77. Nel AE, Mädler L, Velegol D, et al. Understanding biophysicochemical interactions at the nano–bio interface. *Nature materials* 2009; 8(7): 543–557.

78. Kosmalska AJ, Casares L, Elosegui-Artola A, et al. Physical principles of membrane remodelling during cell mechanoadaptation. *Nature communications* 2015; 6(1): 1–11.

79. Decuzzi P, Ferrari M. The receptor-mediated endocytosis of nonspherical particles. *Biophysical journal* 2008; 94(10): 3790–3797.

80. Chen B. Probing the effect of random adhesion energy on receptor-mediated endocytosis with a semistochastic model. *Journal of Applied Mechanics* 2014; 81(8).

81. Nicolson GL. Update of the 1972 Singer-Nicolson fluid-mosaic model of membrane structure. *Discoveries* 2013; 1(1).

82. Yang C, Liu Y, He Y, et al. The use of HA oligosaccharide-loaded nanoparticles to breach the endogenous hyaluronan glycocalyx for breast cancer therapy. *Biomaterials* 2013; 34(28): 6829–6838.

83. Ito S. The enteric surface coat on cat intestinal microvilli. *The Journal of cell biology* 1965; 27(3): 475–491.

84. Allen JA, Halverson-Tamboli RA, Rasenick MM. Lipid raft microdomains and neurotransmitter signalling. *Nature reviews neuroscience* 2007; 8(2): 128–140.
85. Bahrami AH. Orientational changes and impaired internalization of ellipsoidal nanoparticles by vesicle membranes. *Soft Matter* 2013; 9(36): 8642–8646.

86. Dasgupta S, Auth T, Gompper G. Shape and orientation matter for the cellular uptake of nonspherical particles. *Nano letters* 2014; 14(2): 687–693.

87. Hainfeld J, Dilmanian F, Slatkin D, Smilowitz H. Radiotherapy enhancement with gold nanoparticles. *Journal of pharmacy and pharmacology* 2008; 60(8): 977–985.

88. Faucon A, Benhelli-Mokrani H, Córdova L, et al. Are fluorescent organic nanoparticles relevant tools for tracking cancer cells or macrophages?. *Advanced healthcare materials* 2015; 4(17): 2727–2734.

89. Gelperina S, Kisich K, Iseman M, Heifets L. The potential advantages of nanoparticle drug delivery systems in chemotherapy of tuberculosis. *American journal of respiratory and critical care medicine* 2005; 172(12): 1487–1490.

90. Alexis F, P. E, Molnar L, Farokhzad O. Factors affecting the clearance and biodistribution of polymeric nanoparticles. *Molecular pharmaceutics* 2008; 5(4): 505–515.

91. Lin W, Huang Y, Zhou X, Ma Y. In vitro toxicity of silica nanoparticles in human lung cancer cells. *Toxicology and applied pharmacology* 2006; 217(3): 252–259.

## SUPPORTING INFORMATION

The Python code used to display the wrapping configuration, to compute $\overline{\Delta E(f)}$ and to identify the wrapping phases is available at: https://github.com/SarahIaquinta/uptake_of_random_rigid_elliptic.particle