Hengdi Zhang

Chaouqi Misbah
email: chaouqi.misbah@univ-grenoble-alpes.fr

Lattice Boltzmann Simulation of Advection-Diffusion of Chemicals and Applications to Blood Flow

Keywords: lattice boltzmann, blood flow, chemical signaling

Diffusion of solutes is often encountered in many biological processes. In the blood micro circulation system, solutes, such as oxygen and calcium molecules, as well as Adenosine Triphosphate (ATP) and biochemical messengers, are released by cells (like red blood and endothelial cells), and are dispersed via dif-

fusion and advection. Moreover, several targeted drug delivery strategies rely on an encapsulation of chemicals and on their release in the blood stream at specific location. The released chemicals couple to blood flow, in which red blood cells (RBCs) constitute the major component. Thus, the development of numerical methods which take into account both dynamics of RBCs and their coupling with chemicals is of great importance for many biomedical applications. We develop here a lattice-Boltzmann based method that deals with generic moving boundary conditions in an advection-diffusion field representing the chemicals. The boundary condition of the solutes at the cell membrane is based on a modified bounce-back scheme. We prove analytically, and validate numerically, that it enjoys second order precision. The solver is validated with

Introduction

In many biological processes cytoplasmic membranes play an important role regarding molecular and ion (active and passive) transport from the cell interior towards the extra cellular environment and vice versa. Examples are abundant in blood circulation. For instance, red blood cells (RBCs) can release oxygen 5 as well as ATP in the microcirculation zone depending, in particular, on the oxygen pressure as well as on the cell membrane shear stress [START_REF] Forsyth | Multiscale approach to link red blood cell dynamics, shear viscosity, and atp release[END_REF]; [START_REF] Zhang | Atp release by red blood cells under flow: Model and simulations[END_REF]. Other examples are encountered in endothelial signaling pathways modulating vasodilation [START_REF] Davies | Flow-mediated endothelial mechanotransduction[END_REF]; [START_REF] Yamamoto | Fluid shear stress activates ca2+ influx into human endothelial cells via p2x4 purinoceptors[END_REF]; [START_REF] Ando | Flow detection and calcium signalling in vascular endothelial cells[END_REF], or in the lymphatic system where calcium 10 dynamics plays a decisive roleZawieja (2005); [START_REF] Jafarnejad | Measurement of shear stress-mediated intracellular calcium dynamics in human dermal lymphatic endothelial cells[END_REF]. Solute transport is also an active field of research in biotechnologies, such as targeted drug delivery which relies on an encapsulation of chemicals within liposomes which are then released at specific sites; the release is triggered either by intrin-sic properties (e.g. local shear stress) or by means of an external stimulus (e.g. ultrasound) [START_REF] Needham | A new temperaturesensitive liposome for use with mild hyperthermia: characterization and testing in a human tumor xenograft model[END_REF]; Allen and Cullis (2013); Akbarzadeh et al. (2013); [START_REF] Sercombe | Advances and challenges of liposome assisted drug delivery[END_REF]; [START_REF] Kaoui | Computer simulations of drug release from a liposome into the bloodstream[END_REF]. All these examples involve an intimate coupling between blood flow and chemical transport. In other words, the chemical species, besides reactions and diffusion, are advected by the flow.

In addition, the chemical species are bound within moving and deformable domains (e.g. RBCs, drug carriers, such as liposomes, etc.). RBCs constitute the major obstacle against solute dispersion in the blood stream. We have thus to cope the coupling between the moving boundaries, the flow field and chemical transport with specific boundary conditions on the membranes of the suspended entities (e.g. liposomes, RBCs) describing the condition under which the solute is released. This task is, in its full generality, quite complex and presents several numerical and conceptual challenges to be described below.

The main purpose of this paper is to develop an advection-diffusion method in the presence of deformable particles (like RBCs or liposomes). This problem will be formulated and solved by means of a lattice Boltzmann method (LBM). The solution of the pure fluid flow by LBM has now become quite classical [START_REF] Zou | On pressure and velocity boundary conditions for the lattice 785 boltzmann bgk model[END_REF]; [START_REF] He | Theory of the lattice boltzmann method: From the boltzmann equation to the lattice boltzmann equation[END_REF]; [START_REF] Chen | Lattice boltzmann method for fluid flows[END_REF]; [START_REF] Succi | The lattice Boltzmann equation: for fluid dynamics and beyond[END_REF]; [START_REF] Mohamad | Lattice Boltzmann method: fundamentals and engineering applications with computer codes[END_REF]; [START_REF] Krüger | Numerical simulations of complex fluid-fluid interface dynamics[END_REF]; [START_REF] Krüger | The Lattice Boltzmann Method[END_REF]. The diffusion/advection problem is now becoming an emerging field of research from both physical and numerical aspects [START_REF] Mohamad | Lattice Boltzmann method: fundamentals and engineering applications with computer codes[END_REF]; [START_REF] Krüger | The Lattice Boltzmann Method[END_REF]; [START_REF] Huang | An immersed boundary method for restricted diffusion with permeable interfaces[END_REF]; [START_REF] Lee | The immersed boundary method for advectionelectrodiffusion with implicit timestepping and local mesh refinement[END_REF]; [START_REF] Markl | Free surface neumann boundary condition for the advection-diffusion lattice boltzmann method[END_REF]; Chen et al. (2013a); [START_REF] Zu | Phase-field-based lattice boltzmann model for incompressible binary fluid systems with density and viscosity contrasts[END_REF]; [START_REF] Fakhari | Improved locality of the phasefield lattice-boltzmann model for immiscible fluids at high density ratios[END_REF]; [START_REF] Liu | Nanoparticle transport in cellular blood flow[END_REF]. To the best of our knowledge, a LBM with general moving boundary with Dirichlet (the boundary concentration is specified) / Neumann (the normal derivative is specified) / Robin (a linear combination of boundary concentration and its normal derivative is specified) condition has not been addressed yet.

We will develop here a LBM for the diffusion-advection problem and couple it with the fluid solver. The formulation of the coupling of the fluid flow to the chemicals adds an extra layer of complexity. A formidable task is to properly handle the chemical boundary conditions on a curved and moving interface (the cell membrane). We will see that both the curvature as well as the moving boundary character pose a challenge. The difficulty arises from the fact that the chemical concentration can be discontinuous at the boundary, unlike the velocity field which is continuous and where the so-called immersed boundary method (IBM) has been successfully applied. Except for some specific situations (the diffusion profile is smooth; the profile has a finite thickness across the membrane) [START_REF] Peskin | The immersed boundary method[END_REF]; [START_REF] Feng | The immersed boundary-lattice boltzmann method for solving fluid-particles interaction problems[END_REF]; [START_REF] Yang | A smoothing technique for discrete delta functions with application to immersed boundary method in moving boundary simulations[END_REF]; [START_REF] Huang | An immersed boundary method for restricted diffusion with permeable interfaces[END_REF]; [START_REF] Lee | The immersed boundary method for advectionelectrodiffusion with implicit timestepping and local mesh refinement[END_REF]; [START_REF] Chen | A conservative scheme for solving coupled surface-bulk convection-diffusion equations with an application to interfacial flows with soluble surfactant[END_REF], the use of IBM for general problems remains to be shown. Another alternative than LBM for treating diffusion has been dealt with recently [START_REF] Liu | Nanoparticle transport in cellular blood flow[END_REF].

In that work a combination of a Langevin equation (for nanoparticles) and LBM (for RBCs) is adopted. When chemical reactions are present, as well as when general boundary conditions are prescribed at the moving boundaries, a generalization of this spirit remains to be done.

A generic type of boundary conditions at the membrane that we will deal with is the so-called Robin boundary condition. This question has been addressed for a static curved boundaryGebäck and Heintz (2014); [START_REF] Zhang | General bounce-back scheme for concentration boundary condition in the lattice-boltzmann method[END_REF]; Chen et al. (2013b); [START_REF] Li | Boundary conditions for thermal lattice boltzmann equation method[END_REF]; [START_REF] Huang | Second-order curved boundary treatments of the lattice boltzmann method for convection-diffusion equations[END_REF]. This problem was addressed in a way which requires an interpolation strategy and local curvature information, making its parallel implementation for arbitrary moving curved boundaries challenging. Huang et al.Huang and Yong (2015) proposed another simplified scheme for general Robin condition on piecewisely linear segments (which are parallel to the mesh segment) with second order precision. The formulation of the boundary condition is derived via asymptotic analysis [START_REF] Junk | Asymptotic analysis of the lattice boltzmann equation[END_REF]; [START_REF] Yoshida | Multiple-relaxation-time lattice boltzmann model for the convection and anisotropic diffusion equation[END_REF], and requires only local information.

Later on [START_REF] Huang | Second-order curved boundary treatments of the lattice boltzmann method for convection-diffusion equations[END_REF] an extension of this scheme for curved boundaries with second order precision has been proposed. However, this scheme turned out to present an interpolation strategy which is shape-dependent and requires information on the local curvature.

In order to circumvent the problem of inefficiency for parallel implementation and challenges raised by some specific shapes, we have developed, by still adopting the general Robin boundary condition scheme from [START_REF] Huang | Boundary conditions of the lattice boltzmann method for convection-diffusion equations[END_REF], a simplified interpolation scheme. Indeed, our scheme takes into account only a single neighboring lattice point to the membrane, instead of several lattice points [START_REF] Huang | Second-order curved boundary treatments of the lattice boltzmann method for convection-diffusion equations[END_REF]. It will be shown here that the simplified version still enjoys the same precision, but at the same time it offers the possibility of efficient (parallel efficiency) and robust handling of arbitrary and moving boundaries.

In addition to the above complexity due to boundary shape, the motion of the boundary in itself raises another problem. Indeed, in the course of time lattice points belonging to a domain lying on one side of the membrane may shift to the other side, and vice versa. Thus, we must identify at each time step the points which have been swept by the moving boundary. This requires an interpolation and / or extrapolation scheme in order to properly reset the values of the concentration field (in the classical fluid-structure interaction problem , this is the so-called refilling procedure [START_REF] Lallemand | Lattice boltzmann method for moving boundaries[END_REF]). Our method uses a similar spirit as the refilling procedure, but introduces some improvement, as will be described here.

Due to the general complexity of the problem we will focus here on a 2D geometry. The implementation of the chemical problem in LBM is based on the so-called D2Q5-BGK model (2 dimensional, 5 velocity and Bhatnagar-Gross-Krook single relaxation time). To efficiently exploit the parallelization benefit of LBM, the implementation is fully based on a graphic processing unit (GPU) parallel architecture, namely Compute Unified Device Architecture (CUDA).

Several tests are performed in order to demonstrate the precision and validity under static and moving boundary conditions.

After having performed several convergence and validation tests, we integrate this solver into a well-validated immersed-boundary coupled Navier-Stokes LBM solver [START_REF] Shen | Interaction and rheology of vesicle suspensions in confined shear flow[END_REF]. We then present two main applications: (i) we study the liposome drug delivery (caused by an external stimulus) problem and analyze the main outcome. Particularly important contributions have been made by [START_REF] Kaoui | Computer simulations of drug release from a liposome into the bloodstream[END_REF] and others [START_REF] Kabacaoglu | Quantification of mixing in vesicle suspensions using numerical simulations in two dimensions[END_REF] who analyzed the problem of drug release and solute mixing. Another impotant study has analyzed advection-diffusion under steady flow for a given shape (a cylinder) Gekle ( 2017). Here we extend these studies to the case where the suspended entities (e.g. model of RBCs) are both deformable (free moving boundaries) and coupled to the advection-diffusion of the solute. (ii) We shall adapt the method to the case where the membrane boundary condition depends on the 115 local shear stress, with the aim to analyze mechano-involved signaling process in micro-circulation. Of particular interest is the problem of ATP release from RBCs that will be briefly presented.

Methods

Since the problem of solving the Navier-Stokes equations by LBM has by now become quite classical [START_REF] Zou | On pressure and velocity boundary conditions for the lattice 785 boltzmann bgk model[END_REF]; [START_REF] He | Theory of the lattice boltzmann method: From the boltzmann equation to the lattice boltzmann equation[END_REF]; [START_REF] Chen | Lattice boltzmann method for fluid flows[END_REF]; [START_REF] Succi | The lattice Boltzmann equation: for fluid dynamics and beyond[END_REF]; [START_REF] Mohamad | Lattice Boltzmann method: fundamentals and engineering applications with computer codes[END_REF]; [START_REF] Krüger | The Lattice Boltzmann Method[END_REF], we will focus on the convection-diffusion problem, and only briefly recall the LBM for the fluid when needed. The proposed method can handle several solutes which are coupled to each other. However, for ease of presentation we will consider a single solute only. In addition, we will restrict ourselves to a two dimensional domain (denoted as 2D), and thus our model of RBC will be a 2D contour made of an inextensible membrane. It will be recognized that a generalization to many solutes is straightforward. Let c denote the concentration of a given solute that depends on space and time. In its full generality c obeys the following equation

∂c ∂t + u • ∇c = ∇ • [D(t, x)∇c] + R(c, t, x) (1)
u is the velocity field (supposed to be known for the moment), t ∈ R + and

x ∈ R 2 are time and spatial coordinates, respectively. For simplicity, we assume here that the diffusion coefficient D is constant. The reaction term R can be handled within the LBM method (in the same way as an external force in the usual fluid solver). A non-overlapping moving boundary (say a RBC model in

2D) curve is explicitly defined as B(t, s) = [X(t, s), Y (t, s)],
here X and Y are Cartesian components of a given membrane point, s is a scalar parameter, which can typically be chosen as the local arc length. By adopting for the fluid flow the non-slip condition assumption at the membrane, the boundary is advanced by the adjacent fluid velocity

∂B ∂t = u(t, B) (2)
If s is chosen as the local arc-length, the normal vector of the boundary is defined as n(t, s) = [-∂Y /∂s, ∂X/∂s], while tangential vector is defined as t(t, s) =

       α + 1 c + + α + 2 n • ∇c + = α + 3 α - 1 c -+ α - 2 n • ∇c -= α - 3 (3)
Here α i 's are real constant. A schematic of the moving boundary and some definitions are shown in Fig. 1a.

Advection-Diffusion Lattice Boltzmann Method

We adopt the so-called D2Q5 BGK model (2 dimensional 5 velocities and a single relaxation time) to formulate a convection-diffusion lattice Boltzmann scheme. Compared to the classical D2Q9 models (used for the fluid solver [START_REF] Zou | On pressure and velocity boundary conditions for the lattice 785 boltzmann bgk model[END_REF]; [START_REF] He | Theory of the lattice boltzmann method: From the boltzmann equation to the lattice boltzmann equation[END_REF]; [START_REF] Chen | Lattice boltzmann method for fluid flows[END_REF]; [START_REF] Succi | The lattice Boltzmann equation: for fluid dynamics and beyond[END_REF]), D2Q5 requires a smaller memory usage and lends itself to an easier treatment for the geometry of the moving boundary, both of which are important for acceleration in a GPU parallel context. In addition, there are numerical evidences that D2Q5 may enjoy a better stability against D2Q9 for the diffusion problem in some particular situations such as at low or intermediate Peclet numbers [START_REF] Li | Lattice boltzmann models for the convection-diffusion equation: D2q5 vs d2q9[END_REF]; [START_REF] Suga | Numerical schemes obtained from lattice boltzmann equations for advection diffusion equations[END_REF].

Let ∆x and ∆t denote the spatial and temporal mesh sizes, the discrete micro velocities are defined as Boundary conditions are reinterpreted on the zigzag (dash-dot) line

c 0 c 1 c 2 c 3 c 4 =     0 1 -1 0 0 0 0 0 1 -1     ∆x ∆t (
This notation means that the velocity c i has its components in the x -y plane given by the ith column of the matrix. We will define a particular speed as c s = (∆x/∆t)/ √ 3 (which would be called the sound speed in the traditional fluid problem, but here it only has a formal analogy). The direction of c i is defined by the unit vector ĉi = c i /|c i |.

Denoting the micro distribution function as g i (t, x), its temporal evolution follows the two main steps:

(i) the collision step

g * i (t, x) = g i (t, x) + 1 τ (g eq i -g i (t, x)) + w i ∆tR (5)
(ii) and the streaming step

g i (t + ∆t, x) = g * i (t, x -c i ∆t) (6)
Here g * i is known as the post-collision distribution function, w i is the weight factor valued as w 0 = 1/3 and w 1,2,3,4 = 1/6, and τ = 3D • (∆t/∆x 2 ) + 1/2 is the dimensionless relaxation time. The equilibrium distribution function is

g eq i = w i c 1 + u • c i c 2 s (7)
which depends on the macro concentration c and velocity u. The relation between the micro distribution function and the macro concentration is simply given by

c = 4 i=0 g i (8)
It is proven via an asymptotic analysis (see [START_REF] Yoshida | Multiple-relaxation-time lattice boltzmann model for the convection and anisotropic diffusion equation[END_REF]; [START_REF] Huang | Boundary conditions of the lattice boltzmann method for convection-diffusion equations[END_REF] and Appendix Appendix A) that Eqs. ( 5)-( 7) converge to the convection-diffusion Eq. ( 1) with a second-order precision when ∆t/∆x 2 ∼ O(1). A brief derivation is provided in Appendix Appendix A. We 170 define the dimensionless diffusivity and velocity (in the lattice Boltzmann units)

as

D = D • ∆t/∆x 2 u = u • ∆t/∆x (9)
The relaxation time τ and D are related by (see Appendix Appendix A)

τ = 3D + 1/2 (10)

General Moving Boundary Condition Treatment

Since the boundary treatments on both sides of the membrane (Fig. 1) are identical from the technical point of view, we only discuss the handling of c -.

Thus below, we omit "±" sign in (3). We will split the general Robin boundary conditions given by Eq. ( 3) into two pieces (and then combine them in the general case). The first one is the Dirichlet condition written as

c = α 3 (α 1 = 1, α 2 = 0) (11)
and the second one is the Neumann condition written as (Neumann) on the curved boundary value problem. ĉi is the unit vector parallel to

175 ∂c ∂n = α 3 (α 1 = 0, α 2 = 1) (12)
c i , ĉī = -ĉ i .

Treatment for Static Curved Boundary

The moving boundary treatment is composed of two steps. The first one is to deal with the shape itself at a given moment (a static boundary), while the second one consists in reconstructing the boundary lattice points when they flip from one side of the membrane to the other due to the motion of the free 180 boundary.

The static boundary treatment follows closely (with some important modifications; see below) that given in Refs. [START_REF] Huang | Boundary conditions of the lattice boltzmann method for convection-diffusion equations[END_REF]; [START_REF] Huang | Second-order curved boundary treatments of the lattice boltzmann method for convection-diffusion equations[END_REF]. The authors there first extended the halfway bounce-back scheme (known for fluid solvers) to the convection-diffusion LBM. This modified scheme 1 b and 2). For that boundary (which will be called hereafter zig-zag boundary) all its intersecting points with the lattice mesh segments are defined at the middle of the mesh segments (dashed-dotted lines in Figs. 1 b and 2; M is one representative point). Since this scheme requires only local information, it lends naturally itself to an efficient parallel implementation.

Recall that in Fig. 1 a physical boundary is discretized into a series of boundary pairs. Assuming there is a boundary lattice point that resides at position x (see Fig. 2, the hollow circle and hollow square provide an example of a typical boundary pair), due to the existence of boundary B, the distribution function g i (t + ∆t, x) cannot be determined from streaming Eq. ( 6). Indeed, if x is the hollow blue circle, then the streaming operation given by Eq. ( 6)

(where the argument of the right hand side is xc i ∆t) would propagate the information from the hollow red square to the hollow blue circle (see Fig. 2).

The streaming procedure will fail due to the presence of the boundary, and thus the determination of the evolution of g i at a boundary point, designated by x, requires a special treatment.

The concentration field (or its normal derivative) is specified at the real boundary (a representative point is P in Fig. 2) to be equal to α 3 (for Dirichlet boundary condition). The discretized boundary has a representative point denoted as M, at which the concentration (still unknown) is denoted as α (mid) 3

.

Below we show how is this value determined from α 3 and the concentration field at point x. Once this task is performed, we have at our disposal the concentration field inside the domain of interest, satisfying the boundary condition. The idea is then to use the Boltzmann equation in order to determine the soughtafter micro-distribution appearing on the left hand side of Eq.( 6) as a function of the post-collision distribution at point (t, x) and the macro concentration field at the discretized boundary (Eq. ( 13)). The proof is given in the Appendix Appendix B and the result is given by (for Dirchlet and Neumann conditions)

         g i (t + ∆t, x) = -g * ī (t, x)+ 1 3 α (mid) 3 for c = α (mid) 3 g i (t + ∆t, x) = 1 + 3u • ĉi 1 -3u • ĉi g * ī (t, x)+ D ∆x (1 -3u • ĉi ) α (mid) 3 for ĉi • ∇c = α (mid) 3 (13)
The streaming equation given by Eq.( 6) is substituted by the above equation for any lattice point x lying next to the boundary. Our analysis presented in Appendix Appendix B shows that the scheme has a second-order precision. This is consistent with the direct numerical estimates [START_REF] Huang | Boundary conditions of the lattice boltzmann method for convection-diffusion equations[END_REF].

Let us now show how to determine α (mid) 3

. An interpolation procedure was developed in [START_REF] Huang | Second-order curved boundary treatments of the lattice boltzmann method for convection-diffusion equations[END_REF] for such a purpose, but it requires information from several lattice points. In addition, the selection of interpolating points is geometry-dependent. A simpler procedure is required for a practical implementation of moving boundaries and parallel computation. Here, we only use a single neighboring lattice point (x + c i ∆t in Fig. 2).

Interpolate α (mid) 3

for Dirichlet condition.

We perform a simple linear interpolation. Let us define the normalized distance p (Fig. 2) between points x and P. The idea is to use the concentration gradient at x at previous time step to linearly interpolate α

α (mid) 3 =        α 3 + (p - 1 2 )∆xĉ i • ∇c(t, x) if p ≤ 1 2 α 3 2p + (1 - 1 2p )c(t, x) if p > 1 2 (14)
This interpolation has a second-order precision. The separation into two cases (p > 0.5 and p < 0.5) is dictated by numerical stability.

Asymptotic analysis shows that the concentration gradient can be reconstructed locally in terms of the micro-distribution with one-order precision procedure Yoshida and Nagaoka (2010); Huang et al. ( 2016)

∇c = 1 ∆tc 2 s uc - 4 i=0 c i g i + O(∆x) (15)
Interpolate α

(mid) 3

for Neumann condition.

Since the zigzag boundary has its normal vector ĉī which is different from n, the reconstruction of α (mid) 3

involves tangential derivative as well. It is easily seen that

α (mid) 3 = ∂c ∂ĉī = ĉī • n ∂c ∂n M + ĉī • t ∂c ∂t M (16)
The normal derivative at M can be obtained from P and first order extrap-olation from x,

∂c ∂n M =            α 3 + (p - 1 2 ) ∂c ∂n x - ∂c ∂n x+ĉi if p ≤ 1 2 α 3 2p + (1 - 1 2p ) • ∂c ∂n x if p > 1 2 (17)
The tangential derivative of c at point M is unknown, so we will express it in terms of the value at x and at the neighboring point x + ĉi . A linear extrapolation yields

∂c ∂t M = ∂c ∂t x + 1 2 ∂c ∂t x - ∂c ∂t x+ĉi (18)
Finally, we consider the general Robin boundary condition, which is a linear combination of Dirichlet (Eq. ( 11)) and Neumann (Eq. ( 12)) boundary conditions. It reads

α 1 c + α 2 ∂c ∂n = α 3 (19)
Similar to [START_REF] Huang | Second-order curved boundary treatments of the lattice boltzmann method for convection-diffusion equations[END_REF], we reduce the Robin condition back to a Neumann problem by approximating the boundary concentration c in the term α 1 c in Eq. ( 19) with its value at previous time step. Since the boundary position is off-lattice (see point P in Fig. 2), the value of c on the boundary is evaluated with concentration value on the nearest boundary lattice and corresponding gradient information (see Eq. ( 15)). Recalling the diffusive scaling assumption mentioned at the end of section 2.1, namely D = ∆t ∆x 2 D ∼ O(1), it is obvious that the introduced error by this approximation merely introduces a second order error term.

face being advanced in the course of time, and (ii) how is the interface being discretized and the boundary conditions implemented. The first point is quite classical: once the velocity field is known then the interface is simply advected by the local velocity (using Euler scheme). In this case the velocity is defined on the fluid lattice, whereas the interface is off-lattice. In order to transfer the fluid velocity information from the lattice to the interface the immersed boundary method is used [START_REF] Peskin | The immersed boundary method[END_REF]; [START_REF] Feng | The immersed boundary-lattice boltzmann method for solving fluid-particles interaction problems[END_REF]; [START_REF] Yang | A smoothing technique for discrete delta functions with application to immersed boundary method in moving boundary simulations[END_REF]; [START_REF] Shen | Interaction and rheology of vesicle suspensions in confined shear flow[END_REF]. This is what we adopt here for the fluid flow only (but not for the solute). For the diffusion problem, we have to specify how is the boundary condition imposed on the interface. The first step is to determine the discretized interface as described in Fig. 1. This procedure is different from the IBM, since our interface is defined as a geometrical one (a sharp interface description) and not as a thin strip as is the case with the IBM . The reason for this treatment is that, unlike the velocity field, the concentration field (and its derivatives) is generically discontinuous at the interface, so that the IBM is inappropriate for handling this situation with enough precision. Some exceptions are observed for a special type of boundary conditions as described in [START_REF] Huang | An immersed boundary method for restricted diffusion with permeable interfaces[END_REF]; [START_REF] Lee | The immersed boundary method for advectionelectrodiffusion with implicit timestepping and local mesh refinement[END_REF]; [START_REF] Chen | A conservative scheme for solving coupled surface-bulk convection-diffusion equations with an application to interfacial flows with soluble surfactant[END_REF]. Once the interface is discretized the boundary terms are evaluated and then transfered to the concentration lattice points in order to deal with the LBM at the interfacial region (see Eq. 13).

After each interface motion (obtained thanks to the fluid velocity) we have to resolve the problem of points flipping from one side of the interface to the other side. A possible type of treatment is the so-called refilling procedure [START_REF] Lallemand | Lattice boltzmann method for moving boundaries[END_REF]. In this method, a point which passes from one side to the other has to be connected to new neighbors located on the same new side. The new distribution function value of the point having passed the interface is evaluated as an extrapolation from its new neighbors. Here we will use a similar strategy but in a more refined manner, as described below. Our procedure is split into

g k = k ∈S (ref ) w k (2g k (t, x + c k ) -g k (t, x + 2c k )) k ∈S (ref ) w k (20)
Here S (ref ) is a collection of subscripts defined on a lattice point , which flipped from Ω + to Ω -in the new time step (see Fig. 3). An index k belongs to 2), and deal with the new boundary condition with the p value (even if it is larger than 1). (c) Search for neighboring lattice points in the set S (ref ) , and then calculate the new distribution function value via Eq. ( 20) depending on applications. The method developed here lends itself to 320 other boundary conditions. For example, other boundary conditions for the solute at the external boundaries could be used as well, such as the Robin ones, and so on. Here n is the normal vector pointing from Ω -to Ω + . We tested numerically both boundary conditions.

S (
The choice of relaxation time τ in simulation is quite flexible. As we have seen in Eq. ( 9) and Eq. ( 10), once the lattice mesh size ∆x and relaxation time τ are fixed, the total number of the simulation steps is given by

N (step) = T ∆t = T D a ∆x 2 D = 3D a ∆x(τ -0.5) • T
For definiteness, we fixed T = 1 (in an arbitrary unit). Simulations with different values of relaxation time were tested, namely τ = 0.625, 0.75 and 1. The numerical results are then compared to the analytical solution c (ana) at this particular time T = 1. We define the numerical relative error of concentration c (num) as

E(∆x) = ||c (num) -c (ana) || ||c (ana) || (

An advection diffusion problem in an irregular domain

Secondly, we test this boundary condition scheme against an advection diffusion problem in an irregular domain, which has an analytical solution [START_REF] Huang | Second-order curved boundary treatments of the lattice boltzmann method for convection-diffusion equations[END_REF]. The computational domain is defined as Ω = {(x, y) ∈ R 2 |φ(x, y) < 0}, here φ(x, y) is a scalar field which reads φ(x, y) = x 4 -5x 2 -3x + 2y 4 -6y 3 -y -1

The time dependent advection diffusion problem is governed by

∂ ∂t c(t, x, y) + ∇ • (uc) = ∇ • (D∇c) + S, (x, y) ∈ Ω, , t ∈ [0, T ]
The velocity field is set to be a constant, u = (1, 0) and D is a constant scalar.

The time dependent source term S is

S = (x 3 + y 3 )ω cos(ωt) + [3x 2 -D(6x + 2)] sin(ωt)
The initial condition is given as

c| t=0 = 0, (x, y) ∈ Ω
The boundary condition, which is defined on the zeros of φ(x, y), can be defined as either Neumann type or Robin type. For Neumann condition, it reads as ∂c ∂n φ=0 = (3n x x 2 + 2n y y) sin(ωt)

For Robin condition, it is given by c + ∂c ∂n φ=0 = (x 3 + y 2 + 3n x x 2 + 2n y y) sin(ωt)

Here n x and n y are components of the normal vector pointing outward. This several neighboring points and boundary curvature information are employed).

The relative error against the different choices of relaxation time τ (and thus against diffusivity, see eq. 10) is gathered in Fig. 6. This result suggests that a choice for the relaxation time in the range τ ≤ 1 (see Fig. 6) is an appropriate one. Our numerical experiments revealed that this choice favors stability. Indeed, we found that for both Neumann and Robin type boundary condition, simulations with a large relaxation time led to numerical instability (e.g. the choice τ = 3 and N = 16 led to divergence). We identified that the instability results from extrapolation of tangential gradient in Eq. ( 18). Note that omitting the extrapolated term can make the simulation more stable with large relaxation time but as a counterpart the precision of the boundary treatment will degrade down to first-order.

Moving Boundary Validation

Now we perform some validation tests on a problem with moving boundaries. The first example is a cylinder (which is a circle in 2D, which encloses a solute) which is advected by a constant velocity. The second case consists of a deformable opaque surface that encloses a solute and which is distorted by a 360 rotational velocity field.

An Advected Leaking Rigid Reservoir

Consider a cylindrical reservoir which is advected by a constant velocity u 0

in a periodic box Ω = [-a/2, a/2] × [-b/2, b/2].
The boundary of the advected reservoir is described as

B(t, θ) = [r 0 cos(θ) + u 0x t, -r 0 sin(θ) + u 0y t] , θ ∈ [0, 2π)
where r 0 is the radius of the cylinder. A solute concentration field c is defined in Ω and obeys

∂ ∂t c(t, x, y) + ∇ • (u 0 c) = ∇ • (D∇c), (x, y) ∈ Ω, t ∈ [0, T ]
The initial condition is set as The boundary concentration values, c -and c + , are evaluated thanks to the concentration values and their gradients at the lattice points (see Eq. ( 15)).

c ini =        1 inside
Since u 0 is constant, one can infer from Galilean invariance that the actual concentration field in the case where u = [a/T, 0] should be identical to that 370 obtained in the case with zero advection velocity. However, the numerical errors may affect the Galilean invariance. Therefore, it is important to check if this property is preserved after long simulation time (see the comparison of the two results at T = 1 in Fig. 7). τ is fixed to 1 in these simulations. The numerical solution for ∆x = b/1024 and u 0 ≡ 0 is used as a reference solution from which 375 the error is calculated when non zero velocities are considered. A second-order convergence is obtained (Fig. 8).

A Reservoir that undergoes rotation and distortion

In this case we consider a divergence-free rotational velocity field u = [u x , u y ] (see the vector field in Fig. 9 order Runge Kutta method, under which, the numerical error of reservoir area is adequately suppressed. A numerical solution (obtained with via ∆x = 1/512 and τ = 0.75) is presented in Fig. 9. We measure the error at T = 1, when the shape is still easily resolvable for a large enough mesh size ∆x = 1/16.

The finest mesh ∆x = 1/2048 is used to estimate the relative error. (21). A convergence rate between first and second order (around 1.45, see Fig. 10) is observed. We attribute this to the loss of resolving geometry information when the mesh size is not small enough. For example, in the case where ∆x = 1/32, in of precision (from the expected order 2 down to order 1.5) the code robustly handles extreme distortions, such as the situation shown at T = 8 in Fig. 9.

The Arteriole Thrombosis Drug Delivery Simulation

The validations (section 3) have proven that the numerical method's capability of dealing with advection, rotation and distortion of a moving boundary, with even quite ample deformation patterns under flow. We now implement this solver for the study of a first practical example, the arteriole thrombosis drug delivery. This example is inspired by previous numerical studies on simplified models for liposome drug release in channel flow Gekle (2017); [START_REF] Kaoui | Computer simulations of drug release from a liposome into the bloodstream[END_REF].

We consider a liposome that encapsulate some drug to be released in a blood vessel. We take here into account both the plasma fluid and the RBCs. We follow below closely the work of [START_REF] Kaoui | Computer simulations of drug release from a liposome into the bloodstream[END_REF].

Liposomes are closed membranes, and are considered as promising targeted drug carriers. They have typical diameters ranging from 100 nm to 1µm [START_REF] Sercombe | Advances and challenges of liposome assisted drug delivery[END_REF]. They are suitable means for carrying hydrophilic drugs to a particular location in the organism where they may release their content thanks to an external stimulus (for example ultrasound excitation) or when the surrounding fluid shear stress reaches a critical value (shear stress is the highest in arterioles). Although there exist various types of liposomes depending on the precise purpose, we consider it here to be composed of a single lipid bilayer. Actually in 2D and for an incompressible membrane there is only a single mode of deformation (bending), and there is no real distinction between different models (vesicle, capsule, red blood cell, etc...).

We consider a straight long channel for the arteriole model, with a Gaussian hill shaped obstacle on one side of the vessel wall to represent the thrombosis (Fig. 11). We adopt a 2D vesicle as a crude model for RBCs (we will nev- We would like to evaluate how does the presence of RBCs impact the drug delivery process. The role of RBCs is two-fold, they affect the flow field, and they serve as solute obstacles. To elucidate the role of each effect we consider 3 different situations, one without RBCs at all, the second one with passive RBCc (i.e. RBCs are completely transparent to the drug), and finally a realistic case in which the RBCs are opaque to the drug. In all the three cases we consider a When we consider a transparent RBC model, we relax the first of the two above boundary conditions. In that case the RBCs only affect the flow pattern, which 455 in turn affects the advection of the solute, but the RBCs do not affect directly the solute (no solute obstacle).

ertheless

Liposome RBCs

Preliminaries on the fluid-membrane-solute system

We first briefly recall the modeling of 2D membrane problem, which serves to model RBCs and the drug carrier as well. A 2D unstretchable closed membrane is used for the RBC model, and has a bending elastic modulus [START_REF] Betz | Atp-dependent mechanics of red blood cells[END_REF]; [START_REF] Kaoui | Complexity of vesicle microcirculation[END_REF]. The membrane force is obtained from the Helfrich energy

κ b = 3 × 10 -19 J
H(X(t, s)) = κ b 2 κ 2 ds + ζds
Here X(t, s) is the position of RBC membrane, κ is the local curvature, ζ is a Lagrange multiplier that enforces local membrane inextensibility, and s is the curvilinear coordinate. The membrane force (which is obtained as a functional derivative of H with respect to X(t, s)) is applied as a bulk force (albeit localized close to the membrane, in an immersed boundary spirit) to the fluid

ρ ∂u ∂t + u • ∇u = -∇p + ∇ • η(∇u + ∇u T ) + δH δX • δ(x -X)ds,
where the supescript T stands for the transposition operation. The membrane force (which is the functional derivative with respect to the membrane 460 position) is given by [START_REF] Kaoui | Lateral migration of a two-dimensional vesicle in unbounded poiseuille flow[END_REF]

f mem→flu = κ d 2 κ ds 2 + 1 2 κ 3 n -ζκn + dζ ds t (23)
Then, the velocity field u is used to advance the membrane shape

∂X ∂t = δ(x -X)u(x)dxdy
We define the Capillary number as Ca = µ ex γw R 3 0 /κ b , where η is the fluid viscosity (which is position dependent, so that it can be taken to be different inside and outside of the RBC, if need be; see later). γw is the flow shear rate at channel wall in the absence of thrombosis, liposome and RBC, R 0 is a characteristic radius of the RBC, R 0 = A/π. A reduced RBC area is defined as RA = 4πA/P 2 , A and P are constant area and perimeter, respectively.

RA is an index that quantifies the roundness of the RBC, and taken here as RA = 0.7. For the RBC model we take R 0 = 3µm, whereas for the drug carrier (liposome, for which we adopt the same membrane model as the RBC) we take RA (lip) ≈ 0.997 (a quasi-circular shape) and a radius R

y = h t • exp - 1 2 x -L/2 w t 2
Here h t = 0.5 • W and w t = W are the height and the width of the shape.

Several studies [START_REF] Zhao | Shear-induced particle migration and margination in a cellular suspension[END_REF]; [START_REF] Kumar | Mechanism of margination in confined flows of blood and other multicomponent suspensions[END_REF]; Gekle (2016); [START_REF] Müller | Understanding particle margination in blood flow-a step toward optimized drug delivery systems[END_REF]; [START_REF] Krüger | Effect of tube diameter and capillary number on platelet margination and near-wall dynamics[END_REF]; [START_REF] Guckenberger | A boundary integral method with volume-changing 54 objects for ultrasound-triggered margination of microbubbles[END_REF] have shown that due to the small liposomes size, their rigidity and flow conditions, the margination (caused by RBCs) effect will drive them to the so called cell-free layer (CFL) near the vessel wall. Our long term simulation also confirmed this tendency. However, in the absence of RBCs, the liposome remains at its initial lateral position while being advected by the flow along the channel. In order to be able to compare the results with those obtained in the presence of RBCs the CFL (as in the presence of RBCs), 1.25R 0 distant from the upper vessel wall.

We use the following criterion for drug release: the liposome is impermeable to solute until its centroid reaches a given distance from the thrombosis, below which it becomes transparent. In these simulations we have set that distance to be x (permeation) = -17.5R 0 . We have in mind the situation where the liposome develops small pores that allows solute release, whereas the liposome maintain its overall membrane integrity.

We set ∆x = 0.15 µm and P e = max(|u|)R

(ves) 0

/D = 10 and 100, which corresponds to typical drug solute diffusivities D ( 10 -10 and 10 -11 m 2 /s).

Prior to the study of interest (0 value Dirichlet condition in Eq.( 22), we have tested our code with zero-flux boundary condition on vessel walls in order to check numerically the mass conservation. We have found that loss of mass is always smaller than 1.5%.

Results

The solute absorption rate by vessel walls is shown for different cases in Fig.

12 (a). The normalized absorption rate is defined as

R(t) = 1 - c(t)dxdy c(0)dxdy ( 24
)
The integration is treated as a summation among all pixels, which introduce an error term of order 1% only.

It is found that for both P e = 10 and 100, whether we have transparent RBCs or no RBcs at all (solid lines and dashed lines in Fig. 12 (a)) the absorption rate is practically the same. We believe this is due to the liposome's lateral , where the normalized distance between liposome centroid and the wall y (lpsm) (t) is always small enough, so that the absorption curve is diffusion-dominated, rather than advection-dominated.

However, when RBCs are opaque, we find an increased absorption, shown by dotted lines in Fig. 12 (a)). A close inspection suggests that RBCs, as solute barrier, obstruct the diffusion of solute from liposome to the channel center.

Consequently the gradient is increased towards vessel wall; it is proportional to the normalized absorption rate along the upper wall:

f absp (t, x) = - D γR 0 ∂c ∂y y=3.75R0 (25)
Figure 13 shows the absorption rate together with RBCs configuration and the solute pattern at the thrombus location. Here again, we notice the higher absorption rate due to the presence of RBCs playing the role of solute barriers.

490

In conclusion, when having high enough Hematocrit (31.3% in this study)

and high enough rigidity of liposome which lead to margination, the presence of RBCs facilitates the vessel wall absorption process. This is mainly due to the impermeability of RBC membrane to the solute. One can speculate that increasing Hematocrit will strengthen the absorption, although further investigations are needed before reaching a definite answer.

ATP release from Red Blood Cells under Flow

The study of mechanosensing of solute release can be handled by coupling the boundary conditions with other information such as the membrane shear stress. We apply here our solver to the problem of release of ATP (Adenosine

D ∂c ∂n =        γw /P if σ > σ 0 0 if σ ≤ σ 0
Recall that P is the perimeter of the RBC, γw the wall shear rate, defined in section 4. σ is the shear stress along the membrane, defined as σ = η out ∂u t ∂n where u t is the fluid velocity along the membrane tangent.

500

The critical value of the shear stress is set to σ 0 = 0.07P a, and is inspired by reduces the ATP release. In other words, the TB regime triggers less release than the TT one. This finding is consistent with experimental observations [START_REF] Forsyth | Multiscale approach to link red blood cell dynamics, shear viscosity, and atp release[END_REF]. Actually the problem is more subtle, since not only the local shear stress matters, but also the cell deformation amplitude. By taking into account both effects, we have been able to reproduce both qualitatively and 510 quantitatively the experimental results. An extensive study is devoted to this question in a different publication [START_REF] Zhang | Atp release by red blood cells under flow: Model and simulations[END_REF].

Conclusion

In this work, we developed a 2D lattice-Boltzmann based advection-diffusion solver for curved moving boundaries for both Dirichlet, Neumann and linear 515 Robbin boundary conditions. In most cases, a second order convergence is achieved. For highly distorted boundaries, a convergence between first and second order is observed. The boundary scheme requires only one neighboring than 3 we obtain g i (t + a 2 , x + ĉi ) -g i = =w i c (2) + 1 2 aĉ i • uc (1) -τ (ĉ i • ∇)g

(1)

i -τ a

∂ ∂t + 1 2 (ĉ i • ∇) 2 g (0) i g (3) i =w i c (3) + 1 2 aĉ i • uc (2) -τ (ĉ i • ∇)g (2) i -τ a ∂ ∂t + 1 2 (ĉ i • ∇) 2 g (1)
i -τ a(ĉ i • ∇)

∂ ∂t + 1 6 (ĉ i • ∇) 3 g (0) i (A.9)
By summing over i (from 0 to 4) the first equation in (A.9), and using Eq.(A.8) and the fact that the weight factors obey 4 i=0 w i = 1, we find that the resulting equation is automatically satisfied. Performing the same operation with the second equation yields the same conclusion, by using Eq.(A.8) and the first equation in (A.9), and by virtue of the fact that tion function (defined by Eq. ( 5)), but defined for opposite micro velocities; the subscript ī means that the corresponding micro velocity direction is pointing outward (see ĉī in Fig. 2). The classical bounce-back boundary condition is given by g i (t + ∆t, x) = g * ī (t, x) for zero-flux boundary condition in the pure diffusion problem (the relation is analogous to the non-slip boundary condition in Navier-Stokes LBM). This is quite intuitive, since from time t to t + 1 the number of particle crossing the boundary from each side is equal and opposite.

The question is how to extend this relation to the present problem.

To deal with this question, we use the same asymptotic expansion in powers if , as performed in the preceding appendix, for g * ī (t, x) and g i (t + ∆t, x). With the help of the post-collision and equilibrium distribution function definition ( Eq. ( 5) and ( 7)) and the asymptotic expansion (equation (A.9) ), we find for g * ī (t, x) and g i (t + ∆t, x) the following expressions up second order in g * ī (t, x) =gī(t, x) + 1 τ (g eq ī -gī)

=w i c (0) + w i -(3aĉ i • u -(τ -1)(ĉ i • ∇))c (0) + c (1) + 2 w i (τ -1)[(τ - 1 2 )(ĉ i • ∇) 2 -3aĉ i • u(ĉ i • ∇) -a ∂ ∂t ]c (0) -(3aĉ i • u -(τ -1)(ĉ i • ∇))c (1) + c (2 O( 3 ) g i (t + ∆t, x) =g i (t, x) + a 2 ∂ ∂t g i + O( 3 ) =w i c (0) + w i (3aĉ i • u -(τ -1)(ĉ i • ∇))c (0) + c (1) + 2 w i [τ (τ - 1 2 )(ĉ i • ∇) 2 -τ 3aĉ i • u(ĉ i • ∇) -(τ -1)a ∂ ∂t ]c (0) + (3aĉ i • u -τ (ĉ i • ∇))c (1) + c (2) + O( 3 ) (B.2)
Here the symmetrical structure of the micro velocities ĉi + ĉī ≡ 0 has been exploited. A close inspection of the above equations and the boundary con-620 ditions (B.1) allow one to infer the proper writing of the streaming process.

Indeed, summation and subtraction of g i (t + ∆t, x) and g * ī (t, x) provide interesting results if one focuses on zeroth and first order terms only. The results are given by g * ī (t, x) + g i (t + ∆t, x)

2w i = c (0) -( 1 2 (ĉ i • ∇)c (0) -c (1) ) + O( 2 )
g * ī (t, x) -g i (t + ∆t, x)

2w i =3aĉ i • u c (0) -( 1 2 (ĉ i • ∇)c (0) -c (1) ) - (τ - 1 2 )(ĉ i • ∇) c (0) -( 1 2 (ĉ i • ∇)c (0) -c (

[

  -∂X/∂s, -∂Y /∂s] The solute concentration c and the corresponding flux J = uc -D∇c are, in general, discontinuous at the boundary. By defining c ± (B) = lim →0 ± c(B + n), (c ± in short), the general Robin boundary condition (an equation which combines both Neumann and Dirichlet conditions) along the two sides of the moving boundary B can be written as

Figure 1 :

 1 Figure 1: A schematic of the convection-diffusion problem in the presence of a moving boundary: a) physical boundary B(t, s) is piecewisely continuous with countable discontinuities in its derivative. Boundary conditions expressed by Eq. (3) are imposed along each side of the curve. Arrows along the curve show the monotonously increasing direction of s. b) Geometrical information of the physical boundary which is discretized into a series of boundary pairs (a square and a circle which are located on the nearest lattice points from the boundary).

Figure 2 :

 2 Figure2: A typical boundary lattice is highlighted as the hollow circle. g i (t, x) is to be calculated from boundary condition (see text). The curved boundary needs to be interpolated onto the zigzag boundary (dashed-doted line). This process is done via finding a value α (mid) 3

  concentration field (or its derivative) at the real boundary (where c = α 3 for a Dirichlet condition) and at x. The concentration at point M and at time t + ∆t is given by

  Fig.3with a red dashed line) and evaluate the boundary conditions according

  ref ) only if both x + ĉk and x + 2ĉ k belongs to Ω -(t) and Ω -(t + ∆t) in both previous and present time step. We have postulated equation (20) from linear extrapolations and a set of intuitive weight factors as defined in section 2.1. We have successfully tested these weight factors by a number of numerical experiments (examples of validation will follow). Similar extrapolations (for rigid moving boundaries) were adopted for fluid and solute problems in Yin et al. (2012); Chen et al. (2013a).

Finally, note thatFigure 3 :

 3 Figure 3: The procedures dealing with the moving boundary problem: (a) at that time step, we treat the boundary condition as a static boundary. (b) at intermediate time step, keep the boundary grids fixed, calculate the intersection length p due to boundary movement (defined in Fig. 2), and deal with the new boundary condition with the p value (even if it is larger

  Validations on time-dependent problems with static boundaries 3253.1.1. A pure diffusion problem in an irregular domainWe first validate our boundary condition scheme by treating a time-dependent diffusion problem in the presence of a static boundary. This problem has an analytical solution[START_REF] Hu | A coupled immersed boundary and immersed interface method for interfacial flows with soluble surfactant[END_REF]. A periodic square domain is defined asΩ = [-1, 1] × [-1, 1],with the presence of an internal boundaryB = [r(θ) cos(θ), r(θ) sin(θ)], θ ∈ [0, 2π)here r(θ) = 0.4 + 0.1 cos(3θ). As the internal boundary is closed, it divides Ω into Ω -and Ω + , which are the external and internal domains respectively. We focus here on the solution in Ω -. The analytical solution is given by Hu et al. (2018) c (ana) (t, x, y) = 1 + 0.5 exp -2π 2 tD a cos(πx) cos(πy) for (x, y) ∈ Ω - Here D a is a dimensionless diffusion coefficient. For definitness, its value is set here to 0.1. The above solution is valid for both Dirichlet and Neumann conditions The above concentration field satisfies the pure diffusion equation with initial condition c| t=0 = c (ana) (0, x, y) and a corresponding Dirichlet type boundary condition c| B = c (ana) or a Neumann type boundary condition ∂c ∂n B = ∂ n c (ana)

Figure 4 :

 4 Figure 4: The relative error against the mesh size in the static three-leaves problem from Hu et al. (2018).

  problem has an analytical solution c = (x 3 + y 2 ) sin(ωt) By following[START_REF] Huang | Second-order curved boundary treatments of the lattice boltzmann method for convection-diffusion equations[END_REF], we set D = 1 and ω = 1. A computational domain[-3, 3] × [-2, 4] is uniformly meshed into N × N lattices. For the Robin condition, the value of c on the boundary is obtained by the values at the boundary lattices and the gradient information (see Eq. (15)) at previous time step. By choosing different τ and mesh size ∆x = 6/N , the numerical results at T = 1 are compared with the analytical solution. A second order convergence is observed (see Fig.5). Despite the fact that there is only a single neighboring point to a given boundary lattice point and no curvature information (recall Fig.

  340

2Figure 5 :

 5 Figure 5: The relative error against the mesh size in the advection diffusion problem in an irregular domain from Huang et al. (2016).

Figure 6 :

 6 Figure 6: The relative error against relaxation time τ in the advection diffusion problem in an irregular domain from Huang et al. (2016). The lattice mesh size is fixed as N = 128.

  the internal and external sides of B are denoted as c -(t, θ) and c + (t, θ) respectively. The normal vector points into the outward direction. The solute is leaking from the reservoir into the external zone, and this leakage is described by the following boundary condition at the cylinder surface D ∂c - ∂n = -D ∂c + ∂n = k(c --c + ) Parameters are fixed as r 0 = 0.8, a = 4, b = 2, D = 1, and k = 1. The concentration fields computed with two different constant velocities u 0 = [a/T, 0] or [0, 0] are compared at T = 1. The contour plots of the numerical 365 solution with ∆x = b/512 is shown in Fig.7

Figure 7 :Figure 8 :

 78 Figure 7: Snapshots for the concentration field for u = [a/T, 0] (left side) and u = [0, 0] (right side) respectively. In both simulations, we have set ∆x = b/512 and τ = 1. Since the only difference between the left side simulation and the right side one is that the concentration field is advected by a constant velocity, it should be identical to the concentration field on the right side up to a translation (Galilean invariance). This property is accurately reproduced in this simulation, although the LBM is based on an Euclidian mesh. The preservation of Galilean invariance is crucial to scenarios with long term advection such as cell membrane flowing in channels of realistic length.

  colored in gray) in a simulation box Ω = [-1cos(πy)] sin(πx) A localized circular reservoir is initially defined with its centroid at x 0 = [0.5, 0] with radius r 0 = 0.3. Its boundary is denoted as B. The reservoir is then subjected to a rotating velocity field, thus, the motion of the reservoir boundary is expressed as ∂B/∂t = u. A concentration field c is defined on this simulation box, and is governed by the advection-diffusion equation A zeroflux Neumann condition is imposed on both sides of the reservoir boundary B. A representative Dirac delta function is used as an initial condition c 0 = δ(x -x 0 ). More precisely, we consider a regularized Dirac function in the form of Gaussian function with a width of order of two lattice points. We set the Peclet number to P e = max (|u|)/D = 50. The imposed velocity field is expected to extremely elongate and distort the reservoir. This extreme distortion will constitute an interesting test of the code robustness. The advancing of the distorted reservoir boundary is conducted by means of a 4th

Fig. 10 ,Figure 9 :Figure 10 :

 10910 Fig.10, the maximum width of the shape is represented by only 7 lattice points,

  use the abbreviation RBC in what follows). An immersed boundary coupled lattice Boltzmann method is employed to resolve the fluid-membrane coupled system. This means, that the motion of the RBCs is performed by the classical immersed boundary method (unlike the treatment of the diffusion problem)[START_REF] Peskin | The immersed boundary method[END_REF];[START_REF] Feng | The immersed boundary-lattice boltzmann method for solving fluid-particles interaction problems[END_REF];[START_REF] Yang | A smoothing technique for discrete delta functions with application to immersed boundary method in moving boundary simulations[END_REF]; see[START_REF] Shen | Interaction and rheology of vesicle suspensions in confined shear flow[END_REF] for details of our fluid solver. A liposome is initially located at the upstream of a thrombosis site, surrounded by RBCs. The liposome encapsulates a water-soluble drug with concentration c 0 = 1. Its membrane is initially impermeable to the drug. Then the liposome is advected by blood flow until it reaches a particular distance from the thrombosis, where we assume that the liposome membrane becomes suddenly absolutely permeable to solute (mimicking the effect of an external stimulus). A schematic representation is given in Fig.11to show the simulation layout.

Figure 11 :

 11 Figure 11: A schematic representation of liposome drug delivery process near an arteriole thrombosis. The top panel represents some time before drug release, and the bottom one represents the situation after drug release. The colormap represents the concentration of drug solute. The liposome becomes transparent after it reaches a particular distance from the thrombosis

  8µm. We set Ca = 10, meaning that the flow strength is large enough to strongly deform RBCs. This will allow us to test the robustness of the code. In the whole simulation we keep the Reynolds number Re = max(|u|)R (ves) 0 /η ≈ 0.2. The channel width and length are W = 7.5R (ves) 0 and L = 75R 0 , thus the simulation box is Ω (vessel) = [-37.5R 0 , 37.5R 0 ] × [-3.75R 0 , 3.75R 0 ]. The Hematocrit value 31.3% is taken for RBCs, which corresponds to 56 representing cells in the channel. The thrombosis is shaped as

Figure 12 :b

 12 Figure 12: (a) The solute absorption rate as a function of time for P e = 10 and 100. (b)The normalized lateral position of the liposome as a function of time with or without RBCs' presence. In both cases the liposome is close to the wall. The vessel wall is at position 0.5 and the center-line at position 0.

  Triphosphate) from RBCs under shear flow. Experimental studies have shown that when RBCs undergo flow originated shear stress, they can release ATP, a universal energy carrier and messenger molecule that plays important role in vessel dilatation[START_REF] Wan | Dynamics of shear-induced atp release from red blood cells[END_REF];[START_REF] Forsyth | Multiscale approach to link red blood cell dynamics, shear viscosity, and atp release[END_REF]. It has been shown in[START_REF] Forsyth | Multiscale approach to link red blood cell dynamics, shear viscosity, and atp release[END_REF] that the amount of ATP release depends on whether RBC undergo tank-treading (TT) or tumbling (TB). In our 2D model a transition from TT to TB can be achieved by increasing the viscosity contrast between the internal and the external fluids[START_REF] Vlahovska | Vesicles and red blood cells: from individual dynamics to rheology[END_REF]. We define the viscosity contrast as λ = η in /η out , here η in and η out are the viscosity of internal and external domain of the RBC. A shear flow is generated in a simulation box [W, L] = [12.5R 0 , 25R 0 ] by moving the upper and lower walls with the same speed but opposite direction. The solute (ATP) boundary condition on the upper and the lower walls correspond to zero-flux condition, with periodic boundary conditions along the flow direction. The RBC is initially located at the center of the shear flow. Experiments have shown that viscosity (as well as the shear rate) contrast is a key parameter to the motion of RBCs under shear flow Fischer et al. (1978); Tsubota and Wada (2010); Fischer and Korzeniewski (2013). For low enough λ (to be specified below), a vesicle (or 2D RBC) tends to exhibit TT motion, whereas for higher λ TB prevails. Two typical values of λ are shown in Fig. 14. We adopt a simplified model of the ATP release, which is a Heaviside step function for shear stress condition (namely, the release takes place only if the local membrane shear stress exceeds a certain value). Thus, the Neumann boundary condition along the membrane for the solute release from RBC towards the external fluid reads

Figure 14 :

 14 Figure 14: ATP fields under different viscosity contrast values. For λ = 1, the vesicles undergo tank-treading motion, while for λ = 8, the tumbling motion prevails. In the tumbling regime the release pattern along the membrane is more inhomogeneous, but the overall released concentration is lesser than in the tank-treading regime (λ = 1).

  right hand sides of Eqs. (B.3) are linear combination of the boundary condition (Eq.(B.1) ), and in addition they both have the same magnitude regarding the error term O( 2 ). By substituting α Eq. (B.1)) into Eq.(B.3), the boundary scheme (Eq. (13)) can be straightforwardly extracted. Note that in the special case when u ≡ 0 and α the boundary scheme will coincide with the classical bounce-back one. An important point is worth of mention. The choice of the zigzag boundary passing precisely at the middle of the lattice segments is not innocuous. Had we chosen another definition, then the precision would have degraded down to lower order, namely O( ) instead of enjoying O( 2 ) precision. Akbarzadeh A, Rezaei-Sadabady R, Davaran S, Joo SW, Zarghami N, Hanifehpour Y, Samiei M, Kouhi M, Nejati-Koshki K. Liposome: classification, preparation, and applications. Nanoscale Res Lett 2013;8(1):102. Allen TM, Cullis PR. Liposomal drug delivery systems: from concept to clinical applications. Adv Drug Deliv Rev 2013;65(1):36-48.

is suitable for the type of the boundary shown by the dashed-dotted lines in Figs.

(where margination prevails) we selected the initial liposome position to be in

i=0 ĉi = 0, and that

Treatment of the Moving Boundary

When dealing with an interface two questions arise: (i) how is the inter-lattice point with simple strategy, allowing for efficient GPU parallelization.

Since the advection-diffusion solver is designed for generic usage, its implementation for other problems, such as two-sided coupling cases is straightforward. Moreover, the scheme is based on a simple pixelization of curved boundary, with the help of parallel voxelizer, implying that the extension to 3D is also feasible.

We coupled this solver to a well-validated immersed boundary lattice-Boltzmann fluid-vesicle solver and implemented the resulting code for the study of liposome based hydrophilic drug delivery problem. We confirmed that with an assumption of instant absorption on vessel wall, when the liposome margination effect is strong, the presence of red blood cells facilitates the absorption. We have also demonstrated the potential of using this solver for shear stress induced cell-signaling process, and have exemplified it for the problem of ATP release by RBCs under flow. bf Acknowledgements C.M. and H.Z. acknowledge financial support from CNES (Centre National d'Etudes Spatiales), the French-German University Programme "Living Fluids" (Grant CFDA-Q1-14), and the China Scholarship Council (CSC).

Appendix A. Asymptotic Analysis of the Bulk Equations and Boundary Conditions

In this appendix, we show on one hand that the lattice-Boltzmann equations are equivalent in the asymptotic limit to the advection-diffusion equation and on the other hand they allow to derive the boundary conditions (Eq. ( 13)); see [START_REF] Huang | Boundary conditions of the lattice boltzmann method for convection-diffusion equations[END_REF]. Similar procedures were used in [START_REF] Junk | Asymptotic analysis of the lattice boltzmann equation[END_REF]; [START_REF] Yoshida | Multiple-relaxation-time lattice boltzmann model for the convection and anisotropic diffusion equation[END_REF] for other strategy of handling boundary conditions, different from that given by Eq. ( 13) (which we may view as a modified bounce-back scheme). By using our strategy, we will prove analytically here that this boundary condition enjoys a second order precision. We deal here with the zigzag boundary which is defined to pass through the middle of the mesh segment (dashed line in Fig. 1 (b)). It is this choice that allowed us to reach the second-order precision.

Let us first introduce the diffusive scaling which is based on the idea of finding a suitable pair of scaling size (which is also the numerical mesh size) (∆t, ∆x) that makes the diffusivity in numerical simulation D (as it was defined in Eq. ( 9)) close to O(1) in magnitude. This conditions reads

The following choice of scales satisfies our constraint

Here a = ∆t/∆x 2 = D /D is a constant value representing the ratio between numerical diffusivity D and physical diffusivity D. The factor is a small parameter and is often introduced in this way in the context of asymptotic analysis. Accordingly the scaling of velocity is given by

Our asymptotic analysis below obtained in the limit → 0 will show that (i) the Boltzmann equation recovers the advection-diffusion equation, and (ii)

the convergence to the advection-diffusion equation is of order O( 2 ) for both the bulk equations and the boundary conditions.

The Boltzmann equation (combining collision and streaming process (equation (5 -7)) reads, by omitting the reaction term, as

The unspecified argument in g i is (t, x) and is omitted here, ĉi = c i /(∆x/∆t) is the normalized unit vector of micro velocities defined in Eq. ( 4).

We will expand g i and c in powers of up to third order. It turns out that the second order provides the leading order contribution for the advectiondiffusion problem, whereas the third order is necessary for obtaining the desired precision. The expansion reads

Expanding Eq. (A.4) in power series of and ignoring terms of higher order w i enjoys a symmetry property (w 1 = w 2 and w 3 = w 4 ). Performing again the same operation with the third equation, and using the previous orders results leads (after simple algebraic manipulations) to an advection diffusion equation

Finally, performing the same operation with the last equation yields .11) This proves that the lattice Boltzmann scheme (A.4) converges to the advection diffusion equation with an error term O(∆x 2 ). The relation between diffusivity D and relaxation time τ is given by (2τ -1)/(6a) = D, which also

Appendix B. Boundary Condition: A Modified Half-Way Bounce-Back Scheme

In this section, we attempt to explain the boundary scheme for the static zigzag boundary, which is given by Eq. ( 13). This relation has been originally proposed in Ref. [START_REF] Huang | Boundary conditions of the lattice boltzmann method for convection-diffusion equations[END_REF], in which, a second-order convergence is observed in numerical experimentsHuang et al. (2016). We present here a simplified derivation for the particular zigzag boundary which intersects with mesh segments only at middle points (dashed-dotted line in Fig. 1 b). Our derivation will prove analytically the second order precision.

A representative point of the discretized boundary is designated as M in

For brevity, all unspecified arguments on the right hand side are understood to be (x, t).

The main question now is how to substitute the streaming step (6) at the boundary, in a such a way to respect the above boundary conditions. g i (t+∆t, x)

is the unknown incoming distribution function that needs to be determined. Inspired by the traditional bounce back condition (for Navier-stokes equation) we introduce the distribution g * ī (t, x), which is the known post-collision distribu-