## Abstract

We describe a new boundary-integral representation for biphasic mixture theory, which allows us to efficiently solve certain elastohydrodynamic–mobility problems using boundary element methods. We apply this formulation to model the motion of a rigid particle through a microtube which has non-uniform wall shape, is filled with a viscous Newtonian fluid, and is lined with a thin poroelastic layer. This is relevant to scenarios such as the transport of small rigid cells (such as neutrophils) through microvessels that are lined with an endothelial glycocalyx layer (EGL). In this context, we examine the impact of geometry upon some recently reported phenomena, including the creation of viscous eddies, fluid flux into the EGL, as well as the role of the EGL in transmitting mechanical signals to the underlying endothelial cells.

## 1. Introduction

Biphasic mixture theory is commonly used to model linearly poroelastic materials, including cartilage [1], and the brush-like layer that lines the microvasculature, namely the endothelial glycocalyx layer (EGL) [2]. It has long been known how to represent the fluid phase of these biphasic mixture theory models in boundary-integral form [3]. This is a convenient formulation as it describes the flow entirely in terms of quantities defined on flow surfaces, thereby effectively reducing the dimensionality of the problem. However, reformulating the solid phase into boundary-integral form is hampered by the appearance of volume integrals that arise from the coupling between the fluid and solid phases. Here, we demonstrate how these volume integrals can be recast into surface integrals, to yield a boundary-integral formulation for the full poroelastic dynamics. In doing so, we open up biphasic mixture theory models to a boundary element method (BEM) numerical treatment.

Although we present a general purpose boundary-integral formulation for linearly poroelastic materials, our immediate interest in this problem is motivated by the microvasculature, and our choice of physical parameters is informed by this physical regime. The walls of blood vessels are composed of endothelial cells which have non-uniform shape, and which can be coated with a layer consisting of a mixture of macro-molecules. In the literature, this brush-like layer is referred to as the EGL, or sometimes the endothelial surface layer. This layer has a gel-like structure and comprises proteoglycans, glycosaminoglycans, glycoproteins and absorbed plasma proteins [4]. The EGL is believed to serve a number of functions, including mechanotransmission of fluid shear stress (FSS) to the endothelial actin cortical cytoskeleton (degradation of the EGL is seen to be correlated with the endothelial cells becoming less likely to align with the flow, i.e. remodel [5]), a modulator of permeability in the transcapillary exchange of water, and as a regulator in the inflammatory response where it is believed to play a role in the leucocyte adhesion cascade [6]. In addition, the results obtained by Vink *et al.* [7] suggest that, along with other functions, the EGL plays a major role in providing vessels with an anti-adhesive inner lining. As such, the interplay between the EGL, the flow of blood plasma within the lumen and cells within the plasma is expected to be an important one for cardiovascular health.

The motion of cells through the microvasculature is a topic that has received much previous theoretical treatment. Some models have considered rigid cells passing through a straight vessel in the absence of an EGL [8,9]. When the cell is allowed to deform, as do red blood cells (RBCs), cell-depleted regions form adjacent to the vessel walls, and an accompanying drop in the apparent viscosity of the fluid is observed [10]. The role of a rigid, but porous, EGL in the migration of deformable cells from rigid surfaces has recently been considered, where it is suggested that the EGL may act to reduce the thickness of the depletion layer [11]. Boundary-integral formulations already exist that can describe these EGL-free or rigid-EGL scenarios [12,13].

However, the EGL is deformable, and some previous studies have accounted for its poroelastic behaviour using biphasic mixture theory. In a wavy-walled vessel, and in the absence of any cells, it has been demonstrated how the flow can separate and form a recirculation region which may, in turn, influence molecular transport and cellular response [14]. Other studies have adopted a similar approach to consider the passage through the lumen of rigid cells [2], EGL-coated cells [15] and cells which can undergo large deformations (e.g. RBCs) [16–19]. These studies have shown that the poroelastic layer can significantly affect the apparent viscosity of the flow, and that hydrodynamic forces can perhaps explain an observed exclusion of RBCs from the EGL under flow [18,20]. There is also the suggestion that the presence of the compliant EGL can reduce the FSSs exerted upon a RBC as it passes through non-uniformly shaped microvessels [19]. The EGL is also predicted to suppress FSSs on the endothelial wall, which seems at odds with its role as a transducer of mechanical stresses. It is now hypothesized that a significant portion of the stress is actually carried through the solid phase of the EGL.

These previous biphasic mixture theory models of EGL-lined microvessels have typically adopted a lubrication theory approximation, which places certain geometrical constraints upon the model (i.e. a long-wavelength analysis). By implementing our boundary-integral representation for biphasic mixture theory, we are now able to consider the aforementioned effects in a more general setting. Special attention is paid to the effect of geometry upon the system as a whole and the interaction between the flow, particle and poroelastic layer.

In §2, we show how the governing equations for biphasic mixture theory can be recast into boundary-integral form. The results we present in §3 for the two-dimensional case show how wall shape and the presence of a rigid cell affect the flow and solid displacements of the EGL. We draw some conclusions from these observations in §4.

## 2. Formulation

### (a) Geometry

Our geometry consists of a lumen of radius *H**, and volume *Ω*_{l}, through which an incompressible Newtonian fluid can flow unhindered, as well as a poroelastic layer, *Ω*_{m}, which is attached to the vessel walls. The lumen contains a rigid particle, *Ω*_{p}, with surface

The intrinsic masses and the volumes occupied by the fluid and solid phases in the poroelastic layer are *m*_{f},*V*_{f} and *m*_{s},*V*_{s}, respectively. The total volume of the poroelastic layer is therefore *V* =*V*_{f}+*V*_{s}. Introducing the volume fractions for each phase as *ϕ*_{f}=*V*_{f}/*V* and *ϕ*_{s}=*V*_{s}/*V* , we can define two partial densities for each phase: *ρ*_{f}=*ϕ*_{f}*m*_{f}/*V*_{f} and *ρ*_{s}=*ϕ*_{s}*m*_{s}/*V*_{s}. The requirement that there are no voids in the poroelastic layer necessitates that *ϕ*_{f}+*ϕ*_{s}=1.

There is a flow with characteristic speed, *V* *, through the lumen, which in turn drives flows and elastic deformations within the poroelastic layer, as well as transporting a rigid particle through the lumen.

### (b) Particle motility

We assume that the particle moves as an impermeable rigid object (as might be largely the case for white blood cells, for example), and hence the instantaneous particle velocity, ** V***

_{p}, can be written as

*** is a point on the surface of the particle,**

*x******

*x*_{c}is its centre and

**is the axis of rotation (asterisks denote dimensional quantities). The particle translational velocity,**

*k****, and angular velocity,**

*W**ω**

_{p}, are initially unknown. We non-dimensionalize using

*f*_{p}is the traction on the particle surface (see §2c).

### (c) Hydrodynamics in the lumen

As we consider steady flows through very small tubes, such as those that might typify post-capillary venules or capillaries, the governing flow equations in the lumen are the incompressible Stokes flow equations
*P** are the flow velocities and pressure, and *μ*_{f} is the dynamic viscosity of the fluid. The Cauchy stress tensor for the fluid has the form ** σ***=−

*P**

**+(∇**

*I****+(∇**

*v****)**

*v*^{T}) (superscript T denotes transpose). Non-dimensionalizing according to

*c*

_{f}=2

*π*(

*i*,

*j*,

*k*=1,2) for the two-dimensional case, and

*c*

_{f}=4

*π*(

*i*,

*j*,

*k*=1,2,3) for three-dimensional flow. The integrals involving the Stokeslet

*G*

_{ij}are referred to as the single-layer potentials, whereas the second, third and fourth integrals involving

*T*

_{ijk}, the associated stress tensor, are only defined in a Cauchy principal value sense, and referred to as the double-layer potentials (similar terminology applies for the boundary-integral representations which follow). Details of the exact form of these tensors are given in the electronic supplementary material, S1. Here,

**=**

*f***⋅**

*σ***is the surface traction on the Stokes flow surfaces and**

*n***is the inward normal vector.**

*n*### (d) Elastohydrodynamics in the poroelastic layer

The poroelastic layer consists of a solid phase which is linearly elastic, and a fluid phase that obeys the porous flow equations. In what follows, the subscript *s* refers to the solid phase, while subscript *f* refers to the fluid phase. Following earlier work [2,14], we model the poroelastic layer using biphasic mixture theory [21–23] and make the assumption that each phase is incompressible, has the same density, is homogeneous and has negligible inertia.

From incompressibility, we have the mass conservation equation
** w*** is the velocity of the fluid phase and

*****

*v*_{s}is the velocity vector of the solid phase.

The momentum equations then take the form [22]
*ϕ*_{f}** σ***

_{f}+

*ϕ*

_{s}

*****

*σ*_{s})=0) where

*****

*σ*_{f}and

*****

*σ*_{s}are the intrinsic Cauchy stress tensors for the fluid and solid phases, respectively, and

*** is the momentum transfer tensor which expresses the force coupling due to the interaction between the two phases, which we define below.**

*π*For small strains (linear elastic theory), the constitutive equations for each phase are [22,24]
*p** is the flow pressure, ** u*** is the displacement vector of the solid phase,

*μ*

_{f}is the dynamic viscosity of the fluid (which we assume to be the same as that in the lumen), and λ

_{s}and

*μ*

_{s}are the Lame constants of the solid phase. The total stress is then

Following [2,14], momentum transfer is given by
*K** is the hydraulic resistivity of the biphasic mixture. As the solid phase is attached to an immovable solid wall, at steady state, in the Cartesian coordinates fixed on the tube, the solid velocities are zero. Hence,

#### (i) Fluid phase

Conservation of mass (assuming ** v***

_{s}=

**) dictates that**

*0******

*v*_{s}=

**, we can solve (2.18) for the fluid phase before determining the elastic deformation of the solid phase. Non-dimensionalizing according to**

*0**χ*=

*K**

*H**

^{2}/(

*ϕ*

_{f}

*μ*

_{f}) is a measure of porous resistance. Brinkman flow also has a well-known boundary-integral representation (equivalent to that for oscillatory Stokes flow [3])

*M*

_{ij}and

*R*

_{ijk}are the free-space singularity solutions to Brinkman’s equation (see the electronic supplementary material, S1). The coefficient

*c*

_{f}here is as defined for the boundary-integral representation for Stokes flow within the lumen (2.8), and

#### (ii) Solid phase

Substituting (2.13) and (2.16) into (2.11) yields the momentum conservation equations for the solid phase

Non-dimensionalizing using
*ϕ*=*ϕ*_{s}/*ϕ*_{f} and *ν*=λ_{s}/2(λ_{s}+*μ*_{s}) is Poisson’s ratio. Following the usual procedure for writing the Navier equation in boundary-integral form, we can use the Maxwell–Betti reciprocal relation to obtain the following integral form for the behaviour of the solid phase [26]:
*c*_{s}=4*π*(1−*ν*) or *c*_{s}=8*π*(1−*ν*) in two and three dimensions, respectively. Here, *p*_{,j} stands for the partial derivative with respect to the *x*_{j} coordinate, and *S*_{ij} and *K*_{ijk} are the Green’s function and fundamental stress tensor for isotropic linear elasticity (i.e. Kelvin solutions; see the electronic supplementary material, S1). This is not yet in boundary-integral form, due to the volume integrals of the two forcing terms. We shall now show how these can also be converted into surface integrals. We consider each forcing term in turn.

*Pressure forcing*. To convert the volume integral involving pressure gradients into a surface integral, we use divergence theorem and Green’s identities along with the property of fundamental solutions (see details in appendix Aa),
*β*=−1/*r* in two and three dimensions, respectively. In order to use identity (2.26), we need to know both the pressure and its normal derivative, ∂*p*/∂*n*, on the boundary. The boundary-integral flow representation (2.21) yields only flow velocities and tractions, from which surface pressures can be determined from the following boundary-integral relation [27]:
*Q*_{i} and *L*_{ik} are given in the electronic supplementary material, S1). The first integral is a Cauchy principal value integral, and the second a hyper-singular integral that must be regularized for numerical treatment (see appendix Aa for details). Consequently,

*Momentum transfer*. We now turn to the volume integral in (2.25) that stems from the momentum transfer between phases. Our approach involves consideration of the flow within the fluid phase (2.21)
*F*_{j}=*χw*_{j}, alongside a second complementary flow defined by
*S*_{ij} is Green’s function for linear elasticity, and *σ*^{Bi}, *v*^{Bi} are the Cauchy stress tensor and flow velocity for the complementary flow.

We then apply the Lorentz reciprocal relation to flows (*σ*_{f},** w**) and (

*σ*^{Bi},

*v*^{Bi}), integrate over the flow domain and apply the divergence theorem to obtain [3]

*g*^{Bi}=

*σ*^{Bi}⋅

**. It can be seen that, after substitution of the expressions for**

*n***and**

*F*

*F*^{Bi}into (2.32), we obtain the volume integral in terms of boundary integrals

Hence, the problem now reduces to finding a flow which satisfies (2.31), and which is given by

*Boundary-integral representation*. We arrive at the final boundary-integral representation for the solid phase
*q*(** x**)=∂

*p*/∂

*n*.

### (e) Boundary conditions

As we assume that the particle is impermeable, with a no-slip surface, *U*_{0},*W*_{0},*V*_{0} (see the electronic supplementary material, S2), i.e.

As we assume small-strain elasticity, in keeping with earlier models, all boundary conditions on the interfaces can be applied at their undeformed locations, *v*_{s}=** 0**, i.e. no elastic velocities)

The second interface boundary condition is continuity of traction [28]
** n** is a unit normal to the boundary and

**=**

*Γ**ϕ*

_{f}

*σ*_{f}+

*ϕ*

_{s}(

*σ*_{s}/

*ϕ*) is the total stress in the poroelastic material. The proportion of the total stress in the porous medium borne by each phase is proportional to its volume fraction, hence

*ϕ*

_{s}+

*ϕ*

_{f}=1).

## 3. Results

The boundary-integral formulation derived in §2 is very general, applicable in both two and three dimensions. Due to the computational expense, however, we follow [14] and consider a two-dimensional regime, where the vessel is modelled as a channel. Notation specific to this geometry is given in figure 2. We solve the governing integral equations using a BEM scheme, the particulars of which can be found in the electronic supplementary material, S4, alongside validation details.

### (a) Parameter values

We model the elastohydrodynamics of an EGL-coated microvessel containing a rigid cell (such as a white blood cell) with the new boundary-integral formulation. We consider a wavy-walled channel, with top and the bottom channel walls prescribed (in dimensional form) by
*η**_{±}(*x*_{1}), *ζ**_{±}(*x*_{1}) are spline interpolations that guarantee a smooth transition from a straight inlet/outlet to the non-uniform wavy topology (see the electronic supplementary material, S3). Hence, the channel has a mean width of 2*H** with wall undulations of amplitude *a** and wavelength of *Λ**_{e}, and we allow for a phase difference between the top and bottom walls of *Φ*_{+}=0, *Φ*_{−}=*Φ*, similar to [14]. The interfaces between the lumen and the poroelastic layers are located at

Parameter values were chosen to be broadly representative of the movement of a cell in a capillary. Hence, the vessel radius was chosen to be *H**=5 μm, representative of a capillary. We consider a spherical particle having radius *R**=2.5 μm, which is characteristic of a small lymphocyte. The fluid viscosity is assumed to be that of water *μ*_{f}=10^{−3} Pa s. The EGL thickness *ε** varies from 0.2 to 0.4 μm up to 1 μm [5]. Although there are currently no direct measurements of hydraulic resistivity within the EGL, estimates are in the range *K**=10^{10}−10^{11} N s m^{−4} [17]. The shear modulus of the EGL is calculated in [29] to be *ϕ*_{s}*μ*_{s}=3.5−10 Pa, and it is generally assumed [2] that the EGL has a small solid fraction *ϕ*_{s}=0.01. Its Poisson ratio is assumed to be *ν*=0.3. The mean blood velocity in capillaries is *V* *=0.8−1 mm s^{−1} [30]. An endothelial cell has length of approximately *Λ**_{e}=20−50 μm and height of *a**=1−2 μm [14]. Table 1 summarizes these parameters, including values for the non-dimensional quantities which define the dynamics, namely *χ*, *ε*=*ε**/*H**, *δ*=*H**/*Λ**_{e}, *a*=*a**/*H** and *R*=*R**/*H**.

Accordingly, we consider 10 cases corresponding to different combinations of the position of the particle’s centre, *x*_{c}, and shape of the walls. For each cell position, we also consider two distinct channel shapes: varicose (*Φ*=0) and sinuous (*Φ*=*π*/2) (see table 2 for a summary). The variation of channel height against *x*_{1} for both geometries is shown in the electronic supplementary material, S3.

In what follows, we present flow fields and associated flow shear stresses, *Γ*_{f}=*ϕ*_{f}** g**(

**)⋅**

*x***on the channel walls and**

*τ**Γ*

_{c}=

**(**

*f***)⋅**

*x***on the cell, corresponding to a varicose geometry,**

*τ**Φ*=0, and sinuous geometry,

*Φ*=

*π*/2. In addition, we present predictions for the elastic displacements in the EGL and associated elastic shear stresses on the wall,

*Γ*

_{s}=(

*ϕ*

_{s}(

**(**

*h***)−**

*x**ϕp*

**⋅**

*I***)/**

*n**ϕ*)⋅

**(here the direction of the tangential vector**

*τ***always coincides with the positive**

*τ**x*

_{1}direction). The total shear stress is

*Γ*=

*Γ*

_{f}+

*Γ*

_{s}.

### (b) No cell

Firstly, in figure 3, we compare the stresses on the interface with those on the wall and find that the EGL acts to reduce the FSS (as previously reported), but increases the stress in the solid phase. Combining these two contributions, we observe that in general the total shear stress on the wall is in fact greater than that on the interface. The exception is at the widest section of the vessel, although this is only non-negligible for a varicose vessel (figure 4).

Below, we show the flow fields, stresses (both fluid and elastic) and elastic displacements for both sinuous and varicose microvessels in the absence of the cell. These act as base cases, against which we can compare the cases where a cell is present. We observe that elevated stresses and elastic displacements occur at the geometric constrictions, as expected. Moreover, we note that the magnitude of the shear stress on the wall due to the solid phase dominates over that due to the fluid phase for all cases considered. This appears to support current ideas around the solid phase being the main transducer of mechanical stresses to the underlying endothelial cells [31].

### (c) Varicose geometry with the cell

In figure 5, we examine the flow fields and FSSs for the varicose geometry in the presence of a cell. When the cell is located in the geometric constriction (case I), in figure 5*a*, we observe a local amplification of stresses and flow velocities, over that seen when the cell is absent (figure 5*e*). However, immediately above the cell we observe a reduction in the shear stress. This is highlighted further in figure 6, which shows that the presence of the cell leads to increased wall stress (as compared with the cell-free vessel) immediately upstream and downstream of the cell, but decreased stress directly above the cell (i.e. *x*=−2). When the particle is located on the centreline of the vessel, and in the widest part of the vessel (figure 5*c*), we observe that the influence of the cell upon the FSSs on the wall is fairly minimal.

However, when the cell is positioned in the expanding section of the vessel (case II, as shown in figure 5*b*), we again see elevated levels of wall shear in the vicinity of the particle, and a local reduction in shear immediately above and below the cell. In fact, in this scenario, the cell-induced local stress reduction occurs to such an extent that we see a region of negative FSS on the wall. A similar situation occurs when the cell is placed in the widest section of the channel, but close to the upper wall (case IV; figure 6, dashed line). Upon closer examination of the associated flow fields for these two cases (figure 7), we see that this is associated with the presence of a vortex. It was shown in [14] that vortices can appear in a varicose vessel in the absence of a cell, although these were noted to appear in the widest part of the vessel, and for values of *χ*≥1600 which is greater than those considered here. From a physiological standpoint, these flow features are important as they have the capacity to increase the residence time of circulating substances within the EGL. Moreover, the accompanying variations of the shear stress profile could have important implications for mechanotransduction in the microvasculature, as mechanoreceptors located on the surface of the endothelial cells are liable to experience shear stresses exerted by flow within the EGL (fluid component of the total stress).

We also analyse the influence of the cell on the normal velocity through the interface, i.e. *et al.* [14] also previously reported a net flux of fluid into the EGL in the absence of a cell. We see in figure 8 that the presence of the particle leads to greater amplitude EGL fluxes (both positive and negative, as compared with the cell-free vessel) in the vicinity of the particle. The highest value of the normal flux occurs for case IV close to the particle (

When we examine the FSSs exerted upon the cell itself in this varicose geometry, we observe that it experiences a much larger shear stress than the solid walls. We also note the degree to which the distribution of stress over the surface changes with cell location in the channel. The cell takes its greatest translational velocity in the constriction, and lower velocities in the widest section of the vessel. As would be expected, for the case of a varicose channel, the angular velocity of the cell is negligible unless moved substantially from the vessel centreline (case IV), where its (non-dimensional) angular velocity is 0.1. In this case, the cell also experiences an up–down asymmetry in shear stress.

In figure 9, we examine the corresponding displacements and elastic shear stresses. Firstly, as in the no-cell scenario, we see that the elastic stresses are larger in magnitude than the fluid stresses, again offering some support to the notion that a significant proportion of the mechanotransduction in the EGL is performed by the solid phase. We also note that, unlike the FSSs, the elastic shear stresses always remain positive. However, for cases II and IV, where vortex structures were observed in the flow field, we note a corresponding distortion in the elastic displacement field. This is highlighted in figure 10, which compares case IV with the vortex-free case IX. Note that the maximum magnitude of the displacement observed is about 7% of the EGL thickness, and so within the 10% usually accepted for the small strain approximation to be held valid for biological tissue (although a more comprehensive examination of the linear elasticity assumption requires either experimental data or full nonlinear computations; certainly for smaller cells, we expect the linear elasticity assumption to become a better approximation, with the converse being true for larger cells—see the Discussion section for further comments). We can also examine the (non-dimensional) elastic deformations at the interface, and these are shown in figure 11 for the case where no cell is present (case IX), as well as when the cell is close to the upper wall (case IV). We see that the presence of the cell leads to a 65% increase in the displacements induced at the interface.

We also briefly investigate the influence of EGL thickness on the stresses exerted on the solid wall. Figures 12 and 13 compare the elastic and fluid stress distributions in case II with *ε*=0.1, *ε*=0.125, *ε*=0.15, *ε*=0.175 and *ε*=0.2. We observe that a reduction in EGL thickness leads to approximately proportional reductions in the elastic stresses exerted upon the wall (e.g. a 50% reduction in EGL thickness leads to an approximately 50% drop in elastic shear). FSSs, however, increase nonlinearly with decreasing layer thickness. Decreasing EGL thickness by 25% leads to an approximately 20% increase in stress while changing it by 50% increases the stress by 60%. In addition, FSS is always positive for *ε*≤0.15, which indicates that for these EGL thicknesses a vortex is not generated.

### (d) Sinuous geometry with the cell

Let us now briefly examine the changes which occur when the vessel is sinuous (*Φ*=*π*/2). We first note from figures 14 and 15 that the magnitudes of the stresses on the walls and the cell, from both the fluid and solid phases, are comparable with those observed in the varicose case. As the vessel does not expand and contract with downstream distance in the manner of the varicose vessel, the FSSs on the wall have a greater tendency to stay positive. For example, see figure 16, which compares the flow shear stress on the upper wall when *x*_{c}=(−1,0), for both *Φ*=0 (case II) and *Φ*=*π*/2 (case VI). We see that the region of negative FSS disappears.

However, if we move the cell sufficiently close to the upper layer (case VIII), figure 14*d* demonstrates that negative FSS on the wall can still occur in the vicinity of the cell and is comparable in magnitude to that observed for the varicose vessel (figure 17). If we examine this region more closely, we see from figure 18 that there is an associated recirculation in the flow field, and corresponding modification of the elastic displacement field. However, this is slightly offset from the cell centre, as compared with that observed in a varicose vessel (cf. figure 7).

In terms of cell mobility, in the sinuous geometry, the cell tends to rotate more and the translational velocity varies within a smaller range (0.68≤** W**≤0.92) than in the varicose geometry (0.6≤

**≤0.98).**

*W*## 4. Discussion

Biphasic mixture theory is widely used to model linearly poroelastic materials. The linearity of the governing equations suggests that a boundary-integral representation is possible. However, the traditional derivations lead to volume integrals of the terms which couple together the solid and fluid phases. By demonstrating how these can be converted into surface integrals, we have been able to derive a true boundary-integral representation for biphasic mixture theory, which we believe will have wide use and applicability.

This has allowed us to examine the elastohydrodynamics of a rigid particle negotiating a poroelastic-lined channel of general shape, as a model for a small lymphocyte negotiating a microvessel coated with an EGL. The EGL is believed to play an important role in transmitting mechanical stresses to the underlying endothelial cells. Indeed, when the EGL is compromised, it has previously been reported that the endothelial cells are less able to align with the flow direction [5]. Numerically solving the boundary-integral representation of the biphasic mixture theory equations, we have been able to examine the effect of vessel shape, and particle position, upon the poroelastic dynamics of the EGL, and motility of the cell. Table 3 summarizes the main findings. Our computations suggest that the shear stresses on the vessel wall consist of significant contributions from the solid phase, and that this contribution increases with increased EGL thickness. This seems to support current theories that suggest that a significant proportion of the mechanical stress applied at the EGL interface is transferred through the solid phase (which is consistent with another one of the EGL’s hypothesized roles in protecting the endothelial cells from excessive FSSs) [5]. We noted that the magnitude of the wall shear stress was a function of EGL thickness. The contribution from the solid phase was seen to decrease (almost linearly) with decreased EGL thickness, in contrast to the contribution from the fluid phase, which increased (nonlinearly) as the EGL became thinner. These findings have potential implications for mechanotransduction through compromised EGLs, and the damage to the underlying endothelial cells which might ensue as a result.

In the presence of the rigid cell, we see heightened levels of wall shear from both the fluid and elastic phases. We also find that regions of negative shear stress on the vessel wall can develop in the immediate vicinity of the cell, when it is located relatively close to the EGL. These, in turn, are associated with the presence of recirculating flow regions in the EGL. These recirculating flow regions are seen to be produced more readily in the varicose vessel, where the vessel diameter varies more rapidly. These eddies have previously been reported in the absence of a cell, however only for much higher EGL hydraulic resistivities, and at different locations from those reported here (specifically at the widest section of the microvessel [14]). The regions of recirculating flow could be physiologically important, as they have the potential to increase the residence time of blood components.

Following Wei *et al.* [14], we also examined the flux of fluid into the EGL from the lumen, which potentially has implications for the hypothesized filtering function of the EGL [32]. As in their work, we note the sensitivity of this flux to wall shape, which we observe to be greater in magnitude for a sinuous vessel than for a varicose one. However, we also observe the impact of the cell, the presence of which can lead to an almost fivefold increase in the local flux into the EGL.

In terms of transport of the cell itself, we observe that it generally travels more rapidly through the sinuous vessel than through the varicose one. However, the maximum speed we observe occurs within the constrictions of the varicose vessel. It also exhibits non-negligible angular velocity when close to the EGL, a scenario which occurs more often in the sinuous vessel, where it consequently rotates at more locations along the length of the vessel. These rolling motions are potentially material to lymphocyte recruitment, which forms part of the body’s immune response. We also note that the cell experiences greater shear stresses than the vessel walls, perhaps unsurprisingly given the expected role of the EGL in protecting the endothelial cells from potentially harmful levels of FSS. The cell shear stresses, however, appear to be comparable in magnitude for both shapes considered here, i.e. varicose and sinuous. For computational convenience, we have examined the instantaneous motions of the cells at manually specified positions in the geometry. It would, however, be interesting to track the cell’s trajectory from some starting position. Such simulations would necessarily be very computationally expensive, although they would enable us to determine whether certain configurations are stable or unstable, i.e. determine the extent to which a cell close to the vessel wall remains at this distance.

In keeping with earlier studies, we have used a relatively simple poroelastic, continuum model (i.e. biphasic mixture theory) to encapsulate the elastohydrodynamics of the EGL, and all of its structural complexity. It would be interesting to see whether this particular division of stress from the fluid and solid phases holds under more complex representations of the EGL. It would also be useful to apply our new boundary-integral representation to see how the aforementioned elastohydrodynamics are modified in a three-dimensional setting, using biologically informed vessel geometries, and this is work that we are currently undertaking [33]. We have also assumed small-strain deformations of the EGL, which is legitimate provided that the cells are sufficiently small. However, in order to gauge the range of validity of this linear elasticity assumption, and to model situations where cells are large enough to generate displacements comparable to the EGL thickness, it would be useful to extend the model to incorporate nonlinear elastic effects. On a related theme, another useful extension of the model would be to consider the motion of deformable cells, such as RBCs, through the lumen by adopting a capsule model to capture the finite-strain mechanics of the cell [34].

Finally, there are also additional mechanisms besides elastic forces that are believed to restore the EGL to an equilibrium configuration, following its deformation due to the passage of a cell through a microvessel. These include oncotic processes [17], whereby a difference in the concentration of the plasma proteins in the EGL and lumen plasma generates an oncotic pressure. This pressure difference leads to relaxation of the EGL back to its equilibrium profile. It is also hypothesized that mechano-electrochemical effects can play a similar role [35]. This comes about because the EGL is believed to be hydrated by an electrolytic solution, which contains electrostatically charged macromolecules. The changes in the charge density, which occur when the EGL is compressed, thereby provide another restoring force. It would also be valuable to include both of these effects into our EGL model.

## Data accessibility

The data supporting this article have been uploaded as part of the electronic supplementary material.

## Authors’ contributions

P.P.S. developed the numerical code and prepared all the figures. All authors contributed equally to all other aspects of the work.

## Competing interests

We declare we have no competing interests.

## Funding

P.P.S. is funded by a University of Auckland Doctoral Scholarship.

## Acknowledgements

We are grateful to Tet Chuan Lee for useful suggestions regarding the boundary-integral representation of the poroelastic layer. Numerical computations were run on the NZ eScience Infrastructure (NeSI) facility.

## Appendix A

**(a) Pressure forcing terms**

Let us consider the integral corresponding to the *i*th component
^{2}*p*=0),

**(i) Pressure regularization**

In order to regularize the pressure equation (2.27), we first rewrite it as
** V**=

**(**

*w*

*x*_{0}). The corresponding pressure is given by (2.20) as

*p*=−

*χ*

**⋅**

*V***and**

*x***=−**

*g**p*

**. From the boundary-integral representation for pressure (2.27), we obtain**

*n*

*e*_{i}is the unit vector in the

*x*

_{i}direction. Substituting (A 14) into (A 13) gives us

**(b) Complementary flow problem**

**(i) Two-dimensional case**

Let us consider the complementary problem
*σ*^{Bi} is the Cauchy stress tensors for the fluid phase. Let us consider the case *i*=1, as the solution when *i*=2 is derived similarly. Rewriting (A 17) in terms of components, we have
*A*_{1},*A*_{2} are constants and ^{2}*p*=0. There is freedom in the choice of values *A*_{1},*A*_{2} and *ξ*, as, for our purposes, we can choose any boundary conditions to (A 18) and (A 19). Choosing *I*_{α}(*η*) and *K*_{α}(*η*) are the modified Bessel functions of the first and second kind, respectively. The solution consists of two parts, namely a general solution of the homogeneous equation and a particular solution. The homogeneous part contains constants that need to be chosen to remove any algebraic singularities at *r*=0. Analysing the small-scale behaviour of *K*_{α}(*η*) gives [36]
*γ*=0.57721… is the Euler–Maschelroni constant. As a consequence, choosing the constant

The final solutions for both terms, *F*^{B1} and *F*^{B2}, can be written as
*B*_{i} relates to the solution corresponding to the *F*^{Bi} term. This has an associated stress tensor (*i*,*j*=1…2, *i*≠*j*)

For small *r*, we have
*r*=0. Hence, flow velocities and stresses are non-singular and the integrals (2.33) are proper.

**(ii) Three-dimensional case**

For the complementary problem (A 17) corresponding to the three-dimensional case, we have *F*^{B1}. The resulting equations are as follows (*i*=1,2,3):
*i*=1,2,3)
*i*=1. To solve this equation, we consider the solution in the form
*x*_{0}, and seeking a solution for *r* alone we come to the equations
*f*(*r*)
*f*′=d*f*(*r*)/d*r*. Solving equation (A 43), we find an expression for the function *f*(*r*)
*r* or *χ* and it does not contain a strong singularity as

Similarly, we can find the solution to problems (A 36) when *i*=1,2,

- Received December 14, 2014.
- Accepted May 11, 2015.

© 2015 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.