A boundary-integral representation for biphasic mixture theory, with application to the post-capillary glycocalyx.

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.

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.

Formulation

Geometry

Our geometry consists of a lumen of radius H*, and volume Ω, through which an incompressible Newtonian fluid can flow unhindered, as well as a poroelastic layer, Ω, which is attached to the vessel walls. The lumen contains a rigid particle, Ω, with surface . The lumen and the poroleastic layer are bounded by the surfaces, and , respectively. The tube has a rigid wall denoted by . Surface forms the interface between the lumen and the poroelastic layer. Denoting inlet and outlet surfaces as , we have and . See figure 1 for a diagram of the geometry.

Figure 1.

A diagram of the geometry of the microtube containing a rigid cell, and lined with a poroelastic layer.

The intrinsic masses and the volumes occupied by the fluid and solid phases in the poroelastic layer are mf,Vf and ms,Vs, respectively. The total volume of the poroelastic layer is therefore V =Vf+Vs. Introducing the volume fractions for each phase as ϕf=Vf/V and ϕs=Vs/V , we can define two partial densities for each phase: ρf=ϕfmf/Vf and ρs=ϕsms/Vs. 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.

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, *, can be written as where * is a point on the surface of the particle, *c is its centre and is the axis of rotation (asterisks denote dimensional quantities). The particle translational velocity, *, and angular velocity, ω*, are initially unknown. We non-dimensionalize using which gives In the absence of inertia (due to the small scale of the flow) and external forces or torques, the particle must be both force and torque free, i.e. where p is the traction on the particle surface (see §2c).

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 where and 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*+(∇*+(∇*)T) (superscript T denotes transpose). Non-dimensionalizing according to yields This can be recast into the following boundary-integral form [3]: c=2π (i,j,k=1,2) for the two-dimensional case, and c=4π (i,j,k=1,2,3) for three-dimensional flow. The integrals involving the Stokeslet G are referred to as the single-layer potentials, whereas the second, third and fourth integrals involving T, 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, =⋅ is the surface traction on the Stokes flow surfaces and is the inward normal vector.

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 where * is the velocity of the fluid phase and *s is the velocity vector of the solid phase.

The momentum equations then take the form [22] and (note that ∇⋅(ϕ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] and where p* is the flow pressure, * 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 It can be seen that the total stress tensor represents the sum of fluid and elastic stress tensors with rescaled material constants, corresponding to the partial stress tensors.

Following [2,14], momentum transfer is given by where 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,

Fluid phase

Conservation of mass (assuming *s=) dictates that With regards to momentum conservation, upon substituting (2.12) and (2.16) into (2.10), we obtain the Brinkman-type equation for the fluid phase [25], We note that, owing to *s=, we can solve (2.18) for the fluid phase before determining the elastic deformation of the solid phase. Non-dimensionalizing according to gives where χ=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]) where M and R are the free-space singularity solutions to Brinkman’s equation (see the electronic supplementary material, S1). The coefficient c here is as defined for the boundary-integral representation for Stokes flow within the lumen (2.8), and is the surface traction due to flow in the fluid phase.

Solid phase

Substituting (2.13) and (2.16) into (2.11) yields the momentum conservation equations for the solid phase Hence, the solid phase is governed by the steady Navier equation with two forcing terms.

Non-dimensionalizing using gives where ϕ=ϕ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]: with cs=4π(1−ν) or cs=8π(1−ν) in two and three dimensions, respectively. Here, p, stands for the partial derivative with respect to the x coordinate, and is the traction vector for the elastic part of the solid phase. Also S and K 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), where , and or β=−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]: (singularity solutions Qi and L 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, Once surface pressure is known, its normal gradient can be determined through Green’s third identity,

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) where F=χw, alongside a second complementary flow defined by where , S is Green’s function for linear elasticity, and , are the Cauchy stress tensor and flow velocity for the complementary flow.

We then apply the Lorentz reciprocal relation to flows (f,) and (,), integrate over the flow domain and apply the divergence theorem to obtain [3] where =⋅. It can be seen that, after substitution of the expressions for and 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 where and (see appendix Ab for details).

Boundary-integral representation. We arrive at the final boundary-integral representation for the solid phase where q()=∂p/∂n.

Boundary conditions

As we assume that the particle is impermeable, with a no-slip surface, , we have The vessel wall, , is assumed to be rigid and impermeable, hence we have On the inlet and outlet surfaces, we prescribe the analytical solutions for displacement and velocity corresponding to flow through a straight-walled, poroelastic-lined tube, 0,0,0 (see the electronic supplementary material, S2), i.e. As there is no inertia in the flow, we expect entry and exit development effects to be localized, typically scaling with tube radius.

As we assume small-strain elasticity, in keeping with earlier models, all boundary conditions on the interfaces can be applied at their undeformed locations, . The first condition on the interfaces is continuity of flow velocity (accounting for the fact that we assume s=, i.e. no elastic velocities)

The second interface boundary condition is continuity of traction [28] where is a unit normal to the boundary and =ϕff+ϕ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 and (recalling that ϕs+ϕf=1).

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.

Figure 2.

A diagram illustrating two-dimensional geometry of the model, detailing surface labels.

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 where the positive and negative signs correspond to and , respectively. The functions η*±(x1), ζ*±(x1) 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 2H* 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*=1010−1011 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*.

Table 1.

Typical non-dimensional parameter values for a small lymphocyte negotiating an EGL-lined capillary.

R	χ	a	ϕs	ε	δ	ν
0.5	250	0.2	0.01	0.2	0.25	0.3

Accordingly, we consider 10 cases corresponding to different combinations of the position of the particle’s centre, 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 x1 for both geometries is shown in the electronic supplementary material, S3.

Table 2.

Wall shapes and cell positions for the various cases considered.

case	I	II	III	IV	V	VI	VII	VIII	IX	X
xc	(−2,0)	(−1,0)	(0, 0)	(0, 0.4)	(−2,0)	(−1,0)	(0, 0)	(0, 0.4)	—	—
Φ	0	0	0	0	π/2	π/2	π/2	π/2	0	π/2

In what follows, we present flow fields and associated flow shear stresses, Γf=ϕf()⋅ on the channel walls and Γc=()⋅ 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(()−ϕp⋅)/ϕ)⋅ (here the direction of the tangential vector always coincides with the positive x1 direction). The total shear stress is Γ=Γf+Γs.

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).

Figure 3.

Shear stress distribution on the interface and on the solid wall in the absence of a cell (case IX): (black solid line) ; (black dashed line) (elastic stresses); (grey solid line) ; (grey dashed line) (fluid stresses).

Figure 4.

Ratio of the total shear stresses Γ exerted on the interface , compared with that on the wall , in the absence of a cell. Here, ; (solid line) case IX; (dashed line) case X.

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].

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 5a, we observe a local amplification of stresses and flow velocities, over that seen when the cell is absent (figure 5e). 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 5c), we observe that the influence of the cell upon the FSSs on the wall is fairly minimal.

Figure 5.

Varicose vessels ((a–d) cases I–IV and (e) IX, respectively) showing flow fields and shear stresses exerted by the flow. The first colour bar indicates the magnitude of the shear stresses on the vessel walls, whereas the scale beneath corresponds to the stresses on the cell. Corresponding translational and angular velocities are: (a) =(0.98,0), ωp=0; (b) =(0.88,0), ωp=0; (c) = (0.7,0), ωp=0; (d) =(0.6,0), ωp=−0.1. Regions inside dashed boxes in (b,d) are shown magnified in figure 7.

Figure 6.

FSS distributions, showing regions of rapid change: (thin black solid line) case IX; (thick grey solid line) case I; (dashed line) case IV.

Figure 7.

Flow fields in regions of the EGL close to the cell, for (a) case II and (b) case IV. In both cases, the presence of a vortex is evident. These regions of magnification correspond to the dashed boxes shown in figure 5.

However, when the cell is positioned in the expanding section of the vessel (case II, as shown in figure 5b), 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. , since Wei 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 () and the lowest for case X, where there is no cell present (). Although not shown in figure 8, we also note that, for each cell position, phase shifting of the wall does not significantly affect the normal flux distribution.

Figure 8.

Normal fluid fluxes, wn, at the interface between core flow and EGL: (thin black solid line) case IX; (thick grey solid line) case I; (dashed line) case IV.

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.

Figure 9.

Varicose vessels ((a–d) cases I–IV and (e) IX, respectively) showing elastic displacement vectors and shear stresses exerted by the solid phase, the magnitudes of which are indicated by the colour bars beneath. Regions inside dashed boxes in (d,e) are shown magnified in figure 10.

Figure 10.

Elastic displacements in the EGL interior: (a) case IV and (b) case IX. These regions of magnification correspond to the dashed boxes shown in figure 9.

Figure 11.

(Non-dimensional) horizontal u1 (a) and vertical u2 (b) elastic displacements on the EGL interface: (solid lines) case IV ; (dashed lines) case IX.

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.

Figure 12.

Solid shear stress distribution on the rigid wall corresponding to case II: (black solid line) ε=0.2; (dashed line) ε=0.175; (dotted line) ε=0.15; (dashed-dotted line) ε=0.125; (grey solid line) ε=0.1.

Figure 13.

FSS distribution on the rigid wall corresponding to case II: (black solid line) ε=0.2; (dashed line) ε=0.175; (dotted line) ε=0.15; (dashed-dotted line) ε=0.125; (grey solid line) ε=0.1.

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 c=(−1,0), for both Φ=0 (case II) and Φ=π/2 (case VI). We see that the region of negative FSS disappears.

Figure 14.

Sinuous vessels ((a–d) cases V–VIII and (e) X, respectively): flow fields and shear stresses exerted by the flow. The top colour bar indicates the magnitude of shear stresses on the wall, whereas the bottom colour bar indicates the magnitude of the cell shear stresses. Corresponding translational and angular velocities are: (a) =(0.92,0.1), ωp=0; (b) =(0.92,0.1), ωp=0.1; (c) =(0.82,−0.1), ωp=0.1; (d) =(0.68,0), ωp=−0.1. The region inside a dashed box in (d) is shown magnified in figure 18.

Figure 15.

Sinuous vessels ((a–d) cases V–VIII and (e) X, respectively): elastic displacements and shear stresses exerted by the solid phase, the magnitudes of which are indicated by the colour bar. The region inside the dashed box in (d) is shown magnified in figure 18.

Figure 16.

FSS distribution at particle location c=(−1,0), demonstrating the disappearance of negative values when the vessel moves from a varicose shape (case II, solid line) to a sinuous shape (case VI, dashed line).

Figure 18.

Case VIII: (a) flow field, showing the presence of a vortex, and (b) associated elastic displacement.

However, if we move the cell sufficiently close to the upper layer (case VIII), figure 14d 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).

Figure 17.

FSS distribution at particle location c=(0,0.4), comparing the differences between the varicose (solid line) and sinuous (dashed line) geometries, i.e. cases IV and VIII.

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≤≤0.92) than in the varicose geometry (0.6≤≤0.98).

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.

Table 3.

Summary of findings for the different scenarios considered, including maximum and minimum values of FSSs on the wall (), interface () and cell surface (). Similarly, elastic stresses on the wall and interface ( and , respectively). We also summarize the maximum flux of fluid into the EGL (wn).

	I	II	III	IV	V	VI	VII	VIII	IX	X
max{Γ_f(S_10)}	1.1	1.4	0.8	1.1	1.6	1.3	0.8	1.4	0.9	0.7
max{Γ_f(S_4)}	5.2	4.4	3.6	3.7	5.3	4.2	3.2	4.4	3.6	3.2
min{Γ_f(S_10)}	0.1	−0.1	0.02	−1.1	0.1	0.1	−0.02	−1.3	0.1	0.1
min{Γ_f(S_4)}	1.2	1.2	1.1	1.2	1.4	1.4	1.0	1.4	1.2	1.2
max{Γ_s(S_10)}	5.9	6.0	4.8	4.8	6.6	5.5	4.2	4.7	4.8	4.2
max{Γ_s(S_4)}	0.05	0.04	0.03	0.04	0.05	0.04	0.03	0.04	0.03	0.03
min{Γ_s(S_10)}	1.1	0.9	0.9	0.8	1.2	1.2	0.9	0.9	1.1	1.2
min{Γ_s(S_4)}	0.01	0.01	0.01	0.01	0.01	0.01	0.01	0.01	0.01	0.01
max{Γ_c}	7.2	4.2	2.3	7.9	7.3	7.5	3.9	8.2	—	—
max{w_n}	0.026	0.05	0.013	0.07	0.05	0.04	0.013	0.069	0.013	0.012

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.