Hydrodynamics and rheology of a vesicle doublet suspension

TL;DR

This study numerically investigates the hydrodynamics and rheology of vesicle doublets under various flow conditions, revealing how adhesion influences their formation, dynamics, and suspension viscosity.

Contribution

It introduces a high-order boundary integral method to analyze vesicle doublet behavior and identifies flow conditions affecting their stability and rheological properties.

Findings

Vesicle doublets form at equilibrium under quiescent flow due to adhesion.

A critical extensional flow rate can measure vesicle adhesion strength.

Vesicle doublet suspension viscosity depends on vesicle reduced area.

Abstract

The dynamics of an adhesive two-dimensional vesicle doublet under various flow conditions is investigated numerically using a high-order, adaptive-in-time boundary integral method. In a quiescent flow, two nearby vesicles move slowly towards each other under the adhesive potential, pushing out fluid between them to form a vesicle doublet at equilibrium. A lubrication analysis on such draining of a thin film gives the dependencies of draining time on adhesion strength and separation distance that are in good agreement with numerical results. In a planar extensional flow we find a stable vesicle doublet forms only when two vesicles collide head-on around the stagnation point. In a microfluid trap where the stagnation of an extensional flow is dynamically placed in the middle of a vesicle doublet through an active control loop, novel dynamics of a vesicle doublet are observed. Numerical simulations show that there exists a critical extensional flow rate above which adhesive interaction is overcome by the diverging stream, thus providing a simple method to measure the adhesion strength between two vesicle membranes. In a planar shear flow, numerical simulations reveal that a vesicle doublet may form provided that the adhesion strength is sufficiently large at a given vesicle reduced area. Once a doublet is formed, its oscillatory dynamics is found to depend on the adhesion strength and their reduced area. Furthermore the effective shear viscosity of a dilute suspension of vesicle doublets is found to be a function of the reduced area. Results from these numerical studies and analysis shed light on the hydrodynamic and rheological consequences of adhesive interactions between vesicles in a viscous fluid.

Figures

Tables

Table 1. Table 1: The shear rate, period, and effective aspect ratio of a vesicle doublet with reduced area Δ ​ A = 0.9 Δ 𝐴 0.9 \Delta A=0.9 , separation distance δ = 0.4 𝛿 0.4 \delta=0.4 , and Hamaker constants ℋ = 1.5 ℋ 1.5 \mathcal{H}=1.5 and ℋ = 3.5 ℋ 3.5 \mathcal{H}=3.5 . The effective aspect ratios depend on both the Hamaker constant and the shear rate, indicating that for the reported Hamaker constants, the dynamics of the doublet differ from that of a Jeffery orbit.

$\dot{\gamma }$	$p$	${\lambda }_{\mathrm{eff}}$
0.50	36.34	2.49
0.75	25.86	2.72
1.00	20.73	2.96

Equations

$$\mu \Delta \mathbf{u}=\nabla p,$$

$$\nabla \cdot \mathbf{u}=0,$$

$$\mathbf{u}\to {\mathbf{u}}_{\infty },$$

$$\mathbf{u}(\mathbf{x},t)=\dot{\mathbf{x}},$$

$${\mathbf{x}}_s\cdot {\mathbf{u}}_s=0,$$

$$[[T]]\mathbf{n}=\mathbf{\xi },$$

$$\varphi (z)=\mathcal{H}\left[{\left(\frac{\delta }{z}\right)}^m-\frac{m}{n}{\left(\frac{\delta }{z}\right)}^n\right],$$

$$\mathcal{A}\mathbf{x}:=-\mathcal{H}m{\delta }^n\sum _{\begin{array}{c}k=1\\ k\ne j\end{array}}^M{\int }_{{\gamma }_k}\frac{\mathbf{x}-\mathbf{y}}{\Vert \mathbf{x}-\mathbf{y}{\Vert }^{m+2}}\left({\delta }^{m-n}-\Vert \mathbf{x}-\mathbf{y}{\Vert }^{m-n}\right)d{s}_{\mathbf{y}}.$$

$$\frac{dh}{dt}\sim \frac{K_a^{2/3}{F}^{1/3}}{\mu R_0^{10/3}}{h}^3,$$

$$\displaystyle t\sim\frac{\mu R_{0}^{10/3}}{K_{a}^{2/3}F^{1/3}}\frac{1}{2h^{2}}\bigg{|}^{h(t)}_{h(0)},$$

$$F(h)\sim \mathcal{H}\frac{m{\delta }^m}{h^{m+1}}\left(1-{\left(\frac{\delta }{h}\right)}^{n-m}\right)$$

$$t\frac{K_a^{2/3}}{\mu R_0^{10/3}}$$

$$t\sim \frac{\mu R_0^{10/3}}{K_a^{2/3}}\frac{1}{m^{1/3}{\mathcal{H}}^{1/3}{\delta }^{5/3}},$$

$$\frac{t_{\delta =0.2}}{t_{\delta =0.4}}\sim \frac{96}{30}=3.2$$

$$\frac{t_{\delta =0.2}}{t_{\delta =0.3}}\sim \frac{96}{45}=2.1$$

$$\frac{t_{\mathcal{H}=0.1}}{t_{\mathcal{H}=0.5}}\sim \frac{34}{20}=1.71$$

$$p=\frac{2\pi }{\dot{\gamma }}\left(\lambda +{\lambda }^{-1}\right).$$

$$[\mu ]:=\frac{{\mu }_{\mathrm{eff}}-{\mu }_0}{\varphi {\mu }_0}=\frac{1}{\dot{\gamma }{\mu }_0(T_e-T_i)}{\int }_{T_i}^{T_e}⟨{\tau }_{12}⟩dt,$$

$$⟨\tau ⟩=\frac{1}{|\omega |}{\int }_{\gamma }\mathbf{\xi }\otimes \mathbf{x}ds,$$

$$\dot{\mathbf{x}}={\mathbf{u}}_{\infty }(\mathbf{x})+\text{SS}[\mathbf{\xi }](\mathbf{x}),$$

$${\mathbf{x}}_s\cdot {\dot{\mathbf{x}}}_s=0,$$

$$\text{SS}[\mathbf{\xi }](\mathbf{x})$$

$$\mathbf{\xi }=-{\kappa }_b{\mathbf{x}}_{ssss}+(\sigma {\mathbf{x}}_s{)}_s+\mathcal{A}\mathbf{x}.$$

$${\mathbf{x}}^{N+1}-\Delta t{\text{SS}}^N{\mathcal{B}}^N{\mathbf{x}}^{N+1}-\Delta t{\text{SS}}^N{\mathcal{T}}^N{\sigma }^{N+1}={\mathbf{x}}^N+\Delta t{\text{SS}}^N{\mathcal{A}}^N{\mathbf{x}}^N,$$

$${\mathbf{x}}_s^N\cdot {\mathbf{x}}_s^{N+1}=1.$$

$$\varphi (z)=\mathcal{H}\left[{\left(\frac{\delta }{z}\right)}^m-\frac{m}{n}{\left(\frac{\delta }{z}\right)}^n\right],$$

$$U_1={\int }_{{\gamma }_1}{\int }_{{\gamma }_2}\varphi (\Vert {\mathbf{x}}_1-{\mathbf{x}}_2\Vert )d{s}_{{\mathbf{x}}_2}d{s}_{{\mathbf{x}}_1}.$$

$${\tilde{U}}_1={\int }_{{\tilde{\gamma }}_1}{\int }_{{\gamma }_2}\varphi (\Vert {\tilde{\mathbf{x}}}_1-{\mathbf{x}}_2\Vert )d{s}_{{\mathbf{x}}_2}d{s}_{{\tilde{\mathbf{x}}}_1},$$

$$\delta U_1={\int }_{{\tilde{\gamma }}_1}{\int }_{{\gamma }_2}\varphi (\Vert {\tilde{\mathbf{x}}}_1-{\mathbf{x}}_2\Vert )d{s}_{{\mathbf{x}}_2}d{s}_{{\tilde{\mathbf{x}}}_1}-{\int }_{{\gamma }_1}{\int }_{{\gamma }_2}\varphi (\Vert {\mathbf{x}}_1-{\mathbf{x}}_2\Vert )d{s}_{{\mathbf{x}}_2}d{s}_{{\mathbf{x}}_1}.$$

$$\Vert {\tilde{\mathbf{x}}}_1^{\prime }\Vert \approx 1+ϵ(u_s+\nu \kappa ),$$

$$\delta U_1$$

$$\approx ϵ{\int }_{{\gamma }_1}{\int }_{{\gamma }_2}\nabla \varphi (\Vert {\mathbf{x}}_1-{\mathbf{x}}_2\Vert )\cdot \mathbf{y}d{s}_{{\mathbf{x}}_2}d{s}_{{\mathbf{x}}_1}$$

$${\int }_{{\gamma }_2}\nabla \varphi (\Vert {\mathbf{x}}_1-{\mathbf{x}}_2\Vert )d{s}_{{\mathbf{x}}_2}=-m\mathcal{H}{\delta }^n{\int }_{{\gamma }_k}\frac{\mathbf{x}-\mathbf{y}}{\Vert \mathbf{x}-\mathbf{y}{\Vert }^{m+2}}\left({\delta }^{m-n}-\Vert \mathbf{x}-\mathbf{y}{\Vert }^{m-n}\right)d{s}_{\mathbf{y}}.$$

$${\int }_{{\gamma }_2}\nabla \varphi (\Vert {\mathbf{x}}_1-{\mathbf{x}}_2\Vert )d{s}_{{\mathbf{x}}_2}=-4\mathcal{H}{\delta }^2{\int }_{{\gamma }_2}\frac{{\mathbf{x}}_1-{\mathbf{x}}_2}{\Vert {\mathbf{x}}_1-{\mathbf{x}}_2{\Vert }^6}\left({\delta }^2-\Vert {\mathbf{x}}_1-{\mathbf{x}}_2{\Vert }^2\right)d{s}_{{\mathbf{x}}_2}.$$

Full text

Vesicle hydrodynamics, membrane adhesion, drop and bubble phenomena/drop interactions, interfacial flows/thin fluid films, rheology/flow confinement, biological fluid dynamics/collective behavior/clustering, emergence of patterns

††preprint: APS/123-QED

I Introduction

Vesicles (closed fluid-filled phospholipid bilayer membranes) have been widely utilized as biological cell mimics in biophysics and material engineering Sackmann (1996); Fenz and Sengupta (2012); Barthes-Biesel (2016). The ever expanding applications of vesicles have encouraged detailed experimental, theoretical, and numerical investigations of vesicle dynamics in applied flow, electric fields, and magnetic fields. Experimental investigations of vesicles have uncovered vesicle properties and various novel applications Dobereiner (2000); Evans et al. (2013); Sugiyama and Toyota (2018). Theoretical investigations are often limited to small-deformation analysis on a nearly spherical vesicle or spheroidal analysis on a spheroidal vesicle Barthes-Biesel and Rallison (1981); Misbah (2006); Vlahovska and Gracia (2007); Finken et al. (2008); Zhang et al. (2013); Nganguia and Young (2013). On the other hand, various numerical methods have been developed for simulating the transient hydrodynamics of vesicle suspensions Bagchi et al. (2005); Biben (2005); Veerapaneni et al. (2009); Seol et al. (2016).

Hydrodynamics of a single vesicle in Stokes flow has been extensively investigated. In a planar shear flow, the vesicle hydrodynamics is characterized by the reduced volume (reduced area in two dimensions), viscosity contrast between interior and exterior fluids, and shear rate of the imposed far-field fluid flow. In addition, a vesicle with a rigid particle inside is also investigated as a biological mimic of a eukaryotic cell with a nucleus that occupies nearly half of the intracellular volume Veerapaneni et al. (2011). Small-deformation analysis shows that a vesicle tank-treads in a planar shear flow for low viscosity contrast and shear rate. At high viscosity contrast the tank-treading dynamics transitions to tumbling dynamics Misbah (2006); Vlahovska and Gracia (2007), and this leads to a transition in the effective shear viscosity of the vesicle suspension Misbah (2006); Vitkova et al. (2008) that is also validated by direct numerical simulations Ghigliotti et al. (2010) and experiments Vitkova et al. (2008); Deschamps et al. (2009); Kantsler et al. (2008a); Zabusky et al. (2011). Between tank-treading and tumbling vesicle dynamics, a breathing (tremble) mode is also observed Misbah (2006); Kantsler et al. (2008b); Zhao and Shaqfeh (2011); Spann et al. (2014) to alter the effective shear viscosity. In an extensional flow (planar or uniaxial), vesicle shape dynamics depends sensitively on the vesicle reduced volume and the elastic capillary number Kantsler et al. (2008b); Zhao and Shaqfeh (2013); Narsimhan et al. (2014); Dahl et al. (2016): Asymmetric shape and oscillatory undulation of the vesicle membrane are two examples of the complex vesicle hydrodynamics in an extensional flow.

Collective hydrodynamics of vesicles is dictated by the vesicle-vesicle interactions. In a quiescent flow, vesicle-vesicle adhesion leads to the formation of vesicle doublets Ziherl (2007); Ziherl and Svetina (2007) or clusters Svetina and Ziherl (2008); Gu et al. (2016); Flormann et al. (2017). As a model for red blood cell (RBC) aggregates, a simplified model for adhesive vesicle-vesicle interactions can reproduce vesicle shapes similar to those observed in experiments of fibrinogen-induced RBC aggregates Svetina and Ziherl (2008); Gu et al. (2016); Flormann et al. (2017); Hoore et al. (2018); Brust et al. (2014); Clavería et al. (2017). Using the Lennard-Jones (L.-J.) potential for point-point interaction between two vesicle membranes without the molecular details of adhesive interactions between RBCs, Flormann et al. Flormann et al. (2017) found that, under strong adhesion, the membrane may buckle to a sigmoidal shape in the contact region. Such a sigmoidal vesicle shape is also observed in RBC doublets and explained in a slightly different model Ziherl and Svetina (2007). It remains unknown how adhesive interactions between vesicles might lead to different vesicle hydrodynamics that results in different rheological properties of a vesicle suspension. Such studies will provide useful insight into the rheology of blood flow Chien et al. (1967) and how to use nano solutes to control the hypertension by adjusting the adhesive interactions between RBCs.

In numerical simulations of a vesicle suspension, the adhesive interaction between two vesicles is often ignored Veerapaneni et al. (2009); Rahimian et al. (2010) for numerical convenience, mostly because a small time step is often needed to resolve the dynamics of the lubrication thin film between two membranes with adhesion. However, adhesive interactions between RBCs lead to RBC clusters that are expected to alter the rheological properties of the RBC suspension Neu and Meiselman (2002); Svetina and Ziherl (2008). Hydrodynamic simulations of adhesive vesicles in a quiescent flow have revealed physical insights to the equilibrium shapes of RBC aggregates in both experiments Flormann et al. (2017) and theory Ziherl and Svetina (2007). The main goal of this work is to investigate the hydrodynamics of adhesive vesicles in flow conditions common in microfluidics such as a planar shear flow and an extensional flow.

The electrostatic nature of lipid molecules (a hydrophobic tail and a hydrophilic head with an electric dipole) complicates the interactions between a lipid bilayer membrane and another bilayer membrane Evans and Metcalfe (1984); Evans (1988); Israelachvili (1991); Perutkova et al. (2013) or a solid (such as a glass substrate or a nano particle). Adhesion between a vesicle and a solid has been extensively studied with more focus on the static equilibrium shapes Seifert and Lipowsky (1990); Bernard et al. (2000); Shi et al. (2006); Lin and Freund (2007); Gruhn et al. (2007); Das and Du (2008); Keh et al. (2014); Zhang et al. (2009); Agudo-Canalejo and Lipowsky (2015a, b); Steinkuhler et al. (2016); Keh and Leal (2016); Agudo-Canalejo and Lipowsky (2017) than the transient adhesion process Cantat and Misbah (1999); Sukumaran and Seifert (2001); Blount et al. (2013). Adhesive interactions between lipid membranes are essential in many biomedical, biological, and biophysical processes. For example, vesicle adhesion is crucial to initiate membrane fusion and fission in endocytosis, exocytosis, and the transport of small vesicles through membrane surfaces. In the absence of an external electric field and ions in the solvents, it is reported that the adhesive interactions between two membranes can be well approximated by the L.-J.-type potential Flormann et al. (2017), that consists of a long-range attraction component (truncated at some finite distance beyond which membrane interaction becomes negligible) and a short-range repulsion component Israelachvili (1991), and the combination of the two gives rise to an equilibrium distance where the interaction potential is at its minimum Israelachvili (1991).

The strength of membrane-membrane adhesion can be estimated through the membrane contact angle at the edge of contact zone Ramachandran et al. (2010); Mares et al. (2012). Alternatively, researchers have used the micropipette to measure the adhesive force between two vesicles bound by strong adhesion Evans and Metcalfe (1984); Frostad et al. (2014): Evans and Metcalfe Evans and Metcalfe (1984) were able to measure the reduction in the free energy per unit area of membrane-membrane contact formation due to the van der Waals’ attraction. For the adhesive interaction between two lipid bilayer membranes in buffer solutions, a typical range for the energy density is between $1$ to $10$ $\mu J/m^{2}$, similar to the adhesion energy density of $3$ $\mu J/m^{2}$ between two RBCs Flormann et al. (2017). Frostad et al. Frostad et al. (2014) investigated the depletion-attraction induced adhesion between two vesicles, and identified the dynamic role of membrane tension during membrane adhesion/peeling.

For adhesion interaction between a vesicle and a substrate Gruhn et al. (2007), the boundary between weak and strong adhesion is around $1$ $\mu J/m^{2}$: the vesicle-substrate interaction is “strong” when the adhesion energy density is larger than $1$ $\mu J/m^{2}$ and “weak” when the energy density is less than $1$ $\mu J/m^{2}$. The ratio of total adhesion energy to the bending energy (bending modulus of lipid bilayer membranes $\sim 10^{-19}$ $J$) gives a measure of the vesicle deformation in the presence of adhesion Ramachandran et al. (2010).

In this work we focus on regimes where such ratio is of order one, between the weak adhesion (adhesion energy/bending energy $\ll 1$) and strong adhesion (adhesion energy/bending energy $\gg 1$) regimes. In this regime the membrane deformation may increase the “contact area” between two vesicles and enhance the adhesion effects. The equilibration of vesicle membranes in the strong adhesion regime has been well-studied and documented (see Agrawal (2011); Ramachandran et al. (2010); Steinkuhler et al. (2016); Flormann et al. (2017) and references therein). However it is unclear how the adhesion couples to the vesicle hydrodynamics in this regime. This paper aims to address this question with quantitative characterization in terms of physical parameters.

In a quiescent environment, sub-micron size vesicles are found to form a doublet due to their van der Waals’ attractive interactions Ramachandran et al. (2010). Due to the strong van der Waals’ adhesive force, the vesicles are far from spherical shape and the membrane is almost flat in the “contact” region. Gires et al. used small-deformation analysis to investigate the hydrodynamic interactions between two vesicles in a planar shear flow with a long separation distance Gires et al. (2012). They found that the vesicle interaction could be either attractive or repulsive depending on the organization of the two vesicles relative to the shear flow Gires et al. (2014). To the best of authors’ knowledge, the effects of close-range vesicle adhesion on their hydrodynamics under an external flow have not been studied and quantified. The goal of this work is to use state-of-the-art boundary integral simulations to numerically investigate the dynamics of two vesicles in both planar shear flow and extensional flow for a wide range of vesicle shapes and adhesive strength and distance.

Boundary integral equation (BIE) approaches are well-suited for solving the low Reynolds flow problems considered here as they lead to reduction in dimensionality and achieve high-order accuracy even for moving geometry problems. When the vesicles adhere, one major issue for BIE solvers is to resolve the vesicle-vesicle hydrodynamic forces which become nearly singular. We use an interpolation-based quadrature rule Quaife and Biros (2014) to maintain high-order accuracy for all vesicle-vesicle separation distances. To overcome the numerical stiffness induced by the membrane bending forces and to control the global error, we employ a second-order spectral deferred implicit-explicit adaptive time stepping scheme Quaife and Biros (2016).

This paper is organized as follows. In Section II we formulate our model for two-dimensional vesicle hydrodynamics with adhesive interactions between membranes of distinct vesicles. We simulate the adhesion process of two vesicles in a quiescent flow in Section III, where we also present a simple lubrication model to estimate how long it takes for two nearby vesicles to reach the separation distance set by the adhesion potential. In Section IV we first study the hydrodynamics of two vesicles in a planar extensional flow. We then numerically investigate the hydrodynamics of a vesicle doublet in a fluid trap where the stagnation point is actively controlled to be at the center of the vesicle doublet. From these results we propose a novel application of the microfluidic fluid trap to probe the adhesion strength between lipid bilayer membranes in solutions. In Section V we investigate how two vesicles may form a doublet as they move toward each other under a planar shear flow. We map out the parameter regions for bound/unbound vesicles under a planar shear flow, and we also investigate the effects of adhesive interactions on the rheological properties of a dilute suspension of vesicle doublets. Finally in Section VI we discuss the implications of our results and point out a few potential future directions.

II Governing Equations

We consider a suspension of locally inextensible vesicles in an unbounded two-dimensional viscous fluid. For simplicity, we assume the fluid viscosity both inside and outside the vesicles is the same; however, incorporating viscosity contrast is rather straightforward in our numerical algorithm. While we focus on a vesicle doublet suspension for the rest of the paper, here we provide the most general formulation for multiple vesicles interacting with each other through both hydrodynamic and adhesive forces. Individual vesicles are denoted as $\gamma_{j}$, $j=1,\ldots,M$, and they are parameterized in arclength as $\mathbf{x}_{j}(s,t)$. The union of all vesicles is denoted by $\gamma$. Given a background velocity $\mathbf{u}_{\infty}$, the governing (dimensionless) equations are

[TABLE]

where $\mathbf{u}$ is the fluid velocity field, $p$ is the pressure, and the scaled viscosity $\mu=1$ inside the vesicles. Outside the vesicle $\mu=\mu_{e}/\mu_{i}$ with $\mu_{i}$ ($\mu_{e}$) the viscosity of the interior (exterior) fluid. We set the viscosity ratio $\mu_{e}/\mu_{i}=1$ for the rest of the paper. $T=-pI+\mu\left(\nabla\mathbf{u}+\nabla\mathbf{u}^{T}\right)$ is the stress tensor, $\llbracket T\rrbracket$ is the jump across the membrane, and $\bm{\xi}$ is the traction that is the sum of the bending, stretching, and adhesion forces defined in equation (8) of Appendix A. The bending force, arising from the Helfrich energy model, is $\mathcal{B}\mathbf{x}=-\kappa_{b}\mathbf{x}_{ssss}$, where the subscript $s$ denotes derivative with respect to arclength $s$. $\kappa_{b}$ is the bending modulus which we set to be 1 for all examples. The stretching force is $\mathcal{T}\sigma=(\sigma\mathbf{x}_{s})_{s}$, where the tension, $\sigma$, acts as a Lagrange multiplier to satisfy the local extensibility constraint. The resistance to bending and stretching are standard assumptions on vesicles. In this work, we include an adhesive force, $\mathcal{A}\mathbf{x}$, that we now describe.

The adhesion potential for a vesicle membrane has been modeled by a L.-J. type potential

[TABLE]

where $\mathcal{H}$ is the Hamaker constant, $\delta$ is the adhesion length scale, and $z$ is the distance between a patch on the vesicle membrane and a point on the other object, which could be another vesicle membrane or a flat solid wall. The characteristics of the L.-J. potential are summarized in Figure 1(a): The adhesive force is zero at the equilibrium distance $z=\delta$. For large distances ($z>\delta$) the adhesive potential is attractive while for small distances ($z<\delta$) the potential is repulsive to prevent physical contact. The long-range attraction component should decay sufficiently fast to mimic the finite-range attraction between RBCs. In this work we are mainly interested in effects of adhesive interactions between two vesicles drawn together by the surrounding linear flow. Therefore we do not introduce any truncation in the long-range attraction.

The exponents $(m,n)$ in equation (2) depend on the geometry and the molecular details of the two objects under adhesion Israelachvili (1991): $(m,n)=(4,2)$ corresponds to two flat, planar surfaces interacting with each other, and has been a common choice for studying membrane-solid adhesion Seifert and Lipowsky (1990); Sukumaran and Seifert (2001); Shi et al. (2006); Lin and Freund (2007); Blount et al. (2013); Young and Stone (2017). On the other hand, $(m,n)=(12,6)$ corresponds to the L.-J. potential between two molecules, and has been used to model membrane-membrane adhesion Flormann et al. (2017). In our case, the adhesion potential is between two small patches of lipid bilayer membrane (because the lipid molecules are coarse-grained in the continuum modeling). Thus a reasonable choice for $(m,n)$ between two coarse-grained membrane patches would be between those for two planar surfaces and two point molecules.

A large value of $m$ corresponds to a sharp increase in the repulsion force as two objects are within the separation distance $\delta$. This poses a numerical challenge since the problem becomes stiff (i.e., requires a very small time step) for large $m$. The adaptive time-stepping BIE scheme makes it possible to simulate vesicle adhesive dynamics with specified numerical precision for reasonable computation time. We explore several combinations of $(m,n)$ in the simulations of two vesicles forming a doublet in a quiescent flow and in a shear flow. We found that, as long as the close-range interaction is well resolved both in space (using the near-singular evaluations) and time (using the high-order adaptive time integration), there is very little difference in both the dynamic evolution and equilibrium configuration between $(m,n)=(8,6)$ and $(m,n)=(4,2)$. Therefore, in this work, we use $(m,n)=(4,2)$ to regulate the numerical stiffness introduced by the adhesive force.

Focusing on intermediate adhesive strengths, we use Hamaker constants ranging from $\mathcal{O}(10^{-1})$ to $\mathcal{O}(10^{0})$ times the bending modulus. We assume that the adhesion force from equation (2) applies between all pairs of points on different vesicles, and in Appendix B the net adhesive force at a point $\mathbf{x}$ on vesicle $j$ is shown to be

[TABLE]

This summation of adhesive forces is illustrated in Figure 1(b). Figure 1(c) is an example of the calculated adhesive force projected onto their center-of-mass vector. We notice that the adhesive force is repulsive in the contact region, while for the rest of the membranes there is an attractive force that acts to pull the two vesicles together. The summary of boundary integral formulation can be found in Appendix A. For the rest of the paper, the dimensionless bending modulus is set to be one ($\kappa_{b}=1$) unless otherwise specified.

III Adhesion of two vesicles in a quiescent flow

We consider two identical vesicles suspended in a quiescent flow. Without any external forcing such as an imposed electric field or a fluid flow in the far-field, the long-range attraction pulls the vesicles together until their separation distance is close to $\delta$. The L.-J. type potential prohibits physical contact between the two membranes and instead keeps them close to the separation distance $\delta$. In Section III.1, we compute the expected time for the vesicles to reach an equilibrium configuration and compare the analysis with numerical simulations. In Section III.2, we characterize the effects of the adhesion parameters and the vesicles’ reduced area on the shape of the contact region and the bending and adhesive energies.

III.1 Effects of the adhesion parameters on draining times

When two vesicles move towards (or away from) each other under a constant force $F$ without any imposed external flow, the height, $h$, of the thin liquid film between two vesicles follows the draining dynamics Ramachandran and Leal (2010)

[TABLE]

where $K_{a}$ is the area expansion modulus, $\mu$ is the viscosity of the exterior fluid, and $R_{0}$ is the radius of the undeformed vesicle. For the case of a constant forcing $F$ (independent of separation distance $h$ and time $t$), equation (4) can be easily integrated to relate the film thickness $h$ and time $t$:

[TABLE]

where $h(0)$ is the vesicle separation at $t=0$. We note that $F<(>)0$ for an attractive (repulsive) interaction between two vesicles, and consequently $h(t)$ decreases (increases) from the initial separation distance $h(0)$. When the force is attractive, $h$ decreases monotonically, and $t$ in the above equation is the draining time that diverges as $h\rightarrow 0$.

When the force on each vesicle is a function of the separation distance $h$ of the form

[TABLE]

with integers $m>n\geq 2$, attraction is dominant at “large” distances ($h>\delta$) while repulsion is dominant at “small” distances ($h<\delta$). Integration of equation (4) with the $F(h)$ as defined in equation (5) gives the relationship between $t$ and $h$. In this case, the relationship involves an integral of a function of the dimensionless variable $\bar{h}\equiv h/\delta$:

[TABLE]

Solving for $t$,

[TABLE]

and equation (6) tells us that both the adhesion strength $\mathcal{H}$ and the separation distance $\delta$ affect the time it takes for a pair of vesicles to reach the separation distance $\delta$ under adhesion force $F$ in equation (5). In particular, the duration $t$ is proportional to (i) the separation distance $\delta^{-5/3}$ (thus $t\rightarrow\infty$ as $\delta\rightarrow 0$), and (ii) the adhesion strength $\mathcal{H}^{-1/3}$ (as in the constant forcing case Ramachandran and Leal (2010)). From the above analysis we also expect that the scaling of $t$ with respect to $\delta$ depends on the adhesion potential, while the scaling with respect to adhesion strength is independent of the specific form of the potential.

To test the scaling of the duration $t$ with respect to $\mathcal{H}$ and $\delta$, we simulate the vesicle adhesion dynamics in a quiescent flow. Starting with two vesicles at a distance of twice the vesicle radius, the long-range attraction pulls the vesicles together. Figure 2(a) shows the scaled adhesion energy versus time with $\mathcal{H}=0.1$ and three values of $\delta$ as labeled. Adhesion energy reaches minimum at equilibrium, and the scaled $E_{adh}$ evolves towards zero at equilibrium. The arrows indicate the times when the scaled adhesion energy reaches within $1\%$ of equilibrium: $t\sim 96$ for $\delta=0.2$, $t\sim 45$ for $\delta=0.3$, and $t\sim 30$ for $\delta=0.4$:

[TABLE]

The scaling with respect to adhesion strength $\mathcal{H}$ is illustrated in Figure 2(b), where $\delta=0.4$ and $\mathcal{H}$ varies from $0.1$ to $0.5$ as labeled. Again the arrows indicate the times when the equilibrium is reached within $1\%$:

[TABLE]

From the above results, we conclude that the scaling in equation (6) captures the adhesion dynamics of two vesicles that interact with each other via the adhesion force in equation (5).

III.2 Effects of the adhesion parameters on equilibrium

configuration of a vesicle pair

We again simulate two identical adhering vesicles with the adhesion force in equation (5). The initial vesicle separation is smaller than in the previous simulations so that the equilibrium configuration is achieved in a shorter time horizon. Once a vesicle doublet is formed, the membrane shape in the contact region depends on the adhesion strength relative to the membrane bending modulus. For weak to moderate adhesion strength, vesicle membranes are flattened in the contact region while the rest of vesicle maintains a nearly spherical shape Evans and Metcalfe (1984); Evans (1988); Israelachvili (1991); Ramachandran et al. (2010). Under a strong adhesion, however, the vesicle membranes in the contact region buckle and form a sigmoidal shape that has also been observed in RBC doublets Ziherl (2007); Ziherl and Svetina (2007); Flormann et al. (2017). An external electric field is also able to buckle a vesicle membrane that is adhered to a solid substrate Steinkuhler et al. (2016).

In this work we focus on weakly adhesive vesicles whose equilibrium shapes in a quiescent flow are shown in the left plot of Figure 3 for three reduced areas: $\Delta A=0.95$, $0.75$, and $0.6$ with $\mathcal{H}=5$ and $\delta=0.2$.

The equilibrium vesicle shape for $\Delta A=0.95$ is the circular cap with a flat contact region. This is similar to the observed shapes of two vesicles under strong adhesive interaction in Ramachandran et al. (2010). When vesicles are more deflated with a reduced area $\Delta A=0.75$, the equilibrium vesicle shape is elongated with a bigger contact region. This is consistent with the equilibrium shapes of a vesicle doublet under a $(m,n)=(12,6)$ L.-J. potential Flormann et al. (2017). As the reduced area decreases further, we observe undulation of the vesicle membrane on the non-contact side while the contact region remains flat. The color coding along each curve is the tension of the vesicle membrane. We observe that the membrane tension in the contact region is very negative, indicating a dominant compression of membrane when the adhesion force is strong to keep the vesicles bound together.

The right plot of Figure 3 is a zoom (not to-scale) of the membranes in the contact region, where the membranes are not perfectly flat at all three values of reduced area. We observe that the membranes are slightly curved with a dip at the edge. Such membrane undulation in the contact region is predicted by lubrication analyses on an elastic membrane under adhesion with a solid substrate Blount et al. (2013); Young and Stone (2017). Results from the lubrication analysis show that this dip and slightly curved shape in the contact region are independent of the adhesion strength. The high membrane curvature at the edge may pose a problem for using the contact angle there to estimate the adhesion strength. This inspires us to investigate the possibility of using a dynamic fluid trap to measure the adhesion strength (see Section IV).

Also the dip at the edge of the contact region is related to (but different from) the buckling of membrane under strong adhesion: A closer inspection on the dip at the edge shows that the membrane distance is smaller than the neutral separation distance $\delta$ there for more deflated vesicles. For $\Delta A=0.6$ and $0.75$ the adhesion force is mostly attractive in the contact region except at the edge, where the adhesive force turns repulsive. For $\Delta A=0.95$, however, the adhesion force is mostly repulsive, leading to large tension in the rest of the vesicle membrane.

Figure 4(a) shows the total bending energy and adhesion energy as a function of the reduced area. The smaller the reduced area, the more vesicle area is available for deformation and thus the bending energy is higher. In contrast, the total adhesion energy becomes more negative as the reduced area decreases. The sum of the two energies is plotted in Figure 4(b), where a local minimum in the total energy is found around $\Delta A=0.85$.

Figure 5 demonstrates how the adhesion strength $\mathcal{H}$ and separation distance $\delta$ affect the equilibrium configuration of two vesicles under adhesive interactions. Both vesicles have a length of $2\pi$ and a reduced area of $\Delta A=0.95$. The total bending (a) and adhesion (b) energies at equilibrium are plotted against $\mathcal{H}$ for three values of the separation distance $\delta$. We observe that for $\Delta A=0.95$ the equilibrium vesicle shape does not vary much with the adhesion strength $\mathcal{H}$, while the total adhesion energy varies linearly with $\mathcal{H}$.

IV Adhesion of two vesicles in a planar extensional flow

The hydrodynamics of a single vesicle in an extensional flow has revealed novel nonlinear vesicle dynamics not found for a viscous drop Kantsler et al. (2007, 2008b); Zhao and Shaqfeh (2013); Narsimhan et al. (2014); Dahl et al. (2016). Here we focus on a planar extensional flow where the fluid velocity field $\mathbf{u}=\dot{\epsilon}(-x,y)$ is centered at the origin ($(x,y)=(0,0)$), the stagnation point where the fluid flow converges horizontally and diverges vertically. Placed symmetrically at the stagnation point of an extensional flow, a single vesicle may remain steady and symmetric for low extension rate $\dot{\epsilon}$. For sufficiently large extension rate $\dot{\epsilon}$ the vesicle may undergo asymmetric deformation Kantsler et al. (2008b); Narsimhan et al. (2014); Dahl et al. (2016), similar to a surfactant-laden viscous drop under a planar extensional flow Jansen et al. (1997); Hu et al. (2000). Frostad et al. investigated the draining of thin film between two identical vesicles as they collide head-on in an extensional flow Frostad et al. (2013), however they did not include adhesive interaction between vesicles. In the following we first illustrate how two vesicles interact with each other under an adhesive potential in a planar extensional flow. We then illustrate how the adhesion strength between membranes can be measured non-intrusively by using the fluid trap.

At the beginning of the first set of simulations (left panels in Figure 6) we place a vesicle at $(x,y)=(0,0)$, and the other vesicle at $(x,y)=(3,\Delta y)$ where the initial vertical displacement $\Delta y$ varies from $\Delta y=0$ for the top row, $\Delta y=0.02$ for the second row, to $\Delta y=0.1$ for the bottom two rows in Figure 6. The Hamaker constant is ${\cal H}=0.7$ for the top three rows, and is ${\cal H}=0.6$ for the bottom row. The color coding is the membrane tension with the color bar on the right of each row. We note that the negative tension in an extensional flow is consistent with the destabilizing tension for a vesicle with finite bending forces Narsimhan et al. (2014).

Once the planar extensional flow is turned on at $t=0^{+}$, two vesicles move towards each other as shown in Figure 6. For $\Delta y=0$ in the top row, the two vesicles stay on the $x$-axis as they move closer to form a doublet. We observe that the left vesicle (initially centered at $(0,0)$) moves to the left due to the impinging vesicle from the right. As time progresses they form a doublet and reach a steady configuration that is symmetric with respect to the stagnation point (square marks in the top row of Figure 6).

For $\Delta y>0$ we observe similar vesicle dynamics only at the beginning when the right vesicle moves towards the stagnation point and the left vesicle is pushed to the left of the stagnation point. As both vesicles come to close vicinity of the stagnation point, the elevated right vesicle ($\Delta y>0$) is pulled up by the diverging flow in the $y$ direction, causing rotation of both vesicles. We also observe that the larger the initial vertical displacement, the larger the rotation and the faster both vesicles move away from the stagnation point. For sufficiently strong adhesion (${\cal H}=0.7$ for the first three rows in Figure 6), the two vesicles rotate and move away from the stagnation point as a doublet until the end of simulations ($t=100$). As the two vesicles move farther away from the stagnation point, the extensional flow becomes more effective in pulling the two drops apart as shown in the third row of Figure 6. For ${\cal H}=0.6$ (bottom row in Figure 6), the two vesicles are separated by the diverging flow as early as $t=75$, and are clearly not a doublet at the end of simulation.

Here we provide two main observations drawn from the results in Figure 6: (1) The membrane tension is highly non-uniform on the vesicles, leading to a significant Marangoni stress along the membrane as in the case of a surfactant-laden drop at the center of an extensional flow. Thus it is reasonable to expect that the flow around the vesicles is significantly altered. (2) In an extensional flow, vesicles are prone to move away from the stagnation point. A small initial vertical displacement $\Delta y$ in one of the vesicles causes both vesicles to not only rotate but also move away from the center. Once vesicles move off the horizontal axis, their hydrodynamics is dominated by the diverging flow as they move farther away. Thus in the microfluidic experiment it is often more desirable to trap the vesicle around the center of an extensional flow Spjut (2010). With an active control algorithm to place the stagnation point at a desirable location by adjusting the flow with a feedback loop Bentley and Leal (1986); Johnson-Chavarria et al. (2011), a particle can be trapped at the stagnation point for long time scales to facilitate image acquisition or other detailed measures such as particle image velocimetry of flow inside and around the particle.

To explore the application of such a fluid trap to measure the adhesion strength between two bound vesicles, we propose the following experiment. Beginning with two identical vesicles in an equilibrium configuration that form a doublet with a flat contact region (Figure 3), we turn on the fluid trap with the stagnation point placed at the center of the vesicle doublet. Depending on the Hamaker constant, we expect that either the flow overcomes the adhesive force and the doublet is broken, or the adhesive force is sufficiently strong and the doublet reaches a stable stationary configuration.

Based on how the contact region aligns with the extensional flow, we can define the inclination angle $\theta$ of the vesicle doublet as illustrated in Figure 7(b). When $\theta=\pi$ (Figure 7(a)) the diverging flow pulls the vesicles apart from each other and the attractive adhesion force is essential to keep the two vesicles from separating. In contrast, when $\theta=\pi/2$ (Figure 7(c)), the converging flow pushes the vesicles towards each other and the repulsive force is essential to keep the two vesicles at the separation distance. With the converging flow pushing the two vesicles towards the stagnation point dynamically placed at the doublet center, we expect the $\theta=\pi/2$ configuration to be more stable than the $\theta=\pi$ configuration.

In the following numerical experiments we initially place a vesicle doublet with $\theta=\pi$, so that the diverging flow may be strong enough to pull vesicles away from the doublet. At low extension rates, we expect the doublet to stay bound at the fluid trap stagnation point. On the other hand, the vesicle doublet may become unstable and eventually separate at higher extension rates. Thus we expect there to exist a critical extension rate $\dot{\epsilon}_{c}$ above which the vesicle doublet cannot stay bound under a given adhesion potential. Therefore, the dependence of the critical extension rate $\dot{\epsilon}_{c}$ on the adhesion potential and the mechanical vesicle properties provides a means to probe the physics of membrane adhesion.

We consider vesicle doublets with reduced areas 0.70, 0.75, 0.80, 0.85, 0.90, and 0.95, all with a length of $2\pi$, and we vary the extension rate ($\dot{\epsilon}$) between $10^{-2}$ and $10^{-1}$. Since the stagnation point can be controlled in an experimental setting Bentley and Leal (1986); Johnson-Chavarria et al. (2011), we mimic the active control of the microfluidic experiments by moving the center of the doublet at each time step so that the stagnation point occurs exactly in the middle of the doublet. With this adjustment, the vesicles either remain as a doublet in the fluid trap centered around the stagnation point, or the doublet is broken and the vesicles separate from one another.

Figure 8 shows snapshots from simulating a vesicle doublet of reduced area $\Delta A=0.7$ at three different extension rates. A general feature of the doublet dynamics is that the doublet rotates from $\theta=\pi$ configuration towards the more stable $\theta=\pi/2$ configuration. We also observe that with moderate extension rates, the doublet falls short of aligning their long axis with the diverging direction. Moreover, the final inclination angle is closest to the stable $\theta=\pi/2$ for the smallest extension rates. While this doublet remains bound for all the considered extension rates, a doublet with a reduced area of $\Delta A=0.8$ is split with a critical extension rate $\dot{\epsilon}_{c}=0.1$ (Figure 9).

The equilibrium configuration for low extension rate $\dot{\epsilon}=0.02$ (top row in Figure 8) is almost identical to the equilibrium configuration for $\Delta y=0$ (top row of Figure 6). The rotation of the vesicle doublet from $\theta=\pi$ to $\theta=\pi/2$ can be understood as the transition from an unstable configuration $\theta=\pi$ to a stable configuration $\theta=\pi/2$. For higher extension rates (second and third rows in Figure 8) the partially rotated vesicle pairs at equilibrium are also very similar to the case of $\Delta y=0.02$ (second row of Figure 6). Such asymmetric vesicle configuration is reminiscent of the asymmetric deformation of a single vesicle in an extensional flow Kantsler et al. (2008b); Narsimhan et al. (2014); Dahl et al. (2016). To gain more insight into such nonlinear dynamics, in the following we explore more physical characterization of the transition from $\theta=\pi$ to $\theta=\pi/2$ at different reduced area and extension rate.

In Figure 10(a), we plot the inclination angle of the doublet as a function of time for four different extension rates, and the final inclination angle for all doublets that reach an equilibrium state are in Figure 10(b). We observe that not only do smaller extension rates result in smaller inclination angles, but smaller reduced areas also result in smaller inclination angles. We summarize the final inclination angle of the doublet in Figure 10(c). The size of the round dots are scaled to the final inclination angle as defined in Figure 7. At smaller extension rates, the vesicles come closer to aligning their long axis with the diverging direction. When the doublet is broken at reduced area $\Delta A=0.80$, the extension rate is sufficiently large to align the long axis of the vesicle with the diverging direction and then the vesicles separate (Figure 9). For very low extension rates, the time horizon is insufficient for the doublet to reach an equilibrium state, and these simulations are marked with a square. Figure 10(c) also summarizes the reduced areas and extension rates that result in a bound vesicle doublet in a fluid trap: Parameter values with a blue mark result in a fluid trap and parameter values with a red mark result in vesicle separation.

V Adhesion of two vesicles in a planar shear flow

We consider two vesicles under a planar shear flow $\mathbf{u}=\dot{\gamma}(y,0)$, where $\dot{\gamma}$ is the shear rate. To be consistent with simulations in previous sections, the vesicles have a length of $2\pi$ for all reduced areas. The initial placement of the vesicles is chosen so that the flow dives them towards one another in a head-on collision. In the absence of adhesion, two vesicles deform significantly as they collide head-on, deflect to opposite sides of the $y$-axis as they pass one another, and then separate. This head-on collision process of two vesicles or capsules in a planar shear flow is found to be similar between two- Breyiannis and Pozrikidis (2000); Rahimian et al. (2010) and three-dimensions Lac et al. (2007); Lac and Barthes-Biesel (2008); Omori et al. (2013); Rahimian et al. (2015). However, in the presence of adhesion, the vesicles can form a doublet for certain values of the separation distance $\delta$, Hamaker constant $\mathcal{H}$, shear rate $\dot{\gamma}$, and reduced area $\Delta A$. If a doublet does form, then it undergoes a periodic motion with a period that depends on the same parameters.

Figures 11 and 12 show snapshots of two vesicles that have formed a doublet and the color coding is the tension along the vesicles. Each of the vesicles have reduced area $\Delta A=0.9$, shear rate $\dot{\gamma}=0.5$, separation distance $\delta=0.4$, and Hamaker constant $\mathcal{H}=0.7$ (Figure 11) and $\mathcal{H}=2.1$ (Figure 12). Similar to the quiescent example, the membrane tension is negative in the contact region indicating that the membrane is being compressed when the adhesive force is strongest. The two values for the Hamaker constant are chosen since the oscillation period of the vesicle doublet dynamics is both $t^{*}\approx 42$. For both cases the vesicles in the doublet move in tandem with a tank-reading motion that is much slower than the tank-treading motion that occurs in the absence of adhesion. However, the dynamics of a single period for these two Hamaker constants differ significantly. For $\mathcal{H}=0.7$ the individual vesicles undergo a sliding motion as they pass one another (Figure 11). However, for $\mathcal{H}=2.1$, once the contact region is formed, the doublet maintains the same structure and undergoes a tumbling dynamic similar to a Jeffery orbit.

To further understand the formation of the doublet, we compute the minimum distance between two vesicles with reduced area $\Delta A=0.9$, separation distance $\delta=0.5$, and shear rate $\dot{\gamma}=0.5$. Figure 13(a) shows this distance as a function of time for varying Hamaker constants. For Hamaker constants less than a critical value, $\mathcal{H}_{c}$, the adhesive force is not strong enough to form a doublet, and the distance between the vesicles grows with time (red curves). However, when $\mathcal{H}>\mathcal{H}_{c}$, the adhesive force is sufficiently strong to bind the vesicles into a doublet, and then the doublet undergoes a periodic motion (blue curves). In Figure 13(b), we plot the time required for the doublet to make a complete revolution for three different shear rates and two different separation distances. Initially, as the shear rate increases, the period decreases. However, the transition from the weakly adhesive doublet (Figure 11) to the stronger adhesive doublet (Figure 12) results in an increase in the period. Once the Hamaker constant is sufficiently strong that the contact region remains fixed throughout the dynamics, the period of the doublet again decreases as the shear rate increases.

Especially at the large Hamaker constants, the dynamics of the doublet resemble those of a Jeffery orbit. An ellipse with aspect ratio $\lambda$ undergoes a periodic tumbling motion with period Jeffery (1922)

[TABLE]

Instead of measuring the aspect ratio of a doublet in a flow with shear rate $\dot{\gamma}$, we use the computed period to define an effective aspect ratio $\lambda_{\mathrm{eff}}$ that satisfies (7). Therefore, the effective aspect ratio is the aspect ratio of an ellipse in a shear flow that is undergoing a Jeffery orbit with period $p$. For the Hamaker constants $\mathcal{H}=1.5$ and $\mathcal{H}=3.5$, the shear rates, periods, and effective aspect ratio are summarized in Table 1, which illustrates that the effective aspect ratio of the doublet depends on both the Hamaker constant and the shear rate.

Furthermore the oscillatory dynamics in Figure 11 and Figure 12 illustrate that for similar oscillation periods, the actual aspect ratio of the vesicle doublet may be a dynamical variable as well. This indicates that the vesicle doublet’s oscillatory dynamics in a shear flow is more complicated than the Jeffery orbit of a rigid body.

Next, for several shear rates, we determine the critical Hamaker constant, $\mathcal{H}_{c}$, that determines if the vesicles form a doublet or separate. Figure 14 indicates whether a doublet with $\Delta A=0.9$ is formed (blue) or not (red) for two different separation distances. The plot includes a line of best fit between the stable and unstable regions. Since the transition is linear for both separation distances, this indicates that a dimensionless number involving the ratio of the Hamaker constant and shear rate could be used to determine if a doublet forms or not. We note that the critical Hamaker constant also depends on the reduced area. For smaller reduced areas, we expect a smaller contact region to form and this will lower the Hamaker constant, and this behavior will be investigated in future work.

Finally, we investigate the rheology of a suspension of a doublet by computing the effective viscosity of a doublet and compare it to the effective viscosity of a single tank-treading vesicle. The effective viscosity is defined as the viscosity of a homogeneous Newtonian fluid with the same energy dissipation per macroscopic element of fluid. In a simple shear flow, the intrinsic viscosity, $[\mu]$ is

[TABLE]

where

[TABLE]

$\phi$ is the area fraction of vesicles, $\tau$ is the stress due to the vesicles, $\langle\cdot\rangle$ is the spatial average, $|\omega|$ is the total area of both vesicles, and $\bm{\xi}$ is the traction as defined in equation (8).

In Figure 15(a), we compare the intrinsic viscosity of a single tank-treading vesicle to a doublet with Hamaker constant $\mathcal{H}=0.7$, separation distance $\delta=0.4$, and shear rate $\dot{\gamma}=0.5$. To validate our simulations, we superimpose (black marks) the intrinsic viscosity calculated by Ghigliotti et al. Ghigliotti et al. (2010). The presence of the doublet slightly increases the intrinsic viscosity at all the reduced areas. We also compute the intrinsic viscosity of a doublet with $\Delta A=0.9$, $\dot{\gamma}=0.5$, and $\mathcal{H}=2.1$. Again, this doublet has a similar period to the case $\mathcal{H}=0.7$. At this larger Hamaker constant, the effective viscosity increases slightly from $1.97$ to $2.06$. To further characterize the effect of adhesion, in Figure 15(b) we decompose the intrinsic viscosity into the contributions from the bending and tension (black), and the contribution from the adhesion (green). We also superimpose the intrinsic viscosity of a single tank-treading vesicle (blue) to demonstrate that the bending and tension of the doublet behave similarly, but not identically, to a dilute suspensions of non-adhering tank-treading vesicles. We see that the effect of the adhesion on the intrinsic viscosity is largest for vesicles with small reduced areas. We also decompose the intrinsic viscosity contributions of the doublet with $\mathcal{H}=2.1$ and $\Delta A=0.9$. When compared with the smaller Hamaker constant, the intrinsic viscosity contribution due to bending and tension is reduced from $1.85$ to $1.71$, and the contribution due to adhesion is increased from $0.12$ to $0.36$. Therefore, as the Hamaker constant increases, the increase in the intrinsic viscosity is due to the adhesive force rather than the bending and tension.

VI Conclusions

In this work we use a boundary integral formulation with adaptive time-stepping to simulate hydrodynamics of two vesicles with adhesive interactions. In a quiescent flow, two vesicles that are initially sufficiently far apart move towards each other under a long-range attraction. We use a lubrication theory to estimate the time required to reach the separation distance $\delta$, and the theoretical scaling is in good agreement with numerical results. Once two membranes are within separation distance $\delta$, the adhesive force turns repulsive and the membranes flatten to form a contact region. Our simulations show that the membranes in the contact region are actually curved with end points at the shortest distance, consistent with results for moderate adhesion strength in Flormann et al. (2017). Once a vesicle doublet forms, we examine the dependence of membrane bending and adhesion energies on the reduced area, the Hamaker constant, and the separation distance.

Next we conduct numerical simulations of two vesicles in a planar extensional flow. With initial configurations where two vesicles collide around the stagnation point (one vesicle at the stagnation point and the other vesicle placed close to the converging axis), the two vesicles form a doublet and stay around the flow center only when two vesicles collide head-on with no displacement $\Delta y=0$. For $\Delta y>0$ we find that the vesicles rotate and move away from the stagnation point, and eventually move farther from each other due to the strong diverging flow in the far-field.

This inspires us to conduct a numerical experiment where a vesicle doublet is placed at the center of a fluid trap, which can be actively controlled in microfluidic channel so that fluid trap center is effectively the stagnation point of an extensional flow characterized by a extensional flow rate. At low flow rate, we find the doublet to rotate nearly ninety degrees to align with the flow such that the long axis of the doublet is parallel to the divergent axis and the convergent stream is pushing the two vesicles together. As the flow rate increases, the vesicle doublet rotates less, and when the flow rate exceeds the critical value, the diverging flow breaks the doublet structure by pulling the vesicles apart. These results indicate that it is possible to use the fluid trap to separate a vesicle doublet under adhesion, and thus provide a means to probe the adhesion strength between membranes. For a pair of $\mu m$-sized vesicles with a bending modulus of $\sim 10^{-19}$ $J$, an extensional flow rate of $0.5$ $s^{-1}$ is expected to separate a vesicle doublet with reduced area of $\Delta A=0.8$ and a Hamaker constant $\mathcal{H}=0.7$, which corresponds to $\sim 1$ $\mu J/m^{2}$. These conditions are quite realizable in microfluidic experiments, and we hope that our simulations will inspire microfluidic experiments to be designed to measure membrane-membrane adhesion.

We also examine how adhesive interaction may dynamically lead to the formation of a vesicle doublet. We simulate two vesicles approaching each other in a planar shear flow, and examine how their adhesive interactions lead to doublet formation. Once a doublet forms, the two vesicle membranes rotate around each other as they deform dynamically. The usual tank-treading motion of a vesicle under shear flow is not observed in each of the two vesicles. We compute the effective shear viscosity of a dilute suspension of vesicle doublets, and found it to be more than twice the effective shear viscosity of a dilute suspension of single vesicles. Furthermore, we find that the membrane adhesion contribution to the shear viscosity increases with decreasing reduced area while the bending/tension contribution increases with the reduced area.

In our formulation we did not include any electrostatic interactions between membranes under adhesion. When the electrostatic interaction is important, an electro-osmotic pressure in the thin film is found to be responsible for the observed membrane undulation Steinkuhler et al. (2016). In addition, simulations presented in this work are for identical vesicles (same reduced area, length, and bending modulus) in the doublet. Flormann et al. demonstrated that asymmetric vesicle reduced area may lead to various equilibrium doublet shape such as male-female, asymmetric S-shape, and parachute shape. It is also possible that viscosity contrast may also lead to different equilibrium doublet shapes. Future work includes three dimensional simulations, dispersive vesicle properties (such as reduced area and bending modulus), viscosity contrast, and effects of confinement on adhesive interactions. Another future direction is to consider the clustering and packing of an unbounded suspension of vesicles Flormann et al. (2017).

Finally, it is not clear how thermal fluctuations may affect the hydrodynamics of vesicles under adhesion. For example, does the fluctuating hydrodynamics in the thin film between two vesicles enhance adhesion to keep vesicles bound under linear flows as speculated in cells Fenz et al. (2017)? Recently Liu et al. Liu et al. (2019) used immersed boundary simulations to show that, at a separation distance of tens of nanometers, the thin film between the two membranes facilitates the coupling between membranes via strong hydrodynamic interactions. In particular, they demonstrate numerically that the fluctuation in one membrane is highly correlated to the other membrane without any physical contact. We are actively pursuing this direction with hydrodynamic modeling and simulations of adhesive membranes with thermal fluctuations.

Acknowledgements.

BQ acknowledges support from Florida State University startup funds and Simons Foundation Mathematics and Physical Sciences-Collaboration Grants for Mathematicians 527139. SV acknowledges support from NSF under grants DMS-1719834 and DMS-1454010. YNY acknowledges support from NSF-DMS 1614863 and NSF-DMS 1412789. Both SV and YNY were also supported by the Flatiron Institute, a division of Simons Foundation. BQ and YNY contributed equally to model development and simulations. BQ developed the code. All authors contributed equally to writing the manuscript.

Appendix A Integral Equation Formulation

Using potential theory and following Veerapaneni et al. (2009), we recast the governing equations (1) as integro-differential equations for the evolution of membrane positions:

[TABLE]

where the single-layer potential $\SS[\cdot]$ is defined by

[TABLE]

and the membrane force $\bm{\xi}$ is a sum of the bending, tension and, adhesion forces:

[TABLE]

Defining the bending operator as $\mathcal{B}[\mathbf{f}](\mathbf{x})=-\kappa_{b}\mathbf{f}_{ssss}$, the tension operator $\mathcal{T}[\sigma](\mathbf{x})=(\sigma\mathbf{x}_{s})_{s}$, and using IMEX-Euler, the no-slip boundary results in the time stepping method

[TABLE]

and the inextensibility constraint that is discretized as

[TABLE]

We discretize the vesicles at a set of collocation points, compute the bending and tension terms with Fourier differentiation, and apply Alpert quadrature Alpert (1999) to the weakly-singular single-layer potential $\SS$. The source and target points of the adhesion force never coincide since they are always on different vesicles, so the adhesion force (3) is computed with the spectrally accurate trapezoid rule Trefethen and Weideman (2014).

The dynamics of a doublet undergoes many different time scales over time horizons that are sufficiently large to characterize the formation of a doublet and its rheological properties. Therefore, time adaptivity is crucial so that a user-specified tolerance is achieved without using a guess-and-check procedure to find an appropriately small fixed time step size. To control the error and achieve second-order accuracy in time, we use a time adaptive spectral deferred correction method that applies IMEX-Euler twice per time step Quaife and Biros (2016).

Appendix B Adhesion Force

Consider a suspension of two vesicles $\gamma_{1}$ and $\gamma_{2}$ parameterized as $\mathbf{x}_{1}(s)$ and $\mathbf{x}_{2}(s)$, respectively, with $s\in[0,1]$. Here $s$ is the arclength, and we have assumed, without loss of generality, that both vesicles have length one. We use the L.-J. type potential

[TABLE]

where $z$ is the distance between two points on a pair of vesicles. Then, we define the total adhesive energy on $\gamma_{1}$ to be

[TABLE]

Perturbing $\mathbf{x}_{1}$ to $\tilde{\mathbf{x}}_{1}=\mathbf{x}_{1}+\delta\mathbf{x}_{1}$ results in a new vesicle $\tilde{\gamma}_{1}$, and the perturbed adhesive energy is

[TABLE]

and the change in the energy is

[TABLE]

We now decompose the perturbation into normal and tangential components as $\delta\mathbf{x}_{1}=\epsilon\mathbf{y}(s)=\epsilon(u\mathbf{t}+v\mathbf{n})$. The perturbed arclength term, to leading order, is

[TABLE]

where $\kappa$ is the curvature. To leading order, inextensible perturbations satisfy $u_{s}+\kappa v=0$, so the arclength term of $\gamma_{1}$ and $\tilde{\gamma}_{1}$ are identical to leading order. Therefore,

[TABLE]

and the adhesive force applied by vesicle 2 on vesicle 1 is

[TABLE]

When $(m,n)=(4,2)$, the above expression becomes

[TABLE]

A similar expression holds for the adhesive force applied by vesicle 1 on vesicle 2, and equation (3) gives the adhesive force for a suspension of $M$ vesicles.