Wrinkling instability of vesicles in steady linear flow

Wrinkling is a well-known phenomenon frequently observed in nature and everyday life on widely different spatial scales [1,2]. It bridges mechanics, geometry, physics, and biology. It usually occurs either in extensible materials due to external stretching, or as a result of compression of inextensible materials, in which case it reduces to the Euler buckling instability [3]. In contrast to the wrinkling instability in films and structures such as micro-capsules, where it occurs due to stretching elasticity [4], in vesicles with unstretchable lipid membranes it takes place under compression due to negative surface tension [5,6].

Vesicles are fluid droplets encapsulated by a phospholipid-bilayer membrane and suspended in either the same or another fluid [7]. They are considered as simple model objects to simulate red-blood-cell dynamics in a first step towards understanding blood rheology. Two constraints, namely the conservation of both vesicle volume V and surface area A, result in rather involved non-linear dynamics in flow (for a review, see [8]). In transient plane elongation (hyperbolic) flow, the wrinkling instability happening during the relaxation of a vesicle towards a new stationary state was studied experimentally [5]. There, the flow direction was reversed faster than the vesicle's relaxation time scale $\tau\sim\eta_{\text{out}} R_0^3/\kappa$, where κ is the bending elasticity, $\eta_{\text{out}}$ outer fluid viscosity, and $R_0=(3V/4\pi)^{1/3}$ the effective vesicle radius defined via the volume V. This sudden reversal causes a switching from vesicle stretching to vesicle compression, see fig. 1, top. The next reversal then takes place after full relaxation. For strain rates χ lower than a critical value $\chi_c$, this flow reversal leads to an almost smooth transition from an elliptical vesicle with long axis in the direction of compression to one with long axis perpendicular to it. During this transient transformation, only deformations with wave number k = 2 are induced by the flow. Higher-order deformations are only excited thermally, and the deformation power spectrum $P_k\sim k^{-4}$ is observed since bending deformations enter as $\kappa k^4$ in the membrane energy [7]. For $\chi>\chi_c$, higher-order deformation modes, coined wrinkles, are generated on top of an almost elliptical vesicle. A spectral analysis of the shape distortions easily identifies the wrinkling onset due to a sharp growth of higher-order harmonics (see fig. 3 in [5]).

Zoom In Zoom Out Reset image size

In this letter, we explore another mechanism to observe wrinkling caused by the interplay of thermal fluctuations with the trembling (TR) motion. The TR regime of vesicle dynamics in steady linear flow is an intermediate regime between tank-treading (TT) and tumbling (TU) [9]. It is characterized by quasi-periodic noisy oscillations of the vesicle main axis around the shear direction, i.e., oscillations of the inclination angle θ around zero. The lower images of fig. 1 illustrate schematically the flow inversion and vesicle compression in linear flow. The flow direction makes an angle of $\pi/4$ with the direction of the extensional part of the flow. At negative inclination angles $\theta<0$, a reversal of the elongation part of the flow occurs in the vesicle frame of reference. Then, the vesicle is compressed for some part of its oscillation period, leading to wrinkling at sufficiently large values of the normalized strain rate $S = 14\pi s\eta_{\text{out}}R_0^3/(3\sqrt{3}\kappa\Delta)$ [10], where s is the strain rate of the flow, A the vesicle area, and $\Delta=A/R_0^2-4\pi$ the excess area.

In spite of the fact that a possible connection between vesicle shape deformations and a wrinkling-like instability in TR was discussed already in refs. [11–13], only the recent observation of the amplification of thermal noise by TR dynamics [14,15] provides a real basis to this assumption. Indeed, in spite of a small ratio of thermal energy to bending rigidity $k_BT/\kappa \approx 0.05$, where k_B is the Boltzmann constant and T the temperature, the dramatic effect of thermal noise on vesicle dynamics subjected to a linear flow was recently observed in experiment [14]. Numerical simulations including thermal fluctuations reproduce the experimental observations on a quantitative level [15]. Both works such as previous stochastic simulations [16,17] show that thermal noise is the main factor for the strong discrepancies in vesicle dynamics found for the last several years between experiment and various numerical simulations based on deterministic models. The most prominent discrepancies are the difference in vesicle dynamics between TR in experiment [9] and vacillating-breathing (VB) in deterministic simulations [18], and the number of the control parameters to present the phase diagram of the vesicle dynamical states in linear flow [11,13,19]. The experimental results show that two parameters, which represent the competition between relevant time scales, are enough to describe the transition between TT, TR, and TU, namely S and $\Lambda=4(1+23\lambda/32)\sqrt{\Delta}\omega/s$, where $\lambda=\eta_{\text{in}}/\eta_{\text{out}}$ is the viscosity contrast ($\lambda = 1$ in this article) and ω the flow vorticity [10].

The inclusion of thermal noise in the model leads to strongly distorted vesicle shapes due to appearance of higher-order harmonics, in particularly odd ones, which become larger than the second harmonic during the compression at negative θ (figs. 2 in [14,15]). As a result, at large values of S a wrinkling instability similar to that observed in a transient elongation flow could show up [5]. Since the TR period is not well defined and vesicle compression and stretching occur continuously during each cycle, the vesicle shape perturbations can propagate between each flow reversal into the next cycle, additionally contributing to the irregularities of the TR motion and the wrinkling instability. Thus, the quantitative analysis of wrinkling in TR is much less straightforward than in elongation flow, since wrinkles appear on the top of several even and odd modes with sufficiently large amplitudes that continuously vary during a cycle. Furthermore, the power spectrum below the wrinkling instability is very different from the thermal noise spectrum. All these factors hinder an unambiguous identification of the wrinkling onset and scaling above it in TR.

In this letter we present experimental results that clarify the following questions: i) Is the wrinkling instability possible in TR? ii) What is the control parameter? iii) Is it possible to quantitatively characterize the wrinkling transition and scaling above it in the TR regime? Moreover, we present analytical and numerical results which reproduce the main characteristics of the wrinkling instability and show that thermal noise plays a crucial role.

In order to quantitatively investigate the wrinkling instability and the vesicle shape deformations in TR in a wide range of $S\in[10;2400]$, which was chosen as the control parameter, $\Lambda\in[1.3;2.42]$, and $\Delta\in[0.43;1.65]$, we modified our experimental setup (see, e.g., [14]), making it suitable to higher-viscosity solvents by decreasing the resistance of the inlets, and increased the accessible flow rates. As a result, S was increased by more than one order of magnitude while the Reynolds number remained of the order of $10^{-4}$.

We observed 14 vesicles, using $\eta_{\text{out}} = \eta_w$, $5.3\eta_w$, and $10.66\eta_w$, where $\eta_w=1 \text{mPa s}$ is the viscosity of water. Observations and measurements of vesicle dynamics were conducted inside a micro-fluidic 4-roll mill, implemented in silicone elastomer (Sylgard 184, Dow Corning), fabricated via soft lithography [12,14,20]. The key component of this device is a dynamical trap, that allows long time observation of vesicles in a planar linear flow (observation times were long compared to the time scales of the flow). The flow inside the trap is generated using hydrostatic pressure. The physical control parameter of the flow inside the trap, namely the ratio of vorticity to strain rate $\omega/s$ inside the trap, is controlled by the single experimental control parameter, the pressure drop across the dynamical trap that can be varied continuously. The flow rates were measured by particle tracking velocimetry [21]. Vesicles were prepared in water, $38{\%}w/w$ sucrose/water and $60{\%}w/w$ glycerol/water solutions (1,5.3 and $10.66 \text{mPa s}$, respectively, $\lambda = 1$), using the electroformation method [22], as described in [14]. The dynamics of the vesicles was monitored using inverted fluorescence microscope (IMT-2, Olympus). The images of the vesicle at its largest cross-section were collected with a Prosilica EC1380 CCD camera, aligned with the shear axis and synchronized with a mechanical chopper on the path of the excitation beam to reduce unnecessary exposure time. The frame rate/shutter time were adjusted to match the characteristic time scales that vary with the flow rates. Further experimental details were reported elsewhere [14].

The relative position of the membrane $r(\phi,t)$ was determined in the frame of reference of the vesicle using intensity variations along the radial directions. For each image we obtained up to 500 discrete positions sampled along the contour. To quantitatively analyze the TR dynamics the dimensionless shape of the vesicle $u(\phi,t)=r(\phi,t)/R(t)-1$, $0\leq\phi\leq2\pi$, $R(t)=\langle r(\phi,t)\rangle_\phi$, was Fourier-decomposed via fast Fourier transform, i.e. $u(\phi,t)$ was expressed as $u(\phi,t)=\sum_k{u_k(t)e^{-ik\phi}}$. Then the vesicle shape perturbations were quantified via a mean wave number $k^*=\sqrt{\sum_{k=3}^{16}{k^2|u_k|^2}/\sum_k{|u_k|^2}}$, similarly to the analysis done in transient elongation flow [5].

Figure 2 presents snapshots of two time series of TR at two values of $S= 421$ and $S = 1147$, much higher than reported before. One can easily identify, besides the previously observed concavities, wrinkles with significantly higher wave number modes. On the right-hand side, the time series of the second mode $|u_2|$ and $k^*$ corresponding to the same time intervals show their anticorrelation —the attenuation of $|u_2|$ due to compression is compensated by the amplification of higher-order modes. However, during much longer time intervals $k^*$ exhibits irregular temporal dynamics, weakly correlated with the dynamics of the second mode. Therefore, we defined $k^*_{\text{even}}$ and $k^*_{\text{odd}}$, where summation over the even and the odd modes in the expression for $k^*$ was applied. We found that $k^*_{\text{even}}$ shows quasi-regular dynamics with clearly visible peaks anticorrelated with $|u_4|$ and with about the same period, which is about half that of $|u_2|$, whereas the dynamics of $k^*_{\text{odd}}$ is rather erratic (fig. 3). This observation can be explained by the fact that the symmetric even modes are associated with the dynamics, whereas the odd modes are associated with the amplification of thermal noise, as discussed in [14].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Another way to analyze the shape perturbations is to present the probability distribution function (PDF) of $k^*$ and time averaged power spectra of the vesicle shape at different S. It is obvious from both plots in fig. 4 that the higher S, the higher the probability to excite higher harmonics and the larger the spectral amplitude of the higher-order modes. The quantity $k^*$ is a reliable quantitative indication of the wrinkle generation at larger S. Due to strong interference of the amplified thermal noise with the wrinkling, we use further a probabilistic approach to characterize the wrinkling onset. We calculate the probabilities of $k^*$ to be above a certain mode number n as a function of S that is calculated as the ratio of the PDF area at $k^*>n$ normalized by the total PDF area according to ${\cal P}(n,S)\equiv P(k^*>n,S)/P(3<k^*<16,S)$. As one can see in fig. 5, ${\cal P}(n,S)$ increases up to four times between $S\simeq 80$ and $S\simeq 150$ for n = 5 and n = 7 for all Λ and Δ used in the experiment.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

An approach similar to one used in [5,11] can be also applied. In this case, we modify the control parameter to be $S\rightarrow S_{\text{eff}} = S\cos(\theta_{\textit{rms}}+\pi/4)$ that is rotated by $\pi/4$ with respect to the compression axis (see fig. 1). Then the mean value $\langle k^*_{\text{even}}\rangle$ shows a transition as a function of $S_{\text{eff}}$ independent of Λ and Δ used in the experiment. It scales as $S_{\text{eff}}^{0.10}$ in the interval $S_{\text{eff}}\in[40;120]$, and is independent of $S_{\text{eff}}$ outside of this interval (fig. 6). This behavior is in qualitative agreement with the behavior of $k^*$ as a function of the normalized strain in the transient elongation flow reported in [5].

Zoom In Zoom Out Reset image size

From a theoretical perspective, the main mechanism triggering the wrinkling instability can be understood from the model for quasi-spherical vesicles in flow presented in [23]. The three-dimensional radius $r(\theta,\phi)$ of a vesicle is expanded in spherical harmonics ${\cal Y}_{lm}$ as

$$\begin{equation} r(\theta,\phi)=R_0(1+\sqrt{\Delta}\sum u_{lm}{\cal Y}_{lm}(\theta,\phi))\equiv R_0(1+u), \end{equation} \tag{ 1 }$$

where the sum goes over l > 1 and $|m|\leq l$, Δ is the excess area, and R₀ the effective radius. This vesicle is subject to the linear flow

$$\begin{equation} \textbf{v}^{\infty} = \omega(y\textbf{e}_{\textbf{x}}-x\textbf{e}_{\textbf{y}})+s(y\textbf{e}_{\textbf{x}}+x\textbf{e}_{\textbf{y}}) \end{equation} \tag{ 2 }$$

of vorticity ω and strain rate s. Then the equation of motion for the single modes u_lm is

$$\begin{equation} \tau \partial_t u_{lm}=i\frac{S\Lambda}{2}u_{lm}-\mu_l(l(l+1)+\sigma)u_{lm}- i\frac{S}{2}f_{lm} + \zeta_{lm} \end{equation} \tag{ 3 }$$

with $\tau\equiv(385/72)\sqrt{\pi/10\Delta}\eta_{\text{out}}R_0^3/\kappa$, $f_{lm}\equiv \pm\delta_{l,2}\delta_{m,\pm2}$, and $\mu_l\equiv(l-1)l(l+1)(l+2)\tau\kappa/((4l^3+6l^2-1)\eta_oR_0^3)$. The homogeneous surface tension σ is a dynamical parameter which is adjusted such that the vesicle area is kept constant. The correlations of the Gaussian white noise $\zeta_{lm}$ are given by the fluctuation-dissipation theorem in equilibrium.

We first neglect thermal noise and decompose the expansion coefficients as $u_{lm} \equiv r_{lm}\exp(-im\theta_{lm})$. We obtain from (3) the equation of motion

$$\begin{equation} \partial_t r_{lm}=-\mu_l(l(l+1)+\sigma)r_{lm}+ \frac{S}{2}r_{22}\sin(2\theta_{22})f_{lm} \end{equation} \tag{ 4 }$$

for the amplitudes r_lm. For $(l,m)\neq(2,\pm 2)$, $f_{lm}=0$ and eq. (4) describes an exponential decay or growth, depending on the sign of $l(l+1)+\sigma$. The surface tension σ can be obtained from the condition that the excess area Δ remains constant over time, i.e., that the equality

$$\begin{equation} 1=\sum_{l,m}\frac{(l-1)(l+2)}{2}r_{lm}^2 \end{equation} \tag{ 5 }$$

is always obeyed [23]. Taking the time derivative of eq. (5) leads to

$$\begin{eqnarray} &\hspace{-1pc}l(l+1)+\sigma=\nonumber\\ &\hspace{-1pc}\frac{4Sr_{22}\sin(2\theta_{22})+\sum_{l^{\prime},m^{\prime}}[l(l+1)-l^{\prime}(l^{\prime}+1)] \nu_{l^{\prime}}r_{l^{\prime}m^{\prime}}^2}{\sum_{l^{\prime},m^{\prime}}\nu_{l^{\prime}}r_{l^{\prime}m^{\prime}}^2}\quad \end{eqnarray} \tag{ 6 }$$

with $\nu_l\equiv(l-1)(l+2)\mu_l$. Whether the mode $(l,m)$ decays or grows depends on the sign of (6). We first see that positive contributions come from all modes with $l^{\prime} < l$: the excess area is transferred from high-order modes to low-order ones. Moreover, the external flow has a contribution which has the sign of the angle $\theta_{22}$, which is approximately the inclination angle, and a strength proportional to S. In TR and TU motions, the flow compresses the vesicle when the inclination angle takes negative values, leading to negative surface tension and amplification of the modes with $l(l+1)+\sigma < 0$. With increasing S, higher-order modes will be excited during compression phases since the surface tension will take larger negative values. However, harmonics with l > 2 obeying eq. (4) will decay to 0 in the long run due to the stretching at positive inclination angles and the transfer of excess area from high to low-order harmonics.

In eq. (3), the term $iSf_{lm}/2$ is derived from the quantities $X^{\infty}\equiv\textbf{v}^{\infty}(r=R_0)\cdot\textbf{e}_{\textbf{r}}$ and $Y^{\infty}\equiv-R_0\nabla\cdot\textbf{v}^{\infty}(r=R_0)$ which describe the flow at the virtual sphere $r=R_0$ [23]. If we instead define $X^{\infty}$ and $Y^{\infty}$ at $r=R_0(1+u)$, correction terms of the form Su_lm will induce a coupling between $u_{2,\pm2}$ and higher-order modes with even l and m. However, modes with odd l and/or m will still vanish. Therefore, we have to take thermal noise into account to explain the presence of odd modes in experiments. In the following, we present simulation results of eq. (3) with thermal noise ($T=293\text{K}$, $\kappa=25k_{B}T$) and the mentioned additional flow correction terms. To enforce the constraint (5), we replace σ by a quadratic energy term with large spring constant [15]. We simulate trajectories of several thousand periods such that initial transient effects vanish. In order to compare the simulations with the experiments, the shape in the middle flow plane is expanded in Fourier coefficients u_k like in the experiments.

All simulations with the same physical parameters as in fig. 6 predict a TR motion in the steady state, in agreement with experiments. For instance, fig. 7 shows a simulated trajectory for S = 1147, $\Lambda=1.78$, and $\Delta=1.27$ corresponding to fig. 2, bottom. The first panel shows the almost regular oscillation of the inclination angle θ around 0. The amplitude of the second mode $|u_2|$ oscillates with the same period, which is comparable to the experimental one (see fig. 2), although it lags slightly behind. The surface tension σ shown in the second panel, which corresponds to eq. (6) with flow correction terms, is strongly correlated with θ, reaching minima and maxima almost at the same times. This strong correlation is the basis of the wrinkling phenomenon. The third and fourth panels show $k^*_{\text{even}}$ and $k^*_{\text{odd}}$, respectively. As in the experimental data of fig. 3, $k^*_{\text{odd}}$ exhibits a rather erratic evolution while $k^*_{\text{even}}$ is regular since even modes are coupled by the flow. Moreover, like in the experiments, $k^*_{\text{even}}$ has a period approximately equal to half the period of θ. This fact is due to thermal fluctuations, as explained below.

Zoom In Zoom Out Reset image size

Figure 8 compares the evolution of the zeroth, second, and fourth modes with that of $k^*_{\text{even}}$. The top panel shows the deterministic case. It corresponds to a vacillating-breathing motion [18]: when the inclination angle is negative, the amplitude of the second mode $|u_2|$ decreases while $|u_0|$ and modes out of plane increase (see also data in ref. [13]). The dynamics of all higher-order even modes, e.g., $|u_4|$, is perfectly synchronized with $|u_2|$, as is the evolution of $k^*_{\text{even}}$. In addition, $k^*_{\text{even}}$ varies between 4.002 and 4.005, meaning that the fourth mode is always much larger than the higher ones. With thermal fluctuations, the picture changes dramatically. The zeroth and second modes do not show qualitative changes, but $|u_4|$ has a totally different behavior. It now exhibits a maximum not only when $|u_2|$ is maximal but also when it is minimal. This new peak is due to the amplification of thermal noise at negative inclination angles [14,15] and corresponds to the observed wrinkles. Moreover, the mean wave number $k^*_{\text{even}}$ reaches large values, comparable with the experimental ones, and is strongly anticorrelated with $|u_4|$. The explanation is that when the inclination angle reaches its minimum, high-order deformations grow fast, thus increasing $k^*_{\text{even}}$. However, as eq. (6) shows, this excess area is transferred to low-order deformations, leading to a decrease of $k^*_{\text{even}}$ and a growth of $|u_4|$, which promptly decreases again in favor of $|u_2|$. Then the inclination angle becomes positive, and $|u_4|$ increases again due to the deterministic coupling with $|u_2|$. Therefore, thermal fluctuations not only excite odd modes, but change drastically the dynamics of the even modes as well.

Zoom In Zoom Out Reset image size

Figure 9 shows simulation results that can be compared to experimental ones. The probability distribution function of $k^*$ of fig. 9(a) shows the same main features as fig. 5: for S = 1147, it has a larger tail and a higher average than for $S=60.7$, even though the difference is not as large as in the experiments. Furthermore, the simulations could not reproduce the transition shown in fig. 6 due to a marked dependence not only on S but also on Λ and Δ. Figure 9(b) shows, however, that for fixed Δ and Λ, the expected increase in the average $k^*_{\text{even}}$ is observed, although with a smaller magnitude than in experiments. The snapshots at S = 421 of fig. 9(c) show clearly visible wrinkles which are, however, not as pronounced as the ones of fig. 2. The numerical results seem to underestimate the experimental ones. This discrepancy is probably due to the quasi-spherical approximation which does not describe faithfully vesicles with large Δ. Higher-order curvature terms in the bending energy might also trigger short-wavelength instabilities of the membrane [24]. However, in spite of these approximations, the model still captures most experimental features, mainly because it includes thermal fluctuations consistently.

Zoom In Zoom Out Reset image size

We have shown that a wrinkling instability occurs for vesicles in the TR regime at large enough strain rates. Even though the dynamics is much more complex than in pure elongation flow, the experiments show that the normalized strain rate S is the main control parameter triggering the instability. The onset of wrinkling occurs approximately above $S\simeq60$, which comes from a statistical analysis of the power spectrum. The main experimental features are reproduced by a simple numerical model, which shows that thermal fluctuations are crucial to describe wrinkling correctly. This model does not, however, provide a full quantitative agreement with the experiments due to the quasi-spherical approximation. In order to compare theory and experiment on a more quantitative level, it would be interesting to include thermal noise in more complex numerical schemes such as [25,26], which show a strong coupling between even modes at large S, or to analyze single TR vesicles with the algorithm of [27], which already takes thermal fluctuations into account.

This work is partially supported by a grant from the German-Israeli Foundation.