Peripheral bypasses are amongst the most common vascular interventions; however, the reality is that millions of these bypasses fail due to vascular remodeling and this is a real burden for National Health Systems. In the UK alone, costs for surgery exceed 200 M (Peach et al., 2012). Why bypasses fail is a critical issue in vascular surgery today, however, traditional approaches have not provided answers to this problem. It is essential to mention that as of today, animal experiments to study peripheral grafts have failed dramatically. To put it simply, there are no animal models which would provide useful data to understand lower extremity venous bypass failure in humans and other, novel approaches are urgently required. The ideal vascular conduit remains the “Holy Grail” of vascular surgery (Byrom et al., 2010). Variability and uncertainty of the outcome are serious issues and the consequences of a failed graft for every patient, just dire. When a bypass graft blocks, blood supply is usually worse than before bypass surgery. In these circumstances amputation can be inevitable unless the graft can be salvaged and the blood supply restored. Recent randomized controlled trials showed that 40% of lower extremity vein grafts occlude or develop significant stenosis within the first year after implantation (Owens et al., 2015). Results for more complex procedures to the calf vessels have usually slightly worse prognosis, with resultant serious morbidity and mortality. These figures have largely remained unchanged for the past several decades. One can read (Owens et al., 2015), “…it is discouraging to consider that 5 decades of high-powered science has not effectively changed bypass graft outcomes.” Improving results of long-term limb salvage remain critically dependent on understanding the mechanisms of successful vein graft adaptation to the arterial environment. In a clinical landscape with ever-increasing and more aggressive bypass procedures, the use of novel mathematical modeling and simulation tools to understand venous adaptation to the arterial environment would help preventing the significant numbers of excess complications, mortality and cost of re-interventions and alternative therapies. This is key in order to devise a personalized “best management plan” and to deliver the best treatment for a specific patient, at that crucial time.

This paper is focused on the development of a simulation framework to elaborate the mathematical tools needed to understand and model key properties of a multi-scale, clinical problem, i.e., bypass failure, at different levels. The vision is that we can use multi-scale mechanistic models to understand and disentangle the complexity of lower limb graft failure and to start making some headway toward patient stratification in this regard.

Neointimal hyperplasia (NIH) is one of the processes leading to restenosis (and ultimately graft failure). NIH is, in simple terms, the re-narrowing of a blood vessel after a stenting or vein grafting surgical procedure, due to tissue growing at the site of injury (Murphy and Boyle, 2010). Upon the start of the formation of NIH, smooth muscle cells change their phenotype from contractile to synthetic, making them more proliferative and resulting in thickening of the arterial tissue. The new conditions in the vascular tissue trigger a reaction which also leads to the release of growth factors and cytokines, including transforming growth factor beta (TGF-β) (Guerri-Guttenberg et al., 2013), platelet derived growth factor (PDGF) (Huang et al., 2002), fibroblast growth factor (FGF-2) (Nabel et al., 1993) and a group of inflammatory cytokines (Collins et al., 2012). Furthermore, other hypotheses have been studied such as that the fibroblasts in the adventitia might move to the media in the form of smooth muscle cells or that bone marrow-derived progenitor cells coming from the bloodstream might also turn into smooth muscle cells to form part of the vascular tissue (Collins et al., 2012). In addition to smooth muscle cells, other types of cells, such as monocytes, are also involved in intimal volume growth (Stark et al., 1997), in response to the inflammatory process triggered by the surgical cut.

However, studies have shown that biological mechanisms alone are not the sole process leading to changes in the morphology and geometry of the artery. A number of mechanical forces in the arterial wall, and shear stress also play a role (Owens, 2010). As opposed to artificial stents, in the case of vein grafts it is the whole conduit, and not just the arterial tissue around it, to be made of living tissue and thus to respond to the surrounding environment. As a result, it is particularly important to consider both the mechanical and the biological mechanisms behind remodeling. Given the complex interplay between biological mechanisms and mechanical stimuli, mathematical modeling can offer much needed help.

Previous research has shown the potential of computational modeling to describe the connection between hemodynamic factors and NIH. One of the first studies to model the relation between blood flow dynamics and NIH appeared in 2001 (Hill and Spendiff, 2001). In this model, the relationship between flow and tissue growth was described by modeling the permeability of the endothelium as a function of wall shear stress (WSS). Tran-Son-Tay et al. (2007) and Tran-Son-Tay et al. (2008) investigated the problem by setting up an experiment using a rabbit model. Subsequently, they developed a mathematical model and compared it to the results obtained from the animal model, using a differential equation where the rate of change of the intima thickness is proportional to the shear stress. The model by Budu-Grajdeanu et al. (2008) includes cellular and chemical mechanisms and takes the change in lumen radius into account. Other studies (Dexter et al., 2009; Boyle et al., 2010; Hwang et al., 2011, 2013; Garbey and Berceli, 2013) developed agent-based models of the disease. Wu and Cassel (2013) modeled NIH as a diffusion process of smooth muscle cells using a feedback-control system, and simulated a reduction in NIH. Finally, in the model by Fok (2012), intimal thickening was modeled as a free boundary problem, which was an accurate description of cell and chemical dynamics, but did not include the flow characteristics which are an essential aspect in this study. A key feature of the work shown in this paper is that the simulations presented here are compared against patient-specific clinical data, which to the authors' knowledge, has not been done before in the context of NIH and bypass failure. Deidentified patient-specific data for this study was obtained with approval of the institutional human investigation committee (approval AD0009, Veterans Affairs Connecticut Healthcare System, West Haven, CT, USA) as part of an ongoing collaboration between UCL and Yale University.

This paper is organized as follows: the second section describes the methods used to develop the computational framework, including specific details about the different biological scales considered in the model, the mechanisms through which mechanical stimuli combine with biological changes, a description of how blood flow is modeled and how the two aspects of the model are combined in a framework. The third section presents specific simulation results for two vein grafts of a human patient and validation with clinical data. Finally, the fourth and fifth sections show a discussion of these initial simulation results and the conclusions, respectively.

This study is part of a broader research activity aimed at developing multi-scale, patient-specific models of cardiovascular disease (Alimohammadi et al., 2015; Di Tomaso et al., 2015), and in particular, vascular remodeling problems. A similar approach has been successfully implemented for the case of atherosclerotic disease (Díaz-Zuccarini et al., 2014). The patient-specific, multi-scale modeling framework proposed here uses a feedback approach to connect biological and mechanical mechanisms in vascular remodeling and details are provided below. In addition, this study seeks to establish the groundwork toward translational studies to understand the relationship between different multi-scale mechanisms and factors in patient-specific studies of vein graft failure, and to develop predictive tools for NIH-prone areas in vein grafts.

This paper brings together mathematical techniques and biological data to produce a model of NIH that describes the link between hemodynamic forces and significant biological mechanisms, according to the literature to date. The mathematical framework presented here relies on patient-specific imaging and hemodynamic measurements, which enables the possibility to study in detail the characteristics of the blood flow in patient-specific geometries. These data, together with literature findings (Humphrey, 2002; Budu-Grajdeanu et al., 2008; Duru et al., 2015), enable to identify areas most prone to NIH. In a nutshell, the process is as follows: patient-specific blood flow simulations on the vein graft prior to the formation of NIH (right after the procedure) are performed by means of a computational fluid dynamics software (ANSYS CFX). Subsequently, time averaged wall shear stress (TAWSS) results from the simulations are used as an input to a biochemical model of cell growth. This is mainly based on the association between WSS and the production of nitric oxide which has been established in previous research (Andrews et al., 2010), in addition to various other biological mechanisms, which will be explained in this section. A remodeling process occurs and a new geometry is then obtained. Finally, the steps are repeated until the time for comparison with clinical data is reached. Details are provided in the following sections.

From the fluid dynamics simulations, the variable of interest for the model is TAWSS, as this is the mechanical factor that has an influence on the turnover of cells. Results of TAWSS were extracted after each step of the remodeling cycle, in order to be used as input data for the biochemical model. Figure 6A shows TAWSS in the right side graft at the beginning of the cycle, and after 2 and 4 months, while contour maps of the TAWSS values in the left side graft at time zero and at 8 months are shown in Figure 6B. It can be noticed that in both cases WSS distribution changes as the remodeling occurs, and that low TAWSS areas tend to be located in the proximity of bifurcations or turns (as shown by the zoomed in areas in the figure). For instance, in the right side graft case, the areas where TAWSS is lowest are at the proximal and distal anastomosis, with the latter showing a much more extended area where TAWSS is less than 0.5 Pa. Zones with low TAWSS are also present along the graft. However it is interesting to note that due to the remodeling process these areas tend to decrease when not situated near parts of the graft where the flow is likely to be deflected. The same is also valid for the left side graft, in which low TAWSS is located where the diameter of the graft is exposed to sudden changes, or before bifurcations. To illustrate the importance of the remodeling of the graft and how it changes the WSS distribution, a comparison between TAWSS contour plots and hyperplasia growth was made (Figure 7). After the remodeling, the changes in the geometry of the graft have a significant impact on the hemodynamics, which in some cases reduced the areas of low TAWSS, while in others it further reduced TAWSS thus causing the area to become even more prone to growth. By affecting the geometry, the locations of NIH growth have a significant impact on the next remodeling cycle. To illustrate the importance of remodeling in the graft, the areas of critical TAWSS (<0.5 Pa) were measured before and after the full remodeling cycle (Figure 7). In the right side graft, this showed that the remodeling led to a decrease in the areas of critical TAWSS with the total area dropping from 17.79 to 4.75 cm². However, the area close to the distal anastomosis changed more slowly from 7.2 to 3.5 cm², so the extent of low TAWSS remained significant, causing a critical amount of growth in the area near the distal part of the graft. This part of the graft was also where the lowest values of TAWSS were found, with TAWSS as low as 0.2 Pa before remodeling took place. On the other hand, the proximal part of the anastomosis was subject to a change from 1.72 to 0.35 cm² of the low TAWSS area, which, while still causing some degree of growth, did not have such an impact as in the distal segment. A similar behavior was observed in the left graft, with the total area of low TAWSS decreasing from 105.82 to 72.13 cm², the area near the proximal anastomosis from 33 to 21.71 cm², near the distal anastomosis from 77.45 to 44.26 cm². Both grafts present a curvature toward the end, which has an effect on the velocity profile. When going through a turn, the flow at the center for inertia moves slower than the flow near the wall, causing it to move away from the center of the curvature (Giordana et al., 2005). This results in a different flow pattern compared to the case for a straight artery, leaving space for the formation of secondary flows and also slowing down the main flow, both of which are factors causing lower levels of shear stress. For further analysis, the relative residence time (RRT) index, which has previously been used for other cardiovascular diseases to identify areas where TAWSS was low and oscillatory shear index was high, for instance in the study of atherosclerosis (VanderLaan, 2004), was also considered. As shown in Figure 8, this index also indicates the locations of the most severe hyperplasia progression, which suggest high oscillatory shear stress might also play an important role in the development of the disease.

As can be seen from the streamline plots in Figure 9, low TAWSS is located mostly in the proximity of areas where the flow is disturbed due to the presence of bifurcations. In the right side graft, at the onset of systole (point 1), flow was steady and some recirculation zones were observed only after the proximal anastomosis, due to the presence of the occluded vessel. At peak systole (point 2), the velocity was uniform although it decreased at the proximal anastomosis and through the slightly larger diameter of the cross section near the distal anastomosis. As the flow slows down and goes into the diastolic phase (points 3, 4, 5), a different distribution can be seen, with a disturbance of the flow at the distal anastomosis, which will likely cause lower values of WSS, especially where changes in diameter occur. This agrees with previous research showing changes in the intensity of vortices after the flow passes from the graft to the artery (Doorly et al., 2002). A similar behavior can also be noticed in the left side graft, where the flow is more streamlined when at peak systole (point 1), except after the narrowing distal to the bifurcation at the lower end, where some reversed flow can be noticed even at systole. When the flow slows down and goes into diastole (points 2, 3), in addition to separation and reversed flow especially before and after bifurcations and where the graft joins the artery, some helicoidal flow can be seen right after the very narrow section at the distal location.

As expected, the disease progression had a different impact at different locations along the graft. Figure 10 shows five different cross sections of the right hand side bypass with their respective lumen occlusion as obtained from the radiological images. The lumen geometry at time zero is shown in pale yellow, whereas the geometry at 4 months post-operation is marked in dark red. As can be seen, the disease progression was different for the cross sections selected, with most hyperplasia forming in proximity of the anastomoses. Figure 11A shows the locations where the most critical values of NIH development were found in the CT scans in the right bypass, while the critical locations on the left bypass are shown in Figure 11B. As restenosis is defined as the development of 50% or more luminal narrowing in the graft (Hill, 2004), this was chosen as the threshold to select the locations of critical NIH development. The measurements show that the locations of growth were reproduced by the model.

As mentioned in the introduction, several computational approaches have been proposed in the literature to model NIH or restenosis. While some studies focused on developing agent-based models (Dexter et al., 2009; Boyle et al., 2010; Hwang et al., 2011, 2013; Garbey and Berceli, 2013), others (Budu-Grajdeanu et al., 2008) used a deterministic approach based on ordinary differential equations. More recently, the increasing use of computational modeling in medical applications and the latest advancements in these techniques have allowed to study the disease from different angles. Some studies have been using CFD as a tool for analyzing the behavior of the flow in pathological cases. In many cases, studies have been conducted to analyze the flow in stents and grafts and to build hypotheses on its effect on NIH and restenosis. This has confirmed the hypothesis of low shear stress and abnormal flow patterns leading to NIH in the pulmonary artery (Berdajs et al., 2015), carotid artery (Harrison et al., 2014), coronary arteries (Guerciotti et al., 2016) peripheral arteries (McGah et al., 2012; Rivera et al., 2014). In addition, some studies have focused on specific applications, such as vascular access for hemodialysis, looking at minimizing NIH by finding the optimum conditions for performing arteriovenous fistulae operations. Research has been carried out in this area on hemodynamic factors such as WSS (Jia et al., 2015) and oscillatory shear index (Ene-Iordache et al., 2015), but also on variables such as blood flow rate and needle tip position (Fulker et al., 2016). CFD has also been used to analyze possible design optimization routes for grafts for both hemodialysis (Canneyt et al., 2013), coronary artery stents looking at the flow within stent struts (Gundert et al., 2012), comparing different stent geometries (Gogas et al., 2014) and the effect of deformable stents (Martin et al., 2014), peripheral bypasses (Grus et al., 2016), and novel grafting systems such as an endograft for aneurysm repair (Aristokleous et al., 2015) and a mechanism of external vein graft support (Meirson et al., 2015). Moreover, new studies looking at the interaction between mechanical forces and cell response for instance used agent based (Ziraldo et al., 2013), hybrid agent based-continuum (Garbey et al., 2015), coupled agent-based and finite element (Zahedmanesh and Lally, 2011), and mechanistic approaches (Goodman et al., 2016).

The model of NIH presented in this paper shows an intra-patient comparison of two different bypass geometries and the applicability of a patient-specific, multi-scale approach to each case. The physiological parameters used by this model have been tested in previous studies to quantify the behavior of SMC and collagen and its influence on NIH progression. The ability of the model to produce results and to formulate and test hypothesis on a macroscopic scale using data from experiments conducted on a cellular level shows the applicability of a quantitative approach to the interpretation of biological data and processes. The model confirmed that the behavior of the biochemical species play an important role in the disease progression. The model validation was carried out against a patient-specific, unique clinical dataset by detecting the locations of restenosis in the CT scans from the patient, and then comparing them against the simulated results. As can be seen in Figure 11, the locations of restenosis in the patient corresponded to those estimated by the model.

In the right side graft, growth estimated for each point ranges from 0.85 mm (region 3) to 1.2 mm (region 1) whilst in the CT scans, an approximate growth of 1.5 mm (region 3) to 2.6 mm (region 1) is observed. In the left side graft, values of growth at the critical locations measured 0.3–0.4 mm in the simulated geometry, while in the real case were estimated to be around 0.8 mm. It is however difficult to correctly measure NIH development due to the quality of the CT scans and uncertainty in the measurements. Additionally, some sections of the artery occluded in the real case and although the location was relatively well predicted, growth was underestimated by the model. There is a complex interplay between the remodeling/growth in the model and the hemodynamic parameters and this warrants further investigation.

This shows that combining a model describing biochemical interactions and a mechanical model describing the hemodynamics is a promising approach to describe the behavior of the patient's disease. This type of model has a potential use as a clinical tool to inform clinicians on the progression of the disease, and future modifications will allow not only to locate NIH but also to inform on the severity of the growth. The applicability of the model to patient-specific geometries was also evaluated on two different types of vein graft, for both of which the model located areas with the most NIH development. In addition, the study also lead to a further analysis of the importance of including collateral vessels in the model, which were found to have an influence on the fluid dynamics results and consequently on the biochemical processes as well.

This work has shown the potential of using experimental data to develop integrative in silico models to study NIH progression. A combination of data coming from in vivo and in vitro experiments was particularly useful to get a better understanding of the disease and design a computational model able to produce quantitative information related to relevant dynamical processes of NIH progression. In particular, from the fluid dynamics aspect it was important to be able to apply an appropriate type of boundary conditions. Although for the right side bypass this was relatively straightforward due to the presence of only one inlet and one outlet and the validity of the stress-free boundary condition at the outlet, the conditions were more complicated on the left graft. As stated above, boundary conditions for this model were provided by a lumped parameter model. This allowed to account for the effects of resistance of the peripheral vessels, through the resistance parameter, as well as vessel elasticity through the compliance parameter (Shi et al., 2011). Applying a zero-pressure boundary condition in this case would have led to a non-physiological flow split, with most of the flow being directed toward the path of least resistance. In addition, accounting for these aspects of the flow allows to obtain physiological pressure curves (Kim et al., 2010).

Finally in order to study the performance of the whole multi-scale modeling approach proposed in this paper, the sensitivity of NIH development with respect to key parameters in the biochemical model was evaluated. In this initial analysis the kinetic rates related to proliferation and apoptosis of smooth muscle cells were selected (given the key role of these cells in NIH progression).

As seen in Figure 12, NIH growth was sensitive to changes in both the apoptosis and proliferation coefficients. In particular the growth rate was more sensitive to the increase in the apoptosis parameter, with values increasing by up to three folds when the parameter was increased by 20%.

These results show how important is the influence of the biochemical parameters on the final clinical endpoint simulated by this model. So far, with these results it is not possible to conclude about how different smooth muscle cell behaviors could have a different impact on disease progression. However, if a more detailed biochemical model is proposed, the same kind of sensitivity analysis might enable us to draw conclusions about the influence of specific biological/physiological mechanisms (e.g., intra/extracellular mechanisms or processes related to inflammation) that play a key role in NIH.

The multi-scale modeling approach proposed here, has shown how using a mathematical description capturing different biological and physiological scales related to NIH can be used to have a better holistic understanding of the disease allowing the use of advanced computational tools to generate dynamic results related to disease progression. This kind of approach is not only useful for simulation (through validation with clinical data) but also for further analyses (e.g., sensitivity analyses) that can be linked to testing specific biological hypotheses.

The study presented a multi-scale model of the behavior of peripheral vein grafts affected by NIH when subject to WSS, the first to date to simulate the disease progression based on the dynamics of the flow in patient-specific geometries. A comparison intra-patient between a femoro-popliteal and femoro-distal graft was carried out which showed the potential but also, limitations of the applicability of the model to different cases. By means of this model it has been possible to analyze NIH progression and its mechanisms from a hemodynamics and a biochemical standpoint. In the future, the model will be used for the analysis of more patients in order to test its applicability on different cases of vein graft NIH. Moreover, the remodeling framework produced could be expanded in future work, for instance by extending the model to other components of vascular tissue. Mathematical models already exist that describe the cellular behavior in atherosclerosis (Cilla et al., 2013; Díaz-Zuccarini et al., 2014; Pichardo-Almarza et al., 2014). These could be modified to describe behavior of additional cells such as macrophages in NIH, in order to include a better description of the inflammation processes triggered in the disease. In addition, the endothelial layer could be added as already seen in previous literature (Pichardo-Almarza et al., 2014; Alimohammadi et al., 2016) to describe the transport of biochemical species into the vascular tissue. Modeling the endothelium would also enable to take into account the effect of WSS on the behavior of monocytes, as previously done for atherosclerosis (Díaz-Zuccarini et al., 2014), thus considering a further aspect of the mechanical response of cells. Further analysis could be done by considering other shear stress indices, as some of these have been shown to give better predictions in models of other vascular diseases such as atherosclerosis (Alimohammadi et al., 2016), which is also a disease affected by mechanical forces similarly to NIH, although through different cell species and processes. More flow characteristics could be included to describe recirculation and separation zones, all of which have been shown to affect the development of NIH (Murphy and Boyle, 2010). Finally, future work should investigate the use of fluid structure interaction in order to capture the movement of the arterial wall when subject to hemodynamic forces, which might have an influence in the results.

This study was carried out in accordance with the recommendations of the Institutional Human Investigation Committee (Veterans Affairs Connecticut Healthcare System, West Haven, CT, USA) with written informed consent from all subjects. All subjects gave written informed consent in accordance with the Declaration of Helsinki. The protocol was approved by the Institutional Human Investigation Committee (Veterans Affairs Connecticut Healthcare System, West Haven, CT, USA).

FD, CP, and VD conceived the study and carried out simulations and wrote the manuscript. MB created the geometries from CT scans, AD acquired the clinical data and provided medical input. SH conceived the study and provided medical input.

This project is supported by the Engineering and Physical Sciences Research Council (EPSRC) through the Doctoral Training Programme “Healthcare Engineering for an aging population.” The authors are grateful for the support of the EPSRC Network “POEMS” (EP/L001101/1) for facilitating the visit of AD to UCL. The authors gratefully acknowledge support by the EPSRC grant “Personalized Medicine Through Learning in the Model Space” (grant number EP/L000296/1) and the Leverhulme Trust Senior Research Fellowship “Exploring the Unknowable Using Simulation: Structural Uncertainty in Multiscale Models” (Fellowship number RF-2015-482). This work was supported by the resources and the use of facilities at the Veterans Affairs Connecticut Healthcare System (West Haven, CT).