Hydrodynamic Coupling of Particle Inclusions Embedded in Curved Lipid Bilayer Membranes

We develop theory and computational methods to investigate particle inclusions embedded within curved lipid bilayer membranes. We consider the case of spherical lipid vesicles where inclusion particles are coupled through (i) intramembrane hydrodynamics, (ii) traction stresses with the external and trapped solvent fluid, and (iii) intermonolayer slip between the two leaflets of the bilayer. We investigate relative to flat membranes how the membrane curvature and topology augment hydrodynamic responses. We show how both the translational and rotational mobility of protein inclusions are effected by the membrane curvature, ratio of intramembrane viscosity to solvent viscosity, and inter-monolayer slip. For general investigations of many-particle dynamics, we also discuss how our approaches can be used to treat the collective diffusion and hydrodynamic coupling within spherical bilayers.


Introduction
Cellular membranes are complex heterogeneous materials consisting of mixtures of lipids, proteins, and other small molecules [3]. The individual and collective dynamics of these species are fine-tuned to carry out complex cellular processes ranging from cell signalling to shape regulation of organelles [3,18,33,44,46,57]. The effective two dimensional fluid-elastic nature of cell membranes results in interfacial phenomena and interesting geometric shapes effecting both molecular interactions and dynamics that can be very distinct from their bulk counter-parts. To gain a deeper understanding of cellular processes requires insights into the fundamental mechanics of fluid-elastic bilayer membranes.
Early theoretical investigations of the hydrodynamics of flat lipid bilayer membranes include [51,52] and more recently the related works [13,14,33,37,39,43,46]. In the now classic papers of Saffman and Delbrück [51,52], the bilayer is treated as a two dimensional fluid.
The two dimensional fluid is coupled to a bulk three dimensional fluid accounting for the solvent surrounding the membrane on both sides. This description of the hydrodynamics is used to model a protein inclusion within a flat infinite membrane to derive the self-mobility M SD = (1/4πµ m ) (log(2L SD /a) − γ). This asymptotic result assumes a L SD , where a is the protein size, γ ∼ 0.577 is the Euler-Mascheroni constant. The L SD = µ m /2µ f is the Saffman-Delbrück length associated with how dissipation within the entrained bulk solvent fluid of viscosity µ f regularizes the long-range two dimensional flow of viscosity µ m . These results highlight the importance of dissipation in the bulk solvent fluid that if neglected would otherwise lead to the well-known Stokes paradox [27,34,51]. This shows that particle motions even within a flat interface has a very different character than its counter-part in a bulk fluid. From Stokes theory the bulk self-mobility of a particle scales like M ∼ 1/6πµ f a [2,27]. For curved membranes the topology and geometry can result in even more significant differences. This includes providing a finite closed membrane surface and trapped solvent fluid in a bounded interior domain augmenting the hydrodynamics and coupling.
More recent works explore the mechanics of membranes both through coarse-grained molecular models [10,11,16,17,19,20,24,50] and through continuum mechanics approaches [4,10,13,15,20,23,25,28,29,31,33,36,37,41,43,46,48,49,[53][54][55][56]. The particular works [18,28,29,37,40,46] introduce a continuum mechanics description for the hydrodynamics of spherical vesicles and tubules. The self-mobility of an embedded particle is computed as the curvature is varied using a truncation of the series representation of the hydrodynamic flow. In [4] an exterior calculus description of the continuum mechanics of a fluid-elastic membrane sheet is introduced and used to investigate lipid flow during processes such as budding with an asymptotic model for the ambient fluid. The prior work in this area primarily has focused on single particle mobility and transport by hydrodynamics averaged over the two bilayer leaflets.
We introduce here further approaches to investigate the collective hydrodynamic coupling of multiple particle inclusions within leaflets of curved fluid lipid bilayer membranes. We consider the case of spherical bilayer membranes where inclusion particles are coupled through (i) intramembrane hydrodynamics, (ii) traction stresses with the external and trapped solvent fluid, and (iii) intermonolayer slip between the two leaflets of the bilayer. We formulate tractable descriptions of the continuum mechanics of curved fluid bilayers drawing on results from the exterior calculus of differential geometry. We formulate a tractable description for the collective hydrodynamic coupling of the inclusion particles on curved manifolds building on our prior work on immersed boundary approximations [6,8,9,56]. We compute the translational and rotational mobilities of inclusion particles. Relative to infinite flat membranes, we show that spherical vesicles exhibit significant differences arising from the curvature and finite domain size.
In Section 2 we introduce our continuum mechanics description of the bilayer hydrodynamics expressed in terms of the operators of exterior calculus of differential geometry. We use exterior calculus to help take a less coordinate-centric approach in our derivations and to obtain more concise expressions that often have a more clear geometric interpretation. We also show how the exterior calculus can be used to generalize many of the techniques used in fluid mechanics to the context of curved surfaces. In Section 2.2, we use Lamb's solution for the fluid flow exterior and interior to a spherical shell to obtain the traction stresses arising from the surrounding solvent fluid and the trapped solvent fluid. In Section 2.5, we consider the hydrodynamic flow within the lipid bilayer membrane. We use a spherical harmonics representation to derive analytic results for the solutions of the coupled hydrodynamic equations. In section 2.7, we discuss some roles played by curvature in hydrodynamic flows within surfaces.
In Section 3, we introduce immersed boundary approximations on manifolds to account for the coupling between the lipid flow and inclusion particles. We discuss some particular properties of this type of approximation. We then derive mobility tensors for the translational and rotational motions of inclusion particles within curved membranes.
In Section 4, we investigate the self mobility and the coupled mobility of inclusion particles when varying (i) vesicle curvature, (ii) membrane viscosity vs solvent viscosity, and (iii) intermonolayer slip. In Section 4.3, we consider approaches for the collective dynamics of many coupled inclusion particles within spherical vesicles. We consider the collective mobility associated with an attracting cluster of particles and briefly discuss some of the interesting dynamics that can arise from the collective hydrodynamic coupling. In Appendix A, we discuss briefly how we have addressed some of the issues that arise for spher-ical surfaces in practical numerical calculations to obtain our results.
In summary, the work presented here is meant as a starting point to understanding the basic features of the mobility of inclusion particles within curved bilayers. Overall, we expect our approaches introduced here to provide ways to investigate the general collective dynamics of inclusion particles within spherical lipid bilayers relevant in many applications.

Continuum Mechanics of the Vesicle
We formulate a continuum mechanics description of (i) the hydrodynamic flow of lipids within the two bilayer leaflets, (ii) intermonolayer slip, and (iii) coupling to the surrounding solvent fluid, see Figure 1. We derive a set of conservation laws on manifolds using tensor calculus and results from differential geometry similar to Marsden [38]. We then use identities as in Arroyo and Disomone [4] to express our equations in a convenient covariant form that is geometrically invariant. To obtain analytic results for hydrodynamic flows on the curved surface, we use exterior calculus to generalize techniques often employed in fluid mechanics to 2-manifolds. We then use these exterior calculus approaches to perform numerical calculations. While our approaches provide rather general methods for working with hydrodynamics within manifolds, we focus in this paper on the sphere which is relevant to flow within the lipid bilayers of vesicles.

Hydrodynamics of Bilayer Leaflets
We first consider the hydrodynamics within a single bilayer leaflet of the membrane. We treat the membrane as a two-dimensional embedded continuum in the case that the surface velocity V = v + v n n has zero velocity in the direction of the surface normal v n = 0. The conservation of momentum and mass of such a deforming two-dimensional continuum is given by [38] ρ The ∇ denotes the covariant derivative which when expressed in terms of tensor components is where Γ a bc denotes the Christoffel symbols [1,47]. In the notation div(·) and grad(·) the corresponding covariant operations for divergence div(w) = w a |a and gradient grad(w) a b = w a |b . The ρ denotes the mass density per unit surface area, the v the velocity components tangent to the surface, b the body force per unit surface area, and σ the surface stress tensor. We remark that while these equations look superficially similar to the Euclidean case owing to the convenient covariant derivative notation, as we shall discuss the geometry introduces important differences and additional terms. Figure 1: Vesicle Hydrodynamics. We take into account the hydrodynamics within each leaflet of the bilayer, the intermonolayer slip between leaflets, and the traction stresses for both the solvent fluid trapped interior to the vesicle and the solvent fluid exterior to the vesicle. We use a covariant formulation of the continuum mechanics.
The constitutive law for an incompressible Newtonian fluid can be expressed in terms of the rate of deformation tensor of the surface which in terms of tensor components is The Newtonian stress is given by The µ m and µ m are the first and second viscosities of the membrane. The I is the (1,1)identity tensor with (I) a b = δ a b where δ a b denotes the Kronecker delta-function. This has I ab = g ab = g ab , where g is the metric tensor for the surface [38]. We have adopted the notation for raising and lowering indices corresponding to the isomorphisms between the tangent and co-tangent spaces of the surface given by : The ∂ x j denotes the coordinate associated basis vectors of the tangent space and dx j the oneform coordinate associated basis of the co-tangent space. The isomorphisms can also be expressed directly in terms of the components as v i = g ij v j and v i = g ij v j , where we denote the metric tensor as g ij and its inverse as g ij [1]. This extends naturally to tensors. Using these conventions, the steady-state Stokes equations corresponding to equation 1 can be expressed in tensor components as We can express this in a more geometrically transparent manner by considering further the divergence of the rate of deformation tensor We have that where ∆ R v := rough-laplacian(v) = div (grad (v)) and K is the Gaussian curvature of the surface. We have that This follows since div (v) = 0.
It is convenient to express the equations and differential operators in terms of the exterior calculus as follows. Let d denote the usual exterior derivative for a k-form α as The δ denotes the co-differential operator given by δ = d , where for a k-form α = (1/k!) α i 1 ...i k dx i 1 ∧ · · · ∧ dx i k the denotes the Hodge star given by where α i 1 ...i k = g i 1 1 · · · g i k k α 1 ... k , |g| is the square-root of the determinant of the metric tensor, and i 1 ...i k j 1 ...j n−k denotes the Levi-Civita tensor [38]. The generalization of the common differential operators of vector calculus to manifolds can be expressed in terms of exterior calculus as The f is a smooth scalar function and the F is a smooth vector field. There are different types of Laplacians that can be defined for manifolds The ∆ R = div(grad(·)) denotes the rough-Laplacian given by the usual divergence of the gradient. For vector fields, ∆ H (F) denotes the Hodge-de Rham Laplacian, which has similarities to taking the curl of the curl [1]. In fact, in the case that div(F) = −δF = 0 we have Using these conventions, we have div(D) = −δdv + 2Kv (21) where we used that div(v) = −δv = 0. This allows for the steady-state Stokes problem on the surface to be expressed using exterior calculus in the convenient form As we shall discuss, this form provides a very convenient approach for analytic and numerical calculations.

Coupling to External Solvent Fluid
The solvent fluid surrounding the lipid bilayer membrane also exerts a traction stress on the inner and outer leaflets. We account for this using the Stokes equations The Ω = Ω ± denotes either the outside region Ω + of fluid surrounding the vesicle or the domain Ω − of fluid trapped inside the vesicle. The solution to the Stokes equations and traction stress can be conveniently expressed in terms of harmonic functions. This is most immediately seen for the pressure, which when taking the divergence of equation 23, yields For the spherical geometry, the pressure can be expanded as where p n is the solid spherical harmonic of order n where Y n m (θ, φ) = e imφ P m n (cos(θ)).
For the solvent fluid velocity u − in the domain Ω − interior to the vesicle, Lamb showed that the solution can be expressed as [27,34] We shall refer to this as the Lamb's Solution.
The surface flow v determines the solid spherical harmonic functions χ n , Φ n , p n by the following relations The R is the radius of the spherical surface. For the surface fluid velocity V = v + v n n of the membrane, the X n , Y n , Z n are combined surface spherical harmonics of degree n obtained by expanding the following scalar fields on the surface For the solvent fluid velocity u + in the domain Ω + exterior to the vesicle, Lamb's solution is [27,34] where The surface fluid velocity V = v + v n n of the membrane, determines the harmonic functions In the special case of v n = 0 and div(v) = 0 we have that X n = Y n = 0 and that only the Z n term is non-trivial. The Lamb's solution simplifies to In this case, the traction stress of the external fluid on the lipid bilayer membrane is The n ± denotes the unit normal on the surface ∂Ω ± in the direction pointing into the domain.
In these expressions, we emphasize that n ± is to be held fixed during differentiation. From equation 50 and 51 and the properties of solid spherical harmonics, we have that the traction stress on the membrane leaflets can be expressed as

Intermonolayer Slip
We account for the two bilayer leaflets of the membrane by considering two surface velocity fields v + and v − . We model the intermonolayer slip between these two leaflets by the traction term proportional to the difference in lipid velocity

Full Membrane Hydrodynamics
Using an approach similar to Section 2.1, we obtain for a two-leaflet membrane the following hydrodynamic equations.
The M ± denotes the two surfaces representing the inner and outer bilayer leaflets. These equations take into account the internal membrane hydrodynamics of each leaflet of viscosity µ m , the intermonolayer slip γ, and the traction stresses with the bulk solvent fluids of viscosity µ ± trapped within and external to the vesicle. To obtain the coupling in the collective dynamics of inclusions embedded in such bilayer membranes, we must solve these equations for the hydrodynamic flow.

Membrane Hydrodynamics and Modal Responses
We use exterior calculus methods to derive solutions to the membrane hydrodynamic equation 54. For analytic and numerical calculations of flow within surfaces, the exterior calculus provides a number of advantages over more coordinatecentric approaches such as tensor calculus [1]. As already seen in our expressions of the hydrodynamic equations, there are fewer explicit references to the metric tensor with instead more geometrically intrinsic operations appearing such as the exterior derivative and Hodge star [1]. In analytic calculations, we take advantage of this to develop succinct methods for curved manifolds that generalize many of the vector calculus based techniques often employed in fluid mechanics. From the exterior calculus formulation of the Stokes equations we can readily show that the incompressible surface flow can be expressed in terms of a scalar velocity potential Φ as This provides a generalization for the surface geometry of the usual velocity potential used in fluid mechanics. The equation 55 generalizes to surfaces the operation in Euclidean space of taking the curl to obtain an incompressible flow [2,27]. The exterior calculus allows us to readily verify that the generated velocity field on the surface is incompressible This follows since = −1 on a surface (2manifold) and d 2 = 0 holds [1].

Modal Response for Intramembrane Hydrodynamics
To obtain equations for Φ, we use the exterior calculus to determine the eigenfunctions of the operator in the Stokes equations. This can then be used to rigorously derive expressions for the modal responses of the hydrodynamics when acted upon by an applied force in a manner similar to [29]. For this purpose, we consider the eigenproblem This can be satisfied if Φ s is a solution of This can also be expressed as where γ s = (λ s /µ − 2K). For scalar fields, the operator −δd is the Laplace-Beltrami operator of the surface. In the special case of the sphere, the solutions are surface spherical harmonics of the form where s = ( , m) subject to the restriction |m| ≤ . The eigenvalues are γ s = − ( + 1)/R 2 and λ s = µ − ( + 1)/R 2 + 2K . We can express the solution of the Stokes equations 22 by expanding the velocity field as We can also represent the solution with Φ = s a s Φ s . In a similar manner, the applied surface force can be expanded with coefficients c s as c = b − dp = − d s c s Φ s . The problem now becomes to find the coefficients a s for the flow when given an applied force with coefficients c s .
As a demonstration of the utility of this exterior calculus approach, consider the Stokes equations 22 for the surface flow on a sphere associated with a single leaflet, without yet the intermonolayer slip or traction stress. We treat the Stokes problem in equation 22 by taking − d of both sides to eliminate the pressure term. This yields Using the expansion for v in terms of v s = − dΦ s and that Φ s was chosen to solve the eigenproblem in equation 57, we have We remark that we use c as opposed to b throughout our calculations to emphasize that only the solenoidal component of the applied force effects the flow. This is further exhibited in the identity − dc = − db . For mode s, we have λ s a s = c s and K = 1/R 2 which gives This applies for ≥ 2. For the Stokes flow on the membrane surface this gives the modal response to an applied force. We have assumed for this solution that the applied force has net-zero torque. The mode = 1 does not introduce an internal shear stress within the membrane since this mode corresponds to a rigid-body motion of the spherical shell. Since we have not yet included the intermonolayer slip or the external fluid traction stress there would be no stresses to balance a force having non-zero net torque.
We remark that the membrane velocity field is obtained from these calculations by The |g| = det(g) is the determinant of the metric tensor and i is the Levi-Civita tensor (slight abuse of notation). The x denotes the coordinates. We remark if the standard spherical coordinates are used then x 1 = θ and x 2 = φ with the polar angle θ and azimuthal angle φ. However, this has singularities at the north and south poles of the sphere. To use robustly the velocity representation 66 for numerical calculations at each location on the sphere, we need to use different coordinate charts depending on location. The velocity field can be computed robustly by using either the standard spherical coordinates or the spherical coordinates with the poles at the east and west poles, for details see our discussion in Appendix A.

Modal Response when Coupled to External Solvent Fluid and with Intermonolayer Slip
Using this approach, we can also incorporate for a two leaflet lipid bilayer membrane the additional contributions of the traction stress with the external solvent fluid and the intermonolayer slip. We consider the case when the outer bilayer leaflet of the membrane vesicle has radius R + and the inner bilayer leaflet has radius R − . This requires us to derive the modal response to the induced bulk external solvent flow. The solvent flow satisfies the Stokes equations 23-26 with noslip with respect to the flow within the membrane surface. These Stokes equations must be solved twice, once in the domain Ω + exterior to M + and once in the domain Ω − interior to M − . We obtain a representation for these solutions using Lamb's solution [27], see Section 2.2.
We represent the fluid velocity v ± within each leaflet of the membrane using the velocity potential Φ ± . As in Section 2.5.1, we expand the velocity potential in spherical harmonics Φ s as Φ ± = s a ± s Φ s . This allows us to express the membrane velocity as From this representation and equation 52, we can express the traction stress from the external solvent fluid on the membrane leaflet as TheΦ ± denotes the linear combination of modes of degree . In particular,Φ ± = s , = a ± s Φ s where s = (m , ).
By applying − d we have The solution coefficients for v + and v − for the full two-leaflet membrane hydrodynamics in equation 54 can be expressed as with (74) Associated with the inner and outer external fluids, we define the length-scales L − = µ m /µ − and L + = µ m /µ + . The Saffman-Delbrück length-scale [51,52] associated with each leaflet is L − SD = 1 2 L − and L + SD = 1 2 L + and on average L SD = 1 2 L − SD + L + SD . In summary, the equations 72-74 provide the modal responses for the hydrodynamic flow satisfying the two leaflet Stokes problem in equation 54. The model captures the hydrodynamic flow of lipids within the two curved bilayer leaflets that are coupled to one another by intermonolayer slip and that are coupled to the flow of the external solvent fluid. The key parameters are given in Table 1.
We remark that the membrane fluid velocity fields are obtained from these calculations by The |g ± | denotes the determinant of each of the metric tensors g ± associated with the leaflet surfaces M ± and the i denotes the Levi-Civita tensor. We remark that given coordinate singularities on the sphere to use robustly this velocity representation for numerical calculations, we need to use different coordinate charts. For details see our discussion in Appendix A.

Radius of the inner leaflet.
R Average radius of the vesicle.

Characteristic Physical Scales
To characterize the hydrodynamic responses, we discuss a few useful non-dimensional groups. We first consider how the bulk solvent fluid regularizes the two dimensional membrane hydrodynamics. This can be characterized by considering a disk-shaped patch of a flat membrane of radius r. An interesting length-scale is the radius r where the bulk three dimensional traction stress acting on the patch of area A = πr is comparable in magnitude to the intramembrane stresses acting on the perimeter of the patch of length = 2πr. This occurs for the inner and outer leaflets on length-scales scaling respectively like For a vesicle, it is natural to consider these length-scales relative to the radius of the vesicle R. We introduce the non-dimensional groups Π For the intermonolayer slip, we consider for the flow the response of the leading order modes with = 1. These correspond to the rigid rotations of the spherical shell. For a velocity difference between the layers, the drag is given by γ. For the leading order modes with = 1, the traction stresses arising from the entrained surrounding bulk solvent fluid give an effective drag µ + /R + , see equation 73 and 74. To characterize for a lealfet the strength of the intermonolayer slip relative to the traction stress exerted by the surrounding solvent fluid, we introduce the non-dimensional groups Π + 2 = γR + /µ + and Π − 2 = γR − /µ − . For convenience, we also introduce the notation γ ± 0 = µ ± /R ± , so that we can express Π ± 2 = γ/γ ± 0 . We remark that Π ± 2 can be expressed in the more familiar terms of a ratio of rotational drag coefficients. We For a rigid spherical particle subject to torque τ in a fluid with viscosityμ, the angular velocity ω is given by ω = 8πμR 3 + −1 τ , [27]. This shows that the intermonolayer slip contributes similarly to leading order as a bulk solvent fluid of viscosity γR + . We can express similarly Π − 2 . These four non-dimensional groups Π + 1 , Π − 1 , Π + 2 , Π − 2 serve to characterize the relative contributions of the vesicle geometry, shear viscosity within the bilayer leaflets, the shear viscosity of the bulk solvent fluid, and the intermonolayer slip. To simplify our notation, we drop the ± when the same values are used for each leaflet and denote We can non-dimensionalize equations 76-78 by introducing a characteristic velocity v ± 0 and force density f ± 0 . We find it convenient to consider the rigid rotation at angular velocity ω 0 of the spherical shell in the solvent fluid. This motivates the choice of characteristic force density f ± 0 = µ f ω 0 and velocity v ± 0 = Rω 0 . For simplicity, we consider only the case with µ ± = µ f and R ± = R. The non-dimensional velocity is v ± = v ± /v 0 and force densityc ± = c ± /f 0 with coefficientsã andc. With this choice, we can express the non-dimensional problem for the full two-leaflet membrane hydrodynamics in equa- with (78) . As we shall discuss, this analysis will be useful when considering the relative contributions of the solvent traction stress, intramembrane viscosity, and the intermonolayer slip. Other nondimensional scalings can also be considered using a similar approach.

Curvature and Shear
In contrast to a flat membrane, material transported on a curved membrane can undergo shear even by a flow having an effectively constant velocity field on the surface. To investigate the role of intrinsic curvature of the surface, we consider flow on the sphere which has constant positive Gaussian curvature K > 0 and the pseudosphere which has constant negative Gaussian curvature K < 0 [47], see Figure 2.
For concreteness, we parametrize the sphere having Gaussian curvature K = 1 with the coordinates (θ, φ) with x = ψ(θ, φ) = [sin(θ) cos(φ), sin(θ) sin(φ), cos(θ)]. We parametrize the pseudosphere having We first consider a flow having a velocity field v with zero co-variant derivative ∇v = 0 on the surface (constant tangent vector). On both the sphere and pseudo-sphere a velocity having this property is given by v = [− sin(φ), cos(φ), 0]. We remark that it is convenient here to express the velocity in terms of the xyz-components in R 3 given by the embedding from the parametrization above. For a curved surface, this provides the analogue to a flat surface of having a flow with constant velocity. We find that the curvature results in shearing of the transported material. Intuitively, this arises relative to the flat surface by the way intrinsic curvature requires distortion of the distance relationships between points on the surface. More precisely, consider two points located at (θ 1 , φ 0 ) and (θ 2 , φ 0 ) with θ 2 > θ 1 ≥ 0 in the upper hemisphere. While both points travel at exactly the same speed, the point (θ 2 , φ 0 ) which starts closer to the north pole will take less time to traverse fully around the xy-circular cross-section of the surface. This curvature associated distortion of the distances results in shearing of the transported material. This is illustrated in the panel on the left in Figure 2. We can also consider a flow having a velocity field v with dual field v having zero exterior derivative dv = 0 (constant co-tangent vector field). The constant co-tangent case is motivated by the exterior calculus formulation of the fluid equations where for such an incompressible field the flow is determined only from the Gaussian curvature term, see equation 22. We remark that while the co-tangent vector field v = v b dx b is constant on both the sphere and pseudosphere, the velocity field v = v a ∂ x a on each surface is modulated by the local components of the inverse metric factor as v a = g ab v b . For any incompressible velocity field with zero exterior derivative dv = 0, according to equation 22 on any constant Gaussian curvature surface the force density b must also have zero exterior derivative db = 0.
To construct such a flow, we consider for the sphere v = −b /2K = +dφ and for the We see that for both the sphere and pseudosphere the material transported under such a flow is sheared. For the sphere and pseudosphere the shear is expected to be in the opposite direction arising from the difference in sign of the Gaussian curvature K of the surface. This is illustrated in the panel on the right in Figure 2. We remark that similar types of geometry and shear effects can be used for performing rheological experiments as was done in [12].

Particle-Bilayer Coupling : Immersed Boundary Methods for Manifolds
To model the motions of particles within the membrane, we compute a mobility tensor using approximations closely related to the Immersed Boundary Method (IB) [5][6][7][8][9]45]. In IB the fluidstructure interactions are approximated by coupling operators that perform operations on the surrounding flow field to determine the particle velocity and perform operations yielding a force density field to account for particle forces [9,45]. We introduce IB approaches for manifolds to capture both the translational and rotational responses of inclusion particles to applied forces and torques when subject to coupling through the membranes hydrodynamics. We show how our IB approaches can be used to compute an effective mobility tensor for these responses.

Mobility Tensor
We express the mobility tensor M for the velocity response of a collection of particles as The V is the collective vector of velocities and angular velocities of the particles and F is the collective vector of forces and torques applied to the particles. For particle i, the velocity is given by V i = [V] i and the particle force by The response of a single particle to a force applied directly to that particle is given by the diagonal block components M ii . The two-particle response associated with the velocity of particle i in response to a force applied to particle j is given by the off-diagonal block component M ij .

Coupling Operators Γ and Λ for Curved Surfaces
The mobility tensor for the interactions between the i th and j th particle is given by where we have the operators Γ i = Γ X i and Λ j = Λ X j . In the notation, we denote by S the solution operator for the hydrodynamic equations 54. The velocity field for the hydrodynamics v(x) under the applied force density f (x) is given by v = Sf . The operators Γ, Λ approximate the fluid-structure interaction by modelling the velocity response and forces of the particles. The force density generated by an applied force F on particle j is given by f = Λ j F. In response, the velocity V of particle i is given by V = Γ i v. Many choices can be made for the operators Γ and Λ. This can be used for both translational and rotational responses [9]. To ensure that the approximate fluid-structure coupling is non-dissipative, it has been shown the operators must be adjoints [5,9,45]. We require for any choice of field v and vector F that the operators satisfy the adjoint condition The inner-products are defined as where · denotes the dot-procuct in the embedding space R 3 . We use the notation Γ T = Λ to denote succinctly the adjoint condition 82. We remark that from equation 81 this condition has the desirable consequence of yielding a symmetric mobility tensor M when S is symmetric. To obtain the translation and rotational responses of the particles, we introduce the operators The X denotes the collective vector of particle locations. The i th particle is at location [X] i . To obtain the particle velocity in response to the hydrodynamic flow, the tensor W serves to average by sampling and weighing the velocity values on the surface. For a particle force the adjoint tensor W * serves to produce a force density field.
For a curved surface, W must be chosen carefully. A simple form which is widely used in IB is to use a scalar weight W * [F] = η(y − X)F where η is the Peskin δ-function [9,45]. However, for a curved surface this is not a good choice since the force density field produced by ΛF = η(y−X)F is not in the tangent space of a curved surface. Similarly, the averaging procedure Γ may produce a particle velocity which is not in the tangent space.
For curved surfaces, we use a more geometrically motivated operator of the form The sum i runs over the particle indices and the ∂ α | X [i] denotes the tangent basis vector in direction α at location X [i] . The square brackets [i] are used to help distinguish entries not involved in the Einstein conventions of summation. This can be interpreted as the procedure of obtaining the average velocity for particle i by using for each coordinate direction α the inner-product of the velocity field v with the reference vector field q α = w [i],α = w [i],α γ ∂ γ . The adjoint tensor yielding the local force density is given by For translational motions we use the reference vector fields of the form . For rotational responses we use the reference vector field on the surface . We emphasize that we only use these expressions for a coordinate chart chosen so that the particle location X [i] is away from a polar singularity, see Appendix A. The reference velocity fields are illustrated in Figure 3.
In practice for the spatially discretized system, the operators and the associated fields they generate can be expressed conveniently as The index m corresponds to the discrete degrees of freedom, such as the index of a lattice point or harmonic mode, and the i and j index the components of the vector. We have The adjoint condition can be expressed as In practice, we define the operator Λ in numerical calculations using the specified reference velocity fields q α above to generate the force density at lattice sites on the sphere surface. Using the sparse matrix representation of this operation for Λ, equation 94 provides the adjoint velocity averaging operator Γ. This approach for developing consistent operators Λ and Γ on the sphere also extends straight-forwardly to immersed boundary approximations on more general curved surfaces and manifolds.

Dynamics of Inclusion Particles Embedded in Spherical Bilayers
For particle inclusions embedded within spherical lipid bilayer membranes, we investigate their translational and rotational motions in response to applied forces and torques. We consider the case when each embedded inclusion particle only spans one of the fluid bilayer leaflets, see Figure 4. We investigate the mobility of these inclusions when varying the (i) vesicle radius, (ii) membrane viscosity, (iii) solvent viscosity, and (iv) intermonolayer slip. We investigate the four interaction cases (i) outerouter, (ii) outer-inner, (iii) inner-outer, and (iv) inner-inner. We also investigate the coupled motions for the four cases (i) translation-translation, (ii) translation-rotation, (iii) rotation-translation, and (iv) rotation-rotation. We express the translational and rotational responses as where we decompose the mobility tensor into the blocks In the notation, the V denotes the collective translational velocities and ω the collective rotational angular velocities. The F denotes the collective forces applied to particles within the inner and outer leaflets. The τ denotes the collective torques applied to particles within the inner and outer leaflets.
We denote the different ways in which the forces F and torques τ couple to the particle translational motions V and rotational motions ω using the notation M XY . The X denotes the response as either translation (t) or rotational (r). The Y denotes the type of applied force as either standard force (t) or torque (r). The mobility components can be further decomposed into M XY,i 1 , 1 ,i 2 , 2 where i k denotes the location of the i th k particle. The leaflets in which the inclusions are embedded is denoted by k ∈ {inner, outer}. There are a few notable differences between spherical fluid membranes and flat fluid membranes. In the flat case, the membrane domain is often treated as effectively infinite and for theoretical convenience often as having periodic boundary conditions. In the spherical case, the membrane is intrinsically of finite area. Also for a sphere, as consequence of the topology, any in-plane hydrodynamic flow must have a singularity [32]. For the solvent fluid, flat membranes have fluid extending over an infinite domain symmetrically on both sides. In the spherical case, this symmetric is broken with solvent fluid trapped within the interior in a region of finite volume and with solvent fluid extending exterior over an infinite domain. The curvature of the membrane surface can also play an important role in the hydrodynamics. This is particularly apparent from the Gaussian term that appears in equation 54 and the effects we discussed in Section 2.7.
We investigate the mobility of inclusion particles within spherical bilayers in a few different regimes. We consider the characteristic scales introduced in Section 2.6. The regime with Π 1 = L/R 1 and for Π 2 = γR/µ f with Π −1 2 Π 1 = 1 corresponds to the case when the hydrodynamic flow is dominated by the intramembrane viscosity and intermonolayer slip. In this regime, for a force density having a non-zero net torque, the flow is approximated well by the leading order spherical harmonic modes with = 1, see equation 78. The intramembrane viscosity strongly couples the surface fluid resulting in a flow that is a rigid body rotation of the entire spherical shell, see Figure 5. Parameter values are given in Table 2.
Mathematically, this arises from the dominant spherical harmonic modes with degree index = 1 and order index m = −1, 0, 1. Using the exterior calculus formation we apply the generalized curl = − d on the surface to the vector potential Φ given by a linear combination of the harmonic modes of degree = 1. This yields for the velocity field on the surface that of a rigid body rotation, see equation 75. In the case when the surface force has zero net torque in the regime Π 1 1, Π −1 2 Π 1 = 1, the leading order flow is determined by the intramembrane viscosity and intermonolayer slip and depends on the higherorder moments of the torque of degree > 1.
The regime with Π 1 1 and Π 2 1 corresponds to the case when the traction stress from the entrained external solvent fluid dominates the hydrodynamic response relative to the intramembrane viscosity and intermonolayer slip. This results in more localized flow within the surface, see Figure 5.
We remark that the regime when Π 2 1 corresponds to the case when the intermonolayer slip strongly couples the hydrodynamic flow between the two leaflets to make them nearly identical. This effectively doubles the intramembrane viscosity.
We have presented a few regimes indicating the contributing factors in the hydrodynamic responses and the interplay between the entrained solvent fluid, intramembrane viscosity, and intermonolayer slip. We now discuss some features of the hydrodynamic response that arise from the geometry of the spherical membrane. Figure 5: Hydrodynamic flow in response to force acting on an inclusion particle. The L/R = Π 1 is the relative Saffman-Delbrück length-scale scaled by the vesicle radius. For small intramembrane viscosity the force produces a localized hydrodynamic flow on the surface. As the membrane viscosity increases the hydrodynamic flow becomes less localized and eventually approaches the velocity field of a rigid body rotation of the sphere. The flow exhibits two vortices with locations that migrate toward the equatorial poles as the viscosity increases. Parameter values in Table 2. For a force applied to a particle in the outer leaflet located at the north pole, we show as the shear viscosity is varied how the vortex location changes in the outer and inner leaflets. In the nomenclature X-Y in the figure caption, X refers to the leaflet of the applied force and Y the leaflet of the flow response. For low viscosity the vortices are near the north pole θ = 0. As viscosity increases the vortices migrate toward the equator θ = π/2. The inset left to right shows typical progression in flows of the vortex location. The intermonolayer slip Π 2 = γ/γ 0 = 4 moderately couples the inner leaflet to the outer leaflet. We find this results in a flow within the inner leaflet with a vortex location closer toward the equator. Parameter values in Table 2.

Vortices and Membrane Viscosity
As a consequence of the spherical topology of the membrane, any hydrodynamic flow within the surface must contain a singularity [32]. We consider the case of an inclusion particle located at the north pole of the sphere and subjected to a force. These singularities manifest in the flow as two vortices of opposite sign, see Figures 5. The location of these vortices depends on Π 1 = L/R characterizing the relative strength of the intermembrane shear viscosity vs the solvent traction stress. For small Π 1 the vortices start near the north pole and as Π 1 increases they migrate toward the equator, see Figures 5. For a force applied to a particle in either the outer leaflet or inner leaflet we consider how the vortex location changes as the viscosity of the membrane is varied. We show the vortex locations in the outer and inner leaflets in Figure 6. In this case we vary Π 1 = L/R and keep fixed Π 2 = γ/γ 0 where γ 0 = µ f /R with parameters in Table 2. We remark that these results can be used as a reference to estimate the membrane shear viscosity by making observations of the vortex locations of the fluid flow within the leaflets. Some recent experimental work to estimate the membrane viscosity of vesicles using vortex locations can be found in [21,30,58].

Self-Mobility and Coupled-Mobility
We next consider the hydrodynamic responses when a force or torque is applied to an inclusion particle embedded in the outer leaflet when the center of the sphere is held fixed. We take as our convention that this particle is embedded at a pole where we parametrize the sphere with (θ, φ) = 0. We then consider how the resulting hydrodynamic flows within the inner or outer leaflets within the spherical bilayer couple the translational motions and rotational motions of inclusion particles at other locations. Throughout, we use the base-line parameters given in Table 2. These parameters correspond to the nondimensional regime with Π 1 = 0.65 and Π 2 = 4.0.
We investigate the roles played by the bulk solvent fluid, the intramembrane viscosity, and the intermonolayer slip. We use our methods to compute profiles of the mobility responses at different locations when varying the intramembrane viscosity and intermonolayer slip in Figure 8. We show how the mobility varies when changing the intermonolayer slip and membrane viscosity in Figures 9 and 10. For comparison we also compute the mobility responses within a flat membrane shown in Figure 11.
Before discussing in more detail these results, we make a few remarks concerning how the mobility results are reported. The responses are shown along the two great circles on the sphere corresponding to the intersection with the xy-plane and the xz-plane. We consider the hydrodynamic responses when a force or torque is applied to an inclusion particle. For convenience in our calculations, we use by convention the coordinates for the inclusion particle X = (x, y, z) = (1, 0, 0) and we apply force to the inclusion particle in the direction F = (f x , f y , f z ) = (0, 1, 0). To characterize the hydrodynamic responses, we consider the cross-sections of the sphere in the xy-plane and the xz-plane. This gives two great circles of the sphere. We consider the velocity in the directions parallel and perpendicular to the tangents of each of the respective great circles.
We consider the velocity responses in the parallel and perpendicular ⊥ directions along each of these curves. We normalize all of the mobility results by comparing to the case of large intramembrane viscosity Π 1 = L/R = 48 for the leaflet or large intermonolayer slip correspond-ing to Π 2 = γR/µ f = 32. This regime provides a reference case corresponding to the situation when the large intramembrane viscosity yields an effective rigid body rotation of the spherical shell within the bulk solvent fluid or when the two leaflets are tightly coupled.
We remark that this is in contrast to the flat membrane case where the mobility tends to zero as Π 1 = L/R becomes large. In the flat membrane case, we normalize instead our reported results by the self-mobility when Π 1 = L/R = 0.1. Given the mobility model for the hydrodynamic responses discussed in Section 3.1, the self-mobility on the sphere for each type of coupling is given in our model by the results reported at location (θ, φ) = 0.
The mobility profiles reveal a number of interesting aspects of the hydrodynamic coupling between inclusion particles and leaflets. We find that the intermonolayer slip and curvature yield coupling for particles embedded in the inner leaflet significantly different than for particles embedded in the outer leaflet. For a force or torque applied to a particle embedded in the outer leaflet, the intermonolayer slip yields a flow for the inner leaflet with recirculation over a larger scale. This is seen when looking at the vortex locations when applying force at the north pole, where the intermonolayer slip plays a role pushing the vortex location of the inner leaflet closer to the equatorial poles, see Figure 6.
We see this can result in both the translational motions and rotational motions of a particle within the inner leaflet moving in the opposite direction of an inclusion particle within the outer leaflet at the same location. This is seen for the smallest viscosities and intermonolayer slips for the translation-translation responses at location xz with φ = π/4 and for the rotation-rotation responses at location xy with θ = π/4, see Figure 9 and 10. In each case a force or a torque is applied to a single inclusion particle within the outer leaflet located at (θ, φ) = 0. The resulting inclusion particle hydrodynamic response within the outer leaflet or inner leaflet is shown in terms of the mobility M defined in Section 3.1. We use in the nomenclature in the titles of X-Y to indicate a forcing of type X and a response of type Y. We normalize the mobility by the self-mobility response obtained in the case when Π 1 = L/R = 48 and Π 2 = γR/µ f = γ/γ 0 = 32. The intramembrane viscosity or intermonolayer slip is held fixed in panels displaying respectively Π 1 = L/R = 0.13 or Π 2 = γ/γ 0 = 4. For the translational and rotational response to forces in the outer leaflet, we find that the intermonolayer coupling smooths the flow over a larger scale within the inner leaflet.
We next consider for fixed intramembrane viscosity how the intermonolayer slip effects the flow. We see for a force acting on the outer leaflet as the intermonolayer slip becomes small the flow within the inner leaflet approaches a rigid body rotation, see the bottom curve in the upper-right panel of Figure 5.
From an analysis of the hydrodynamic response equation 74, we have two interesting cases for the modes of the inner leaflet: (i) = 1 and (ii) > 1. In the first case, the inner-leaflet rotates as a rigid spherical shell and entrains the fluid trapped within to a rigid body motion. As a consequence there is no traction stress with the external solvent fluid for the inner leaflet and no intramembrane shear stress. This means there are no other stresses acting against the intermonolayer drag so −γ(a − − a + ) = 0 and the inner leaflet matches the outer leaflets rotation with a − = a + for = 1.
In the second case with > 1, the intramembrane stress and traction stress balance the intermonolayer drag. In this case, the hydrodynamic modes of the inner leaflet scale in proportion to the intermonolayer slip and the modes of the outer leaflet. As the intermonolayer slip decreases, the modes a − of the inner leaflet become small for > 1.
This can be seen mathematically from equation 74 where the inner leaflet modes satisfy a − = − A 2 /γ − 1 −1 a + . This can be expressed as where for convenience we denote Π 3 = Π − 2 /Π − 1 . For = 1 this shows a − = a + independent of the magnitude of Π 3 = 0. For > 1, we have as the intermonolayer slip becomes small Π 3 1 the hydrodynamic response for the mode of the inner leaflet with > 1 become small a − 1. This shows that the resulting hydrodynamic responses in the inner leaflet become dominated by the rigid rotation mode = 1 for small intermonolayer slip. This can be seen in the upper-right panel of Figure 8. This has a number of interesting consequences for the motions of inclusion particles embedded within leaflets of spherical bilayers. From the different hydrodynamics of the two spherical shells, we have that for small intermonolayer slip the self-mobility and coupled-mobilities can result in large motions when forces or torques act on an inclusion particle within the inner leaflet. For small intermonolayer slip this arises since the rigid body mode = 1 of the hydrodynamic response for the inner leaflet is not damped by the trapped solvent fluid but only by the weak intermonolayer coupling. This manifests in a near rigid rotation of the inner leaflet and a large self-mobility and coupled-mobility in response to an applied force or torque, see Figure 10.
We remark that it is important to keep in mind this behaviour arises when forces applied to inclusion particles result for the inner leaflet in a force moment with non-zero net torque. This is what drives a significant hydrodynamic response for the rotation mode = 1. In contrast, for the case of a collection of inclusion particles with total force acting on the inner leaflet that yields a zero net torque, super-position of the particle hydrodynamic responses cancel for the = 1 mode and the behaviour of large motions for inclusions from the rigid shell rotation is suppressed. This means for inclusion particles embedded within spherical bilayers it is important to consider carefully the different ways forces and net torque act on the leaflets.
As the intermonolayer slip becomes large, the hydrodynamic flows within each of the two leaflets approach a common velocity. The self-mobility of inclusion particles embedded in the inner and outer leaflet also approach a common value. It is interesting to note that the common value is not simply 1/2 of the self-mobility for the uncou-pled leaflets, see Figure 10. This arises from the asymmetric way in which the leaflets couple to the external solvent. For the outer leaflet the solvent is within an unbounded domain exterior to the spherical shell. For the inner leaflet the solvent is within a bounded domain trapped interior to the spherical shell. As a consequence, we see there are different tractions acting on the inner and outer leaflet, see equation 74. As we saw for the rigid rotation mode = 1, there is no traction stress on the inner leaflet since the solvent fluid rigidly rotates within the spherical shell but there is traction stress from the solvent on the outer leaflet. For the other modes > 1 there continue to be asymmetries in the strength of the traction stress. As a result, the mobility of inclusion particles depend on the particular leaflet in which they are embedded. In the large intermonolayer slip limit, the mobility is determined by a combination of these different solvent tractions from each of the leaflets.
When investigating the mobility of membrane inclusion particles, the finite spatial extent and curved geometry of the bilayer can result in important hydrodynamic effects not captured when treating the membrane as an infinite flat sheet. We remark that the key consideration is how large the spatial extent or curvature is relative to the Saffman-Delbrück length L SD . For a very large vesicle radius or small curvatures, we do expect of course to recover similar behaviours as in the case of an infinite flat sheet. The interesting case is when the vesicle radius or membrane curvature yields a scale comparable or smaller than the Saffman-Delbrück length L SD .
We show the self-mobility of an inclusion particle embedded in a membrane treated as an infinite flat sheet in Figure 11. These results were obtained by solving in Fourier space the hydrodynamic flow in response to an applied force density following closely the analytic approach presented in [43,51] and our method for computing the mobility tensor discussed in Section 3.2. We see significant differences compared to the mobility responses in spherical bilayers.
In the regime of a vesicle radius comparable to the the Saffman-Delbrück length L SD , the finite spatial extent of the membrane and topology can play an important role. For spherical leaflets, it is required that mobility responses result in recirculation flows of the material within the finite leaflet. As we have seen, this can yield non-trivial behaviours in the coupling and provide possibly useful flow features for estimating viscosity as discussed in Section 4.1.
In contrast when treating the membrane as an infinite flat sheet, no vortex arises in the flows generated from single particle responses. The infinite flat sheet also no longer results in trapped fluid within an interior domain. The bulk solvent fluid is treated as occupying an effectively infinite domain on both sides of the bilayer. This results in more traction stress acting on the infinite flat sheet relative to the spherical shell which as a result reduces the self-mobility and strength of the coupled mobilities. In particular, as the intramembrane viscosity increases the rotational mode of the spherical case is no longer available and the self-mobility decays to zero, see Figure 11. Our results show that significant differences can arise with treatment of the bilayer as an infinite flat sheet requiring treatment of the finite domain size and curved geometry of the bilayer when these length scales are comparable to the Saffman-Delbrück length L SD . Figure 9: Membrane Viscosity and Particle Mobility. For a torque applied to a particle embedded within either the inner or outer leaflet, we show as the membrane viscosity is varied the translational and rotational responses of inclusion particles embedded within the inner or outer leaflet. The intermonolayer slip is kept fixed at Π 2 = γ/γ 0 = 4. This is discussed in more detail in Section 4.2.
. Figure 10: Intermonolayer Slip and Particle Mobility. For a torque applied to a particle embedded within either the inner or outer leaflet, we show as the intermonolayer slip is varied the translational and rotational responses of inclusion particles embedded within the inner or outer leaflet. The membrane viscosity is kept fixed at Π 1 = L/R = 0.13. This is discussed in more detail in Section 4.2. Figure 11: Mobility for Flat Membranes. For a force or torque applied to a particle embedded within a large flat membrane, we show as the intramembrane viscosity is varied the translational and rotational responses of inclusion particles. Normalized by the mobility response when Π 1 = L/R = 0.13. These results were obtained from solving in Fourier space the hydrodynamic flow in response to an applied force density following closely the analytic approach presented in [43,51] and our method for computing the mobility tensor discussed in Section 3.2.

Many-Particle Dynamics : Hydrodynamic Coupling and Diffusion
The collective dynamics of multiple particles within a spherical membrane can be modelled as The X denotes the collective particle configuration and F the collective forces acting on the particles. The mobility M is obtained from the hydrodynamic-coupling methods introduced in Section 3.1. The thermal fluctuations driving diffusion are accounted for by the drift k B T ∇ · M and the Gaussian random force F thm (t) which is δ-correlated in time with mean zero and covariance F thm (s)F thm (t) T = 2k B T Mδ(t − s) [5,9,26]. The equations for the particles are to be interpreted in the sense of Ito Calculus [26,42]. The thermal drift term arises from the configuration-dependent correlations of the thermal fluctuations [5,9].
We focus here in this paper on how our approaches can be used for the collective hydrodynamics of particles within the membranes of spherical vesicles. We defer to a future paper the full use of our introduced methods for the entire stochastic Brownian-hydrodynamic model in equation 97. As a demonstration of the introduced methods, we consider the specific case of 4 particles that are actively attracted to a central particle located on the positive x-axis at the east pole. We consider the hydrodynamic flow and particle dynamics within the outer-leaflet of the curved spherical bilayer. In addition to the 4 attracting particles, we also consider 195 passive tracer particles that are advected by the flow, see Figure 12. We see that the hydrodynamic coupling can result in interesting dynamics with the passive particles either moving in the opposite direction of the attracting particles or swept along depending on their relative location. This can be characterized by looking at the coupled mobility M of the passively advected particles defined by M = V /F T . The V is the passive particle velocity, F T the total force acting on the attracted particles. We consider the responses in the circular section in the yz-plane of radius r 0 = 0.5R centred about the x-axis near the east pole and in the circular section in the xz-plane of radius r 0 = R about the center of the sphere, see Figure 12 and Figure 13. The parameters in these calculations are taken to be the same as in Table 2.
We see from the yz-responses M x that for locations half-way between the attracted particles, the passive particle move in the opposite direction to the attracted particles. This change in direction is a consequence of the incompressibility of the fluid which results in an out-flow to compensate for the in-flow toward the east pole generated by the attracted particles, see Figure 14. We consider the hydrodynamic responses of the inclusion particles on two cross-sections of the sphere. The first is the great circle of the sphere when intersected with the xy-plane. The second is the circle of radius r 0 on the sphere surface parallel to the xz-plane. For forces applied to the four attracting inclusion particles, we consider for the motions of the other inclusion particles as characterized by the mobility components parallel and perpendicular to the tangent of the respective cross-section curves. We parametrise the xz-section using angle θ with 0 corresponding to the location (x, y, z) = (1, 0, 0) and the xy-section using angle θ with 0 for location (x, y, z) = (0, 0, r 0 ).
We also see this manifest in the yz-responses M which are out of phase with M x reflecting that the passive particles move laterally toward the out-flow half-way between the attracted particles. The xz-responses correspond to passive particle motions when located on the same great circle in the xz-plane as two of the attracted particles. We see in these responses that the passive particles always move toward the attracting particle at the east pole, see bottom panel of Figure 14.
These results indicate some of the rich dynamics that can arise from hydrodynamic coupling even for relatively simple configurations of particles and force laws. The analytic approaches and computational methods we have introduced for the collective mobility tensor M allow for incorporating such effects into simulations of manyparticle dynamics within spherical lipid bilayers. Many of the approaches we have introduced can also be extended for more general curved bilayers. Figure 14: Multi-particle Mobility. We show the location dependent mobility of the passively advected particles in response to the hydrodynamic coupling to the four attracting particles. We show M = V /F T where V is the particle velocity, F T the total force, M is the mobility tangent along the circular section. The yz indicates the circular section in the yz-plane of radius r 0 = 0.5R about the east pole and xz indicates the circular section in the xz-plane of radius r 0 = R about the sphere center, see Figure 13. In the response, depending on the position, the passive particles either move in the opposite direction or are swept along with the attracting particles. The maximum response M 0 corresponds to the self-mobility of each of the attracting particles.

Conclusions
We have investigated the hydrodynamics of inclusion particles embedded within curved lipid bilayers. We have performed an extensive study of the hydrodynamic flows and mobility responses within spherical bilayers. We have studied both forces and torques applied to inclusion particles embedded within distinct inner and outer leaflets of the bilayer. We have found significant differences relative to the case when the membrane is treated as a single infinite flat sheet. We have found that the interplay between curvature, topology, and the two leaflet bilayer structure can yield interesting effects in the hydrodynamics. We found for spherical bilayers that the difference between the infinite exterior solvent domain and the finite interior solvent domain result in different traction stresses acting on each of the leaflets. We found this can have significant consequences for the mobility of inclusion particles especially when embedded within the inner leaflet. We further found that the intermonolayer slip can play an interesting role with flow regimes where inclusion particles at the same location but in different leaflets can move in opposing directions in response to the lipid flows generated by other inclusion particles. Our results show the rich individual and collective dynamics that can arise for inclusions within bilayers of spherical vesicles. Many of our analytic approaches and computational methods also can be extended to study inclusions embedded within more general curved lipid bilayer membranes.

Acknowledgements
On a spherical surface, coordinate singularities arise, such as exhibited by spherical coordinates at the north-south pole. We handle this issue by using two alternative coordinate charts A or B depending on the location. Chart A corresponds to the spherical coordinates with singularities at the north and south poles. Chart B corresponds to the spherical coordinates with singularities at the west and east poles. A notable feature of the Lebedev sampling is that its symmetry allows us to make use of the same quadrature nodes in the two coordinate charts, see Figure 15. To perform numerical calculations for operations on field values at the Lebedev points, such as divergence, gradient, and curl, we make use whenever possible of the intrinsic geometric meaning of such operations (as opposed to the coordinatecentric formulas). When coordinate-centric formulas are used, we express the operations in one of the two different coordinate charts A or B. The particular chart is chosen to yield significant distance to the coordinate-system singularities. This allows for robust numerical calculations at all points on the sphere surface. This further highlights one of the advantages of our less coordinate-centric exterior calculus approach to the hydrodynamics.
We obtain the membrane velocity field in our calculations by The |g| = det(g) is the determinant of the metric tensor and i is the Levi-Civita tensor (slight abuse of notation). The x denotes the coordinates. For spherical coordinates in chart A, x 1 = θ A , x 2 = φ A and in chart B, To obtain the velocity, we express in each of the charts the coordinate derivatives of the spherical harmonic modes ∂Φ/∂x and the basis vectors ∂ x i in terms of the embedding space R 3 . We then choose at each given location the expression for the chart that has a significant distance to the coordinate-system singularities. In this manner, we compute robustly the velocity field over the entire surface.