## Abstract

Vein graft failure is a prevalent problem in vascular surgeries, including bypass grafting and arteriovenous fistula procedures in which veins are subjected to severe changes in pressure and flow. Animal and clinical studies provide significant insight, but understanding the complex underlying coupled mechanisms can be advanced using computational models. Towards this end, we propose a new model of venous growth and remodelling (G&R) based on a constrained mixture theory. First, we identify constitutive relations and parameters that enable venous adaptations to moderate perturbations in haemodynamics. We then fix these relations and parameters, and subject the vein to a range of combined loads (pressure and flow), from moderate to severe, and identify plausible mechanisms of adaptation versus maladaptation. We also explore the beneficial effects of gradual increases in load on adaptation. A gradual change in flow over 3 days plus an initial step change in pressure results in fewer maladaptations compared with step changes in both flow and pressure, or even a gradual change in pressure and flow over 3 days. A gradual change in flow and pressure over 8 days also enabled a successful venous adaptation for loads as severe as the arterial loads. Optimization is used to accelerate parameter estimation and the proposed framework is general enough to provide a good starting point for parameter estimations in G&R simulations.

## 1. Introduction

Vein graft maladaptation and failure are ubiquitous in vascular surgeries, including coronary artery bypass graft (CABG) surgery, peripheral vascular surgery and arteriovenous fistula (AVF) procedures for haemodialysis access. The long-term failure rate is as high as 50% at 10 years for CABG grafts [1,2], 25–55% at 5 years for infragenicular bypass grafts [3,4] and 40% for the primary failure of AVF [5]. Although the primary mode of long-term failure differs across each of these pathologies, for example atherosclerosis for CABG venous grafts [6] and intimal hyperplasia for AVF [7], the dramatic change in mechanical environment remains a common driver of failure. Pressure- and flow-induced stresses provide critical stimuli to soft tissue remodelling [8] and influence the long-term outcomes of these grafts. There is, therefore, a pressing need to understand how alterations in mechanical environment affect venous (mal)adaptation.

Here, we subject an idealized vein to altered haemodynamics consistent with implantation into an arterial environment and study its response. We adapt the constrained mixture theory of growth and remodelling (G&R), which has been successfully used to model arterial adaptations in health and disease [9–11], to probe mechanisms of venous adaptation in a perturbed haemodynamic environment. This theory incorporates cellular level information phenomenologically within fundamental equations of nonlinear continuum mechanics and can account for different mechanical properties and rates and extents of turnover of individual wall constituents (e.g. elastin, collagen and smooth muscle). Initial efforts in extending the constrained mixture theory to veins aimed to understand adaptation (one to 12 months) by a coronary bypass graft but examined only the altered pressure state [12]. We now extend our previous venous G&R theory to model the adaptive response to altered pressure and flow while accounting for both acute changes in vasoactivity and chronic changes in matrix turnover. As in our prior study, we use optimization tools to accelerate parameter estimation.

Arterial models of G&R have demonstrated that the combined effects of mass turnover, vasoactivity and constituent prestretch affect vessel adaptation in the presence of altered haemodynamics [13]. In this work, we use optimization first to determine constituent prestretch and other G&R parameters that achieve homeostasis at moderate loads. We then fix these values and subject the numerical vein to a series of combined loads, increasing from moderate to severe changes in flow and pressure. Results uncover plausible mechanisms of adaptation versus maladaptation. We then explore potential ways to mitigate maladaptation by subjecting the vein to gradual changes in flow and pressure and quantifying the resulting adaptation.

## 2. Material and methods

The elastodynamics of the vascular wall can be described well quasi-statically; hence one must satisfy the equilibrium equation
2.1where ** σ** is the overall Cauchy stress, which has contributions from all structurally significant passive (e.g. smooth muscle, collagen, elastin) and active (contraction of the smooth muscle) components of the vessel wall. Within the continuum framework, Cauchy stress can be written as
2.2where

**F**is the deformation gradient tensor associated with the passive motions, with ,

*W*is strain energy density function, and

*σ*

^{act}

_{θ}is the active stress in the circumferential direction.

### 2.1. Growth and remodelling theory

Details of the constrained mixture theory of G&R and its application to a range of problems in the arterial circulation can be found elsewhere [9–11,14]. Here we discuss relevant aspects of the model for a thin-walled cylinder and its extension to a vein.

Each constituent within the mixture is endowed with its own material properties, rates of production and degradation, and natural configuration, but is constrained to deform with the bulk wall. While it is difficult to distinguish the effects of individual constituents from the others, for mathematical convenience, and prompted by prior success with arterial models, we model the elastin-dominated contributions using a neo-Hookean model and collagen-dominated contributions along the axial, circumferential and diagonal directions using Fung exponential models (figure 1). Because one cannot delineate the passive contributions of the smooth muscle to the bulk behaviour, it is modelled along with the collagen. Modelling these constituents as hyperelastic, according to a simple rule-of-mixtures, the total energy stored in the wall is given conceptually by *W* = *Σ**ϕ*^{k}*W*^{k}, where *ϕ*^{k} is the mass fraction of constituent *k* and *W*^{k}, its stored energy density. Specifically, for elastin,
2.3where *c*^{e} is a material constant, and for smooth muscle and collagen,
2.4where *c*^{k}_{i} are material constants and the superscript *k* represents the collagen fibre family or smooth muscle. *λ*_{θ}, *λ*_{z} and *λ*_{r} denote stretch along the *θ*, *z* and *r* directions, respectively. *λ*^{k}_{n(τ)}(*s*) = *G*^{k}_{h}*λ*(*s*)/*λ*(*τ*) is the stretch of constituent *k* that was deposited within the mixture with prestretch *G*^{k}_{h} at instant *τ* [9,14], where *n*(*τ*) denotes constituent-specific natural configuration, and ∀*τ* ∈ [0, *s*]. These constitutive models have successfully captured the passive anisotropic, nonlinear biaxial behaviour of a vein [12,15,16].

To capture turnover of the constituents, mass density and strain energy density are allowed to evolve according to a ‘fading memory’ behaviour [14], namely
2.5and
2.6where *M*^{k}(*s*) and *W*^{k}(*s*) are mass per unit area and energy stored in constituent *k* per unit volume, respectively, and *ρ*(*s*) is the current overall (constant) mass density. Here, *Q*^{k} ∈ [0, 1] is the fraction of constituent *k* produced at or before *τ* = 0, while *q*^{k} ∈ [0, 1] is the fraction of constituent produced at time *τ*, both surviving until current time *s*. *m*^{k} > 0 is the true mass density production and is the constituent-specific stored energy per area, which depends on **C**^{k}_{n(τ)} the right Cauchy–Green tensor for constituent *k* defined with respect to its individual natural configuration *n*(*τ*) [9].

#### 2.1.1. Active smooth muscle.

The active stress, *σ*^{act}_{θ}, generated by the circumferentially oriented smooth muscle is modelled as [17]
2.7where (Nm^{−2}) is a (calcium-dependent) scaling parameter, *ϕ*^{m}(*s*) is the mass fraction of active smooth muscle, *λ*_{M} and *λ*_{0} the circumferential stretches at which the active stress is maximum and zero, respectively, and *C*(*s*) is the shear modulated net ratio of constrictors to dilators: *C*(*s*) = *C*_{B} − *C*_{s}Δ*τ*_{w}, where *C*_{B} is the basal constrictor to dilator ratio and *C*_{s} is a scaling factor for shear stress-induced changes in constrictor/dilator concentration. The vasoactive response of a vessel can evolve over a range of diameters, thus shifting the stretch at which maximum and minimum tone occurs [18]. This shifting behaviour motivates an active circumferential stretch *λ*^{m(act)}_{θ}(*s*) = *a*(*t*)/*a*^{m(act)} that evolves via the first-order equation
2.8where *a*^{m(act)}(0) = *a*(0) and *K*^{act} is a rate parameter [9,10]. We note that there are currently no data on active stress generation in the axial direction, and it is the circumferentially oriented smooth muscle that is most likely to affect radius and thickness at fixed axial stretches.

#### 2.1.2. Mass production and degradation functions.

Prior mass production and degradation functions for veins were based on a model for a healthy basilar artery [9]. While those mass production functions were adequate for (relatively) moderate perturbations in pressure (*γ* = *P*/*P*_{h} = 1.5), they were found to be inadequate (i.e. they did not restore the theoretical homeostatic state within numerical tolerance), for a moderate combined load, defined as a moderate increase in both pressure (*γ* = 1.5) and flow (*ε* = 1.1), with *γ* = *P*/*P*_{h} and *ε* = *Q*/*Q*_{h}. Hence, there was a need to modify the mass production constitutive equation for general vein graft responses. Prior applications of this model to pathological arterial adaptations (e.g. abdominal and cerebral aneurysms [11,14]) required a mass production function that accounted for an increase in cell density in response to severe insults. Because the acute insult from increases in haemodynamic loading experienced by a vein subjected to arterial conditions is more severe, the mass production functions for the vein were similarly modified to handle this increase in cell density,
2.9where *m*^{k}_{0} is a basal production, *K*^{k}_{1} and *K*^{k}_{2} are gain-type parameters and Δ*σ* = (*σ*_{θ} − *σ*^{h}_{θ})/(*σ*^{h}_{θ}) and Δ*τ*_{w} = (*τ*_{w} − *τ*^{h}_{w})/(*τ*^{h}_{w}) are normalized differences between the current values of intramural and wall shear stress and their homeostatic targets, and superscript h denotes homeostatic.

The constituent degradation relation was modelled consistent with first-order type kinetics, of the form
2.10where *K*^{k} is a rate-type parameter that can be a constant or depend explicitly on changes in the tension within the fibres; here we consider the specific form [9] , where
2.11is the level of tension on constituent *k* that was produced at time *τ* and is the difference in fibre tension from its homeostatic value. In homeostasis, the vessel wall density is roughly a constant, thus changes in volume must offset changes in mass. This was best achieved when smooth muscle was allowed to degrade at basal rates (*K*^{m} = const.), while collagen degradation had contributions from both basal and tension-dependent degradation [9,11]. Note that these general functional forms, including mass production and degradation, are consistent with biological experiments and prior models of arterial G&R [9,19].

In summary, our implemented G&R formulation solves for equilibrium at each G&R time *s* by accounting for the evolving nature of the wall constituents, including active and passive responses. For simplicity, we keep axial stretch constant. We numerically solve the traction equilibrium equation for radius, at every time step, using the Newton–Raphson method, then time advance the radius and mass of constituents using an explicit Euler scheme.

### 2.2. Parameter estimation

The G&R formulation has many parameters, making brute force estimation of best-fit values intractable. We systematically determined parameter values in serial steps (figure 2). We broadly classified parameters into following categories: passive, gain, half-life, prestretch and prescribed. Passive material properties were determined by curve fitting published biaxial experimental data for a murine vena cava [12,16]. The choice of constituent directions were consistent with our prior arterial and venous G&R models and the resulting material properties consistent with biaxial measurements on veins from the literature [15,20]. These parameters were then held fixed for subsequent parameter estimations and G&R simulations.

We then determined the remaining G&R parameters through use of the surrogate management framework (SMF), a derivative-free optimization method with proved convergence to a local optimizer based on pattern search theory [21,22] (see appendix A for a brief description of the algorithm). We exploited the vessel's inherent character to maintain a homeostatic state of stress to define an intramural- (*σ*_{θ}) and shear stress (*τ*_{w})-based normalized cost function of the form
2.12We minimize in successive iterations of the SMF method to determine the G&R parameters, including gain, prestretch and half-life parameters, that lead to an optimal adaptation for a prescribed moderate haemodynamic load.

Owing to a lack of experimental data, some G&R parameters were chosen to be consistent with reported experiments or previous models of arterial G&R, or fixed to biologically reasonable values when data were not available, as described below.

## 3. Results

Values of the G&R parameters—prescribed from the literature, determined from mechanical data or optimized—are listed in table 1.

### 3.1. Numerical experiment 1: moderate load

We expect a vein to adapt well to modest-to-moderate perturbations in haemodynamic load. Based on pilot simulations, we chose *γ* = 1.5 and *ε* = 1.1 as representative upper range values for moderate perturbations. By minimizing the cost function (equation (2.12)), we identified G&R parameters that led to an optimal adaptation (table 1) under moderate perturbations. We kept these parameters fixed for the remaining numerical experiments.

Initial optimization results showed that values of radius and thickness at 1500 days were suboptimal using mass constitutive relations prescribed for a healthy artery [12], regardless of the choice of parameters for optimization and their bounds. By contrast, optimization identified a near optimal adaptation (figure 3) with a pre-multiplier form of the mass production equation (equation (2.9)) that was previously prescribed for cases of severe insults. Based on these findings, we concluded that an optimal adaptation, even to a moderate perturbation in haemodynamic load, is dominated by the rate of mass turnover of its structural constituents. Thus, the pre-multiplier form of mass production equation (equation (2.9)) was used in all subsequent simulations.

### 3.2. Numerical experiment 2: combined loads (moderate to severe)

To recapitulate a vascular bypass graft scenario, we subjected the numerical vein to a series of increasing combined loads, from venous to arterial conditions: *γ* ∈ [1, 20] and *ε* ∈ [1, 4]. We sampled 300 points from the *γ* − *ε* space, with approximately 75 points in each quadrant (with origin at *ε* = 2 and *γ* = 10), and simulated adaptations for up to 1500 days at each of the load combinations (figure 4*a*). As it can be seen failed simulations (denoted by ‘X’) clustered in the upper right corner of the *γ* − *ε* space, with failure defined as an inability to satisfy the equilibrium equation at every time step in the simulation.

Veins typically failed during the first few days after haemodynamic perturbation, which motivated us to probe vasoactive changes as they can be immediate. While exploring possible mechanisms of failure, an altered rate of tonic adjustment (*K*^{act}) enabled successful adaptations, that is equilibrium was successfully established at every time step in the simulation (figure 4). All simulations on the *γ* − *ε* grid converged for values of . The cost function distribution in *γ* − *ε* space for is depicted in figure 5; it represents the deviation of stress from a homeostatic value. The increase in cost function is more sensitive to a change in flow than a change in pressure from the basal value. Nevertheless, representative evolution curves (figure 6) for an illustrative severe load of *γ* = 18.96 and *ε* = 3.95 show near optimal adaptation with .

Contrasting with values for cerebral arteries to [9,10], we infer that venous adaptation at severe loads can be impaired by an inability to maintain or rapidly evolve the smooth muscle tone.

In summary, results from numerical experiment 2 led us to conclude that endowing a vein with enhanced mass production and vasomotor capabilities can ameliorate maladaptation, even at severe loads.

### 3.3. Numerical experiment 3: step versus gradual loading

To evaluate the potential benefits of a gradual change in load and to understand the interplay between different factors affecting adaptation, we performed simulations at the same 300 sampled (*γ*, *ε*) combinations for the following conditions:

(a) gradual change in pressure over 3 days, but a step change in flow

(b) step change in pressure, but a gradual change in flow over 3 days

(c) gradual change in pressure and flow over 3 days

(d) gradual change in pressure and flow over 8 days.

Note that the pre-multiplier form of mass production (equation (2.9)) but (more conservative) were used for all simulations. Results for cases (a)–(c) are summarized in figure 7, and contrasted with results from a step change in pressure and flow (as used for figures 3–6). Overall, the gradual load cases have fewer failed points than a combined step change in load (Experiment 2). Unexpectedly, gradual flow accompanied by a step change in pressure (case (b)) had the lowest number of failed points, fewer than the gradual change in flow and pressure (case (c)). All simulations fully converged when the gradual load was applied over 8 days, however (case (d), figure 8). Results from this numerical experiment suggested that a gradual change in load can mitigate maladaptation, even with a non-optimal evolution of active properties.

## 4. Discussion

Blood vessels respond to sustained changes in haemodynamic loads by altering vasoactivity and modulating the turnover of wall constituents to restore the homeostatic mechanical environment, often represented in terms of mechanical stresses [28–31]: intramural circumferential stress *σ*_{θ} = *Pa*/*h* and axial stress *σ*_{z} = *f*/*π**h*(2*a* + *h*) as well as wall shear stress *τ*_{w} = 4*μ**Q*/*π**a*^{3}, where *P*, *f* and *Q* are the applied ‘loads’, *μ* is the viscosity, and *a* the radius and *h* the wall thickness. Although there is motivation to build computational models based on the underlying changes in gene expression, phenomenological G&R models focusing on gene products (e.g. collagen production) have proved very useful in describing and predicting changes in geometry and mechanical properties [9,32]. In this study, we extended a prior G&R model to evaluate possible (mal)adaptation of an otherwise healthy vein subjected to increased values of blood pressure and flow, up to arterial levels. We uncovered two possible mechanisms leading to maladaptation: an inability of the smooth muscle cells to maintain or rapidly evolve the vessel-level active tone and an insufficient rate of production of matrix constituents. Together these findings are consistent with prior experimental [33] and computational [9] studies that show complementary roles of vasoactivity and matrix turnover in arteries. Nevertheless, this is the first study to confirm such roles based on vein-specific parameters. Importantly, the finding that the mass production term must be augmented via the pre-multiplier term in equation (2.9) emphasizes the importance of heightened turnover in the severe case of vein grafts, hence suggesting either a marked proliferation of resident intramural cells or recruitment of progenitor or inflammatory cells to promote matrix synthesis in response to the dramatic increase in mechanical loading. The essential role of vasoactivity emphasizes the importance of possible endothelial—smooth muscle cell interactions.

We also found possible ways to mitigate maladaptation. Increasing the flow gradually over 3, and especially 8, days was found to be advantageous relative to an abrupt ‘step’ increase. Although this result was not surprising, an unexpected finding was that a gradual increase in flow over 3 days is particularly beneficial when combined with a step, rather than gradual, increase in pressure. This result emphasizes that it is often difficult to intuit G&R responses because of the complex, coupled mechano-responses exhibited by blood vessels. That a gradual increase in flow plus a step increase in pressure resulted in fewer maladaptations than the gradual increase in flow and pressure (cf. figure 7) can be explained as follows. Increases in pressure-induced intramural stress stimulate increases in matrix production but are modulated by changes in flow-induced wall shear stress. Namely, increases in intramural stress stimulate increases in local angiotensin II production by smooth muscle cells and fibroblasts and associated transforming growth factor-beta induced matrix synthesis [27]. By contrast, increases in wall shear stress increase endothelial production of nitric oxide (NO), which attenuates matrix synthesis by smooth muscle cells and fibroblasts. Hence, a step increase in pressure would rapidly initiate matrix production, whereas gradually increasing the flow would delay the attenuation of this process. Regardless, more gradual pressure and flow loading over a longer period (e.g. 8 days) clearly presents some benefits.

Traditional clinical definitions of the failure of a vein when placed within an arterial environment, as a graft or fistula, include loss of patency, formation of an aneurysm or rupture. In this work, however, we adopted a much more stringent definition of failure—an inability to mechanoadapt properly. In particular, mechanoadaptation of a blood vessel to an altered pressure (parametrized by *γ*, where *γ* = *P*/*P*_{h}) and altered flow (parameterized by *ε*, where *ε* = *Q*/*Q*_{h}) assumes that the circumferential intramural stress and wall shear stress are maintained at, or restored to, normal values. This process requires that the luminal radius and the wall thickness via G&R, where the subscript h denotes the homeostatic value [27]. Newton's second law of motion requires forces to balance, which for a thin-walled quasi-statically loaded cylindrical vessel requires *Pa*/*h* = *σ*_{θ} and *f*/*π**h*(2*a* + *h*) = *σ*_{z}. The novel focus herein was to determine if G&R processes, which can evolve the vasoactivity as well as vessel microstructure, could mechanoadapt the vessel within the prescribed period; failure to do so was indicated by an inability of the numerical code to converge at a prescribed (*γ*, *ε*) pair. Hence, we did not address classical metrics of failure such as loss of patency, rupture or formation of an aneurysm.

Perhaps surprisingly, we found that the modelled veins have a remarkable ability to mechanoadapt to marked increases in pressure and flow, though not uniformly so, at levels expected upon exposure to an arterial environment. We emphasize, however, that this mechanoadaptive capacity was achieved in the assumed absence of load-induced cellular damage, particularly no damage to endothelial cells or smooth muscle cells. Damage to endothelial cells, resulting in a lost ability to produce the vasodilatory molecule nitric oxide (NO), could both limit the pressure-induced distension and allow underlying smooth muscle cells to produce more matrix more rapidly, which could help restore circumferential stress toward normal with less dramatic wall thickening. Yet, loss of endothelial cells could also promote adverse thrombotic or inflammatory processes that could easily outweigh any potential benefits of lower NO. There is a pressing need to delineate diverse effects of endothelial cell damage. Conversely, damage to smooth muscle cells would be thought to be doubly detrimental, resulting in a decreased ability to synthesize the additional matrix needed to thicken the wall and likely inciting a further inflammatory response. Again, there is a need for more detailed data on the effects of cellular damage to enable more detailed modelling, particularly related to the role of inflammation. Related to this, there is also a need for data that enable multilayered models to be constructed. The present model captures transmurally averaged behaviours, both mechanical and mechanobiological. Although data are accumulating on layer-specific differences in venous histology and mechanical properties [15,20,34] and it is straightforward to construct multilayered G&R models [35,36], we must await layer-specific mechanobiological data. That is, there is a need to delineate possible differences in smooth muscle cell and fibroblast homeostatic targets. There is similarly a pressing need to determine if such targets could evolve, particularly when a vein is subjected to an arterial environment. It is known that venous identity (indicated by Eph-B4) is lost, but arterial identity (indicated by Eph-B2) is not gained in venous grafts [37]. Hence, although vein grafts and fistulas do not arterialize, there is a need to determine if they can reset their homeostatic targets. Consistent with our goal, we assumed that these targets were preserved throughout remodelling.

The choice of parameters and cost function for problems with few experimental data is a non-trivial exercise in parameter estimation, which we addressed through a novel application of optimization methods. We chose stress-based cost functions to maintain consistency with the constitutive relations for constituent mass production (appendix A) and because we wanted the model to retain its ‘emergent’ behaviour, where evolving radius and thickness emerge naturally from the equilibrium solution at each timepoint in the G&R simulation.

Associated parameter values and their bounds, while primarily motivated by previous work (table 2 in [12]), were identified through a series of optimization steps. Importantly, the results agree with the few measured values. Monos *et al.* [38] reported the active stress generated by smooth muscle in a canine femoral vein to be 5.1 ± 1.5 kPa at a stretch of 1.36. The active stress for the vein deduced from table 1, for a murine vena cava, is approximately 3.2 kPa at a smooth muscle prestretch of 1.31. By contrast, the active stress of smooth muscle in the arterial circulation is of the order of 100 kPa, and the range on during optimization was chosen to reflect the disparity in active stress between an artery and a vein.

In the past, the constrained mixture theory of vascular G&R has typically been employed for modest-to-moderate perturbations in haemodynamics, making this study perhaps the first when the theory has been applied in severe haemodynamic loading scenarios. The consistent clinical observation of maladaptation in vein grafts would lead one to expect that simulations should fail at higher loads. Because greater perturbations from homeostatic states also pose challenges to numerical schemes, one must take care to differentiate a numerical failure from an adaptation failure. In the G&R framework, a failed simulation could be due to a failure in the root finding Newton–Raphson scheme, the time marching explicit scheme, or a physical maladaptation. In Newton–Raphson schemes, high gradients can lead to a bad initial guess and hence push the algorithm outside the range of guaranteed convergence. This issue was resolved in our simulations using a bound check condition on the roots at every iteration. In addition, values of the cost function were confirmed with 4th order Runge–Kutta time advancement schemes and tested to ensure time step convergence. Thus, we are confident that none of the failed simulations could be attributed to lack of convergence or failures in numerical algorithms.

Some of our findings have translational potential as a pharmacological therapy or medical device. For example, a pharmacological therapy that upregulates constituent mass production, alters active muscle properties or thresholds for phenotypic modulation could hold promise as an avenue to prevent vein maladaptation. Increased matrix production can be achieved with certain antagomirs, for example [39], and various agonists can increase smooth muscle contractility. The challenge with regard to the latter, of course, would be to increase venous or large artery but not resistance vessel contractility so as to not increase systemic blood pressure. There is a need, therefore, to understand what modulates smooth muscle contraction in veins. Also a step change in load, as occurs during surgical procedures when the arterial clamps are released, can lead to a graft failure at arterial conditions. A gradual change in load, possibly achieved using a polymer wrap or an external biodegradable sheath or stent [40], holds promise in preventing venous graft maladaptation. Towards that end, we explored the effect of the duration of the gradual loading, which we found need not be that long (less than or equal to 8 days). We hope this will inform design considerations and specifications, particularly since long-term foreign body responses will drive local inflammation and limit adaptive capacity. Of course, many other issues such as uncertainties associated with parameter values [41], and sensitivity to loads would need to be addressed prior to translation. The proposed framework is, however, general enough to provide a good starting point for each of these considerations.

The absence of simulated changes in axial stretch and cell damage are two major assumptions in our model. Axial stresses play fundamental roles in compensatory adaptations [42]. A vein undergoes changes in axial stretch when transposed from its native environment to an arterial environment. It is not yet understood how and if this change in stretch affects venous adaptation. Adaptive responses to flow and pressure, for the most part, have been achieved without independent altering of vessel length [18,42]. Also, lack of clinical evidence on maladaptation due to axial loading, such as kinking or twisting of grafts, justified our choice to neglect axial stretch changes. Harvesting, overdistending with saline solution prior to implantation, surgical insult and arterial haemodynamic loads are all potential sources of cell damage during a surgical procedure [43,44]. While there is evidence of a preserved endothelium during the first week of AVFs, [45] damage to venous smooth muscle cells is expected within an arterial environment [46]. The data are currently too limited, however, to inform models of damage or the inflammation that accompanies such damage. Changes in surgical techniques, such as atraumatic, ‘no touch’ technique, have shown superior long-term patency by reducing vascular damage [44]. The extension of G&R model to compare different approaches might lead to interesting insights into venous adaptation. However, this type of comparison will need to await additional data, not the least of which is characterization of the baseline (damaged) properties following different methods of preparation. We also neglect complexities associated with pulsatile blood flow. While the theory has been coupled with fluid-structure solvers in the past [47], and an associated G&R theory is available [48], lack of quantitative data on venous cell responses to oscillatory loads makes extensions to pulsatile flow premature and likely non-trivial.

Computational models can not only capture experimental observations and inform the design of experiments, they also hold potential for use in clinical decision making. Computational tools to stratify patients at risk for graft failure would require incorporation of the proposed G&R models into patient-specific simulations. Recently, patient-specific CABG simulations quantified mechanical stimuli in coronary grafts with high spatial and temporal resolution, and revealed significant differences in mechanical stimuli between venous and arterial grafts [49]. Patient-specific simulations are required to capture the complex interactions among local geometry, wall constituents and global haemodynamics [50]. At the same time G&R models, such as the continuum-based constrained mixture models used here, are well suited for simulating vein graft (mal)adaptations, by providing a means to incorporate cellular level mechanisms into a continuum framework. These models have been coupled with idealized patient-specific fluid–structure interaction simulations in the past [47], constituting an essential step towards a coupled framework. Of course, longitudinal clinical data, focused experiments and validation are necessary to create coupled predictive models of the underlying mechanisms of adaptation. *In vivo* models of vein graft failure [51,52] and *ex vivo* set-ups [53], coupled with histochemical and biaxial studies, are natural pathways for model validation. Future efforts should be directed towards coupling patient-specific simulations with G&R models of vascular adaptation, establishing *in vivo* models of vein graft failure, exploring ways to induce a gradual change in flow and pressure in vivo, and elucidating underlying mechanisms of venous maladaptation through focused experiments.

## Data accessibility

Data supporting this publication are available at Dryad (http://dx.doi.org/10.5061/dryad.33th6) [54]. Additional data can also be obtained upon request by emailing: amarsden{at}stanford.edu.

## Author's contributions

A.B.R, J.D.H. and A.L.M. were involved in design of the study and interpretation of results. A.B.R. was involved in computational implementation and simulations.

## Competing interests

We declare we have no competing interests.

## Funding

This work was supported by NIH grant nos. (R01 HL123689 to A.L.M., R01 HL128602 to J.D.H.), a Burroughs Wellcome Fund Career Award at the Scientific Interface and an NSF CAREER award (OCI 1150184) to A.L.M.

## Appendix A. Optimization

**A.1. The surrogate management framework**

The SMF is a derivative-free pattern search method with a well-established convergence theory, employing a surrogate function to accelerate convergence. Our implementation of the algorithm starts by running function evaluations (in this case, each a full G&R simulation) at a set of well-distributed initial points in an *N*-dimensional space using Latin hypercube sampling. Then the algorithm moves to a SEARCH step that accelerates convergence by fitting a Kriging surrogate function through evaluated points and finding a candidate point, usually the minima of a Krigged surface, where a function evaluation is then performed. The surrogate is updated with each function evaluation and this process is repeated until a SEARCH fails to identify an improved design point. When the SEARCH step fails, the algorithm moves to a POLL step. In a POLL step, a positively spanning basis is identified around the current design point and new trial design points are selected by performing function evaluations along these directions with a specified step size. In our implementation, due to its superior convergence properties, we use a Mesh Adaptive Direct Search (MADS) algorithm to generate POLL sets [21]. When the POLL step also fails to generate a better design point, the grid is refined and the algorithm returns to the SEARCH step. Note that all points evaluated are required to lie on a grid, which can be refined or coarsened, and the algorithm stops when the grid is refined to a preset tolerance. The threshold grid size was set to in our optimizations. Thus, the algorithm converges to a local minimum using a combination of SEARCH and POLL steps. Note that the SEARCH step is not strictly required for convergence whereas the POLL step is necessary and sufficient for convergence [21,22]. These methods have been successfully coupled to computational fluid dynamics simulations in the setting of shape optimization for cardiovascular surgery and aerodynamics [55–57].

**A.2. Cost function**

There are many ways to define a cost function for a multi-objective problem such as G&R; we use a weighted sum of multiple objective functions. One could define an illustrative cost function as a weighted sum of individual geometric and stress-based objective functions. For example,
Weighting of terms plays a vital role in estimating the optima, and consequently the evolution and adaptation of a vessel. Experiments are critical in informing the weight on each of these terms, the lack of which led us to explore different weights. For example, setting *w*_{3} = 0 and *w*_{4} = 0 produces a stress based cost function, and further setting *w*_{1} = 1 and *w*_{2} = 1 recovers the cost function used in our work (equation (2.12)). On the other hand, setting *w*_{1} = 0 and *w*_{2} = 0 gives us a pure geometric cost function.

For comparison, figure 9 illustrates the predicted evolution of radius and thickness for the optimal solution for a modest haemodynamic perturbation (*γ* = 1, *ε* = 1.1) with a geometric cost function (*w*_{1} = 0, *w*_{2} = 0, *w*_{3} = 1 and *w*_{4} = 1). A 10% increase in flow, in an incompressible vessel, should lead to a quick increase in radius and a decrease in thickness. The geometric cost function based optimization predicted a non-physiological evolution of thickness, however, in which incompressible behaviour was lost. Thus, a cost function inconsistent with constitutive relations can predict non-physiological behaviour. Our use of a stress-based optimization was, in contrast, consistent with most mass production and survival functions.

- Received December 8, 2016.
- Accepted April 28, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.