## Abstract

Peripheral arterial disease (PAD) is generally attributed to the progressive vascular accumulation of lipoproteins and circulating monocytes in the vessel walls leading to the formation of atherosclerotic plaques. This is known to be regulated by the local vascular geometry, haemodynamics and biophysical conditions. Here, an isogeometric analysis framework is proposed to analyse the blood flow and vascular deposition of circulating nanoparticles (NPs) into the superficial femoral artery (SFA) of a PAD patient. The local geometry of the blood vessel and the haemodynamic conditions are derived from magnetic resonance imaging (MRI), performed at baseline and at 24 months post intervention. A dramatic improvement in blood flow dynamics is observed post intervention. A 500% increase in peak flow rate is measured *in vivo* as a consequence of luminal enlargement. Furthermore, blood flow simulations reveal a 32% drop in the mean oscillatory shear index, indicating reduced disturbed flow post intervention. The same patient information (vascular geometry and blood flow) is used to predict *in silico* in a simulation of the vascular deposition of systemically injected nanomedicines. NPs, targeted to inflammatory vascular molecules including VCAM-1, E-selectin and ICAM-1, are predicted to preferentially accumulate near the stenosis in the baseline configuration, with VCAM-1 providing the highest accumulation (approx. 1.33 and 1.50 times higher concentration than that of ICAM-1 and E-selectin, respectively). Such selective deposition of NPs within the stenosis could be effectively used for the detection and treatment of plaques forming in the SFA. The presented MRI-based computational protocol can be used to analyse data from clinical trials to explore possible correlations between haemodynamics and disease progression in PAD patients, and potentially predict disease occurrence as well as the outcome of an intervention.

## 1. Introduction

Peripheral arterial disease (PAD) is often referred to as a vascular disease inducing the progressive formation of plaques (atherogenesis), specifically in the arteries that carry blood to the limbs, head, kidneys and stomach. PAD has been associated with an increased risk of atherothrombotic cardiovascular events and mortality and, as such, it is considered an early marker of coronary artery disease, strokes and myocardial infarction [1–4]. PAD management can be categorized into pharmacological (lipid-lowering treatments), surgical (deployment of stents, revascularization and vascular surgery) and a combination of the two [1,5]. Following the well-accepted oxidative modification hypothesis, the first step in atherogenesis is the accumulation of circulating low-density lipoprotein particles into the sub-intimal space of the arterial wall [6]. This accumulation is generally higher at sites of disturbed flow: low flow, separation, flow reversal and even turbulence [7–12]. At these sites, typically, larger rates of endothelial cell (EC) proliferation and apoptosis, higher vascular permeability, overexpression of EC adhesion molecules (ICAM1, VCAM1 and E-selectin) and the release of chemo-attractants from ECs are consistently observed [13,14]. This clearly emphasizes the importance of an accurate determination of the local flow field to potentially predict the occurrence of the disease and the outcome of any intervention. The main objective of this work is to bring together a number of available technologies to build a patient-specific analysis framework that can be used to analyse data from clinical trials to provide insights into PAD management. To that end, this paper aims to serve as a proof of concept and presents some relevant results to demonstrate the feasibility of our approach.

Finite-element-based isogeometric analysis is one of the prominent approaches [15] that can be efficiently used to simulate cardiovascular biomechanics problems in patient-specific three-dimensional vascular models [16–20] (for a more comprehensive review, readers are referred to the following, to cite but a few: [9,10,21–23]). The isogeometric framework uses geometric primitives (NURBS: non-uniform rational B-splines) developed in computer-aided design and computational geometry, which allows us to represent geometrical objects with a higher degree of precision and smoothness (higher order basis functions with continuous derivatives) than piecewise polynomials (continuous, but not smooth, basis functions) commonly used in finite-element analysis [24]. The ability of NURBS to accurately represent smooth exact geometries (i.e. arterial systems), which is not attainable in the faceted finite-element representation, renders fluid and structural computations more physiologically realistic. Also, the parametric definition of NURBS meshes allows automatic refinement of the boundary layer region near the arterial wall, which is crucial for the overall numerical stability and for the accurate computation of wall quantities, such as the wall shear rates and stresses, wall permeability and adhesion [16]. Following this approach, the authors have solved several structural, fluid and solid mechanics problems within the biomedical field and others [25–27]. The methodology has been shown to be suitable for cardiovascular applications because of its precise and efficient geometric modelling, capability to appropriately capture both laminar and turbulent flow regimes, and accurate representation of stresses [28].

Magnetic resonance imaging (MRI) provides excellent soft-tissue contrast without ionizing radiation and it has been extensively used to image atherosclerotic plaques and plaque composition [29]. However, the study of atherosclerosis in lower extremity vessels, including the superficial femoral artery (SFA), presents various MRI challenges, due to potential involuntary muscle twitching causing motion artefacts, vessel calibre and tortuosity [30]. Although recent advances in coil technology allow full coverage of the SFA, sampling a large enough field of view (FOV) remains challenging in the domain of high-resolution plaque imaging. New three-dimensional MRI sequences such as three-dimensional motion-sensitized driven equilibrium prepared rapid gradient echo (3D-MERGE) have shown promising results for imaging of the peripheral vessels. Multi-contrast MRI, two-dimensional phase contrast (PC), four-dimensional flow MRI and magnetic resonance (MR) angiography have been described extensively [31–33]. These sequences can be used to obtain patient-specific information including the vascular geometry, inflow/outflow velocity profiles and the velocity field within the vascular tree of interest over time. These data can be readily used to prescribe the physical domain (geometry), the initial and boundary conditions in any patient-specific computational framework.

In this work, the geometry and haemodynamic conditions for an SFA of a PAD patient are extracted from MRI at baseline and 24 months after stent implantation. Then isogeometric analysis is used to analyse the blood flow dynamics within the SFA over a full cardiac cycle. Near-wall quantities such as time-averaged wall shear stress (TAWSS), the oscillatory shear index (OSI) and relative residence time (RRT) are also investigated. Finally, the vascular transport of systemically injected nanoparticles (NPs), targeted to different inflammatory vascular molecules such as VCAM-1, E-selectin and ICAM-1 that are known to overexpress at the diseased site, is analysed. In particular, the surface density of these NPs within the SFA is estimated as a function of time and vascular molecules (receptors) to explore the use of nanomedicine in the non-invasive imaging and treatment of PAD patients. The possible relationship between local haemodynamics and the particle distribution pattern is also studied.

## 2. Material and methods

### 2.1. Imaging data acquisition and magnetic resonance imaging

The MRI data used in this analysis were obtained from the Effect of Lipid Modification on Peripheral Artery Disease after Endovascular Intervention Trial (ELIMIT) [30]. ELIMIT was a double-blind and double-placebo-controlled trial in PAD patients with the goal of studying the efficacy of intensive lipid-modifying triple therapy (simvastatin, ezetimibe and extended-release niacin) versus standard lipid-modifying mono-therapy (simvastatin) on the progression of atherosclerosis in the SFA as quantified by high-resolution MRI. The main findings of ELIMIT have been reported previously [30]. Briefly, patients underwent MRI of the distal SFA at baseline, six months, 12 months and 24 months with a 3.0 T system (Signa Excite; GE Healthcare, Milwaukee, WI, USA) using a unilateral phased array coil with a FOV of 8 cm (along the *z*-axis) and 12 cm (in-plane *x*- and *y*-axes; Pathway Biomedical, Inc.). Lipid data were also recorded. The geometries for the current analysis were extracted from proton-density-weighted (PDW) fast spin-echo scans. Sequence parameters of the PDW scans were described previously [30]. ELIMIT participants were also imaged using gated two-dimensional PC sequences, acquired at select locations within the corresponding PDW volumes—typically proximally and distally to SFA lesions. Gated two-dimensional PC sequences were acquired during the same examination with slice thickness = 4 mm, repetition time (TR) = 10.6 ms, echo time (TE) = 4.97 ms, echo train length = 1, bandwidth = 244 Hz pixel^{−1}, 20 frames per cardiac cycle and PC-encoding velocity (VENC) of 120 cm s^{−1}. Serial MRI scans were carefully co-registered across follow-up visits using anatomical landmarks (artery, vein and muscle) [30]. Velocity profiles over the cardiac cycle were obtained at each location using the luminal cross-sectional areas (region of interest (ROI)), the mean signal intensity of the ROIs and the VENC factor.

### 2.2. Reconstructing the patient-specific vascular geometry

Here, an image-based pipeline, as detailed in [24], is adopted to construct hexahedral solid NURBS meshes for patient-specific vascular geometric models to be used in isogeometric analysis. Briefly, the image quality is first improved by passing the raw imaging data through a preprocessing pipeline for enhancing the contrast, filtering noise, classifying and segmenting regions of various materials. In this work, the segmentation of the lumen and the vessel wall is performed manually with Amira. The target surfaces (the lumen–wall interface) are then extracted by isocontouring the pre-segmented imaging data, followed by the extraction of the vascular skeleton (vessel centrelines) using open-source software [34,35]. Finally, an optimized skeleton-based sweeping method is used to construct hexahedral control meshes employing an in-house numerical code [24] as follows.

The geometry is first decomposed into mapped meshable piecewise linear patches. Each patch is then meshed using a heuristic sweeping method by iteratively running three steps: (i) clustering projected control points to obtain an optimized centreline; (ii) computing the normal of the cross section at each control point on the skeleton (centre points); and (iii) projecting boundary points onto the target surface preserving perpendicularity of cross sections and centreline normal. The first step improves the distribution of centre points on the skeleton, and the second and third steps guarantee good shape of hexahedral elements. The iterative procedure serves as a compromise of these two constraints. Such heuristics could be formulized as an alternating optimization problem [36] whose convergence has been thoroughly studied and thus quality meshes are obtained. Finally, hexahedral solid NURBS are constructed from the control mesh and used in isogeometric analysis to simulate blood flow and particle delivery in the SFA. Figure 1 shows the right SFA geometry generated from MR images of a patient with PAD.

### 2.3. Governing equations for the fluid flow and particle transport

A continuum-based approach was adopted to simulate blood flow and particle transport within the patient-specific SFA. The details of the mathematical equations and the numerical methodology are reported elsewhere [16,19,25,37]. Briefly, considering blood as an incompressible Newtonian fluid with a dynamic viscosity (*μ*) of 0.003 Pa-s and a density (*ρ*) of 1060 kg m^{−3}, blood flow was assumed to be governed by the unsteady Navier–Stokes equations for incompressible flow subjected to appropriate boundary conditions (figure 2). An inflow velocity profile was specified at the inlet, a no-slip boundary condition was prescribed at the rigid and impermeable wall, and a traction-free outflow boundary condition was set at the branch outlet. In order to analyse the flow-field characteristics and its near-wall behaviour, two wall shear-based quantities were calculated over a cardiac cycle with a time period *T*: (i) the TAWSS
2.1and (ii) the OSI [38,39]
2.2The surface traction vector **t**_{s}, defined as the tangential component of the traction vector, is calculated from **t**_{s} = **t** − (**t** · **n**)**n**, where the traction vector **t** is computed from the stress tensor **σ** using the relation **t** = ** σn**. Here, the OSI characterizes the changes in flow direction and the velocity gradient during a cardiac cycle and can range from 0 (unidirectional flow) to 0.5 (oscillatory flow). RRT is associated with the particle residence time near the wall and combines the information provided by the TAWSS and OSI as follows:
2.3

High RRT leads to increased uptake of blood-borne particles, which has implications for atherosclerosis [40].

Treating the NPs as passive scalars, the mass transport of the particles is assumed to be governed by an unsteady scalar advection–diffusion equation. This assumption is deemed reasonable for sufficiently small particles travelling in large vessels [41,42]. A bolus of NPs is located at the inlet of the artery segment for the entire duration of simulation (figure 2). This translates to a Dirichlet boundary condition with a volumetric concentration *C*^{0} prescribed at the inlet, where *C*^{0} denotes NP concentration within the bolus. At the outflow, a homogeneous Neumann boundary condition was specified. A special Robin-type boundary condition developed by Hossain *et al*. [16] is prescribed at the lumen–wall interface to account for NP deposition to the vessel wall. The mass flux of particles diffusing through the lumen–wall interface and adhering to the vessel wall is assumed to be a function of particle size *d*_{p}, local wall shear rate *S* and probability of adhesion *P*_{a}, defined as the probability of having at least one ligand–receptor bond. The parameter *P*_{a} can be considered as a measure of particle avidity towards the vessel wall or the strength of adhesion. That is, the larger the value of *P*_{a}, the greater is the avidity and strength with which the particle firmly adheres to the wall.

The propensity of a particle to adhere stably to the vessel wall withstanding dislodging hydrodynamic forces is strongly influenced by the physico-chemical properties of the NPs. Accordingly, Decuzzi & Ferrari [43] modelled *P*_{a} as a function of physiological parameters such as target receptor density and local wall shear stress (WSS), as well as certain design parameters including particle shape, size, ligand type and density. This particle firm-adhesion model, which has been validated *in vitro* in a previous work [16] for spherical particles in point contact with the vessel wall, is used herein to account for NP wall deposition. A summary of the modelling details, including the mathematical expressions and parameter values used, is provided in appendix A. For a more comprehensive description of the modelling approach and parameter selection, see [16,44] and references therein.

#### 2.3.1. Solution approach

The governing equations were solved by applying finite-element-based isogeometric analysis [37,45] using quadratic NURBS to describe the exact geometry as well as the solution space [17,24]. A residual-based variational multi-scale method [18] was implemented to solve the system of equations, using a Newton–Raphson procedure with a multi-stage predictor–corrector algorithm applied at each time step. The generalized−*α* method [46,47], an implicit second-order time-accurate method that is also unconditionally stable, was used for time advancement. Readers are referred to the numerical procedures described in [17,25,45,48,49] for further details.

A portion of the right SFA in a patient affected by PAD has been reconstructed from MR images using a hexahedral solid NURBS model. A representative case is shown in figure 1 where the three-dimensional MRI reconstruction and three cross sections of the SFA lumen and blood vessel wall are displayed. Similar MRI data were extracted from the PAD patient followed-up longitudinally over a period of 24 months. The baseline geometry and the corresponding problem set-up are presented in figure 2, together with the initial and boundary conditions, and the governing equations. The patient-specific mean inflow velocity waveform extracted via PC MRI with a period *T* of 1 s (heart rate = 60 beats min^{−1}) was imposed at the corresponding SFA inlet assuming a paraboloid-shaped profile with a maximum at the centroid of the vessel inlet cross section and zero values at the boundary. The inflow velocity vectors were oriented normal to the inlet plane. The simulations were run on a computational mesh consisting of 25 992 quadratic NURBS elements (baseline case) for 10 cardiac cycles (10 s total) with a time step of 0.05 s employing the solution strategies described above. Boundary layer meshes were used with the finest boundary element thickness of the order of 10^{−6} cm for more accurate computation of wall quantities of interest, namely the WSS and OSI. The flow remained laminar throughout the cardiac cycle. The upper bound of the Reynolds number was calculated to be 774 and 1642 for the baseline and the 24-month post-intervention cases, respectively.

## 3. Results and discussion

### 3.1. Blood flow analysis

In figure 3, the velocity magnitude is presented for baseline (top) and 24-month post-intervention case at three different cross sections, namely proximal (slice I), distal (slice III) and an intermediate location near the stenosis referred to as ‘mean’ (slice II). In addition, three different times in the cardiac cycle are depicted: (A) end diastole, (B) post-peak systole and (C) end systole. Significant luminal enlargement is observed post intervention because of stent implantation, resulting in a 453%, 313% and 264% increase in cross-sectional area at the stenosis (‘mean’), proximal and distal locations, respectively. It is worth noting that the SFA segment enclosed between the proximal and distal slices represents overlapping regions across follow-up visits. That is, the proximal, mean and distal slices in the baseline case essentially correspond to the same locations in the 24-month post-intervention geometry.

At baseline, the velocity magnitude levels seem to drop notably from the proximal to the distal sections consistently during the cardiac cycle because of flow restriction at the stenosis. However, this phenomenon appears to be less pronounced during peak systole (see electronic supplementary material, figure S1). Mean flow velocity magnitude measured via PC MRI (see electronic supplementary material, figure S2*a*) also shows markedly decreased values distal to the SFA narrowing. While a tri-phasic flow profile typical of healthy subjects is not preserved above and below the lesion due to the diseased condition, there is no peak delay. On the other hand, because of the implanted stent in the 24-month post-intervention SFA geometry blood flow is no longer restricted, resulting in a 500% increase in peak flow rate (see electronic supplementary material, figure S2*b*). Consequently, there is less of a variation in velocity magnitude between proximal and distal slices in contrast with the baseline case. It also appears that a bi-phasic waveform is recovered 24 months post intervention (see electronic supplementary material, figure S2*c*).

In figure 4, the streamlines coloured according to velocity magnitude are depicted for the baseline (top row) and 24-month post-intervention cases. Figure insets represent magnification of the region downstream of the stenosis (‘mean’ slice) at various times during the cardiac cycle as departures from unidirectional flow are known to occur distal to stenosis [50]. The time evolution of the streamline trajectories in figure 4 clearly captures the pulsatile nature of the blood flow. In the baseline case, driven by the inflow condition, the general flow direction reverses between (A) *t* = 0.2 s and (B) *t* = 0.4 s as it enters the systolic phase of the cardiac cycle, only to reverse back again during diastole at (C) *t* = 0.6 s. Slightly after peak systole (B), streamlines downstream of the stenosis appear to be tangled, indicating disturbed flow. Such localized complex flow features are notably absent in the corresponding post-intervention case (see bottom row of figure 4). In fact, streamlines in general appear to be fairly regular post intervention, indicating largely unidirectional flow in the highlighted region. At the advent of the systolic phase (A), there is a hint of change in flow direction as it transitions from the diastolic phase, resulting in a slightly disturbed flow downstream of previous stenosis (the ROI). After peak systole (B), during the downstroke of the systole phase, a reversal of flow occurs. At end systole (C), the streamlines become mildly perturbed as they once again try to adjust to the backflow during the diastolic phase.

Figures 5 and 6 report the spatial distribution of the TAWSS and OSI, respectively. In the baseline case, alternate areas of lower (less than 1 Pa) and higher (more than 5 Pa) TAWSS regions appear near the stenosis (figure 5). Twenty-four months after stent implantation, relatively lower levels of TAWSS magnitude (less than 5 Pa) are observed that also appear to be more uniformly distributed throughout the SFA segment when compared with the baseline case. Similarly, as a result of a more regularized blood flow in the post-intervention geometry, the distribution pattern and magnitude for the OSI is altered appreciably (figure 6). At 24 months post intervention, OSI levels appear to diminish considerably near the original constriction, in some areas by a factor of 2, when compared with the baseline case. It appears that regions of relatively higher OSI generally correspond to lower WSS values (see electronic supplementary material, figure S3). This seems to be largely the case for the 24-month post-intervention geometry where the flow is essentially unidirectional. Interestingly, in the disturbed flow region downstream of the stenosis in the baseline case, this inverse relationship between the TAWSS and OSI does not seem to hold.

Figure 7*a*,*b* shows histograms of TAWSS and OSI distributions, respectively, for both baseline (light grey) and 24-month post-intervention (dark grey) cases. Each bar of the histogram represents the amount of normalized area with a defined range of each quantity of interest. It is worth noting that the *x*-axis labels denote the upper bound of each range. For example, in figure 7*a*, the first bar has an *x*-axis label of ‘1’, indicating that, for the 24-month post-intervention distribution (24 M PI), 12% of the arterial wall area has a TAWSS between the values 0 and 1. The mean, the skewness and the kurtosis [40] of the TAWSS and OSI distributions are presented in the electronic supplementary material, table S1. Mean TAWSS is reduced by 22% from the baseline to the 24-month post-intervention simulations. Baseline TAWSS distribution (figure 7*a*) has a positive skewness (approx. 1.27), indicating a distribution with an asymmetric tail extending towards more positive (higher TAWSS) values from the mean. This can be attributed to the higher TAWSS values present at the stenosed region at baseline (see figures 5 and 7*a*). Stent implantation induced a 40% decrease in the maximum TAWSS (from 7 to 4 Pa) value (along with a drop in mean TAWSS), which resulted in a reduced skewness (approx. 0.2319) for the 24-month post-intervention case. OSI levels also diminished post intervention (32% drop in mean OSI), signifying reduced disturbed flow. Additionally, a relatively smaller proportion of the vascular wall area (80% versus 34%) is subjected to higher OSI (more than 2.5) values post intervention (figure 7*b*). This indicates appreciable improvement in haemodynamics as lower OSI regions are considered less at risk for restenosis [29].

### 3.2. Predicting the vascular deposition of nanoparticles in the superficial femoral artery

Next, the computational tool is used to predict the wall deposition patterns in a simulation of systemically injected NPs in the pre- and post-intervention SFA to explore possible use of nanomedicine in the non-invasive imaging and treatment of PAD patients. It should be emphasized here that, although the same patient-specific data as for ELIMIT have been used, no patients in this clinical trial received an NP injection. In the *in silico* simulations, spherical NPs (*d*_{p} = 100 nm) were introduced directly into the bloodstream from a bolus of NPs located at the inlet of the artery segment. Assuming that the receptors are uniformly distributed throughout the entire vessel wall with a uniform surface density, flow-coupled particle transport simulations were run for 10 cardiac cycles following the methodology outlined in the Material and methods section to quantify particle adsorption to the arterial wall. Figure 8 shows the resulting wall deposition pattern of NPs under such non-specific targeting for both baseline and 24-month post-intervention cases. A distinctively inhomogeneous distribution of NPs is observed in the baseline case (figure 8*a*) where the vast majority of the particles appear to accumulate near the source with NP surface density gradually decreasing in the downstream direction. Because of the restricted blood flow near the constriction and the resulting increased proportion of backflow during the cardiac cycle, the NPs struggle to move past the stenosis (electronic supplementary material, figure S4). The NPs therefore largely remain confined upstream of the stenosis throughout the simulation. Conversely, a much more homogeneous distribution of particles is observed in the 24-month follow-up case. This can be mostly attributed to the significantly greater particle availability downstream compared with the baseline case (as evidenced by the number of NPs leaving the SFA outlet, which increases by a factor of approx. 260). From electronic supplementary material, figure S5, it takes the NPs roughly three cardiac cycles to fill up the entire SFA segment, while failing to do the same (i.e. saturating the blood with particles) in the baseline case, even after 10 cardiac cycles.

In order to gain further insights into disease management, we explored how surface-functionalized NPs would fare when specifically targeting inflammation, a precursor to many diseases including atherosclerosis. To that end, the vessel wall deposition of cell adhesion molecule (CAM)-directed particles was investigated. These CAMs, in particular VCAM-1, are regarded as early markers of atherosclerosis due to their upregulation at the site of inflammation. First, a phenomenological model (see appendix A) correlating CAM expression to local WSS as detailed in [16,44,51] was used to estimate the surface density of the target receptor with respect to their unstimulated expression under static conditions (see electronic supplementary material, figure S6). Spherical NPs (*d*_{p} = 100 nm) uniformly coated with anti-VCAM-1 antibodies (aVCAM-1 NPs) were then introduced into the bloodstream by placing a bolus of NPs at the SFA inlet as before. Carried by the blood flow, these NPs can recognize VCAM-1 expression and firmly adhere to the vessel wall through ligand–receptor bond formation. Once again, flow-coupled particle transport simulations were carried out for 10 cardiac cycles. The resulting spatial distribution of NPs in terms of particle surface density (cm^{−2}) is presented in figure 9*a*,*b* for the baseline and 24-month post-intervention cases, respectively. In addition to local haemodynamics, receptor surface density appears to modulate particle adhesion, resulting in a spatially inhomogeneous particle distribution pattern in both cases. That is, the higher the receptor density and higher the OSI, the greater the particle concentration at a given site. From figure 9, aVCAM-1 NPs appear to preferentially accumulate near the diseased region in the baseline case, which is encouraging from a therapeutic (and also diagnostic) perspective as this may potentially translate to a higher drug/imaging agent concentration in the targeted area (i.e. the stenosis). Conversely, the same NPs appear to adhere less in the 24-month post-intervention case, resulting in a 47% decrease in mean surface density.

In a similar vein, ICAM-1- and E-selectin-targeted particle distributions were investigated (electronic supplementary material, figures S7 and S8). The mean surface densities for the aICAM-1 and aEsel NPs in the baseline case were estimated to be 1.544 × 10^{−6} cm^{−2} and 2.85 × 10^{−6} cm^{−2}, respectively. aICAM-1 NPs accumulate least efficiently near the constriction when compared with the aVCAM-1 and aEsel NPs. The latter seems to achieve a more uniform distribution than the aVCAM-1- and aICAM-1-coated particles. Interestingly, the E-selectin-directed particles appear to have better targeting in the baseline case, that is, they have the highest concentration in the most severely diseased (narrowest) part of the SFA, thereby making them viable candidates for local drug delivery systems that use NPs as drug carriers to diffuse atherosclerosis [37]. Much like the aVCAM-1 NPs, both aICAM-1 and aEsel particles adhere less in the 24-month post-intervention geometry, resulting in a 52% and 50% reduction in mean surface density, respectively. aVCAM-1, aICAM-1 and aEsel NP distributions are presented in the electronic supplementary material, figure S9.

In order to investigate possible spatial correlation between local haemodynamics and NP adhesion, the OSI and aVCAM-1 NP distribution colour maps were merged for the baseline and 24-month post-intervention cases (figure 10). Qualitatively, there appears to be a strong spatial correlation between the two quantities throughout the computational domain, except for small areas very close to the inlet and outlet boundaries. Generally, regions of high OSI are associated with regions of greater aVCAM-1 deposition. In a similar qualitative comparison made between the TAWSS and aVCAM-1 (electronic supplementary material, figure S10), the spatial correlations are not as evident, especially in the baseline case. In regions downstream of the stenosis, the aVCAM-1 NP deposition seems to follow the pattern of high OSI values rather than that of low TAWSS values. RRT values also show strong spatial correlation with aVCAM-1 surface density (electronic supplementary material, figure S11). RRT, which combines both the OSI and WSS, appears to perform better than TAWSS in predicting NP deposition, especially in the baseline case.

There are a few limitations and possible extensions to this work. Although any error is undesirable, in all practical sense they are unavoidable when simulating patient-specific haemodynamics. This is particularly true in the absence of a measured three-component velocity profile prescribed at the inlet [52–55], as realistic three-dimensional inflow velocity boundary conditions can potentially amplify three-dimensional bulk flow structures [53] and promote particle transport [56]. In these limiting cases, different strategies of applying MRI measured flow data as boundary conditions can be implemented. At their simplest, inflow boundary conditions can be imposed by choosing *a priori* velocity profiles under a fully developed flow assumption: (a) Womersley's profile derived from analytical solutions [57], (b) parabolic profile, or (c) blunt profile [58]. In such a case, the accuracy of the numerical solution is impacted by the choice of velocity profile. While mean WSS distributions are expected to be qualitatively similar, this approach in general has been shown to influence oscillating WSS [54]. Alternative strategies include the augmented formulation of the flow problem using Lagrangian multipliers [59,60] and other subsequent improvements [61,62]. However, some of these approaches can be computationally expensive and difficult to discretize, potentially posing additional challenges.

For simplicity, a fully developed inflow profile was assumed in this work. This may not be the case as the truncated segment upstream may be not straight enough, containing undeveloped velocity profiles. However, if placed sufficiently far upstream of the ROI, as it was in this study, the influence (or ‘memory’) of the inlet velocity profile should be essentially lost [63]. The specific choice of a paraboloid-shaped velocity profile, a reasonable assumption that has been widely used in previous computational fluid dynamics studies [40,64,65], was made considering the following observations. Given the alternative Womersley profile's complexity of implementation and minimal, if any, superiority in haemodynamic accuracy over other inlet velocity profiles, its use might not be necessary in many cases [66–68]. Campbell *et al*. [58] have recently shown that the difference in errors produced by the parabolic velocity profile and the Womersley profile was not statistically significant most of the time. In cases where an assumption must be made, a pulsatile parabolic profile was recommended over the blunt or Womersley profiles as it produced the lowest magnitude of error in the WSS and OSI. The Wormersley number in those cases was in a similar range as in this work (approx. 4 in the baseline case for example). It was further contended that both these haemodynamic metrics are affected more by variation in three-dimensional vessel geometry (geometric factors such as curvature of the vessel) among patients [66–68], which is the main objective of this work and a planned future study. In addition, using the subject's own flow waveform over a non-specific (generalized) waveform leads to more accurate results [58]. Hence, prioritizing an accurate reconstruction of vessel anatomy and acquisition of the subject's own flow waveform over selecting the perfect velocity profile appears to be a reasonable approach.

Another limitation of this is the rigid wall approximation. This simplifying assumption can lead to an overestimation of local WSS [69], especially when there is significant vessel curvature or a bifurcation. However, as demonstrated by Kim *et al*. [70] in a previous study on a compliant SFA with bifurcation, bulk flow quantities are generally insensitive to wall motion. While minor changes (less than 25%) in the OSI and TAWSS were reported in select areas like bifurcations, it was concluded that wall compliance did not significantly influence these quantities [70]. In light of the above findings, the rigid wall assumption was deemed reasonable as a first approximation for the single-branch SFA geometry under consideration [63].

Due to lack of suitable data, the stent struts were not explicitly modelled in our simulation of the 24-month post-intervention case, which can alter the local haemodynamics as previously observed by Chiastra *et al*. [40]. This limitation of our approach might have led to an underestimation of per cent area exposed to lower WSS (i.e. atheroprone regions), thus potentially affecting our prediction of NP deposition pattern in the post-intervention case. Additionally, the precise identification of anatomical markers, and consequently the determination of overlapping SFA regions across the two MRI scan time points (baseline and 24 months post intervention) were subject to usual human error [15].

Planned future work includes the application of the presented analysis framework to investigate a larger population of PAD patients. Longitudinal studies averaging multiple patient-specific models over time can better elucidate relationships between haemodynamics and atherosclerosis progression (and regression). Additional topics that might be of interest include geometric characterization of pre- and post-intervention geometries (vascular segment curvature, torsion and so forth) in multiple patients [71] to investigate, and potentially establish, relationships between various geometric factors and disturbed flow to gain further insight in PAD management. Extending the current model to include vessel compliance along with more realistic outflow boundary conditions is also planned [72,73]. Furthermore, as vascular expression of CAMs is known to be sensitive to the temporal gradient of WSS [74], which can potentially influence CAM-targeted NP deposition patterns, it would be interesting to carry out a systematic analysis of this important quantity at critical geometric locations (e.g. recirculation zones distal to the stenosis) [50]. Additionally, extending the current analysis to the intermediate time points (six and 12 months) to investigate if any of the detected haemodynamic improvements occurred earlier in the post-operative process could provide important insights into PAD management.

## 4. Conclusion

An isogeometric analysis framework has been used to predict the blood flow and vascular deposition in a simulation of systemically injected nanomedicines into the SFA of a PAD patient. Upon intervention, a dramatic improvement in blood flow dynamics has been documented by the *in silico* simulations showing a 32% reduction in OSI and an approximate 500% increase in overall peak blood flow rate. Also, it has been observed that systemically injected nanomedicines, targeted to inflammatory vascular molecules such as VCAM-1, E-selectin and ICAM-1, would preferentially accumulate near the stenosis in the baseline configuration. In particular, for the VCAM-1-targeted 100 nm NPs, the vascular accumulation has shown a maximum at the stenosis with a surface density up to 1.9 times higher than the mean. This high, specific vascular accumulation could be effectively used for the detection and treatment of plaques forming in the SFA. The development of accurate image-guided patient-specific predictive tools, such as the isogeometric analysis framework presented here, can help in assessing disease burden and monitor progression/regression of atherosclerosis in response to therapeutic and interventional treatment.

## Funding statement

This work was supported by the Cancer Prevention Research Institute of Texas through grant no. CPRIT RP110262, and through grants from the National Institutes of Health (USA) (NIH) U54CA143837 and U54CA151668. The additional support for S.S.H. and T.J.R.H. from the Portuguese CoLab grant no. UTA06-894 is also gratefully acknowledged.

## Acknowledgements

The authors thank Travis Sanders for his assistance with data post-processing and Benjamin Urick for generating the artwork in figure 2.

## Appendix A

**A.1. The flow problem**

The strong form of the continuity and momentum balance equations in the fluid domain *Ω* with its boundary *Γ* divided into three non-overlapping parts, the inflow (*Γ*_{in}) and outflow (*Γ*_{out}) boundaries and the vascular wall (*Γ*_{s}), are as follows:
A 1
A 2
A 3
A 4
A 5
A 6where **x** is a point in the spatial domain ** Ω** and

*t*is a point in the time domain [0,

*T*]. In these equations,

**u**(

**x**;

*t*) represents the fluid velocity vector,

*p*(

**x**;

*t*) the pressure,

**f**the external body force,

**n**the unit outward normal to the surface, denoted

**, and**

*Γ***u**

_{0}(

**x**; 0) the initial velocity vector. An inflow velocity vector

**g**(

**x**;

*t*) was specified at the inlet (

*Γ*

_{in}), a no-slip boundary condition (equation (A 4)) was prescribed at the rigid and impermeable vascular wall (

*Γ*

_{s}) and a traction-free outflow boundary condition (equation (A 5)) was implemented at the branch outlet (

*Γ*

_{out}).

**A.2. The mass-transport problem**

In the strong form, the transport problem can be stated as
A 7
A 8
A 9
A 10
A 11where *C*(x; *t*) is the particle concentration, **K** is the diffusivity tensor and *σ* is the reaction coefficient. At the inlet (particle injection site) *Γ*_{in}, a Dirichlet boundary condition is prescribed (equation (A 8)) where *C*^{0} is the particle concentration at the inlet. At the outflow *Γ*_{out}, a homogeneous Neumann boundary condition was specified (equation (A 9)); and a Robin-type boundary condition (equation (A 10)) was prescribed at the rigid wall interface (*Γ*_{s}), where the normal component of the flux of particles ‘diffusing’ out of the fluid domain *Ω* and adhering to the vessel wall *Γ*_{s} is assumed to be directly influenced by the local wall shear rate *S*, the particle diameter *d*_{p}, and the probability of adhesion *P*_{a}, as follows (figure 2):
A 12Here, d*φ*/d*t* is the rate of NP accumulation within the vessel wall with *φ* denoting the NP surface density.

The parameter *P*_{a} for a spherical particle in point contact with the vessel wall can be expressed as
A 13Here, *m*_{l} is the uniform surface density of ligand molecules decorating the NP surface that can specifically interact with counter molecules (receptors) expressed on the vessel wall with a surface density *m*_{r} (figure 11). is an affinity constant characterizing the molecular interaction between ligands and receptors at zero mechanical load, *λ* is a characteristic length of the ligand–receptor bond typically of the order of 0.1 to 1 nm; *k*_{B}*T* is the Boltzmann thermal energy ( = 4.142 × 10^{−21} J); *F*^{s} ( = 1.668) is the coefficient of hydrodynamic drag force on the spherical particle; *μS* is the WSS; and *r*_{0} is the radius of adhesion defined as
A 14with Δ being the separation distance between the particle and the substrate, at equilibrium (i.e. firm particle adhesion). The selected parameters are listed in table 1. For a more comprehensive discussion, readers are referred to [16].

## Footnotes

↵† Present address: Department of Molecular Cardiology, Texas Heart Institute, Houston, TX, USA.

- Received January 1, 2015.
- Accepted March 24, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.