## Abstract

The delivery of endothelial ligands and macromolecules, such as lipoproteins, from the circulation to the arterial wall is of central importance in modulating endothelial cell function and physiology and, consequently, in the onset of vascular disease. Given the strong spatial correlation between areas of disturbed blood flow and occurrence of atherosclerotic plaque, a detailed understanding of the effects of different fluid flow characteristics on the delivery of factors in the bloodstream to the endothelium is an essential step towards understanding the observed localization of vascular disease to certain focal sites within the vasculature. In this paper, a model of biochemical mass transport in a two-dimensional flow chamber with spatially varying wall shear stress is presented. The advantage of the relatively simple arterial geometry is that all the essential features of blood flow *in vivo* are captured, but the underlying effects on mass transport, and, hence, on endothelial cell function, are not masked by complex three-dimensional flow. Indeed, it is demonstrated that a previously derived similarity solution, in terms of the shear stress on the endothelial wall, is an asymptotically close approximation to the exact solution to the advection–diffusion equation. The mathematical analysis is then used to identify fundamental links between the wall shear stress and the distribution of chemical species along the endothelium. The physiological implications of these links are discussed, in particular the location of maximum chemical concentration on the endothelium, which is of great significance to the regulation of intracellular signalling and, consequently, endothelial cell function.

## 1. Introduction

The fluid mechanical characteristics of arterial blood flow are known to have a significant impact on the delivery of blood-soluble biochemical species to the endothelial wall of the vessel (Giddens *et al*. 1993). For instance, arterial bifurcations and bends frequently exhibit flow detachment and the formation of recirculation zones (Malek *et al*. 1999), resulting in a spatially varying chemical concentration along the vessel wall. Physiologically important biochemical factors known to interact with the artery wall via endothelial cell (EC) receptor-mediated reactions include (but are not restricted to) adenosine nucleotides, thrombin, histamine and bradykinin (Jaffe *et al*. 1987; Tran *et al*. 2000; Yamamoto *et al*. 2000). In addition, blood-borne macromolecules such as lipoproteins are transported from the bloodstream to the arterial intima via trans-endothelial filtration (Phelps & DePaola 2000; Stangeby & Ethier 2002). The spatially varying concentrations of all such chemical species are determined by the relative rates of advection, diffusion and reaction, which are profoundly affected by the local fluid flow characteristics.

The binding of ligands such as adenosine triphosphate (ATP) and thrombin to EC surface receptors initiates a complex intracellular signalling cascade involving G-proteins, phospholipases, inositol trisphosphate (IP_{3}) and calcium (Tran *et al*. 2000). Thus, the delivery of ligands to the endothelium affects EC phenotype and function, as well as regulating vascular tone and platelet–endothelium adhesion (Shen *et al*. 1993). The delivery rate of these species is, therefore, an important factor in the onset of atherosclerosis, which is closely linked to endothelial dysfunction (Schachter 1997; Traub & Berk 1998). Similarly, the rate of mass transport of cholesterol-carrying macromolecules, such as low-density lipoprotein (LDL), to the endothelium, and subsequent transmission into the intima, is clearly a highly important factor in atherogenesis (Libby 2003). Therefore, spatial variations in the concentrations of endothelial ligands and of LDL will play a key role in the observed localized occurrence of atherosclerotic plaques in regions of disturbed flow and low wall shear stress (Gimbrone 1999).

The strong spatial correlation between disturbed flow characteristics (such as flow separation and recirculation, oscillatory flow and low time-averaged wall shear stress) and incidence of atherosclerosis is well known, but poorly understood (Malek *et al*. 1999; Shaaban & Duerinckx 2000; Resnick *et al*. 2003). A detailed understanding of the effects of different fluid flow characteristics on the mass transport of biochemicals in the bloodstream, and, in particular, their delivery to the endothelial wall, is an essential step towards understanding the spatial relationship between flow characteristics (in particular, wall shear stress) and atherogenesis. In this paper, we present a model for predicting the concentrations of generic biochemicals in the bloodstream and, importantly, at the endothelial wall, where receptor–ligand binding and/or transmission into the intima occur. The model assumes that the fluid flow is steady: clearly this is not the case *in vivo*. However, it is a reasonable first approximation, since the time-scale over which variations in EC phenotype and in intimal accumulation of LDL affect the progression of atherosclerosis is of the order of several years, whereas the time-scale of the cardiac cycle is approximately 1 s. A possible future model refinement would be to address the effects of pulsatile flow.

The model is quite general in the sense that it is applicable for any chemical species in an arbitrary two-dimensional boundary layer flow. The only restriction is that the Peclet number (i.e. the ratio of advection rate to diffusion rate) is high. This is a valid assumption for a very wide range of physiologically important factors with low diffusion coefficients, including adenosine nucleotides and lipoproteins. Under this restriction, the fluid flow characteristics and the resulting balance between advection, diffusion and reaction can be wholly characterized by the wall shear stress—the tangential force imposed by the viscous fluid flowing past the EC. Given that shear stress is observed to be one of the most important factors in the occurrence of atherosclerosis (Malek *et al*. 1999), this is highly advantageous as it enables the relationship between wall shear stress and the delivery of agonists and lipoproteins to the endothelium to be directly studied.

A number of computational fluid dynamical (CFD) studies that model the concentrations of particular chemical species in specific two- or three-dimensional arterial structures may be found in the existing literature. For example, Wada & Karino (2002) modelled LDL concentration in the human right coronary artery; Stangeby & Ethier (2002) modelled LDL concentration in a stenosed two-dimensional arterial geometry; Shen *et al*. (1993) and John & Barakat (2001) modelled adenosine nucleotide concentrations in a parallel-plate flow chamber. Other studies have looked at blood cell distribution and adhesion to the endothelium by treating the cells at a continuum level. For instance, Longest & Kleinstreuer (2003*b*) modelled the near-wall residence time of platelets in an end-to-side arterial anastomosis.

The model used in this study focuses on relatively simple two-dimensional flow geometries with spatially varying wall shear stress. These are not intended to replicate accurately the complex situation *in vivo*, which is possible only in numerically intensive CFD studies, but rather to identify the most important features of the generic behaviour. The advantage of this approach is that it captures the important flow characteristics present in realistic arterial geometries (flow detachment and regions of recirculating and fully developed flow separated by a stagnation point), yet is sufficiently simple that the model equations are mathematically tractable and the cause-and-effect relationship between shear stress, advection, diffusion and reaction is not masked by complex three-dimensional flow. This enables the model analysis to identify the most important flow characteristics affecting the delivery of ligands to the endothelium and to elucidate the underlying relationship between mass transport in the bloodstream and chemical concentrations and reaction rates at the endothelial wall.

In §2, the mathematical model for mass transport in a vessel with spatially varying wall shear stress is introduced. A pseudo-similarity solution to the model equations was originally proposed in David *et al*. (2001) for the distribution and flux of platelets and in David (2003) for adenosine nucleotide concentrations. In §3, it is demonstrated that this solution is a close approximation to the exact solution, provided certain conditions are met. In §4, some important physiological implications of the mathematical analysis are then discussed in the context of adenosine nucleotides and their interaction with EC. Finally, our conclusions are presented in §5.

## 2. The mathematical model

Under the assumption of steady, two-dimensional flow, the concentration of a chemical species in the bloodstream satisfies the following advection–diffusion equation:(2.1)where is the fluid flow vector, is the species diffusion coefficient and the reaction rate is permitted to be non-zero only on the boundaries. The diffusion coefficient is, in general, dependent on shear rate due to the enhancement of diffusive motion by red blood cell mixing in shear flow (Sorensen *et al*. 1999). However, as a first approximation, we use a constant value of ; the extension of the model to include shear-dependent diffusivity is left for future work.

A reactive surface is positioned at , representing the endothelium. On this surface, there is uptake (at rate ) and production (at rate ) of species via the following boundary condition:(2.2)

We assume that is small relative to (i.e. the Peclet number is large), in which case there is a thin mass transfer boundary layer at the reactive surface. The boundary layer thickness is *O*(*Pe*^{−1/3}), so, provided *Pe*≫1, this implies that the upper and lower boundaries will not interact, and, therefore, that the upper boundary condition can effectively be applied at infinity. The following boundary condition for the concentration in the bulk fluid is, therefore, adopted:(2.3)

We now define the following dimensionless variables and constants:(2.4)where is the wall shear stress on the surface , is a characteristic length-scale, is a characteristic fluid velocity and *ρ* and *ν* are, respectively, the density and kinematic viscosity of blood. Note that the *y*-coordinate has been scaled by *Pe*^{−1/3}, so that the change in chemical concentration with respect to *y* (i.e. the mass transfer boundary layer) occurs for values of *y* that are *O*(1).

The system (2.1)–(2.3) becomes(2.5)

(2.6)

(2.7)

Again, the thinness of the mass transfer boundary layer means that, in this region, the component of fluid velocity tangential to the wall increases almost linearly with *y*. This is valid provided that, in the dimensional coordinates, the thickness *y*_{BL} of the mass transfer boundary layer satisfies , where *a* is the radius of the vessel. This is true for , where *Q* is the flow rate through the vessel, since grows linearly with *x* at a rate proportional to . Then it can be shown by use of an appropriate stream function (David *et al*. 2001) that

Substituting into (2.5) gives

Finally, since the Peclet number is large, the term of order *Pe*^{−2/3} may be neglected to obtain(2.8)

Note that the transport of large particles, such as platelets and leukocytes, cannot be characterized solely by wall shear stress, due to inertia-driven contact between finite-size particles and the vessel wall (Longest & Kleinstreuer 2003*a*). However, for the transport of dilute chemical solutes (such as adenosine nucleotides), the problem is specified entirely in terms of a single function—the wall shear stress *τ*—of the stream-wise coordinate *x*. The wall shear stress may be chosen as a simple continuous function for an idealized fluid flow, such as planar or axisymmetric stagnation point flow or flow past a wedge, or defined in terms of CFD-calculated data for a particular two-dimensional geometry, such as a backward-facing step or an arterial bifurcation. The effects of different wall shear stress functions may be readily investigated, thus revealing fundamental aspects of the link between flow characteristics and transport of blood-soluble factors to the endothelium, and consequent blood-to-wall flux and EC signalling. In the case of flow with a stagnation point streamline, the *x*-coordinate is chosen so that *τ*(0)=0, i.e. the flow reattachment point corresponds to *x*=0. This is illustrated in figure 1, which shows the flow characteristics in a two-dimensional backward-facing step.

For a given wall shear stress, the problem is specified in terms of just two dimensionless parameters, *K* and *S*, defined by (2.4), which are essentially Sherwood numbers for the system.

## 3. Analysis

Let us first make the change of variables(3.1)This effectively divides the domain into two regions: *x*>0 and *x*<0. The change of variables (3.1) is invertible in each region. On dropping the asterisks, (2.8) becomes(3.2)with boundary conditions(3.3)

(3.4)

We seek a solution of the form *ϕ*=*ϕ*^{sim}+*u*, where *ϕ*^{sim} is the pseudo-similarity solution to (3.2)–(3.4) (David *et al*. 2001) given by(3.5)the similarity variable is(3.6)and(3.7)

(3.8)

This is not a true similarity solution because of the dependence of the coefficient *A* on *z*. Hence, although *ϕ*^{sim}(*z*) satisfies the boundary conditions (3.3) and (3.4), it does not satisfy the partial differential equation (3.2). Nevertheless, extensive numerical simulations have shown that *ϕ*^{sim}(*z*) provides an accurate approximation to the true solution over a wide range of parameter values. For example, figure 2 shows both the pseudo-similarity solution and the numerical solution to the system (3.2)–(3.4) in a backward-facing step (see figure 1). The wall shear stress on the bottom surface of the backward-facing step was calculated using the CFD package FLUENT.

We seek to show that the pseudo-similarity solution *ϕ*^{sim}(*z*) is a close approximation to the exact solution to the system (3.2)–(3.4) by showing that the error term *u* is small; *u* is the solution to the corresponding system(3.9)

(3.10)

(3.11)

(3.12)

Differentiating (3.7) with respect to *z* giveswhere *r*=*S*/*K*.

If *τ*=*ax*+*O*(*x*^{2}) for some *a*>0 and sufficiently small |*x*| (which is true for any physiologically realistic stagnation point flow), then and henceIt follows that lim_{z→0}*A*′(*z*)=0. Also, if lim_{z→∞}*τ*(*z*) is a finite positive constant, then lim_{z→∞}*A*′(*z*)=0. Let us now consider the dependence of the magnitude of *A*′(*z*) on the parameter *K* for *z* bounded away from zero and infinity:

Thus, (d/d*K*)|*A*′(*z*)|=0 whenand this value of *K* maximizes |*A*′(*z*)| at the particular value of *z*. Hence(3.13)where , provided the ratio *r* of production rate to uptake rate satisfies 0≤*r*≤2. Therefore, *ϵ* is a small parameter and (3.9) may be written as(3.14)where is the parabolic operator *ζ*∂_{z}−∂_{ζζ} andNow, seeking a solution in the form of a power series in the small parameter *ϵ*,it is apparent that *u*_{0}(*z*, *ζ*)≡0 is a solution to the partial differential equation (3.14) at order *O*(*ϵ*^{0}), satisfying the boundary conditions (3.10)–(3.12).

At *O*(*ϵ*^{1}), we haveUsing (3.13), we have an upper bound on |*G*(*z*, *ζ*)|:

Hence, if we can find a function *v* such thatand *v* satisfies the same boundary conditions (3.10)–(3.12) as *u*_{1}, then and it follows by use of a comparison principle (e.g. Jones & Sleeman 2003) that either *v*≡|*u*_{1}| or *v*>|*u*_{1}| for all *z* and *ζ*. Either way, *v* is an upper bound for the magnitude of the first-order correction term *u*_{1}.

We know that *H*(*z*)→0 as *z*→∞, so *v*(*z*, *ζ*)→0 as *z*→∞ because of the reactive boundary condition of the form (3.11) on *v*. Thus, *v*(*z*,0) will attain a maximum, *M*, at some *z*=*z*_{m}>0, which implies that the error term *u*(*z*, *ζ*) is not greater than *O*(*ϵM*). We have thus demonstrated that *ϕ*^{sim}(*z*, *ζ*) given by (3.5) is a leading-order solution to the system (3.2)–(3.4), which is asymptotically accurate in the limit *r*→1 (i.e. as the ratio of production to uptake tends to one). The pseudo-similarity solution is also valid for the case, where the uptake and production rates, *K* and *S*, are functions of *z*, provided |*K*′(*z*)| and |*S*′(*z*)| are sufficiently small and *S*(*z*)/*K*(*z*)≤*O*(1).

In summary, *ϕ*^{sim}(*z*, *ζ*) is an exact solution to the partial differential equation (3.2) for an arbitrary wall shear stress function, *τ*(*x*), only when *A*(*z*) is independent of *z*, which occurs when *r*=1 or, in dimensional terms, . The latter condition means that the concentration in the bulk fluid is equal to the equilibrium concentration on the reactive boundary, and hence there is no change in the concentration as the fluid is advected downstream because it is already at equilibrium. Intuitively then, the reason that *ϕ*^{sim}(*z*, *ζ*) is a leading-order solution is that *A*(*z*) varies relatively slowly with respect to *z*, no matter how small or large *K* and *S* are, provided that their ratio *r* satisfies *r*≤*O*(1). In this case, the pseudo-similarity solution is ‘close’ to being an exact similarity solution to (3.2)–(3.4) and is asymptotically accurate in the limit *r*→1.

## 4. Physiological implications

In this section, some physiological implications of the pseudo-similarity solution are discussed. The first thing to note is that, downstream of the flow disturbance as the flow approaches a fully developed parabolic profile, the concentration on the endothelium approaches an equilibrium value *ϕ*_{s}. The dimensional value of this equilibrium is the ratio of the rate of production to uptakeIf the bulk concentration is greater than this equilibrium value, then the concentration will be high near the stagnation point, where there is rapid advection from the bulk fluid to the mass transfer boundary layer, and will gradually decay towards the equilibrium value away from the stagnation point. Conversely, if the bulk concentration is less than the equilibrium, the concentration will be low near the stagnation point and will gradually increase towards the equilibrium value.

We now address more specific questions regarding the occurrence of maximum chemical concentration on the endothelium. Attention is focused on adenosine nucleotides, adenosine triphosphate (ATP), diphosphate (ADP) and monophosphate (AMP), as these are known to be important regulators of endothelial function and vascular tone (Shen *et al*. 1993), although the analysis is applicable to any chemical species transported in the bloodstream. Adenosine nucleotides are sequentially dephosphorylated by EC-derived nucleotidases in the reaction sequence (Dull & Davies 1991; Shen *et al*. 1993):(4.1)AMP and adenosine have little or no effect on EC signalling or the regulation of vascular tone (Shen *et al*. 1993), so we focus on ATP and ADP. In this section, we consider the effects of the fluid transport, via the wall shear stress function, on the distribution of these nucleotides along the vessel wall and, in particular, on the locations of maximal chemical concentration. This is of importance in the localization of endothelial dysfunction and consequent atherogenesis to certain focal sites within the vasculature.

In §4*a*, a fundamental relationship is identified between the ‘shape’ of the wall shear stress function at the stagnation point and the position of maximum ATP concentration relative to the stagnation point. An expression for the value of the ATP concentration at this point is also found in terms of the spatial gradient of the wall shear stress. In §4*b*, the pseudo-similarity solution is extended to give the concentration of ADP, which is produced by EC-mediated dephosphorylation of ATP (i.e. as a result of ATP hydrolysis on the reactive surface). It is shown that maximum ADP concentration can either occur at the same point on the endothelium as the maximum ATP concentration, or at a separate point on the endothelium. A condition is derived for which of these cases occur, again in terms of the wall shear stress gradient. Finally, an expression is also found for the maximum concentration of ADP.

### (a) Location of maximum ATP concentration

Returning to the original coordinate system (*x*, *y*), the concentration *ϕ*_{1}(*x*, 0) of ATP on the endothelium, as given by (3.5)–(3.8), is(4.2)where *S*_{1} and *K*_{1} are the production rate of ATP and the rate of the ATP→ADP reaction in (4.1), respectively, and(4.3)

In general, there will be a relatively high concentration of ATP in the bulk fluid, and a thin boundary layer at the endothelium in which ATP is dephosphorylated by nucleotidases (Shen *et al*. 1993). Thus, (or in non-dimensional terms *K*_{1}>*S*_{1}) and the ATP concentration will be high near the stagnation point and will decay towards the equilibrium value (which will be zero if there is no production of ATP at the endothelium) downstream of the flow disturbance.

As shown in §3, the pseudo-similarity solution given by (4.2) and (4.3) is valid for all *x*, provided that *τ*(*x*)=*ax*+*O*(*x*^{2}) for sufficiently small |*x*|. However, the shape of the concentration profile, and, in particular, the location of the maximum concentration, depends on the higher-order terms in the Taylor series for *τ*(*x*) at *x*=0. Since *K*_{1}>*S*_{1}, *ϕ*(*x*, 0) is an increasing function of *β*(*x*) and hence the maximum of *ϕ*(*x*, 0) occurs at the maximum of *β*(*x*). Suppose that, for sufficiently small |*x*|,(4.4)Provided that |*b*/*a*|≪1, it follows that the Taylor series for *β*(*x*) at *x*=0 is given byand hence that the gradient of *β*(*x*) at *x*=0 isThus, if the Taylor series (4.4) for *τ*(*x*) about *x*=0 contains an *O*(*x*^{2}) term (i.e. the wall shear stress function has non-zero curvature at the stagnation point), then the position, say *x*=*x*_{m}, of the maximum ATP concentration will be offset from the stagnation point. If the wall shear stress has negative curvature at the stagnation point (*b*<0), then the maximum will occur to the left of the stagnation point (*x*_{m}<0); if it has positive curvature (*b*>0), it will occur to the right of the stagnation point (*x*_{m}>0); if the wall shear stress has zero curvature (*b*=0), the maximum will occur at the stagnation point (*x*_{m}=0).

The derivative of *β*(*x*) with respect to *x* is given bywhere *z* is given by (3.1). It follows that the maximum value of *β*(*x*) is(4.5)which gives the corresponding maximum ATP concentration according to (4.2).

Computational solution of the advection–diffusion equation (2.1) is now used to verify the above results. Figure 3 shows the wall shear stress *τ*(*x*) and the ATP concentration *ϕ*_{1}(*x*, 0) (calculated using the CFD package FLUENT) on the endothelium in a backward-facing step for four different Reynolds numbers. We are interested in the position of maximum ATP concentration relative to the flow reattachment point, so the *x*-coordinate has been shifted such that the reattachment point (i.e. the point where *τ*=0) occurs at *x*=0 in each case, and *τ*(*x*) and *ϕ*_{1}(*x*, 0) are plotted in a small region around *x*=0. Interestingly, the maximum concentration occurs on different sides of *x*=0 for different Reynolds numbers. In order to confirm the theoretical relationship between the curvature of the wall shear stress at *x*=0 and the position of the maximum concentration, the second derivative of *τ*(*x*) with respect to *x* was calculated at *x*=0 for a range of Reynolds numbers. The results are summarized in table 1 and the sign of *x*_{m} is the same as the sign of *τ*″(0) in each case.

During the cardiac cycle, the pulsatile nature of the pressure gradient causes a variation in the effective Reynolds number of the flow through a given vessel. In a region exhibiting a stagnation streamline, this will result not only in movement of the point of flow reattachment, but, as our analysis and numerical simulations show, in movement of the position of maximum ATP concentration relative to the reattachment point. In an area exhibiting flow recirculation, it appears that, at relatively low Reynolds number, the ATP maximum occurs downstream of the reattachment point, whereas at high Reynolds number, the ATP maximum occurs in the recirculation zone.

Finally, note that in the case of chemical species, where *K*<*S* (i.e. where there is a net production on the endothelium), the opposite situation occurs: *ϕ*(*x*, 0) is a decreasing function of *β*(*x*). Therefore, the point *x*=*x*_{m} corresponds to a minimum in *ϕ*(*x*, 0).

### (b) Location of maximum ADP concentration

Unlike ATP, the bulk concentration of ADP is very low (or zero) and the situation in the boundary layer is reversed as ADP is produced as a result of ATP hydrolysis. We show that the position of maximum ADP concentration is determined by the reaction rates *K*_{1} and *K*_{2} in (4.1). There is now an additional dimensionless parameter, namely the ratio of the ADP and ATP diffusion coefficients. Assuming that there is no production of ATP by the EC (i.e. *S*_{1}=0) and that the concentration of ADP in the bulk fluid (i.e. as *y*→∞) is zero, the concentrations on the endothelium are given by

According to (4.1), the (dimensional) rate of ADP production is equal to the (dimensional) rate of hydrolysis of ATP. The corresponding dimensionless value isNote that *S*_{2} is a function of *x*, but by the arguments of §3, the pseudo-similarity solution is still valid. The ADP concentration on the endothelium is, therefore,and differentiating with respect to *x* we find that

Hence *ϕ*_{2}(*x*, 0) has stationary points where d*β*/d*x*=0 (i.e. at the maximum ATP concentration) and at points whereThere are, thus, two possibilities: (a) *β*_{max}>*β*_{c}; (b) *β*_{max}<*β*_{c}.

In case (a), the maximum in *β*(*x*) at *x*=*x*_{m} corresponds to the maximum ATP concentration (as shown in §4*a*) and a local minimum in the ADP concentration. The point at which *β*(*x*)=*β*_{c} corresponds to a maximum in the ADP concentration and, at this point, the concentration is given byThis maximum occurs at a different position on the endothelium to the maximum ATP concentration. This is illustrated in figure 4*a*, which shows the ATP and ADP concentration profiles in a backward-facing step for a Reynolds number of 5 (note that the reaction rates *K*_{1} and *K*_{2} in figure 4 are higher than the values shown in table 2 for illustrative purposes). The maximum ATP concentration occurs at , as calculated in §4*a*. The maximum ADP concentration occurs where *β*(*x*)=*β*_{c}=0.14, which is at . The reason that the ADP maximum occurs away from the ATP maximum is that, near the stagnation point, there is a relatively high rate of advection of ATP to the endothelium, resulting in net production of ADP via the ATP→ADP part of (4.1). However, the ATP concentration, and hence the rate of ADP production, gradually decays as one moves downstream. Where production of ADP from ATP is balanced by ADP hydrolysis, the ADP concentration attains a maximum. Further downstream, the ADP→AMP part of (4.1) dominates the ADP kinetics and the concentration of ADP thereafter declines.

In case (b), the maximum in *β*(*x*) corresponds to maxima in the concentrations of both ATP and ADP, which hence occur at the same point on the endothelium. This is illustrated in figure 4*b*, which shows the nucleotide concentration profiles for very high reaction rates *K*_{1} and *K*_{2}.

Case (a) occurs if *β*_{max}>*β*_{c}. Using the value of *β*_{max} given by (4.5), it follows that case (a) occurs ifIn dimensional terms, this condition is equivalent to(4.6)So if the spatial gradient of the wall shear stress at the point, *x*=*x*_{m}, of maximum ATP concentration is less than this threshold value, then the maximum ADP concentration will occur at the same point *x*=*x*_{m}; otherwise, it will occur away from *x*=*x*_{m}. Thus, the behaviour depends on the ‘strength’ of the stagnation point flow relative to the effective Sherwood numbers of the reaction (4.1). The reason for this is that if the wall shear stress gradient near the stagnation point is high, then there is rapid advection of ATP from the bulk fluid into the mass transfer boundary layer. As a consequence, there is rapid production of ADP near the stagnation point and ADP concentration initially increases as it is advected downstream. However, if the wall shear stress gradient is low, advection of ATP from the bulk fluid to the boundary layer is weak and there is always a net loss of ADP via hydrolysis to AMP, resulting in a monotonic decreasing concentration profile.

For the parameter values shown in table 2, the threshold value (4.6) of the wall shear stress gradient is approximately 0.3 Pa m^{−1}. Physiologically, this is extremely low and one would, therefore, expect case (a) to be the dominant phenomenon. Furthermore, the low value of *β*_{c} (approximately 0.01 for the parameter values in table 2), coupled with the fact that *β*(*x*) decays slowly in the region of fully developed flow (*β*(*x*)∼*x*^{−1/3}), means that the theoretical maximum of ADP concentration will occur an extremely long way (of the order of several metres) downstream of the flow disturbance. Hence, on a realistic length-scale, the ADP concentration *ϕ*_{2}(*x*, 0) is expected to be monotonic increasing with *x*.

## 5. Discussion

There are several examples in the literature of large-scale computational fluid dynamic models of the mass transport of chemical species in particular vascular geometries (Stangeby & Ethier 2002; Wada & Karino 2002). However, such studies, which are computationally intensive by nature, cannot reveal underlying links between macroscopic fluid flow characteristics and the balance between advection, diffusion and reaction of chemical species, and their physiological consequences, such as EC signalling and intimal accumulation of lipoproteins. The model presented in this paper focuses on a relatively simple two-dimensional geometry, which, although not intended as a complete model of complex three-dimensional arterial structures, does capture the essential macroscopic characteristics of physiological blood flow. These consist of flow detachment and a resulting region of secondary recirculation and point of flow stagnation. Downstream of the stagnation point, the flow steadily approaches a fully developed parabolic profile. This results in spatially varying wall shear stress; a major advantage of the model is that the wall shear stress function, *τ*(*x*), captures all the essential flow characteristics. In this way, the model is capable of addressing important questions regarding the relationship between wall shear stress and the physiology of the arterial wall. For instance, blood-to-tissue transport of LDL is a key factor in the onset of endothelial dysfunction and atherogenesis (Schachter 1997). Spatial variation in the concentration of LDL is at least partly responsible for the localization of atherogenesis to certain focal sites within the vasculature (Gimbrone 1999), and predicting the location on the vessel wall of the maximum concentration of LDL is of central importance in understanding this phenomenon. The model is applicable for all chemical species whose associated Peclet number (i.e. the ratio of advection rate to diffusion rate) is sufficiently large, which includes many blood-soluble species of physiological importance, such as adenosine nucleotides and lipoproteins.

The pseudo-similarity solution proposed by David *et al*. (2001) provides an exact solution to the advection–diffusion equation in both the diffusion-limited (i.e. large *K*) and reaction-limited (i.e. small *K*) cases. In this paper, it has been shown that this solution is a close approximation to the concentration profile in the vessel when the transport is intermediate between these two limiting cases. Furthermore, the approximation is valid even when the reaction rates are spatially varying. This is advantageous as it allows situations in which endothelial function is regulated by external stimuli to be modelled. For instance, it is known that the permeability of the endothelium to lipoproteins is modulated by local wall shear stress (Davies 2000), and this can be modelled by taking the uptake rate *K* to be a function of *τ*(*x*). Another possibility is that shear stress- or ligand-induced intracellular signalling results in spatially varying reaction rates. Endothelial calcium signalling in response to ATP and shear stress has been modelled in Plank *et al*. (in press). Calcium signalling is known to affect synthesis of nitric oxide (Lin *et al*. 2000), which is an important vasodilator (Shaul 2003), and this could be modelled by taking the nitric oxide production rate to be dependent on the concentration of free calcium within the cell. These processes will be addressed in a future paper.

The pseudo-similarity solution has been used to identify the physiological implications of the relationship between wall shear stress and the spatial distribution of chemical species along the vessel wall. In particular, it was found that, for species such as ATP and LDL that are present in the bulk fluid and are consumed at the endothelium, the location of maximum chemical concentration on the vessel wall can occur on either side of the stagnation point, depending on the sign of the curvature of the wall shear stress function at the stagnation point. Furthermore, the sign of the curvature can change as the effective Reynolds number changes during the cardiac cycle, resulting in oscillation not only in the position of the flow reattachment point, but also in the position of maximum concentration relative to the reattachment point.

In the case of ADP, which is produced as a result of ATP hydrolysis at the endothelium, two cases are possible. If the spatial gradient of wall shear stress at the point of maximum ATP concentration is less than the threshold value given by (4.6), then the maximum ADP concentration occurs at the same point as the maximum ATP concentration. Otherwise, the maximum ADP concentration will occur away from the maximum ATP concentration. For physiological wall shear stress, reaction rates and diffusivities, the model predicts that the latter case will occur, and that, on a realistic length-scale, the ADP concentration is likely to be monotonic increasing with the distance downstream of the region of disturbed flow.

As a final caveat, it should be noted that these results are based on the assumption that the ATP/ADP diffusivity in the bloodstream is constant and equal to the diffusivity in plasma. The presence of red blood cells in the flow can cause a significant shear-dependent enhancement of diffusive motion. For large particles, such as platelets and leukocytes, the effective diffusion coefficient in shear flow containing red blood cells can be 2–3 orders of magnitude greater than the Brownian motion value (Sorensen *et al*. 1999). The effect of red blood cell mixing on smaller particles, such as adenosine nucleotides, is likely to be less pronounced and the effective Peclet number will still be large. The effects of a shear-dependent diffusion coefficient on the model results will be addressed in a future paper.

## Footnotes

- Received April 28, 2005.
- Accepted October 12, 2005.

- © 2005 The Royal Society