A mathematical model of venous neointimal hyperplasia formation.

BACKGROUND: In hemodialysis patients, the most common cause of vascular access failure is neointimal hyperplasia of vascular smooth muscle cells at the venous anastomosis of arteriovenous fistulas and grafts. The release of growth factors due to surgical injury, oxidative stress and turbulent flow has been suggested as a possible mechanism for neointimal hyperplasia.
RESULTS: In this work, we construct a mathematical model which analyzes the role that growth factors might play in the stenosis at the venous anastomosis. The model consists of a system of partial differential equations describing the influence of oxidative stress and turbulent flow on growth factors, the interaction among growth factors, smooth muscle cells, and extracellular matrix, and the subsequent effect on the stenosis at the venous anastomosis, which, in turn, affects the level of oxidative stress and degree of turbulent flow. Computer simulations suggest that our model can be used to predict access stenosis as a function of the initial concentration of the growth factors inside the intimal-luminal space.
CONCLUSION: The proposed model describes the formation of venous neointimal hyperplasia, based on pathogenic mechanisms. The results suggest that interventions aimed at specific growth factors may be successful in prolonging the life of the vascular access, while reducing the costs of vascular access maintenance. The model may also provide indication of when invasive access surveillance to repair stenosis should be undertaken.

Background

Vascular access dysfunction in chronic hemodialysis patients

Healthy kidneys filter wastes from blood and regulate electrolyte, acid-base, and volume homeostasis. When the kidneys fail, one needs treatment to replace the work the kidneys normally perform. One available treatment is hemodialysis, which utilizes an artificial kidney. The patients' blood is pumped into the artificial kidney where metabolic waste products diffuse out of the blood, and the cleansed blood is then returned back to the body. In order to perform hemodialysis, the patient must have suitable vascular access to allow adequate flow of blood to the hemodialysis circuit.

The most common types of vascular access used for hemodialysis are the arteriovenous (AV) fistula and the expanded polytetrafluoroethylene (ePTFE) graft. A surgeon creates an AV fistula by directly connecting an artery to a vein, usually in the forearm. The increased blood flow causes the vein to hypertrophy so that it can be used for repeated needle insertions. A graft connects an artery to a vein by using a synthetic tube of ePTFE, usually in the shape of a loop. It does not require as much time to mature as a fistula, so it can be used soon after placement. The direct purpose of the graft is to provide a vessel that is close to the skin (unlike the arteries) and has a high enough pressure to provide a sustained flow rate over 300 ml/min without collapsing (unlike the veins).

Both types of vascular access can have complications that require further treatment or surgery [1,2]. The data analysis of the Dialysis Outcomes Quality Initiative panel [2,3] suggests a primary patency of 85% for AV fistulas at one year and 75% at two years, whereas the ePTFE graft patency can be as low as 50% after one year and 20% at two years. These data exclude fistulae that did not mature adequately to support hemodialysis.

Over the last thirty years, hemodialysis vascular access dysfunction has been a major cause of morbidity and hospitalization among hemodialysis patients worldwide [4]. In the US alone, it is responsible for the hospitalization of more than 20% of patients with end-stage renal disease, at an annual cost of 1 billion dollars [2]. Novel monitoring and intervention programs, such as balloon angioplasty and surgery to open or bypass the stenosed segment, have improved the patency of native fistulae as well as ePTFE grafts, but at a significant financial cost. The expense of creating and maintaining vascular access for patients on dialysis accounts for a significant portion of any health care system. The intervention rates for ePTFE grafts are currently running six times higher than for fistulae [5]. While infections account for 10–15% of the failure of the ePTFE grafts, the leading cause of access failure is from loss of patency due to venous stenosis. Venous stenosis is the result of neointimal hyperplasia and luminal narrowing or occlusion [6-8], either at the site of venous anastomosis or in the downstream (proximal) vein. We assume that both AV fistulae and ePTFE grafts have similar mechanisms of venous neointimal hyperplasia. However, these accesses are inherently different with different flow characteristics. The model described here is more likely to be applicable to ePTFE grafts, rather than AV fistulae, due to exuberant inflammation produced by synthetic ePTFE graft.

Pathogenesis of venous neointimal hyperplasia (VNH)

The most important events initiating the pathogenesis of VNH are: (a) surgical injury at the time of creation of the vascular access, as the vein is often stretched and manipulated; (b) hemodynamic stress at the graft-vein or artery-vein anastomosis, as a result of a combination of high shear stress and turbulence [2,9,10]; (c) the presence of the ePTFE graft itself, as a foreign body, which can attract macrophages that release cytokines and growth factors [2,11]; and (d) vascular access injury from dialysis needles. Other possible causes for VNH formation are: (e) differences in diameters between arteries and veins and less defined intimal layer may cause harmful fluid ebbs and backflow [2]; and (f) genetic predisposition of veins to vasoconstriction and neointimal hyperplasia after injury to endothelial and smooth muscle cells [12,13]. Treatment of an initial stenosis is often accomplished by balloon angioplasty. However, this treatment may inflict endothelial and smooth muscle cell injury, predisposing the vein to exaggerated VNH and repeated stenosis [2].

All the above stenosis-initiating events result in the activation of the smooth muscle cells and fibroblasts of the vascular media and adventitia, with migration into the intima and proliferation. In addition, there is a significant adventitial angiogenesis and excessive intimal synthesis of collagen [7,11]. This excess extracellular matrix (ECM) creates a neointimal expansion that contributes to access stenosis [14]. Access stenosis predisposes to access thrombosis and subsequently to access failure [15]. Thus, the so-called neo-intima is composed of vascular smooth muscle cells that are derived from all three layers of the vein.

Various groups [11,15-17] have demonstrated the expression of a number of chemical mediators during the pathogenesis of VNH, some of which could be potential therapeutic targets [2]. It has been demonstrated that (i) transforming growth factor-beta (TGF-β) stimulates smooth muscle cell growth and matrix production, and inhibits the degradation of matrix proteins [15,18,19]; (ii) platelet-derived growth factor (PDGF) has strong mitogenic and chemotactic effects on smooth muscle cells [7,20]; and (iii) endothelin-1 (ET-1) is a potent mitogenic peptide, and causes constriction of smooth muscle cells [16,21]. Each of these growth factors has been implicated in the occurrence of neointimal hyperplasia [16]. Several mechanisms have been suggested for enhanced production of these growth factors in neointimal hyperplasia including, in particular, oxidative stress [16] and turbulent flow [7,22].

Oxidative stress is characterized by circulating tissue proteins by oxidative activity [16]. Several studies have shown that increased levels of oxidative stress induce the production of TGF-β [16,23,24]. Other studies have implied that increased oxidative stress levels contribute to the platelet-activated release of PDGF and the production of ET-1 by endothelial cells [16,25,26].

It has also been suggested that turbulent flow of blood stimulates the mechanoreceptors on smooth muscle cells and shear-stress receptors on endothelial cells [27,28]. Turbulent flow might also stimulate the production of TGF-β since it is thought to be produced locally by smooth muscle cells as well as by macrophages and lymphocytes within the lesion created by the intimal hyperplasia [29]. Blood flow rate and the corresponding wall-shear stress can influence platelet aggregation, which, in turn, effects the production of PDGF [7,22,27]. Also, ET-1 levels increase in response to increased blood flow in the AV fistula [16,30].

Present work

Based on the above cited work, a schematic diagram illustrating some causes and effects of VNH formation is represented in Figure 1. For simplicity, some of the intermediate factors are not included in the diagram. For example, we assume that fibroblasts produce basic fibroblast growth factors (bFGF) [31]; in turn, bFGFs stimulate the production of smooth muscle cells [27]. These two facts account for the arrow going from the fibroblast to smooth muscle cells (i.e., the intermediate factor bFGF is dropped out). Also, the fibroblasts contribute to the intimal hyperplasia [2]. The fibroblasts in the neointima may acquire a smooth muscle cell-like phenotype by expressing smooth muscle actin, and thus be called myofibroblasts.

Figure 1

A schematic diagram illustrating some causes and effects of intimal hyperplasia. The red letters represent the variables in our model, while the blue numbers indicate the sources cited.

While the occurrence of VNH is well recognized, the pathogenesis of it is complex and still not well understood. Few studies have attempted to analyze the pathways that lead to dialysis access stenosis and direct attention to potential therapies [2,11]. Computational and mathematical tools have been applied to many areas of biology resulting in descriptive models with predictive capabilities. However, to our knowledge, there is no mathematical model to account for cellular and molecular interactions relevant to hemodialysis vascular access dysfunction. In the present work, we propose such a model for venous neointimal hyperplasia development describing:

• the interaction among growth factors, smooth muscle cells, and fibroblasts;

• the effect of these interactions on the venous stenosis;

• the effect of the stenosis on the level of oxidative stress and degree of turbulent flow;

• the influence of oxidative stress and turbulent flow on growth factors.

In the next section we introduce the mathematical model and illustrate how the model can potentially be used to predict vascular access failure based on the concentration of growth factors. The goal of any surveillance method is to detect access stenosis in a timely manner so that appropriate corrective steps can be undertaken prior to thrombosis. This is of critical importance, since the access survival after an episode of thrombosis is markedly reduced. With this in mind, we discuss possible applications of our results, not only to identify vascular access at the risk of thrombosis, but also for using the model to develop innovative strategies to prevent or delay vascular access failure. We conclude the work with comments on the mathematical model and future directions.

Methods and Results

Model description

The mathematical model that describes the VNH development is based on a simplification of the network diagram of Figure 1. However, we hope that the features retained for discussion are those of greatest importance in the present state of knowledge. The process of developing the model will identify important parameters and relationships that have not yet been investigated and can thus promote refinement in future studies.

To begin with, we identify the model variables and consider their movement, production and death in a radially symmetric control domain, Ω, that represents the intima and the lumen of the blood vessel at a cross-section where a stenosis develops. The geometry of the domain is specified by the radius L = R0 + d, where R0 is the average radius of the lumen before the neointimal layer starts to form, and dstands for the average thickness of the venous intimal layer. In this setting, the boundary of the domain, Γ, corresponds to the interface between the media and the intima.

We now motivate our choice of the variables. For simplicity, we lump together several chemical species elemental to the process of neointima formation, as well as several cells and extracellular matrix components:

• a(x, t), general chemical species (TGF-β, PDGF, ET-1);

• s(x, t), general cellular species (smooth muscle cells, fibroblasts);

• ρ(x, t), extracellular matrix (collagen, fibronectin, elastin).

The quantity a(x, t) represents the concentration (in g/cm3) of growth factors at x ∈ Ω in time t. In the absence of more detailed information on each factor, a accounts for all growth factors that potentially have a chemotactic effect on the cells. However, it is possible to separately describe the mechanism of action of particular growth factors as the model expands.

The quantity s(x, t) represents the density (in g/cm3) of cells at x ∈ Ω in time t. We do not distinguish between various cells that are known to be involved in the formation of the neointimal hyperplasia, assuming instead that they all follow the same process of diffusion, chemotaxis and growth.

The quantity ρ(x, t) represents the density (in g/cm3) of extracellular matrix at x ∈ Ω in time t. Although the matrix ρ and the cellular species s have different geometric features, for the purpose of this paper we assume that they both act as a source of material filling in the intimal-luminal space, and consequently we treat them in the same way.

To study the impact of the chemicals, cells and ECM on stenosis, we chose to monitor the reduction of the luminal volume ω(t), which is initially ω0 (according to clinicians, vascular access needs clinical intervention when the neointimal hyperplasia obstructs more than 50% of the initial luminal space, that is, when ω(t) = ω0/2). As the luminal space gets partially filled with cells s and extracellular matrix ρ, the boundary of the luminal space is not clearly defined. We take the point of view that the more material there is in the intimal-luminal domain, the smaller the luminal space will be, and simply define

where k is a dimensional constant.

Applying the laws of mass conservation to each of our variables we obtain the equations governing the evolution of a, s and ρ.

Chemical species

At the time t > 0 and the position x ∈ Ω, the concentration of chemicals changes according to

We assume that the chemical species undergo random motion (i.e., diffusion). Although the diffusion coefficient Dmay in general depend on position, we take it here to be constant. Due to chemical signaling, the chemical species decrease through uptake by the cellular species. The value of the parameter λ is determined by the receptivity of cells to the growth factors. In the absence of more detailed information, we simply assume that the production rate of all growth factors is proportional to ω0 - ω(t). This term represents the observation that the production of chemical species depends on the oxidative stress and turbulent flow caused by the narrowing of the luminal space. We assume that the smaller the luminal space, the larger the oxidative pressure and shear flow, and also the larger the concentration of growth factors. Thus, the production of chemicals within the lesion is triggered by a large number of factors, which includes inflammation, hemodynamic and mechanical stresses.

Cellular species

The density of cells is assumed to follow the equation

The cellular species undergo random motion, are chemotactically attracted to the chemicals in the presence of extracellular matrix, and grow up to a maximal value S. The chemotactic force is proportional to s∇a. We assume that the movement of cells due to chemotaxis cannot occur without extracellular matrix, which has maximum density P. For simplicity, the diffusion coefficient Dand the chemotactic coefficient χare considered constants. The parameter c2 of the logistic growth term depends on the whole family of growth factors, but for simplicity we have taken it to be constant. We note that in the expression for the chemotaxis we have lumped together all the cells (by s) and all the growth factors (by a). In an extended model one would quantify the effect of each specific growth factor on the proliferation of each cell type when the growth factors are separately modeled.

Extracellular matrix

We assume that extracellular matrix is being produced by cellular species, up to a maximum value P,

We assume that the overproduction of extracellular matrix during the formation of VNH exceeds the degradation of the extracellular matrix, so that there is a total gain of the ECM density at rate c3, as long as the density is not saturated; for simplicity, we assume that c3 is constant.

Boundary and initial conditions

To complete the description of our model, it remains to specify the boundary and initial conditions for each of the variables. To begin with, we denote by a(x, 0) = a0 > 0 the initial concentration of growth factors in the proximal vein, at a cross-section characterized by the radius R(0) = R0. We further assume that no cellular species or extracellular matrix are present in the intimal-luminal space at time t = 0, hence s(x, 0) = 0 and ρ(x, 0) = 0.

If there is an influx of growth factors from the media-adventitia into the intima, we assume it is negligible compared to the production of the growth factors due to oxidative stresses and turbulent flow.

Consequently, we do not model the contribution of any factors from the medial-adventitial layers or nonvascular wall tissues, and therefore take

At low concentrations of chemicals inside the domain, there is no tendency for cells to cross the boundary into the intima. As the concentration of growth factors increase, a threshold concentration (a = A) is reached inside the domain, triggering the migration of cellular species from the medial-adventitial layers into the intima through the media-intima boundary. We assume a constant influx rate, β, and write

although in a more general case, the rate of this influx of cells could depend on the concentration of chemicals. The term H(.) is the Heaviside step function, defined as H(v) = 0 when v < 0 and H(v) = 1 for v ≥ 0, and it is used to represent the chemical signal that switches on as soon as the density arises above a threshold A.

Finally, to account for the inability of extracellular matrix to pass through the boundary, we impose a no-flux condition for ρ, namely

Parameter values

Table 1 gives a summary of the parameters and their numerical values used in the computer simulations to solve the PDE system (2)–(4) with the boundary conditions (5)–(7). The model parameters were obtained from a wide variety of experiments on many different human or animal models. Whenever such data were not available, we estimated the order of magnitude of the parameters and made choices that gave biologically reasonable results.

Table 1

Model parameters and their numerical values. Where no reference is given, the value chosen is our estimate.

	Description	Dimensional values	References
Da	Diffusion coefficient of chemicals	2.6 × 10-7 cm2 s-1	[32-34]
Ds	Diffusion coefficient of cells	8 × 10-10 cm2 s-1	[35]
C1	Production rate of chemicals	0.01 cm-4 s-1 g
C2	Proliferation rate of cells	0.0202 day-1	[36]
C3	Production rate of ECM by cells	0.01 h-1
λ	Removal rate of chemicals	6.5 × 104 cm3 s-1 g-1	[37-39]
χa	Chemotactic coefficient	2.9 × 102 cm5 s-1 g-1	[36],[40]
A	Intimal chemical concentration threshold	10-9 g cm-3	[40]
S	Maximum cells density in O	1 g cm-3
P	Maximum ECM density in O	1 g cm-3
R0	Healthy vein luminal radius	1.35 mm
dINT	Thickness of the venous intimal layer	0.15 mm
βs	Cells movement across the boundary rate	0.2 cm-4 g
k	Neointima formation rate	0.02 g-1 cm3

Model can be used to predict vascular access stenosis

In order to make predictions using the model described in the previous section, we numerically solve Eq. (2)-(4) for a0 = 25 × 10-7 g cm-3 and R0 = 1.35 mm. That is, when the lumen radius is 1.35 mm and the total concentration of growth factors inside the intimal-luminal cross-section of a vein is 25 × 10-7 g cm-3. Given this initial data, we compute how the luminal radius changes over a year of dialysis treatment; for simplicity we do not take account of any effects caused by the actual treatment (i.e., needle punctures).

In the 2-D case the lumen is a disc of radius R(t) and ω(t) = πR2(t). Since some important parameters are still currently unknown, initial understanding of the values of these parameters can be gained by doing the simulation in the simpler 1-D case. Furthermore, the results in the 1-D case already suggest strategies for delaying stenosis.

In the 1-D case, the lumen at each time t occupies an interval 0

The resulting time-dependent graph is the black curve shown in Figure 2. The red horizontal line marks the critical luminal radius associated with stenosis, here considered to be half of the luminal radius of a healthy vein. As time increases, the neointimal hyperplasia forms, decreasing the luminal radius. From this result, we see that the blockage caused by the formation of neointimal hyperplasia in a patient with 25 × 10-7 g cm-3 initial growth factors concentration will reach the critical stenosed state after approximately 8 months of dialysis.

Figure 2

Numerical simulations showing that decrease of initial concentration of growth factors in the proximal vein, a(0), by a factor of 5 delays the onset of stenosis by more than 2 months.

Impact of the growth factors on the vascular access lifespan

To better understand the effect of the concentration of growth factors on the development of VNH, we performed another numerical simulation, with a different input value for a0. We decrease the initial growth factors concentration from a0 = 25 × 10-7 g cm-3 to a0 = 5 × 10-7 g cm-3. As before, we computed the changes in the luminal radius when the patient undergoes dialysis for over a year.

The result is the blue curve shown in Figure 2. We see that a drop in the initial concentration of growth factors delays the access stenosis by more than 2 months, prolonging the lifespan of the vascular access to more than 10 months. This implies that one mechanism by which the functional state of the hemodialysis vascular accesses can be extended is to control the concentration of the growth factors in the proximal vein. Our model and simulations, which build on cellular events leading to VNH formation, suggest that interventions aimed at specific chemical mediators involved in VNH formation may be successful in reducing the human and economic costs of vascular access dysfunction.

Conclusion

The process of VNH formation is complex, involving a number of growth factors, different types of cells, ECM, oxidative stress, and fluid flow. Figure 1 illustrates the main interactions among these players. These interactions can be described in terms of a large system of partial differential equations. In this paper, we have developed a simple model in which we have lumped together all the chemical species into one variable, all the cellular species into one generic cell type, and treated the ECM as one concentration of connective tissue. We have also accounted for oxidative stress by having the growth factors increase as the luminal space decreases. Although our model is relatively simple, it captures some of the main features of VNH formation; in particular, it realistically predicts the stenotic event as a function of the initial concentration of the growth factors inside the intimal-luminal space.

Future modeling opportunities

Our model represents a first step toward the development of a more realistic model that can be used by clinicians to identify vascular access at the risk of thrombosis, and to prevent or delay vascular access failure. Of the future modeling extensions needed to achieve this goal, perhaps the most important is a more accurate numerical approach. Rather than treating the stenotic lesion symmetrically in an 1-dimensional environment, in a future model we plan to develop higher dimensional numerical methods to investigate different geometries consistent with the complex nature of the VNH. Before embarking on a more detailed inspection of VNH formation it seems crucial to have a set of robust parameters. As only limited empirical data for various parameters is available at present, clinical studies need to be conducted in parallel with the development of the model to improve its reliability.

There are other different aspects of this project that can be improved, all aimed at better understanding of the cellular events leading to VNH formation. For simplicity, we have conglomerated all cytokines and cell types into one category, giving equal importance to all cytokines and cell types, which is unlikely to be true. However, with the current state of knowledge, it is not unreasonable to make such an assumption. As the relative importance of such factors is determined by future experiments, the model can be adjusted. Other issues that remain to be investigated concern: the contribution of chemicals from the medial and adventitial layers or from the nonvascular wall tissue; understanding the mechanical properties of the ECM building up the hyperplasia; understanding how the cells interact with the ECM; quantifying individual cell motion and cell-cell/cell-chemical interactions.

Clinical relevance

With cooperative effort (i.e., interplay between computational experiments and data) this model can be (expanded and) used by clinical researchers as a testbed for exploring and evaluating various therapies that can target both the traditional and the alternative pathways that are involved in the pathogenesis of VNH and vascular stenosis. In particular, our model suggests that clinical trials need to be conducted to examine the currently available agents that are known to inhibit the production of growth factors by smooth muscle cells, fibroblasts, or various other cells involved in the process of VNH formation.

Assuming this model is validated clinically, it could be applied in two main ways to address access function. First, because the model is predictive of 50% stenosis, it will provide an indication of when invasive access surveillance with the intention to repair stenosis should be undertaken. The model could be prospectively compared to current indicators of access intervention (declining flow rate, venous pressure) for predictive value, and the efficacy of repair on access lifespan could be determined for the model, and compared to current practices. This could be particularly useful in the case of patients that have undergone repeated access repair, since the current prognosis for such cases is less certain.

Second, because the model is based on pathogenic mechanisms, it can be used to design and test interventions that may prevent access stenosis. For example, methods to decrease oxidative stress could be prospectively tested to determine how they effect time to stenosis. Similarly, when clinically available, agents to reduce cytokine/growth factor expression in the access could be tested to determine how they extend the time to failure.

Competing interests

The author(s) declare that they have no competing interests.

Authors' contributions

PBG, RCS and AF formulated the model equations and wrote the manuscript. PBG performed the numerical calculations. CV, AKA and BHR were consulted on the model during the preparation of the paper, and all authors read and approved the manuscript.