Translocation kinetics of vesicles through narrow pores

The passage of vesicles through a small pore induced by pressure difference or chemical reactions is relevant to biological processes, e.g., transendothelial migration in immune systems [1], transdermal applications [2], or drug delivery [3] in pharmacology. Owing to their flexible surface, lipid bilayer vesicles readily undergo shape deformations enabling them to move through confining tissue. The problem has gained ample attention by researchers and has been treated theoretically both analytically [4,5] within a statistical mechanical approach as well as by means of computer simulations of coarse-grained models [6,7]. Earlier studies have shown [7], that a realistic concentration gradient of osmotically active molecules is suffient to overcome the free energy barrier for liposome translocation. Under the approximations of constant volume and/or surface area of the vesicle as well as assuming sufficiently slow dynamics of the translocation process so that the existence of thermodynamic free energy barrier in equilibrium can be considered [4–9], these studies have provided important insights on the effectiveness of the translocation process for various governing parameters as the driving forces, initial size and elastic moduli of the vesicle as well as geometry and size of the pore. While Shoaei and Muthukumar [4] have considered the effect of various parameters on the translocation time by using the Fokker-Planck formalism to derive the free energy profile for vesicle translocation, the work of Khunpetch et al. [5] employed the Onsager variational principle which relies on the existence of slow variables during the translocation kinetics to study the free energy variation along the actual kinetic path. Yet a number of questions concerning the general dependence of the total translocation time $\tau_{tr}$ on the pore radius R_p , the pressure difference $\Delta P$, or the elasticity of the membrane shell, in particular, in the absence of the volume and area conservation assumption, remain still open. In the present work we report on results based on comprehensive molecular dynamics (MD) computer experiments varying vesicles and pore sizes as well as the bending rigidity of the membrane surface in broad intervals.

We consider polymerized ("tethered") membranes of fixed connectivity comprised of cross-linked vertices as in the case of spectrin network of erythrocyte plasma membranes and of the lamina network of eukariotic nuclear envelopes [10]. The vesicle is modeled by a triangular network of M spherical vertices of diameter σ, fig. 1, and contains a dense gas of athermal beads with diameter σ at volume density $\approx 0.42$. Both vertices and fillers repel one another with a short-range repulsive Weeks-Chandler-Andersen (WCA) (i.e., truncated and shifted Lennard-Jones) type potential [11] between any pair of vertices at distance r,

$$\begin{align} U_{WCA}(r) = \begin{cases} 4 \varepsilon \left[\left(\displaystyle\frac{\sigma}{r}\right)^{12} - \left(\displaystyle\frac{\sigma}{r}\right)^6\right] + \varepsilon, & r \leq 2^{1/6}\sigma,\\ 0, & r > 2^{1/6}\sigma, \end{cases} \end{align} \tag{ 1 }$$

choosing the potential strength $\varepsilon = 1$ and putting the thermal energy k_B T to unity as well. The vertices are connected by flexible bonds, each being represented by the finitely extensible nonlinear elastic (FENE) potential [12]:

$$\begin{align} U_{\textit{FENE}}(r) = \begin{cases} -\displaystyle\frac{k r_0^2}{2} \ln\left(1-\displaystyle\frac{r^2}{r_0^2}\right), & r < r_0,\\ \infty, & r \geq r_0. \end{cases} \end{align} \tag{ 2 }$$

Zoom In Zoom Out Reset image size

The constants are chosen as $r_0 = 1.5\sigma, k = 30 \epsilon / \sigma^2$. Note that eqs. (1), (2) together create a potential for the bond length $\ell_b$ with a sharp minimum at $\ell_b = 0.97 \sigma$. The chosen interactions, eqs. (1), (2), and the values of the parameters $\sigma,\;r_0$ and ensure no leakage of fillers through holes in the membrane shell during the translocation process.

The membrane bending energy may be cast in the convenient form [13],

$$\begin{equation} {\cal H}_{bend} = \kappa_{imp} \sum{i,j} [1 - \cos(\phi_{ij} - \phi_0)], \end{equation} \tag{ 3 }$$

where $\phi_{ij} = \arccos(\hat{\textbf{n}}_i \cdot \hat{\textbf{n}}_j)$ is the improper dihedral angle between the normals $\hat{\textbf{n}}_i$ and $\hat{\textbf{n}}_j$ of each pair of adjacent triangles (i, j) sharing a common edge and the hat denotes a unit vector. Spontaneous spherical curvature of the vesicle is induced by the offset $\phi_0 = 0$. As shown in fig. 2, there exists an optimal value of bending rigidity, $5 \le \kappa_{imp} \le 10$, which ensures the fastest translocation process almost regardless of the pore diameter. Interestingly, for an optimally fast passage through the pore vesicle rigidity becomes less essential with growing membrane size M since for fixed $\kappa_{imp}$ larger shells deform more easily just like semiflexible polymer chains with contour length exceeding significantly their persistence length. The ruggedness of the stiff pore rim, on the other hand, leads to increased friction for too soft smaller membranes which need time to unhook from local bumps in the rim.

Zoom In Zoom Out Reset image size

The parameter $\kappa_{imp}$ may be related to the Young elastic modulus $Y_0 = 2\kappa_{imp}/\sqrt{3}$ [13] and the bending rigidity $\kappa_0 = \sqrt{3} \kappa_{imp}/2$ [14] of the 2D surface material. Recently, Tortora and Jost [15] have demonstrated that these relationships hold remarkably well for triangulated membrane shells even in the presence of topological Stone-Wales 5/7-disclination defects, inevitable on curved surfaces.

We employ the Matlab S²-Sampling Toolbox [16] to generate a nearly uniform triangular tessellation of spheres with $M=2400,\,4800,\,9800,\, 19600,\, 39200$ vertices and mean gyration radii, $\langle R_g^2\rangle^{1/2} / \sigma = 12.45$, 17.58, 25.75, 36.82, 51.73, respectively.

The number of 5/7-defects, n_SW , is found to scale as $n_{SW} = 0.138 M^{0.81}$ with the membrane size M so that, e.g., for M = 9800 one has $n_{SW} \approx 250,$ cf. fig. 1. Each vesicle with M vertices is comprised of $3M-6$ bonds, $3M-6$ dihedrals, and $2M - 4$ angles within the triangles.

In all studies the translocating membrane goes through a cylindrical pore of radius R_p and length $L=4\sigma$ in a closely packed solid wall made of discrete beads of diameter σ. The wall separates the donor compartment, where the vesicle is exposed to pressure P due to highly concentrated solution of beads, from the receiver compartment which is initially empty so that the initial pressure gradient is $\Delta P$, see fig. 1. The translocation time $\tau_{tr}$ is measured from the moment the vesicle touches the pore entrance until it exits the orifice into the received compartment. Our integration step was chosen as $\delta = 0.005 \tau_{MD}$ where the MD time unit, or "mdu" in the following, $\tau_{MD} = \sqrt{(m\sigma^2 / \epsilon)}$ when we choose the vertex mass m = 1 as well. The MD simulations are carried out at constant temperature T with a Nosé-Hoover thermostat using the HOOMD - Blue software on graphics processing units (GPUs) [17].

How do the units chosen in our coarse-grained model relate to those used in typical experiments? Taking our unit length $\sigma \approx 3\ \text{nm}$, the vesicle size for M = 19600 would be of the order of $\approx 200\ \text{nm}$ which is the typical diameter of a DOPC [C₄₄H₈₄NO₈P] liposome GUV, cf. [18]. If we ascribe a mass of 62000 Da to each of our 19600 vertices so as to match that of the DOPC GUV with an energy unit $\epsilon = k_BT = 4.14^{-21}\ \text{J}$ at $T=300\ \text{K}$, then our unit of time (mdu) would be of the order of $10^{-8}\ \text{s}$ which yields typical translocation times $\approx 10^{-6} \div 10^{-5}\ \text{s}$ at pressures $\Delta P = k_BT/\sigma^3 \approx 1.5 \times 10^5\ \text{Pa}= 1.5\ \text{atm}$, cf. [19].

Figure 3 illustrates the changes which the vesicle shape undergoes during the translocation in terms of the eigenvalues e₁, e₂, e₃ of the mean squared gyration radius tensor, $\langle R_g^2\rangle$, with elapsed time since the onset of the process. Due to the circular shape of the pore, the vesicle retains spindle-like axial symmetry, $e_2(t) \approx e_3(t)$, in the direction of motion perpendicular to the confining wall while the largest eigenvalue, e₁(t), undergoes strong variations that differ qualitatively for a narrow, fig. 3(a), and a wider, fig. 3(b), pore.

Zoom In Zoom Out Reset image size

The initial impact of the vesicle on the narrow pore leads to an immediate strong flattening at minimal translocated fraction, $f \approx 0.12$, followed by a series of jerks reflected by shape oscillations while the membrane remains partially stuck in the hole and gradually protrudes forward. During this process the translocated fraction changes negligibly, $f \le 0.25$, indicated by a long plateau until the very end of the process when an almost abrupt transition into the free compartment occurs. In contrast, for the wider pore, fig. 3(b), the initially equal eigenvalues $e_1 \approx e_2 \approx e_3$ split gradually upon entering the pore whereby e₁ goes through two successive maxima at $f=0.25$ and $f=0.78$, separated by a minimum at $f=0.39$, before all three eigenvalues join again upon exiting the pore. In both cases one observes some relatively weak variation of the surface area, A, and the volume, V, in the course of the translocation process. Thus, the width of the yellow-shaded regions in fig. 3, spanning the difference of e₁ and e₂, e₃, displays the asphericity of the protruding membrane at different stages of the passage through the pore, as reflected by the translocated fraction f. We interpret the oscillations of the overall shape of the vesicle as soon as it gets partially stuck in the narrow pore as being caused by a "to and fro" motion of the embedded fillers which, owing to their initial momentum, suffer several collisions with the front and rear part of the membrane shell, pushing thereby the whole vesicle forward. Figure 3(a) indicates that the translocation process takes nearly an order of magnitude longer as compared to the case of wider pore, fig. 3(b), where the vesicle gets briefly arrested when already about half of it has crossed the pore. Then the enclosed fillers overtake and hit the front part of the shell, reflect subsequently from the rear part and hit again, pushing thus the whole vesicle ultimately into the receiver compartment. Concomitantly, one observes some changes in the membrane area A and volume V, and most notably, a slight decline in A and an increase in V at the very end of the process when the vesicle exits the pore and attains again its spherical (i.e., maximal volume, minimal area) shape yet with somewhat larger radius R_g due to the negligible pressure in the receiver compartment.

One might assume that all the aforementioned transient changes are also accompanied by an ongoing variation of the energy of the vesicle whereby an activation barrier is built up from the different energy contribution mentioned in the description of the model. In fig. 4 we demonstrate for the case, shown in fig. 3(b), that this barrier is mainly formed by the elastic energy of the membrane bonds, $U_{\textit{FENE}}$, cf. eq. (2), and by the short-range repulsion of the fillers inside the membrane, eq. (1), while the inputs from vertices repulsion, U_mem , and from surface bending, U_imp , remain negligible during the entire translocation process. The maximum of U_tot is attained at $f\approx 0.60$ which occurs shortly before the second largest deformation of the vesicle membrane takes place, see fig. 3(b). It should be noted that the shape of this energetic barrier resembles strongly that of the analytically predicted Helfrich free barrier [4], or the one derived from Monte Carlo simulations [7] albeit we do not assume here the conservation of area and volume during translocation or a semi-spherical shape of the caps on both sides of the pore, as in the theory.

Zoom In Zoom Out Reset image size

The main focus in this investigation of translocation kinetics is on the dependence of mean translocation time, $\tau_{tr}$, on pore size R_p , and in fig. 5(a) this is shown in log-normal coordinates for a number of vesicles of different size, M, and for two values of the applied pressure, $\Delta P=7.01$ and $\Delta P=10.0$. Evidently, in all cases $\tau_{tr}$ is found to diverge as the pore radius R_p approaches some critical value, R_cr (P, M), for which the vesicle gets stuck in the hole. The critical size R_cr (P, M) naturally increases with growing vesicle size M, and decreases with increasing applied pressure P. As shown in the inset to fig. 5(a), however, with growing size M the vesicles manage to penetrate more and more narrow pores in comparison to their own radius R_g . We attribute this behavior to the observed easier formation of wrinkles on the gradually softening surface and to a decrease of the effective roughness of the pore rim as the membrane size increases.

Zoom In Zoom Out Reset image size

We find that all results, shown in fig. 5(a), can be plotted as a universal master curve in log-log coordinates, see fig. 5(b), whereby the following power law holds: $\tau_{tr} \propto (R_p-R_{cr}(\Delta P,M))^{-\rho}$ with an exponent $\rho = 2/3$. Regarding the power law divergence of $\tau_{tr}$, one should make here the caveat that the value $\rho=2/3$ emerges as a best fit from our data while the error in determination of $\tau_{tr}$ when $R_p \to R_{cr}(P,M)$ (exponentially) rapidly grows due to limited statistics of sampling as manifested by the correspondingly increased scatter in the values in fig. 5(b).

The effect of the driving force on $\tau_{tr}$ is displayed on fig. 6(a) where we plot $\tau_{tr}$ as function of the pressure gradient $P - P_{cr}(R_p)$ for a vesicle with M and several pore sizes R_p . As critical pressure P_cr (R_p ) we use the least pressure still capable of pushing the vesicle through a pore with radius R_p . Apparently, one finds again a power law relationship, $\tau_{tr} \propto (P - P_{cr}(R_p))^{-\omega}$, however, the scaling exponent ω is found to vary with the pore radius, R_p . For comparatively wide pores, which are rapidly passed by the vesicle, the value of ω is $\approx 0.22$ whereas as the pore gets narrower, ω increases, tending to $\omega \approx 0.85 \pm 0.05$. In fig. 6(b) we show the variation of P_cr (R_p ) with the dimensionless relative size of the pore, $R_p/R_f$, where R_f denotes the smallest pore radius still penetrable by a vesicle of size M. It is seen that P_cr (R_p ) increases when $R_p \to R_f$, yet one should note that for relatively wide pores P_cr may become so small that the vesicle fails to enter the pore. The inset in fig. 6(b) indicates the gradual increase of the exponent ω when $R_p/R_f \to 1$. It should be noted that the observed power law variation of $\tau_{tr}$ with $\Delta P$ is at variance with the exponential decay, predicted in [4], while it qualitatively agrees with the results of Khunpetch et al. [5], who obtained an exponent $\omega =1$.

Zoom In Zoom Out Reset image size

Discussion

Support by COST CA17139 (European Cooperation in Science and Technology (see http://www.cost.eu and https://www.fni.bg) and its Bulgarian partner FNI/MON under KOST-11) is gratefully acknowledged. BR acknowledges the financial support from the Centre of Excellence "National center for mechatronics and clean technologies" (Project BG05M2OP001–1.001–0008-C010) supported by the European RD Fund within the Bulgarian OP "Science and Education for Smart Growth 2014–2020".

Data availability statement: All data that support the findings of this study are included within the article (and any supplementary files).