UFDC Home  UF Theses & Dissertations   Help 
Material Information
Thesis/Dissertation Information
Subjects
Notes
Record Information

Full Text 
PAGE 1 1 MODELING OF HEMODYNAMICALLY I NDUCED VEIN GRAFT REMODELING By MINKI HWANG A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORID A IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2011 PAGE 2 2 2011 Minki Hwang PAGE 3 3 To my parents PAGE 4 4 ACKNOWLEDGMENTS I thank my advisor Dr. R oger TranSonTay for giving me this opportunity to study and his guidance throughout this study. I thank Dr. Scott A. Berceli for his supports and helps during the entire period of my study. I thank Dr Marc Garbey for helping me much improve the research resu lts. I thank my committee professors Dr. Balachandar, Dr. Kim, Dr. Sarntinoranont, and Dr Batich for their advisory on this study. I thank the lab members of Dr TranSonTays lab, and the members of Dr. Bercelis lab. I thank my brother and his family for t heir support during this study. I couldnt thank my parents enough for their love and support. My parents made me who I am today with their unbelievable love towards me and my brother. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDG MENTS .................................................................................................. 4 LIST OF TABLES ............................................................................................................ 7 LIST OF FI GURES .......................................................................................................... 8 ABSTRACT ................................................................................................................... 10 CHAPTER 1 INTRODUC TION .................................................................................................... 12 2 BACKGRO UND ...................................................................................................... 15 Hemodynamic Impact on Vein Graft Rem odeling ................................................... 15 Previous Models of Bl ood Vessel Rem odeling ....................................................... 19 Modeling and Role of Leukocyt es in Blood Vessel ................................................. 19 Leukocyte Rolling Mechanism .......................................................................... 22 Leukocyte Rolling Models ................................................................................. 24 RuleBased Simulation of Biological S ystems ........................................................ 26 3 SINGLELAYER MATHEMATICAL MODEL OF VEIN GRAFT REMODELING ...... 41 Model Devel opment ................................................................................................ 42 Intimal Thickening in Rabbit M odel ................................................................... 42 Predictive Model of Intimal Gr owth ................................................................... 43 Model Applic ation ................................................................................................... 47 Focal Vein Graft St enosis Cons truct ................................................................ 47 Simulation of Intimal Thicke ning around a Foca l Stenosis ............................... 47 Parameter Sensitiv ity Analys is ......................................................................... 52 Discussio n .............................................................................................................. 53 4 MULTILAYER MATHEMATICAL MODE L OF VEIN GRAFT REMODELING ........ 56 Model Devel opment ................................................................................................ 57 Remodeling Characteristi cs of Vein Graft Wall ................................................. 57 Intima Mo del ..................................................................................................... 58 Media Model ..................................................................................................... 61 EEL Model ........................................................................................................ 63 Sensitivity and Stocha stic Anal ysis .................................................................. 65 Model Applic ation ................................................................................................... 68 Simulation of Idealized Stenosis ....................................................................... 68 Simulation of Human Vein Graf t ....................................................................... 75 Discussio n .............................................................................................................. 75 PAGE 6 6 5 RULEBASED MODEL OF VEI N GRAFT REMO DELING ...................................... 79 Determination of Probabilities of Cellular Acti vities ................................................. 80 Probability of Smooth Muscl e Cell Divis ion ...................................................... 80 Probability of Smooth Muscle Cell Apoptos is ................................................... 84 Probability of Extracellular Matrix Behav iors .................................................... 87 Effects of M onocytes ............................................................................................... 88 Monocyte Entry Model ...................................................................................... 88 Monocyte Effects on EC M synthesis ................................................................ 88 Algorithm ................................................................................................................. 90 1D Algorit hm ................................................................................................... 90 2D Algorit hm ................................................................................................... 90 Results .................................................................................................................... 92 1D Model Vali dation ........................................................................................ 92 2D Model Vali dation ........................................................................................ 94 Focal Stenosis Simulati on ................................................................................ 97 Discussio n .............................................................................................................. 98 6 CONCLUSION ...................................................................................................... 101 LIST OF REFE RENCES ............................................................................................. 105 BIOGRAPHICAL SK ETCH .......................................................................................... 117 PAGE 7 7 LIST OF TABLES Table page 31 Parameter sensitivity analysis for singlelaye r model ......................................... 53 41 Parameter sensitivity analys is for multila yer model ........................................... 6751 Assumptions used in the rulebas ed model of vein graft remodeling .................. 88 PAGE 8 8 LIST OF FIGURES Figure page 11 Vein bypass surgery in coronar y artery and infrai nguinal artery ......................... 13 21 Hemodynamic impact on ve in graft re modeli ng .................................................. 17 22 Biologic mechanism of intimal hyper plaisa ......................................................... 18 23 Schematic of computational domain for rulebased simulation of multicellular biological system ................................................................................................ 27 24 Schematic of cell behavior impl ementation on t he lattice ................................... 28 31 Schematic of vein graf t ....................................................................................... 42 32 Intimal thickenss as a function of s hear stress at each measurement day. ........ 45 33 Surface plots of intimal thick ness and intimal thi ckening ra te ............................. 46 34 Iteration procedure for simulati on of intimal th ickening. ...................................... 48 35 Hemodynamic characteristi cs inside the foca l stenosis ..................................... 50 36 Comparison of focal stenosis simula tion results with exper imental data. ........... 51 41 Schematic of the crossse ction of venous wall. .................................................. 58 42 Normalized geometical variables of vein graft with respect to shear stress. ....... 59 43 Normalized medial area as a function of time ..................................................... 61 44 Rate of normalized medial area change as a function of wall tensio n. ............... 62 45 Normalized EEL radius as a function of time ..................................................... 64 46 Surface plots of the normalized geomet rical variables of vein graft. ................... 66 47 Stochastic analysis of the vein graft mathematical model. ................................. 68 48 Schematic of stenosis simulation. ....................................................................... 69 49 Simulation results for intimal area around the stenosis under low flow conditi on ............................................................................................................. 72 410 Simulation results for intimal area around the stenosis under high flow condition ............................................................................................................ 73 PAGE 9 9 411 Simulation results for EEL radius and m edial area under low flow condition. ..... 74 412 Simulation results for human vein gr aft. ............................................................. 76 51 Schematic of cells undergoing cell cycles during 1 day prior to harvest. ............ 81 52 % of BrdU positive cells as a f unction of normalized shear stress. ..................... 82 53 % of BrdU positive cells as a functi on of actual distance from endothelium ....... 83 54 % of TUNEL positive cells as a function of time ................................................ 85 55 % of TUNEL positive cells as a func tion of actual dist ance from IEL. ................. 86 56 Simulated intimal area with normaliz ed rate of mono cyte entry .......................... 89 57 Schematic of 1D domain. .................................................................................. 90 58 Flow chart for the rulebased simu lation of vein gr aft remodel ing....................... 91 59 2D algorithm of rulebased simu lation of vein gr aft remodel ing ......................... 92 510 Element redistribution algorit hm for stenosis simulati on. .................................... 93 511 1D simulati on result s. ........................................................................................ 94 512 An example of 2D simulati on. ............................................................................ 95 513 Simulation results and experimental dat a of intimal area and comparison of 1D and 2D simula tions. .................................................................................... 96 514 An example of focal stenosis simu lation ............................................................. 97 515 Comparison between simulation resu lt and experimental data for intimal thickening around focal stenosis ......................................................................... 98 PAGE 10 10 Abstract of Dissertation Pr esented to the Graduate School of the University of Florida in Partial Fulf illment of the Requirements for t he Degree of Doctor of Philosophy MODELING OF HEMODYNAMICALLY I NDUCED VEIN GRAFT REMODELING By Minki Hwang August 2011 Chair: Roger TranSonTay Major: Mechanical Engineering Though vein grafting, which is one of the primary treatment options for arterial occlusive disease, provides satisfactory resu lts at an earlier stage of the treatment, the patency is limited to a few mont hs for many patients. Early adaptation of vein graft to arterial flow environment, which precedes the later vein graft failure, has been found to be influenced by hemodynamic forces. The goal of this study is to develop mathematical and computational models of vein graft re modeling induced by hemodynamic forces, which can potentially be used as predictors of vein graft failu re. This study consists of three components: singlelayer mathematical model, multilayer mathematical model, and rulebased model. The singlelayer model describes intimal layer remodeling as a function of wall shear stress and time. The mathematical form is simple and predicts intimal thickening around focal stenosis, whic h agrees very well with experimental data. The multilayer model describes the remodel ing of three layers: intima, media, and external elastic lamina. T he effect of wall tension is included as an additional hemodynamic factor. All the three layer m odels are based on experimental data, and the models can be applied to the simulation of idealized stenosis as well as irregular geometry of human vein graft. The rulebased model incorporates biologic mechanism PAGE 11 11 of vein graft remodeling. Rules of behav iors are applied to smooth muscle cell, monocyte, and extracellular matrix, and t he wall thickening (global behavior) is observed emerging from the local interact ions among the components. Monocytes have been found to have significant contribution to the wall thickening. The models developed in this study provide detailed insight into the mechanism of vein graft remodeling, and are a foundation for the future development of longterm predictor that can be used in clinical environment. PAGE 12 12 CHAPTER 1 INTRODUCTION Every year since 1900 except 1918, card iovascular disease (CVD) has been the leading cause of death in t he United States (Rosamond et al., 2007; Roger et al., 2011) accounting for 33.6% of all deaths in the Un ited States in 2007. One in three American adults are estimated to have 1 or more types of CVD. The CVD d eath rates are higher than those of cancer. There are three major CVDs: heart attack, stroke, and peripheral vascular occlusive disease (PVOD). Heart a ttack occurs when the coronary artery which supplies oxygenrich blood to the heart muscle is blocked due to atherosclerotic lesion. Stroke occurs when the vasculature in the br ain is blocked by atherosclerotic or blood clot originated from diseased blood vessel in other part of the body. PVOD is caused by atherosclerotic plaque in the ar tery resulting in blockage of the supply of blood to the organs other than heart and brain. All of these major CVDs ar e related to the occlusion of artery. One of the primary treatments for the occluded arte ries in heart attack and PVOD is vein bypass grafting. Autologous vein su ch as saphenous vein in leg is extracted and attached to the diseased artery to make the blood flow bypass the occluded region (Figure 11). Approximately, 400,000 co ronary and 100,000 lowerextremity bypass graftings are performed annually in the United St ates (Jacot et al., 2004). However, the problem is that, although the vein grafting prov ides satisfactory resu lts at an early stage of the treatment, the patency is limited to a few months for many patients (Fitzgibbon et al., 1996). Vein graft failure occurs in 30% to 40% of lower extremity cases at 5 years, and more than 50% of coronary grafts fail wit hin 10 years (Jacot et al., 2004). Although PAGE 13 13 it is due to atherosclerotic lesion when vein graft failure occurs, t he early remodeling of vein graft precedes the at herogenesis (Figure 11). Figure 11. Vein bypass graft surgery in co ronary artery and infraing uinal artery (Jeremy et al., 2007 with permission). When vein graft is implanted in the arteri al system, the vein gr aft is exposed to higher blood pressure and wall shear st ress than those of the venous system. The PAGE 14 14 morphology of the implanted vein graft changes such that the vein graft wall becomes thicker and the graft diameter increases. It is believed that these early stage remodeling processes are adaptive processes to the high pressure and high shear environment of the arterial system (Dobrin et al., 1989; Fillinger et al., 1994; Zwolak et al., 1987; Galt et al., 1993). Though the biologic mechanism under lying the vein gra ft failure are not completely understood, early vein graft remodeling induced by hemodynamic forces such as wall shear stress stands among the pr imary regulators for these events (Varty et al., 1993; Mills et al., 1995). In this dissertation, early vein graft re modeling is studied using mathematical and computational modeling approach. More specifically, mat hematical and computational models of early vein graft wall remodeling are developed in order to gain insight into the exact mechanism of vein gr aft remodeling and to lay the foundation for developing predictors of vein graft failure. PAGE 15 15 CHAPTER 2 BACKGROUND Hemodynamic Impact on Vein Graft Remodeling Venous wall consists of three layers: the intima, media and adventitia. The intima is the innermost layer in direct contact with the blood flow. The intima consists of monolayer of endothelial cell s though there are sometimes other components such as smooth muscle cells. The media, which is t he layer between the intima and adventitia, consists of smooth muscle cells and extracellu lar matrix. The adventitia is the outermost layer which is in contact with surrounding ti ssue. Between the intima and media lies a connective tissue layer known as the internal elastic lamina (IEL), and the external elastic lamina (EEL) between the media and adventitia. Shear stress has been recognized as one of the major hemodynamic factors which modulate vein graft remodeling. It has been found that lower shear stress increases the vein graft wall thickness. One of the earliest observations that the flow rate is related to the vein graft failur e is the one by Grondin et al. (1971). They examined the vein grafts implanted into hum an arteries, and observed that all grafts with flow of 20ml/min or le ss measured at implantation became occluded after 10 to 21 days while all grafts with the implantati on flow greater than 45ml/min remained open. Later, Berguer et al. (1980) hy pothesized that the low shear stress was related to intimal hyperplasia in vein graft, and since then, numerous evidences which support the low shear hypothesis have been reported. Dobrin et al. (1989) separated the effects of various hemodynamic factors in canine ve in graft experiment and showed that the intimal hyperplasia is best associated with low flow velocity. Galt et al. (1993) observed the inverse relationship between flow rate and in timal thickness in their rabbit vein graft PAGE 16 16 experiment. Jiang et al. (2004) observed a se venfold greater intimal thickness 28 days after implantation when the flow rate re mained 15fold lower throughout most of the perfusion period in their rabbit bi lateral vein graft experiment. While lower shear stress increases inti mal thickness, higher shear stress has been known to increase vein graft diameter Fillinger et al. ( 1994) reported a linear relationship between the percent change in ve in graft diameter over 12 months and initial shear stress in human saphenous vein grafts after in frainguinal bypass. Jiang et al. (2004) observed a 24% increase in lumen diameter at 28 days in the vein graft under elevated flow condition in their rabbit bilate ral vein graft experimen t. Owens et al. (2006) reported a positive correlation between the in itial shear stress and the change in lumen diameter during the time interval of 1 month after autogenous ve in implantation in human patients undergoing lower extremity bypass. The ve in graft enlargement under high flow condition can be understood with regard to the Glagov phenomena (Glagov et al., 1987) which is the observati on that the coronary arteries enlarge in the early stages of atherosclerosis to preserve lumen area. Wall tension induced by blood pressure is another important hemodynamic factor in vein graft remodeling. Higher wall tensio n has been known to increase the vein graft wall thickness. Zwolak et al. (1987) observed that the thickening wall in rabbit jugular vein segments transplanted into the carotid ar tery caused the ratio of luminal radius to wall thickness to decrease to a level equal to that of normal artery. Schwartz et al. (1992) reported a linear relationship between myointimal area and wa ll tension in their rabbit vein graft experiment. PAGE 17 17 The effect of each hemodynamic factor on vein graft remodeling is summarized in Figure 21. Figure 21. Hemodynamic impac t on vein graft remodeling. Two of the representat ive processes of the intimal hyperplasia in vein grafts are smooth muscle cell (SMC) proliferation and extr acellular matrix (ECM) accumulation in the intima (Figure 22). The SMCs migrate fr om the media, and t hey proliferate before migrating from the media to the intima. One of the factors th at promote the SMC proliferation in the medi a is the endothelial damage. Damaged endothelium produces the SMC growthstimulating factors, and t he amount of SMC pro liferationinhibiting factors which are produced from normal end othelium decreases due to the endothelial damage (Lemson et al., 2000). SMC in the media starts to migrate into the intima a few days after implantation due to degradation of ECM which normally prevents the SMC migration. The migrated SMCs proliferate in the intima and produce the ECM. Growth PAGE 18 18 factors such as TGFand PDGF are known to stimulate ECM production by SMC (Lemson et al., 2000). Figure 22. Biologic mechanism of intimal hyperplasia. Shear stress has been known to be an important factor in the vein graft remodeling. Altered shear st ress stimulates the expressi on of adhesion molecules and chemokines involved in vein graft wall thicke ning (Lemson et al., 2000). After the vein graft is implanted, leukocytes are recruit ed to the vein graft wall, and they release cytokines which influence the SMC prolif eration and migration (Golledge, 1997). Low shear stress has been observed to promote the leukocyte entry into the vein graft wall (Walpola et al., 1995; Ric hter et al., 2004). Outward remodeling in the early stage of the vein gr aft remodeling compensates the wall thickening resulting in the preser vation of lumen area (Galt et al., 1993). PAGE 19 19 Change in shear stress induces the endothe lial cells to release numerous biologic mediators related to blood vessel diameter change (Filli nger et al., 1994). Recently, Jiang et al. (Jiang et al., 2007) reported that the reorganizati on of adventitial myofibroblasts was the dominant histologi cal events that limit ed the early outward remodeling of vein grafts in their rabbit vein graftfi stula experiment. Previous Models of Bl ood Vessel Remodeling Previous investigators have sought to develop predictive models of intimal growth in diseased arteries or in vein grafts. Among the first were Friedman et al. (1986), who developed a model of sheardependent intimal thickening of human arterial intima based on experimental data. Because their model is based on experimental data which were measured at only one time point, they introduced the concept of a time variable based on measured intimal thick ness to describe the timedependent behavior of intimal growth. Lee and Chiu (1992) used this model to simulate th e growth of intima in a curved artery. Yang et al. (2003) assu med a linear inverse relationship between the intimal growth rate and shear stress, and appli ed the model to the simulation of intimal growth in a rat vein graft. Unique was thei r use of a modified piecewise linear growth function based on their experimental data. Z ohdi (2005) developed a theoretical model of blood vessel lumen reduction based on the concept of feedback control, employing a linear inverse relationship between the lum en radius changing rate and shear stress to examine the lumen reduction as a function of time. Modeling and Role of Leukocytes in Blood Vessels Leukocytes (white blood cells) are infla mmatory cells (diameter on the order of 10 m) circulating in the blood stream. Leukocytes migrate to the inflammation sites by sensing the inflammatory signals originated at the infection sites in the tissue. When the PAGE 20 20 circulating leukocytes migrate to the infl ammation sites, they go through several steps: initial capture, rolling, firm adhesion, and tr ansendothelial migrat ion. Leukocytes are captured from the blood stream through the binding between th e selectins expressed on the leukocytes and the ligands on the surfac e of the endothelial ce lls. The captured leukocyte rolls on the surface of the endot helium under the combined effect of shear stress imposed by the blood flow and t he association/dissociation between the endothelial selectins and the leukocyte ligands. When the chemokine receptors on the rolling leukocyte detect chemokines on the endothelial cells, the leukocyte integrins which were previously low affinity stat e undergo conformational change to become high affinity state. Through the binding between thes e high affinity integrins on the leukocyte and the ligands on the endothelium, the rolli ng leukocytes firmly adhere to the endothelium. The firm adhesion triggers the cytoskeleton di apedesis machinery leading to transendothelial migrati on (Shuhaiber et al., 2002). Monocytes enter the vein graft wall as part of an inflammatory response, and the cytokines released by macrophages are k nown to promote the intimal hyperplasia (Hoch et al., 1999). It has been reported that the monocyte entry into the blood vessel wall is affected by the wall shear stress impo sed by the blood flow. Walpola et al. (1995) decreased shear stress 73% in a rabbit caroti d artery by surgical manipulation, and observed adhered monocytes in 5 days, wh ile they did not detect any adhered monocyte when they increased the shear st ress by 170%. Richter et al. (2004) observed densely adhered mono cytes at the region of lo w shear stress and complete absence of monocytes in high shear region in side an arterial branch in their porcine model. PAGE 21 21 Leukocytes and other blood ce lls are differentiated from the pluripotent hematopoietic stem cell in t he bone marrow. Several hundr ed billion leukocytes are produced daily due to the short halflife which is typically hours to days (Silverthorn, 2001). There are five types of leukocyt es: lymphocytes, monocytes, neutrophils, eosinophils, and basophils. Monocytes develop into macrophages when they enter the tissue from the blood vessel. Neutrophils monocytes, macrophages, and eosinophils are phagocytes because they engulf foreign invaders such as bacteria. Basophils, eosinophils, and neutrophils ar e called granulocytes because granules are observed in their cytoplasm. Neutrophils are also ca lled polymorphonuclear leukocytes because their nucleus consists of several nuclear lobes. Neutrophils, the most abundant type of leukocytes, typically kill five to twenty bac teria during their life s pan of several hours to few days. On the other hand, monocytes constitu te only 1 to 6% of all leukocytes in the blood, but kill about 100 bacteria during their life span of several months (Silverthorn, 2001). Lymphocytes play an im portant role in the acqui red immune response of the body. On the surface of the leukocytes, ther e are projections of plasma membrane, called microvilli. There are 500 to 10,000 microvilli on a leuko cyte surface (Knutton et al. 1975). The length of the microv illi ranges from 50 to 3500 nm (Bruehl et al. 1996), and the surface area of t he tip of a microvillus is approximately 0.01 m2 (Bongrand and Bell 1984). At the tip of the microvilli, there ar e various adhesion molecules which mediate cell to cell interactions. The microvillus elongat es when force is exerted. At small force, the microvillus stretches like an elastic spring, and the spring constant has been measured to be 43 pN/ m in case of neutrophil micr ovilli (Shao et al. 1998). The PAGE 22 22 leukocyte plasma membrane has numerous wrinkles and folds which accommodate surface area increase when the leukocyte deforms from its in itial spherical shape. The surface area increases more than 100% when the leukocyte spreads on the endothelium before extr avasation (Dewitt and Hallett 2007). Among the adhesion molecules on the surf ace of the leukocytes and endothelial cells, selectin and integrin play critical role s in the leukocyte recruitment to inflammation sites. Selectins are transmembrane molecule s composed of Nterminal extracellular domain, epidermal growth factor, cons ensus repeats, tran smembrane domain, and cytoplasmic tail. Lselectin is expressed on leukocytes and binds to the ligands on endothelial cells. Pselectin is express ed on endothelium and platelets, and binds to PSGL1 on leukocyte. Eselectin is expre ssed on endothelial cells and binds to ligands on leukocytes. Integrins are heterodimeric gl ycoproteins which are subdivided according to their subunits. Integrins are also transmembrane molecules composed of extracellular domain and cytoplasmic tail which is linked to the cytoskeleton of the cell. VLA4 integrin on leukocyte binds to VCAM 1 ligand on endothelium upon activated by chemotactic agents. Leukocyte Rolling Mechanism Leukocyte rolling before transendothelial migrat ion in the inflammatory response is a movement coordinated by fluid shear and binding between adhesion molecules on the leukocyte and endothelium. When the bonds at the trailing edge of the contact area dissociate, the leukocyte rota tes forward due to the flui d force, and the probability of bond association at the leading edge increases. Leukocyte rolling is mediated mostly by PAGE 23 23 selectins, and different rolling velocities have been observed on different types of selectins (Jung et al. 1996). It has been observed that the rolling velocity increases as fluid shear stress increases. Lawrence and Springer (1991) repr oduced rolling of neutrophils on artificial lipid bilayers containing purified CD62 in vitro, and observed that the rolling velocity increased with increasing shear stress. They also observed that the rolling velocity increased with decreasing selectin density. Al on et al. (1995) reported that the rolling velocity of blood T lymphocytes increased w hen shear stress was increased in a parallel flow chamber assay. Dong and Lei (2000) observ ed that rolling velocities of HL60 cells which express comparable levels of PSG L1 (Pselectin ligand) as most human leukocytes increased on a surface with Pselectin adsorbed as wall shear stress increased in their in vitr o sideview flow assay. An interesting phenomenon in the leukocyte rolling is the stopandgo motion. The leukocytes frequently st op for some time as they roll on the ligands, and it presumably results from the fact that the binding between the adhesion molecules is stochastic in nature. Goetz et al. (1994) m easured the instantaneous velocity of bovine neutrophils interacting with lipopolysaccharid estimulated bovine aortic endothelium at various shear stresses, and confirmed that the neutrophils translat ed with a nonconstant velocity. The population average variance in t he instantaneous velocity was at least 2 orders of magnitude higher t han the theoretical variance generated from experimental error. Smith et al. (1999) obtai ned the dissociation constants (koff) for P, E, Lselectin from the pause times in the neut rophil rolling in their parallel plate flow chamber. The elasticity of the microvilli on t he surface of the leukocyte also has a role in the leukocyte PAGE 24 24 rolling. Park et al. (2002) compared the rolling characteristi cs of PSGL1 coated microbead and neutrophil, and repor ted that the dissociati on rates for the PSGL1 microbeads were briefer than t hose of neutrophils. They also observed the elongated microvilli at the rear of the rolling neutrophils using sca nning electron microscope. Fluid shear affects the number of tethered and rolling leukocyt es. Finger et al. (1996) reported that the number of te thered and rolling neutrophils on Eselectin and Pselectin decreased as the shear stress increased. However, they observed that, in case of Lselectin, the number of tethered and rolling neutrophils increased as the shear stress increased up to around 1 dynes/cm2, and then the number decreased with increasing shear stress over about 1 dynes/cm2. Later, Lawrence et al. (1997) reported that the threshold shear exists for Eselectin and Psele ctin as well. Dwir et al. (2003) identified shortduration (4 to 20 ms) tethers of Lselectin below 0.3 dynes/cm2 of shear stress, and postulated that shearmediated multivalent tether is required for the Lselectin tethers to become stabilized. An inte resting observation was made by Chen and Springer (2001) that the tether formation is related to shear rate rather than to shear stress. They showed that the tether frequency of neutrophils to Pselectin reached the maximum value at the shear rate of ~100 s1 regardless of viscosity. Leukocyte Rolling Model Hammer and Apte (1992) modeled the leuko cyte as a rigid sphere with microvilli on the leukocyte surface. They used the bond model described in the previous section for the reaction between the adhesion molecule s on the tip of the microvilli and on the surface on which the leukocyte is rolling. The hydrodynamic force was calculated based on Stokes solution (Goldman et al. 1967). Their model enabled them to examine the PAGE 25 25 effects of microvilli density on cell surface, receptor density on microvillus tip, the stiffness of the receptorli gand bond, and the magnit ude of the hydrodynamic force. The simulation results showed the velocity fluctuat ion which is one of the characteristics of leukocyte rolling. NDri et al. (2003) simulated the roll ing of deformable compound drop model which has nucleus inside the cell. They solved the NavierStokes equations to simulate the flow field, and used immersed boundary te chnique (Peskin, 1977) to incorporate the effect of the cell in the simulation. The Eu lerianLagrangian techniq ue (Shyy et al. 1999, 2001) was used to track the moving interfaces between the outer flow and the cell, and between the cytoplasm and nucleus. They simulated the association and dissociation of the adhesion molecules on the leukocyt e and the surface using the bond model described in the previous section. They f ound that the presence of nucleus inside the cell increases the bond life time, and decreases the cell rolling velocity. They also found that a bigger cell rolls faster, and causes bond lifetime decrease, possibly due to the fact that a larger cell exper iences greater hydrodynamic force from the blood stream. Jadhav et al. (2005) extended the model of NDri et al. (2002) to 3dimesional model. They modeled the cell as an elastic capsule, and used the immersed boundary method to simulate the motion of the caps ule in a shear flow. The nucleus was not included in the cell model. T hey simulated the forward a nd reverse reactions between the adhesion molecules on the leukocyte and t he surface based on the probabilities of each reaction (Hammer and Apte 1992). They s howed the rolling velocity fluctuation for 0.7 seconds at 200 s1 of shear rate. They observed that higher membrane stiffness PAGE 26 26 made the cell roll faster, and the number of bonds decreased with increasing membrane stiffness. Bond lifetime also decreased with increasing membrane stiffness. Pappu and Bagchi (2008) examined a sim ilar model as that of Jadhav et al. (2005). They showed the forward movement of the contact area between the cell and the surface, and the displacement of the cell as a function of time for up to 3 seconds, which showed the stopandgo motion. They observed a sideway motion during the rolling due to discrete nature of microvilli distribution. They also found that most of the adhesion force is sustained by only 1 to 3 te thered microvilli in the rearmost part of a cell. Tang et al. (2007) developed an in silico model of leukocyte rolling, activation, and adhesion. They modeled the leukocyte as having rectangular contact area which is discretized such that each element has adhes ion molecules. The surface on which the leukocyte rolls is also discretized and adhesion molecules are contained in the elements. The bindings bet ween the adhesion molecules were simulated using the probabilities of bond formation/disso ciation. They applied a fluid force which affects the dissociation of the bonds, and when all the bonds have been dissociated at the rear row of the contact area, the leukocyte rolls fo rward. Their simulation results of rolling behavior agreed well with in vitro experimenta l data, and the transit ion from rolling to adhesion was also simulated. RuleBased Simulation of Biological Systems Introduction Modeling multicellular biol ogical systems (MCBS) poses challenges in that the global system behaviors result from individual cell behavior and the interactions among PAGE 27 27 the cells. Cells are live creatures which di fferentiate, proliferate, move, and die, and these behaviors of living cells constitute the dynamics of MCBS. Moreover, these cell behaviors are influenced by the environmental factors such as extracellular matrix, chemicals, and forces. Modeling MCBS r equires incorporating these complex interactions among the individual ce lls and the environment (Figure 23). Figure 23. Schematic of co mputational domain for rulebased simulation of multicellular biological system. Modeling approaches for MCBS can be gr ouped into two categories: continuum models and cellbased models (Byrne and Dr asdo, 2009). Continuum models usually take the form of partial differential equations (PDE). One of the advantages of the PDE models is that the model equations provide insight into the relationship among the components in the system. PDE can be used for modeling all levels of biological systems from molecular to organ levels. Howe ver, the continuum models do not catch PAGE 28 28 the discrete nature of MCBS co nsisting of individual cells, and become cumbersome when modeling complex process involving m any variables. On t he other hand, the cellbased models such as cellular automata (CA) and agentbased (or individualbased) models (ABM) simulate each individual ce ll behavior and interactions among them enabling the observation of t he emergent system behavior (F igure 24). Hybrid models combine these two approaches taking the advantages of each method (Gerlee and Anderson, 2008; Guo et al., 2008). This review focuses on the cellbas ed models of MCBS, and especially, the technical aspect of the rulebased simula tion method for MCBS is reviewed. How to implement the cell behaviors and their in teractions with other cells and with the environment into the computational domain is discussed, and application examples from recent publications are introduced. Figure 24. Schematic of cell behavi or implementation on the lattice. PAGE 29 29 Cellular Automata and AgentBased Modeling CA and ABM are two of the widely used methodologies for rulebased simulation of MCBS (Alber et al., 2003; Chavali et al., 2008; Deutsch and Dormann, 2005; Ermentrout and EdelsteinKeshet, 1993; Thorne et al., 2007; Thorne et al., 2007). Both methods are believed to have originated from an idea by John von Neumann in the 1940s (Gilbert, 2008; Wolfram, 2002). Later, John von N eumann with the help of Stanislaw Ulam introduced the concept of CA (Deutsch and Dormann, 2005; Neumann, 1966; Wolfram, 2002), and ABM al so started to be establis hed by other researchers thereafter (Gilbert, 2008). Traditional CA is defined by the following components: a regular discrete lattice, a fi nite set of cell states, a fi nite set of neighboring cells, and rules for the transition of cell states (Deutsch and Dormann, 2005). In reallife application of CA, the traditiona l definition is considered t oo restrictive, and there have been many relaxations such as nonuniform grid, asynchronous update of the cell states, and extension of the cell neighbor hood (OSullivan, 2001). ABM is defined in a similar way. ABM is a computational method in which decisionmaking agents interact with the environment fo llowing a set of rule s (Bonabeau, 2002; Gilbert, 2008). CA and ABM are similar in that the behaviors of cells or agents are governed by the rules in their neighborhood or environment, and gener ate global behavior of the system emergent from the local interactions. CA, however, seems more mathematical in its formulation. Mathematical rigorousness of the CA can be seen in the graphCA proposed by OSullivan (2001), a new CA in troduced to accommodate nonuniform grid whereas the traditional CA is defined on uni form grid. Although the relationship between CA and ABM is not obvious, some researcher s think of ABM encompassing CA (Galle PAGE 30 30 et al., 2009; Zhang et al., 2009) possibly due to the more general definition of ABM (Bonabeau, 2002). Both CA and ABM have been widely used in the modeling of a variety of systems such as social, ec onomic, and biological systems (Bonabeau, 2002; Gilbert, 2008; Wolfram, 2002). Detailed modeling techniques of these methods applied to MCBS are discussed in the following sections. Lattice Rulebased simulations are usually performed on a lattice system. Cells and other components occupy some of the grid elements and can move from one element to another (Figure 23). Regular, irregular, or latticefr ee approach can be used, and each of these grid systems is discussed in this section. In the case of regular lattice, the di stance between two adjacent elements remain constant over the entire simulation domain, and hence the griddependency of the simulation results can be minimized. The latti ce can be triangular, square, or hexagonal (Deutsch and Dormann, 2005; Kim et al., 2009). Engelberg et al. (2008) reported that, in their simulation of tumor spheroid growth square grid required a higher order implementation of discrete diffusion co mpared with hexagona l grid, and generated artifacts that are not present with the hexagonal grid. The size of individual lattice element can be made to be comparable to that of the biological cell (Gerlee and Anderson, 2007; Mallet and De Pillis, 2006), or can be any value (Grant et al., 2006; Simpson et al., 2007). One biological cell (C heng et al., 2006; Simps on et al., 2007) or multiple cells (Kansal et al., 2000; Piotrow ska and Angus, 2009; Ferreira et al., 2002) can occupy one lattice element. Piotrowska and Angus (2009) reported that assigning many biological cells to one lattice site can reduce the total number of lattice sites for a PAGE 31 31 given number of biological cells, and hence the computational time, and provides the flexibility in positioning the newly created cells from cell divisions. Cellular Potts model is a type of CA in which one biological cell occupies more than one lattice sites (Alber et al., 2003; E ngelberg et al., 2011). This model enables the incorporation of the shape chan ge of each cell into the simulation. Jiang et al. (2005) simulated an avascular tumor growth using an extended largeQ Po tts model. LargeQ means that the number of possible cell states is comparable to that of the connected subdomains of different cell types (Al ber et al., 2003). Robertson et al. (2007) segmented one biological cell into nine subco mpartments enabling different portions of a cell to respond to different stimuli. Thes e subcompartments can also incorporate cell polarity and sense the spatial gradient of the environment across the cell. Kansal et al. (2000) used an irregular la ttice for CA simulati on of brain tumor growth. They used Voronoi tessellation to generate the lattice, and each resulting automaton cell takes the form of polyhedra in threedimensional space. They also used varying lattice density with higher density at t he center to allow the tumor to grow to a large size during the simulation. Due to their varying size of the el ements, the number of real cells contained in the automat on cells ranged from roughly 100 to 106 depending on the size of the automaton cell (Kansal et al., 2000). Gevertz and Torquato (2006) adopted a similar irregular lattice to investigat e the effects of vasculature on early brain tumor growth. In the case of latticefree models, cells can be at any location in the computational domain. The positions of t he cells are usually determined by solving equations of motion incorporatin g the forces acting on the cells (Galle et al., 2005; Galle PAGE 32 32 et al., 2009; Schaller and MeyerHermann, 2007). Rheological properties of the cells also can be included in the simulation (Palsson, 2008). Galle et al. (2005; 2009) simulated multicellular syst ems using latticefree models incorporating contactdependent regulation mechanisms. Schalle r and MeyerHermann (2007) simulated steadystate flow equilibrium of skin using an offlattice agentbased model. Cell Behavior Division, migration, apoptosis, necrosi s, and differentiation are among the cell behaviors that are commonly modeled in the rulebased simulations of MCBS. The rules for these behaviors are dependent upon specific biological applications. Some of the rules found in recent publications are discussed in this section. Two of the important deci sions regarding the cell division are determination of the division probability and where to position the two daughter cells. For the cell division probability, cells can be programmed to divi de after cell cycle time (Cheng et al., 2006; Grabe and Neuber, 2005; Piotrowska and An gus, 2009). Experimentally obtained cell cycle time can be applied (Cheng et al., 2006; Piotrowska and Angus, 2009; Walker et al., 2004). Each cell can be assigned diffe rent cell cycle time based on a normal distribution (Cheng et al., 2006; Piotrowska and Angus, 2009), and in that case, the daughter cells can inherit the cell cycle time of the parent cell (Piotrowska and Angus, 2009). When incremental time step t is used for time depend ent simulation, the probability of cell division c an be applied such that Pd= t/tc where Pd is the cell division probability at the time step and tc is the cell cycle time (Schaller and MeyerHermann, 2007). The cell division probability can also be calculated based on other parameters such as local nutrient concentra tion (Ferreira et al., 2002). PAGE 33 33 Cell cycle can also be modeled. In their model for the epidermis, Schaller and MeyerHermann (2007) incorporated a cell cycl e in which cells can enter different phases depending on the local environment. Wa lker and coworkers (2004) modeled a cell cycle in their simulation of epithelial cells. Jiang et al. (2005) used a simplified protein regulatory network to model the tr ansition between different phases of the cell cycle in their simulation of avascular tumor growth. Regarding where to position the two daught er cells, one of the common rules is to position them randomly at the adjacent vacant site s (Checa and Prendergast, 2009; Cheng et al., 2006; Gerlee and An derson, 2007). Contact inhibition is a commonly used rule which prevents cell division when all t he adjacent sites are al ready filled (Bartha and Rieger, 2006; Checa and Prendergast, 2009; Piotrowska and Angus, 2009). Depending on specific applicat ion, the daughter cells can be placed away from each other with a site between them (Simpson et al., 2007) or always adj acent to each other (Perez and Prendergast, 2007). Prez and Prendergast (2007) assigned different probabilities at different available sites in their anisotropic mitosis model. Piotrowska and Angus (2009) applied a proba bilistic overlay to determine the location for the daughter cells to avoid morphological artifact in their simulation of in vitro multicellular spheroid tumour growth. In t heir simulation of threedime nsional brain tumor growth, Kansal et al. (2000) used the intercellula r mechanical stress algorithm in which a daughter cell pushes an adjacent cell outward unt il the cell at the tumor edge fills the adjacent empty space. Ferreira et al. (2002) enabled the cancer cells to pile up in a given lattice site in their simu lation of avascular tumor grow th. If the dividing cancer cell was inside the tumor, the daughter cell piled up at the site, and if the cell is on the tumor PAGE 34 34 border, the daughter cell replac es the normal or necrotic cell at the nearest neighboring site (Ferreira et al., 2002). In the case of cellular Potts model in which one biological cell occupies more than one lattice sites, half of the lattice si tes in the parent cell can become a new cell (Alber et al., 2003; Jiang et al., 2005). In the case of lattic efree models, the orientation of cell division can be determined by the dire ction of the total fo rce the dividing cell experiences (Galle et al., 2005) Galle et al. (2005) included t he effect of the substrate in the determination of the cell division orientation in their si mulation of the growth of epithelial cell populations in vitro Cell growth can be modeled between the cell divisions in the models such as cellular Potts model or latticefree model both of which can accommodate cell shape change. In their Potts model of avascular tu mor growth, Jiang et al. (2005) set a target volume which each cell tries to reach, and r eaching the target volume is a condition for cell division. Galle et al. (2005) modeled the cell growth such that the cell doubles its mass and volume during the interphase in their si mulation of the growth of epithelial cell populations. When a cell is compressed by its neighbor cells and the resulting cell volume is less than a threshold value, the gr owth is inhibited (Galle et al., 2005). For random movement of cell, new loca tion of the cell can be selected randomly from one of the neighboring sites (Checa an d Prendergast, 2009; Gerlee and Anderson, 2007; Perez and Predergast, 2007). The migrati on, however, can occur several times for each proliferation step because the time sca les of the migration and the proliferation are different (Checa and Prendergast, 2009; Perez and Predergast, 2007). Contact inhibition is commonly used as in the case of cell division such that the cell cannot move PAGE 35 35 if all the neighboring sites are occupi ed (Checa and Prendergast, 2009; Gerlee and Anderson, 2007; Perez and Prendergast, 2007). C heng et al. (2006) developed a tissue growth model in which cell migration is m odeled as a persistent random walk. Each cell moves in one direction for a certain period of time until it changes its direction and continues to move. After the cell cycle time, the cell stops to divide. The two daughter cells resume their persistent random movement s. When the two cells collide, they stop for some time and start moving again (Cheng et al., 2006). For directed movement of cell, Deisboeck and coworkers (Mansury and Deisboeck, 2003; Zhang et al., 2007) chose the best location for migration am ong the neighboring sites based on the amount of nutrients, levels of toxicity, and mechanical confinement in their tu mor model In their simulation of cell movement in the prosta te epithelium, Lao and Kamei (2008) tested different movement behaviors of transit amplif ying/intermediate cells and luminal cells in the prostate duct, and compar ed with experimental data. In their model of avascular tumor growth Ferreira et al. (2002) used a probability of migration which increases with the number of tumor cells in the element. This probability also increases with t he level of nutrient. In a similar model by Mallet and De Pillis (2006) for tumorimmune system inte ractions, immune cells move randomly until they encounter a tumor cell. In the case st udy of their hybrid agentbased model for microbiological systems, Guo et al. (2008) modeled the chemotactic displacement of cells such that it is proportional to the di fference in newly bounded receptors at the front and rear of the cell. Rober tson et al. (2007) modeled th e cell migration based on the relative concentrations of fibronectin, in tegrin, and cadherin in their simulation of Xenopus laevis morphogenesis. PAGE 36 36 In the case of latticefree models, t he movements of the cells are usually computed from the equations of motion which incorporate the forces acting on the cells (Galle et al., 2009; Schaller and MeyerHermann, 2007). When apoptosis occurs, cells are usually removed from the lattice, and the site can remain vacant until it is filled with ot her cells (Galvao et al., 2008). In their simulation of tissue differentiation, Checa and Prendergast (2009) implemented apoptosis to the cells differentiated by a type of stimulus when the stimulus changed to other type. In the model of chronic c hagasic cardiomyopathy after stem cell transplantation by Galvo el al.(2008), the apoptos is of inflammatory cell occurs if there is at least one bone marrow stem cell in the neighborhood of the inflammatory cell. In their simulation of the growth of epithelial cell populations in vitro Galle et al. (Galle et al., 2005; Galler et al., 2009) made the cell undergo anoikis when the contact area to the substrate is smaller than a threshold value. One of the common ways to model necrosis is to make it occur when local nutrient concentration is below a threshold value (Jiang et al., 2005; Piotrowska and Angus, 2009). In case of tu mor growth simulations, the necrotic cells usually do not vanish from the lattice but are added to the necrotic material inside the tumor (Dormann and Deutsch, 2002; Jiang et al., 2005). When a cell differentiates into other cell, t he type of the cell in the lattice site can simply be changed. Checa and Prendergast (2009) made a portion of the mesenchymal stem cells differentiate into fibroblasts chondrocytes, or osteoblasts depending on the level of mechanical stimulus and local vascu larity when the cells have reached the maturation age. Grant et al. ( 2006) determined whether a cell would make a transition to PAGE 37 37 a more differentiated form based on the a rrangement of cell, matrix, and free space around the cell in the hexagonal grid system in their simulation of in vitro epithelial cell morphogenesis. They also included the dediffer entiation of the cell in their simulation. Environment Cell behaviors are influenced by their in teraction with the environment. Gerlee and Anderson (2007) linked the environmental factors to the cell behaviors using a response network in their model of tumor growth. They consider the environmental factors such as neighbors, oxygen concentra tion, glucose concentration, and hydrogen ion concentration, and the cell behaviors such as proliferation, quiescence, apoptosis, metabolism, and movement. The cellular res ponses to those environmental factors are determined through the response network that each cell is equipped with (Gerlee and Anderson, 2007). In this section, modeling of extracellular matrix (ECM), chemicals, microvasculature, and forces are reviewed. Chemical concentrations and forces are usually computed in a different spatial scale from that of the cell, and this information from different spatial scale is projected to th e cellular level for the cells to react to those factors. Extracellular matrix (ECM) can be placed in the lattice sites which are available for cells as well (Checa and Prendergast, 2009; Grant et al., 2006). Checa and Prendergast (2009) let their lattice sites av ailable for either cell or ECM in their simulation of tissue differentiation. They made the cells synthesize ECM after cell division so that the number of cells and matrix production can be in a proportional relation. Grant et al. (2006) made the cells produce matrix depending on the arrangement of the surrounding elem ents in their simulation of in vitro epithelial cell PAGE 38 38 morphogenesis. When matrix is generated by a cell, the cell moves to a neighboring site, and the resulting vacant site is filled with the produced matrix. Robertson et al. (2007) enabled each pixel to have different am ount of fibronectin, and affect the cell behavior. One of the popular ways to determine the concentrations of the chemicals such as nutrients is solving the reactiondiffusion par tial differential equations (Jiang et al., 2005; Ferreira et al., 2002). Finite difference method can be used for the computation of the equations (Gerlee and Ander son, 2007; Guo et al., 2008), and the grid system which cells occupy can be used for the numeric al computation of the reactiondiffusion equations as well (Gerlee and Anderson, 2007; Mallet and De Pillis, 2006). Mallet and De Pillis (2006) gained sufficient accuracy in nutrient concentrations when they used the same grid system for cells and for the reacti ondiffusion equations for nutrients in their simulation of tumorimmune system interacti ons. Their equations are based on the ones used by Ferreira et al.(2002), and for the bou ndary conditions, they assumed that the nutrients are constantly supp lied from the blood vessels located at the top and bottom sides of the computational domain. Palss on (2008) used a regular 3D grid for the calculation of cAMP concentrations while letti ng the cells unrestricted to the grid points in his 3D model of multicellular systems. As a result, the cAMP concentrations are interpolated between the cells and the grid system where the cAMP concentrations are calculated. Incorporating microvasculature in the m odel provides the sour ce for nutrients. In their modeling of the effects of vasculatu re evolution on early brain tumor growth, Gervertz and Torquato (2006) modeled the micr ovasculature evolution on a triangular PAGE 39 39 lattice overlaid on top of the lattice for cells. The interaction between the microvasculature and the cells were simulated such that the interaction is mediated by the key proteins involved in the vessel grow th and regression. The concentrations of the proteins were obtained from the numerical solution of the reactiondiffusion equations. Checa and Prendergast (2009) modeled each capillar y as a sequence of endothelial cells in their model for tissue differentiation. The growth of the vascular network on a regular lattice was simulated followi ng their random walk model. High oxygen concentration was assumed within a distance from any blood vessel. Peirce and coworkers (Bailey et al ., 2007; Bailey et al., 2009) modeled the microvasculature consisting of endothelial cells and smooth muscle cells based on the confocal microscopy image of mouse muscle, and simulated circulating inflammatory cell trafficking (Bailey et al., 2007) and adipo sederived stromal cell trafficking during ischemia (Bailey et al., 2009). Qutub and Popel (2009) simulated capillary sprouting by applying local rules to the i ndividual endothelial cells. In their model for tissue differ entiation, Checa and Prendergast (2009) incorporated the interaction between mechan ical stimuli and the cell behavior in the tissue. The mechanical stimuli affect t he differentiation of precursor cells and angiogenesis, and in turn, t he mechanical properties of the tissue change due to the resulting change of the tissue composition. The mechanical stimuli we re calculated from finite element analysis. Ausk et al. (2006) si mulated the realtime signaling induced by mechanical stimuli in osteocytic networks. Bailey et al. (2007; 2009) used the pressure differential calculated from network flow m odel to drive leukocyte movement in their microvasculature model. In latticefree model s, the forces that the cell experiences PAGE 40 40 usually determine the direction and magnitude of the cell disp lacement (Galle et al., 2009; Schaller and MeyerHermann, 2007). PAGE 41 41 CHAPTER 3 SINGLELAYER MATHEMATICAL MODE L OF VEIN GRAFT REMODELING Transposition of a vein segment into t he arterial system where flow rate and intraluminal pressure are higher than thos e of venous system makes the vein segment experience increased wall shear stress and intr amural wall tension. These forces act together with the biologic injury response to modulate the morphologic changes of the implanted vein graft (Dobrin et al., 1989; Fillin ger et al., 1994; Galt et al., 1993; Zwolak et al., 1987). Because of the complexity of the remodeling process, there is a need for a mathematical model to better understand the role of the various factors involved. Toward that goal, a mathematical model is developed in this chapter that will account for the effect of wall shear stre ss. Because the intima is the wall layer that is directly exposed to blood flow, the change of its thi ckness will be monitored as a function of time and shear stress, and the model will be derived for that layer. Reduced shear forces lead to accelerat ed intimal hyperplasia development and reduction in lumen size, while increased shear retards intimal thickening and leads to an expansion in vein graft diameter. A quantitative model of these conceptualized relationships between shear and hyperplasia will help gain insight into the dynamics of the remodeling process. To develop and vali date a mathematical model relating the dynamic changes in intimal thickening to the imposed shear stress during early vein graft adaptation, a set of experimental data obtained from the vascular biology laboratory (Dr. Scott A. Berce li) at the University of Fl orda is used. The experimental data are from a bilateral rabbit carotid vein graft model which was developed and validated in the vascular biology laboratory (Jiang et al., 2004) This model induces a 10fold difference in flow rate between the right and left vein grafts (Figure 31). PAGE 42 42 Morphology data obtained from grafts exposed to these high and low shear conditions were used to create the data matrix (shear st ress, implantation time, graft diameter and intimal thickness) for the model developmen t. Performance of the model and its utility in nonlinear geometries are te sted by applying the model to the complex flow environment of a focal stenosis in a rabbit vein graft. Figure 31. Schematic of vein grafts. A. Schem atic of bilateral vein graft. B. Schematic of focal stenosis. Model Development Intimal Thickening in Rabbit Model Data used in the development of the model were obtained from the work by Jiang et al. (2004). A brief description of the ex periment is provided her e. Thirtyfour New Zealand White male rabbits (3.03.5 kg) under went bilateral carotid vein grafting, and differential flow conditions were created by di stal branch ligation as described previously (Fernandez et al., 2004; Jiang et al., 2004). An anastomotic cuff technique was used to PAGE 43 43 implant bilateral external jugular vein segmen ts into the common carotid arteries. A flow differential between grafts was created by ligati ng the internal carotid artery and three of the four primary branches of t he external carotid artery, re sulting in an approximate 10fold difference in mean flow rates. At 1 (n=8), 3 (n=7), 7 (n=8), 14 (n=5 ), and 28 (n=6) days after implantation, dynamic flow rates (Transonic Systems, mo del T106, probe 2SB) were recorded, and the vein grafts were harvested for morphom etric analysis. In the parentheses above, n represents the number of experiments performed. Exte rnal vein graft radius was measured using a video camera (Sony, m odel DXC151A) attached to a dissecting microscope. Formalinfixed graft segment s (n=34) were embedded in paraffin and 5 m histologic cross sections stained with Ma ssons trichrome and van Giesons elastin stains. Digitized images were collected and analyzed (Zeiss Imaging, AxioVision v4.4) to measure crosssectional area of each layer of the graft wall. The in vivo radius of each layer was then calculated from the me asured crosssectional areas using the measured external vein graft r adius. It was assumed that the vein graft was originally in a cylindrical shape and the ar ea of each layer was conserved. The shear stress, was calculated from the measured flow ra te and assuming a Poiseuille flow ( = 4 Q/ r3; where is the blood viscosity (0.035 poise), Q is the mean flow rate, and r is the lumen radius). Predictive Model of Intimal Growth The intimal thickness data obtained from t he bilateral vein graft experiment as a function of time and shear stress are shown in Figure 32. At Day 1 and 3, minimal intimal thickening is identified beyond the single layer of endothelium and basement PAGE 44 44 membrane that characterizes the normal v enous architecture. Intimal growth is identified at Day 7, 14 and 28, and is consistent with previous experimental observations (Meyerson et al., 2001; Schwartz et al., 1992; Zwolak et al., 1987). The data show an inverse relationship with the imposed shear. As a consequence, the intimal thickness is assumed to be of the form: h h0 Rlumen[1 e A ( t t *)] 1 BC t t (31) where h0 is the intimal thickness at implantation, Rlumen is the initial lumen radius, is shear stress, t is time, t* is the time when the onset of intima thickness change occurs, and A, B and C are constants to be determined from experimental data. The units used for length, shear stress and time are m, dynes/cm2 and day, respectively. An important feature of this model is the use of a decaying ex ponential term to describe the initial burst and t hen the gradual reduction in the rate of intimal growth as a function of time. Noting a period of quiescence prior to t he start of active gr owth, a latency period (t*) is incorporated to model this behavio r. Because a direct inverse sheardependent relationship will lead to an infinite intimal th ickness at zero shear stress, a form (1+ B C) is used. At infinite time and zero shear stress, the maxi mum inward growth is bounded by the value of the initial lumen radius, wh ich corresponds to a complete occlusion of the vessel. Coefficient A is a measure of t he rapidity at which intimal growth reaches steady state. Coefficients B and C r epresent the proportionality and power characteristics to the shear stress term, respectively. PAGE 45 45 Figure 32. Intimal thickness as a function of shear stress at each measurement day. Measured intimal thickness data for both low (filled circle) and high (open circle) flow conditions are shown. So lid curves represent model curves. Experimental data are fr om Jiang et al. (2004). PAGE 46 46 Equation 31 is used to fit all the intima l thickness data points from Day 1 to Day 28 (Figure 32). The coefficients A, B and C are obtained by minimi zing the sum of the squared differences between the data and t he fitted values, using the nonlinear regression function in MATLAB (Version 7.0, The MathWorks Inc.). They are found to be equal to 0.170 0.221, 10.9 5.05, and 0.523 0.351, with the associated 95% confidence intervals, respectively. Fr om the experimental morphology data, no response to shear is seen in the intimal thickness before 3 days so t* was assigned a value of 3 days. The initial lumen radius,lumenR and intimal thickness, oh are found from specimens collected at the time of vein gr aft implantation to have an average value of 1402 m and 4 m, respectively. Figure 33A show s a surface plot of the intimal thickness as a function of shear stress and time as predicted by Equation (31). It is found that 77% of the expe rimental data points are lo cated between the upper and lower bounds of the 95% confidence interval of the model surface. Figure 33. Surface plots of intimal thi ckness (A) and intimal thickening rate (B). PAGE 47 47 The dynamic nature of the intimal growth in response to shear stress can be further explored through differ entiation of Equation (31): C t t A lumenB Ae R dt dh 1*) (, t t (32) Figure 33B provides a surface plot of t he rate of change of the intimal thickness as a function of time and shear. A rapi d intimal growth is seen under low shear conditions during the first week following implantat ion. A significant dr op in the rate of growth as time progresses or when the s hear stress increases is also observed. Model Application Focal Vein Graft Stenosis Construct In order to assess the intimal growth m odel, it is applied to a rabbit focal vein graft stenosis construct. Bilateral carotid inte rposition vein grafts with partial distal branch ligation were performed by the member s of the vascular biol ogy laboratory (Dr. Scott A. Berceli) as described in Model Deve lopment Section. Prio r to restoration of flow, a 3.0 French (1.0 mm) polyethylene mandrel is placed external to the graft, and a Midgraft focal stenosis is created by placement of a single 80 silk ligature, resulting in an approximate 80% stenosis of the lumen crosssectional area. Flow measurement and video morphometry are recorded at the ti me of implantation (day 0) and 28 days later. At 28 days following implantation, vein grafts are surgically exposed and perfusion fixed with 10% buffe red formalin at 80 mm Hg. Grafts are paraffin embedded and longitudinally sectioned at a thickness of 5 m. Every tenth section is collected, stained with Masons trichrome, and analyzed for morphology. Simulation of Intimal Thicken ing around a Focal Stenosis PAGE 48 48 Due to the complex geomet ry around the focal stenosis the shear stress along the length of the graft is det ermined numerically using a commercially available finite element analysis software package ADIN A (Version 8.3, ADINA R&D Inc.). Calculations of the local shear stress and estimation of graft remodeling are performed in an iterative fashion where the shape of the wall is tracked as a function of time, as described in Figure 34. Figure 34. Iteration procedure for si mulation of intimal thickening. The profile of the radius of the internal elastic lamina along the axial direction is reconstructed using a smoothing spline al gorithm (MATLAB) through the measured data points. The radius of the internal elastic lamina at Day 28 is assumed to be the lumen radius at Day 0 where the intimal thickne ss is negligible. Though there must be some degree of diameter change, it is assumed t hat the diameter change is minimal due to PAGE 49 49 the ligature at the focal stenosis. The wall s hear stress along the length of the graft is obtained from the steady state axisymmetric computation of blood flow in the vessel using ADINA. Blood is assumed to be Newt onian, and incompressible, and a no slip boundary condition is applied at the lum en wall. The comput ational domain is constructed such that the total axial length fr om the inlet to outlet is 20 times the inlet diameter with the focal stenos is located in the middle. T he length of the graft from proximal to distal cuffs is appr oximately 3cm. As a first order approximation, the effect of the cuffs is neglected. Althoug h the cuffs are going to affect the flow environment and model coefficient values used, they shoul d not have an effect on the general method presented here. Based on the experimental measurements of the flow rate and lumen radius at the time of implant ation, a parabolic velocity prof ile with a center velocity of 16.2 cm/s and 4.7 cm/s, for high and low flow grafts respectively, is applied at the inlet at time zero. The center velocities are calculated by assuming a Poiseuille flow (ucenter = 2Q/ r2, where Q is flow rate and r is lumen radi us.). The values of the calculated wall shear stress are entered into Equation (32 ) to estimate the intimal growth and determine a new lumen contour. This updated gr aft geometry is then used to calculate the new wall shear stress. With the most dynamic changes in graft geometry occurring early following implantation, a time step interv al of 0.1 day is used from Day 3 to Day 10, followed by use of a 1 day interval between Day 10 and Day 28. As defined in the model (Equation (31)), it is assumed that no change in intimal thickness occurs prior to Day 3, i.e., *3 t Figure 35 shows the pathlines and wall shear stress profiles at Day 28 for stenotic grafts exposed to low and high flow c onditions. It also shows the shear stress PAGE 50 50 gradient along the vein graft. Peak values in shear and gradient are observed at the focal stenosis, and recirculation zones just downstream to t hem. That zone of recirculation is getting larger as the shear is increased. Figure 35. Hemodynamic characteristi cs inside the focal stenosis at Day 28. The corresponding experimentally measur ed and predicted graft geometries at Day 28 for each flow simulation are provided in Figure 36. Figur e 36A provides a schematic of the focal stenosis experiment. It also shows histology images at two different locations, including t he one at the center of the stenosis. The empirical model predicts enhanced intimal growth in the area around the constriction (Figures 36B and 36C). This is observed experimentally for both flow conditions, but is most dramatic under low flow. PAGE 51 51 Figure 36. Comparison of focal stenosis simulation results with experimental data. Schematic of focal stenosis exper iment (A). Measured and simulated morphology data for low (B) and high (C) flow conditions. Experimental data are from the vascular biology labo ratory (Dr. Scott A. Berceli). PAGE 52 52 Of particular importance in validating the presented approach is the relative fidelity between predictive and ex perimental intimal growth at various locations along the graft. The low flow graft data (Figure 36B ) generally demonstrate a good agreement between the experim ental and predictive result s, with a modest trend towards overestimation of the magnitude of intimal thickening in the downstream region. In the high flow graft case (Figur e 36C), reasonable agr eements proximal and distal to the stenosis are observed. Howe ver, marked difference near the stenosis is seen. Significant intimal thickening is exper imentally observed in this high shear region, a site where the model suggests a limited hy perplasia response. This is not too surprising since the model takes into consid eration only the effect of shear stress. Parameter Sensitivity Analysis To examine the sensitivity of the simulated intimal gr owth data within the focal stenosis construct to the empirical parameters in Equation (32), low flow data at day 28 are used. Equation (32) is us ed because it reflects the rate of change of the thickness. It is more practical than Equation (31) sinc e it does not depend on the knowledge of an initial thickness value. Each coefficient (A B, and C) is varied 10%, and its effect on the simulation is analyzed. For axial positions ranging from 500 m to 1500 m with respect to the center of the stenosis, T able 31 shows the percent change of intimal growth associated with the variation of each coe fficient. It is found that coefficient A has a limited influence on the value of the intima l thickness due to the exponential growth and rapid rise to steady state at Day 28. C oefficient B has the greatest overall effect, with average intimal thickness changes of 8.5% and 10.3% following a 10% increase and decrease, respectively. Coefficient C amplifies any local shear stress changes, and PAGE 53 53 provides the largest increase in intimal th ickness (24.1%) at 0. 06mm upstream to the stenosis when C is decreased by 10%. Table 31. Parameter Sensitivity Analysis % Change of each coefficient in Eq. (32) % Change of intimal growth (MeanSD) Coefficient A Coefficient B Coefficient C +10 10 +10 10 +10 10 0.50.05 0.80.07 8.50.6 10.30.9 3.64.0 3.94.8 Discussion In this chapter, a mathematical model of vein graft intimal hyperplasia has been developed based on experimental data obtained fr om a rabbit vein graft construct exposed to a wide range of physiologic shea r stresses. It is understood that shear stress is not the only factor leading to intima l hyperplasia. Intimal hyperplasia is a very complex process involving a large number of biologic and hemody namic factors. Nevertheless, a general approach has been dev eloped to model vessel thickening due to shear. Additional factors can be included following a similar approach in order to get a more comprehensive vessel remodeling model. The model utilizes a nonlinear inverse relationship between the intimal thickening rate and shear stress, which is consistent with the ex perimental data of our group (Fernandez et al., 2004) and others (Mey erson et al., 2001; Schwartz et al., 1992). It successfully describes the time dependent behavior of early vein graft remodeling (Zwolak et al., 1987). Zhang et al. (2002) found t hat the neointimal hyperplasia develops rapidly, and stabilizes by about 4 weeks postoperatively in a PAGE 54 54 murine vein graft model. Zou et al. (1998) observed increasing intimal thickness up to 16 weeks in a mouse model, with the rate of thickening decreasing after 4 weeks. The present study confirms that the rate of change of intimal thickness depends not only on shear stress, but also on time. Also unique to the model is the use of a time delay prior to the initiation of the hyperplasia response. From the experimental data, this time delay is found to be 3 days. For other cl inical conditions or diseases, this time delay is expected to be different. This time delay is not known a priori and needs to be determined experimentally or treated in t he model as a variable. While the underlying mechanisms for this time delay may include mu ltiple factors, it appears to be related to the development of microvasculature gr owth from the surrounding tissue (Westerband et al., 2001). For example, it has been report ed that the smooth mu scle cell migration from media to intima starts 4 days after injury (Lemson et al., 2000). Also critical may be the initial disequilibrium between cell pr oliferation and cell death, with several investigations confirming an early burst of apoptosis following vein graft implantation (Jiang et al., 2007; Westerband et al., 2001). As such, the inco rporation of biologicallydefined parameters into the pr edictive model may provide a critical insight into the remodeling process. Among the strength of the current model is the us e of an initial experimental construct to define the model parameters, and application of the model to the simulation of a second, independent experimental cons truct, the focal stenosis geometry. The secondary analysis revealed a differenc e between experimental data and model predictions at the throat of the stenosis where shear stress is higher (Figure 36C). It suggests that the elevations in shear do not universally yield reductions in the rate of PAGE 55 55 intimal growth, instead suggesting elevat ed shear leads to enhancement of the hyperplasia response. It is not clear either if it is shear stress or its gradient that is responsible for this enhanced intimal thickeni ng at the maximum constriction (Figure 35). Also a limitation of the cu rrent model is the focus on intimal growth as the dominant mechanisms for changes in lumen geometry. Work by our group (Jiang et al., 2007) and others (Glagov et al., 1987; Korshunov and Berk, 2004; Owens et al., 2006) have established the early proc ess of outward remodeling is critical in maintaining lumen patency. The shear stress data used in this study (Figure 32) also show that there was outward remodeling at Day 28 in t he high flow construct case. Experimental constructs designed specifically to eval uate outward remodeling and incorporation to our predictive model could be a futu re direction of this work. The process of early adaptation of vein graf t to arterial environment also involves medial thickening, though it has been reported t hat the medial thickening has little or no relation to shear stress, and correlates more with circumferential strain (Dobrin et al., 1989; Fernandez et al., 2004; Galt et al., 1993). It is evident that other factors than shear (e.g. tension and biological environment ) need to be taken into account in order to fully describe the thickening process. PAGE 56 56 CHAPTER 4 MULTILAYER MATHEMATICAL MODE L OF VEIN GRAFT REMODELING Though the relationship between wall shear stress and the remodeling of vein graft intimal layer is described in the previous chapter, the lumen reduction in vein graft adaptation is the result of t he combination of multiple la yers remodeling. Also, as described in chapter 2, wall tension as well as wall shear stress is another important hemodynamic factor influencing vein graft remodeling. Although numerous experimental evidences relating hemodynamic forces to vascular remodeling have been reported, a comprehensive mathematical model describ ing the relationship among the wall shear stress, wall tension, and the remodeling of eac h individual wall layer of vein graft is lacking. One of the characteristics of vein graft implantation is the injury that occurs due to the altered hemodynamic environment the vein graft experiences. In response to the damage, repair events such as leukocyte recruitment, reendothel ialization, smooth muscle cell proliferation, and matrix deposition occur in a time frame ranging from hours to months. The imposed hemodynamics has diffe rent effects on vein graft adaptation depending on the time frame. This temporal dependence between the local hemodynamic environment and the response of the vascular wall is frequently lacking from models that have been developed to qu antitatively examine the morphologic changes of vein graft. In this chapter, a multilayer mathem atical model which describes timedependent behaviors of multiple layers of vein graft is developed based on experimental data from the wellcha racterized rabbit model. Founded on the model developed in the previous chapter, which focused on the infl uence of shear stress on intimal adaptation, PAGE 57 57 the current model examines the changes in in timal and medial areas and the radius of external elastic lamina (EEL). Also unique in the current model is the integration of a logistic relationship between independent and dependent variables to describe the early acceleration and later reduction in the rate of adaptation. The det ailed understanding of the temporal changes in vein graft morphology that can be extract ed from the current model is critical in identifying the dominant contributors to vein graft failure and developing strategies to improve their longevity. Model Development Remodeling Characteristics of Vein Graft Wall Experimental data of rabbit bilateral ve in graft construct used in the model development were obtained from the vascular biology laborator y (Dr. Scott A. Berceli), as described in Chapter 3 (Jiang et al., 2004; Jiang et al., 2009). In the previous chapter, the mathematic al model of intimal layer remodeling (TranSonTay et al., 2008) was based on ex perimental data up to 28 days. In the present development, to accommodate additiona l wall layers and additional data at 3 month and 6 month time point, a logistic m odel is used to describe the timedependent adaptation of the intimal ar ea, medial area, and EEL radius (Figure 41). The experimental data of wall thickening (Figures 42, 43, 45) show an initial accelerated growth and then damping once the wall thickne ss has reached a critical value, which is the hallmark of most growth models. T he Sshaped behavior para llels the population growth of smooth muscle cells. For that r eason, a logistic model, which is a common and one of the simplest models of population gr owth, is used to describe the observed Sshaped behavior. A logistic model with t he following general fo rm was used (Britton, 2003): PAGE 58 58 b N aN dt dN 1 (41) where N is the variable of interest, t is time, and a and b are constants representing the growth rate and carrying capacity of N, respectively. The logistic relationship was appropriately modified to model the timedependent changes in intimal area, medial area, and EEL radius. To minimize the influence of inherent variations in initia l vein geometries, each outcome variable was normalized to the corresponding value at the time of implantation. The exception is intimal area, which was normalized with respect to the medi al area due to the inherent inaccuracies in measuring the crosssectional area of this singlecell layer. Figure 41. Schematic of the crosssection of venous wall. Intima Model Experimental data detailing the effect of shear stress and implantation time on the normalized intimal area is shown in Figure 42. PAGE 59 59 Figure 42. Normalized intimal area, medial area, and EEL radius with respect to shear stress at each time point. Experimental data are from Jiang et al. (2004). PAGE 60 60 Consistent with the finding of other investigators (Meyerson et al., 2001; Schwartz et al., 1992), an inverse rela tionship between the intimal area and shear stress is observed. Due to outward remodel ing of the graft and an increase in lumen diameter, a reduction in shear stress is init ially observed at 28 days after implantation and continues throughout the 180 day observa tion period. To incorporate these characteristics, the logistic model for t he growth of the intima is as follows. ) ( 1 ) ( ) (* * i I I I Ib t A t A a t t A x x x, with Id I I I it D c b b )) ( ( 1 x (42) where I A is the normalized intimal area, is the shear stress, DI =(R/R0)3, R is the lumen radius at time t, R0 is lumen radius at implantation, x is a position vector, and Ia Ib Ic and Id are constants. The coefficient Ia represents the slope of the rate of intimal area growth. The coefficient Ib provides the limit of intimal area growth for the shear stress. The exponents Ic and Id represent the dependence on the shear stress. Note that ib represents the upper limit that can be obtained by the intimal area, and incorporates the inverse relationship with sh ear stress that was pr eviously described. Integration of Equation (42) yields: t a I I i t a I i II Ie A A b e A b t A ) ( ) ( ) ( ) (* 0 0 0 *x x x x (43) The coefficients Ia Ib Ic and Id were obtained through nonlinear regression (MATLAB v7.7, The MathWorks Inc., Natick, MA) against the experimental data (Figure 42) and were found to be 0.350.06, 8. 00.9, 0.0130.035, and 2.61.4 with 95% confidence intervals, respectively. 10 A is the normalized intimal area at implantation and PAGE 61 61 was experimentally measured to be 0.186. T he units used for shear stress and time are dynes/cm2 and day, respectively. Wall thickness depends on both cell prolifer ation and vein graft diameter change while wall area changes only by cell proliferat ion. By using wall area as the modeling variable, the effects of cell proliferati on and graft diameter change can be separated. Media Model Experimentally measured medial areas, norma lized to the medial area at the time of implantation, are illustrat ed in Figures 42 and 43. Time (months) 0123456 Normalized Medial Area 0 2 4 6 8 10 Exp. (Error Bar=SD) Model Figure 43. Normalized medial area as a f unction of time. Experim ental data are from Jiang et al. (2004). Unlike the notable influence of shear stress on intimal hyperplasia, the process of medial thickening appears independent of the imposed wall shear stress. Instead, previous investigators have shown intramural wall tens ion to have an important dependence on medial thickening (Schwartz et al., 1992; Zwolak et al., 1987). Although the experimental design was not developed to directly evaluate the influence of wall PAGE 62 62 tension on vascular remodeling, examination of our data s upports a positive correlation between the rate of medial hypertrophy and the imposed wall tensi on (Figure 44). Figure 44. Rate of normalized medial ar ea change as a function of wall tension. Experimental data are fr om Jiang et al. (2004). Using this as a starting point, a linear correlation is assumed between the rate of medial area growth and wall tension, m odifying the logistic model as follows. ) ( ) ( ) ( ) (* *t A b t T D t A a t t AM M c M M M MMx x x x (44) where M A is the normalized medial area, DM=R0/R, R is the lumen radius at time t, R0 is the lumen radius at impl antation, T is the circumfe rential wall tensile stress, x is a position vector, and Ma, Mb, and Mc are constants. 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.0E+002.0E+064.0E+066.0E+068.0E+06 Tension (dynes/cm2)Rate of Normalized Medial Area Change (/day) PAGE 63 63 The coefficient Ma represents the slope of the rate of medial area growth. The coefficient Mb gives the limit of medial area grow th for the wall tensile stress. The exponent Mc is a function of the tensile stress. In order to obtain the coe fficients in Equation (44), the rate of change of the normalized medial area (dA*M/dt) is needed. The equation MA=AB eCt (t is time, and A, B, and C are constants) is fitted to the exper imental data in Figure 43 to obtain the rate of change (dA*M/dt=BCeCt). Knowing that rate, the values for Ma, Mb, and Mc in Equation (44) are obtained thro ugh a nonlinear regression of the logistic model using the experimental data for MA, DM, and T where the wall tension stress is calculated from the relation T=PR/h (T: tension per unit area, P: blood pressure, R:lumen radius, and h: medial thickness) Using this approach, the coefficients aM, bM, and cM in Equation (44) were found to be 0.0140.175, 0.798.33, and 0.302.03 with 95% confidence intervals, respectively. The units used for blood pressure and time are 104 dynes/cm2 and day, respectively. Figure 43 shows the model curve obtained by numerically integrating Equation (44). EEL Model EEL radii were determined experimentally over a 6 month period and are shown in Figures 42 and 45. Simila r to our observations related to adaptation of the media, changes in the graft radius were relatively independent of wall shear stress (Figure 42). Hence, the logistic model fo r the normalized EEL radius (R* EEL) utilizes time as the only dependent variable and is provided by the following expression PAGE 64 64 E EEL EEL E EELb t R t R a t t R ) ( 1 1 ) ( ) (* *x x x (45) where t is time, x is a position vector, and aE and bE are constants. The coefficient aE represents the slope of the rate of EEL radius increase. The coefficient bE is the limit of EEL radius increase. Figure 45. Normalized EEL radius as a func tion of time. Experimental data are from Jiang et al. (2004) The analytical solution to Equation (45) is t b b a EEL EEL E t b b a EEL E EEL E EELE E E E E Ee R R b e R b R b t R) 1 ( 0 0 ) 1 ( 0 0 *) 1 ) ( ( ) ( ) 1 ) ( ( ) ( ) ( x x x x x (46) and the coefficients Ea, Eb, and 0 EELR are found to be 0.1330.092 and 1.940.19, and 1.020.04 with 95% confidence intervals, respectively, through nonlinear regression using the shearaveraged experimental data in Figure 45. Note that coefficients Ea PAGE 65 65 and Eb are proportionality constants and t hey define the maximum limit of EELR. The unit used for time is day. The model curv e describing the normalized EEL radius as a function of time is shown in Figure 45. Figure 42 shows that there is no c onsistent correlation between WSS and EEL remodeling. This is somewhat surprising co nsidering that it has been suspected that vein graft expands to reduce the WSS imposed by faster arterial flow. Wall tension does not seem to be the driving force of EEL remodeling either. EEL expansion tends to increase lumen radius which increases wall tension. EEL remodeling seems to be related to the biologic environm ent of the arterial system wh ich is represented by time dependent terms in t he model equation. Figure 46 shows the surface plots of the normalized EEL radius as a function of time, the normalized intimal area inversely re lated to shear stress, and the normalized medial area positively related to tensi on and reaching plateau after about 3 months. Sensitivity and Stochastic Analysis To examine the effects of the variation of each coeffi cient on the calculation of the change of the area or radius in Equations (42), (4 4), and (45), each coefficient was varied 100% and 50%, and the model predi ctions at 6 months were obtained (Table 41). In the intima model (Equation (42)), bI has the greatest impact where the % change of the coefficient was directly refl ected on the model prediction. The other coefficients have minimal impacts. In the medi a model (Equation (44)) the coefficients Mb and Mc have the greatest impacts on the m odel prediction. In the EEL model (Equation (45)), Eb has direct impact. However, for a ll the coefficients, the % change of the model prediction was within the range of the % change of the coefficient. PAGE 66 66 Figure 46. Surface plots of normalized EEL r adius (A), normalized intimal area (B), and normalized medial area (C). PAGE 67 67 Table 41. Parameter Sensitivity Analysis Coefficient % Change of Coefficient % Change of Model Prediction at 6 months aI 100 50 0 0 bI 100 50 100 50 cI 100 50 6 4 dI 100 50 24 4 aM 100 50 0 0 bM 100 50 80 44 cM 100 50 102 36 aE 100 50 0 6 bE 100 50 100 48 In order to examine the robustness of the model to variations in the experimental data, a stochastic analysis was performed. Using the mean and st andard deviation (SD) of the experimental data at each time point, 104 sets of new data (consisting of intima, media, and EEL) were generated assuming a no rmal distribution at each data point. In the absence of a standard deviation for t he shearand timedependent intimal area experimental data (Figure 42), these values were assumed to have a coefficient of variation (CV) equal to that of the medial area (CV=0.39) Using the experimentally measured flow rate, blood pressure and initial vein graft geometry, Equations (42), (44), and (45) were used to update the vein graft geometry at every 0.5 day. Because each experimental data point represents an indi vidual rabbit, separate simulations were PAGE 68 68 performed for each harvest day using the av erage of the initial conditions for that specific day. To parallel the experim ental design, simulations were performed separately for the low and high flow groups of the vein grafts. Figure 47 shows both experimental data and model predictions for these stochastic simulations. Figure 47. Stochastic analysis of the vein graft mathematic al model. Experimental data are from Jiang et al. (2004). In the cases of the EEL and lumen radii, the mean va lues of the experimental data are located within one standar d deviation of the model predictions after 28 days. In the cases of the intima and media, some in creased deviations at the 6 month time point are observed. However, for all the data point s, the SDs of the m odel predictions are comparable to those of the ex perimental data, confirming t he robustness of the model. Model Application Simulation of Idealized Stenosis The model described in the previous sectio ns is now applied to the simulation of wall remodeling in an idealized stenosis to test the applicability of the model. Figure 48 shows a schematic of the stenosis geometry and the equations describing dynamic PAGE 69 69 remodeling within each layer of the vein gra ft wall. A smooth geometry within the throat region was modeled using a sinusoidal curve. Figure 48. Schematic of stenosis simulation. The initial wall shear stress was esti mated using the ADINA finite element analysis software package (Version 8.5, AD INA R&D Inc., Watertown, MA), assuming axial symmetry, a Newtonian fluid, and a nonslip boundary conditi on at the wall. Supported by our histological examination of the developing intima, which demonstrates a disorganized extracellular ma trix structure, we have assumed this layer to be a nonload bearing structure in t he analysis schema, and the re sulting wall tension is determined by the relationship T=PR/h (T: tensile stress, P: blood pressure, R:lumen radius, and h: medial thickness). Using wa ll shear and tension as input parameters for Equations (42), (44), and (4 5), a new wall geometry (i.e intimal and medial thickness and EEL radius) was determined. Updated wa ll geometries were determined at 1 day PAGE 70 70 intervals for the duration of the simulation as follows. At the time interval, ti, the geometry of the vessel, WSS and pressure were computed using the mesh and values generated at ti1. This coupling between computational fluid dynamics (CFD) and growth model is in line with the one used by Figueroa et al. (2009). In order to evaluate the influence of the 1 day time interval on the numer ical integration of E quations (42), (44), and (45), an analytical solution to these equations was obtained. Comparison between the forward Euler time integr ation using a 1 day interval and the analytic solutions demonstrated a less than 0.002% difference for all expressions at Day 180, using shear and tension values in the physiologic range. Analyses we re performed assuming an initial 50%, 70%, and 80% stenosis (diame ter) subjected to low and high flow conditions. To parallel our experimental rabbi t vein graft model, the initial inlet lumen radius was set to 1.5 mm. Inlet parabolic pr ofiles with a centerline ve locity of 2.45 cm/s and 12.6 cm/s were used for the low and high flow conditions, respectively. The total axial length of the geometry wa s set to be equal to 80 times the inlet diameter to ensure no reverse flow at the exit. Figure 49 shows the intimal area along t he graft at 1 week, 2 weeks, 3 weeks, and 6 months after implantation under low fl ow conditions. Intimal thickening is observed along the graft except in the regi on of the stenosis, where shear stress is maximal. The intimal area increases rapi dly and reaches an asymptote at about 1 month after implantation. In association with an initial 50 % stenosis (Figure 49A), a localized increase in intimal thickening is observed between axial positions 0.1 and 0.3 cm and corresponds to the flow separation and reattachment distal to the stenosis. In association with the 70% and 80% stenoses (Fi gures 49B and 49C), the site of flow PAGE 71 71 reattachment is located at 1.1 and 1.5 cm respectively, leading to a broad region of recirculation and a mild increase in intimal area in these regions. Figure 410 shows the intimal area al ong the graft under high flow conditions. Within regions of uniform/unidirectional flow proximal and distal to the stenosis, high flow conditions result in a reduction in intima l thickening. In contra st to the low flow simulations, the marked reductions in shear distal to the stenosis lead to a marked augmentation in intimal area in this region. Notable are the 70% and 80% high flow simulations where the reduced shear within the region of flow separation has shifted distally, resulting in enhanced intimal thickening that is several centimeters distal to the site of maximum stenosis. Figure 411 illustrates the change in the EEL radius and medial area associated with a 50% stenosis exposed to low flow co nditions. A uniform increase in both graft radius and medial area is observed pr oximal and distal to the stenosis. In contrast to the effect on intimal th ickening, the perturbations in the flow patterns induced by the stenosis have a limit ed effect on these geom etric parameters. The patterns of graft expansion and medi al growth associated with 70% and 80% stenoses demonstrated a similar asymptot ic behavior and were not notably influenced by the magnitude of the fl ow (data not shown). One of the limitations of the stenosis study is that the shear stress was obtained from steady flow simulations. However, it is expected that under pulsatile flow the recirculation zone behind the stenosis would o scillate, which would result in different prediction of wall thickening in those region s. Thus, the intimal thickening around the reattachment point shown in Figur e 410 could be an overprediction. PAGE 72 72 A B C Figure 49. Simulation results for intimal ar ea around the stenosis of 50% (A), 70% (B), 80% (C) reduction in stenosis (diame ter) under low flow condition. (50% Stenosis, Low Flow) 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 10.500.511.522.53 Axial Position (cm)Intimal Area (cm^2) (70% Stenosis, Low Flow) 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 10.500.511.522.53 Axial Position (cm)Intimal Area (cm^2) (80% Stenosis, Low Flow) 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 10.500.511.522.53 Axial Position (cm)Intimal Area (cm^2) 6 months 6 months 6 months time time time Re = 53 Re = 82 Re = 116 PAGE 73 73 A B C Figure 410. Simulation resu lts for intimal area around the stenosis for 50% (A), 70% (B), 80% (C) reduction in stenosis (d iameter) under high flow condition. (50% Stenosis, High Flow) 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 3113579 Axial Position (cm)Intimal Area (cm^2) (70% Stenosis, High Flow) 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 3113579 Axial Position (cm)Intimal Area (cm^2) (80% Stenosis, High Flow) 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 3113579 Axial Position (cm)Intimal Area (cm^2) 6 months 6 months 6 months Re = 238 Re = 366 Re = 528 PAGE 74 A B Figure 4 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 EEL Radius (cm) 0 0.002 0.004 0.006 0.008 0.01 0.012Medial Area (cm^2) 6 months 6 months 4 11. Simul a under lo w 10.8 10.8 a tion result w flow con d 0.6 0.6 s for EEL r d ition. (50 % 0.40.2 A (50 % 0.40. 2 74 adius (A) a % Stenosis, Lo w 0 A xial Position ( c % Stenosis, Low 2 0 A xial Position ( c Re = 53 Re = 53 nd medial a w Flow) 0.20 c m) Flow) 0.20 c m) a rea (B) fo r .40.6 .40.6 r 50% sten o 0.8 0.8 time time o sis 1 1 PAGE 75 75 Simulation of Human Vein Graft The present mathematical model was appl ied to the simulation of human vein graft remodeling, encompassing the geometric complexities wh ich are inherent in these conduits. Figure 412A illustrates a 14 cm se gment of human saphenous vein whose structure was extracted from a computer tomographic (CT) image obtain several days after implantation. An approach similar to that described above for the idealized stenosis was used to simulate the human vein graft wall remodeling. Assuming a blood flow rate of 100 cm/sec and an in traluminal pressure of 1.3x104 Pa, the resulting wall shear and tension served as the input parameter s for Equations (42), (44) and (45) to determine the interval change in intimal area, medial area, and EEL radius. Alternately calculating the local hemody namics and graft geometry at 1 day intervals yielded a temporal understanding of vein graft remodeling through the in itial 6 months following implantation. Figure 412B shows the intimal area, medi al area, and EEL radius along the graft at 1 week, 1 month, 3 months and 6 months. The remodeli ng patterns observed in the human vein graft approximate those observ ed in the idealized stenosis, with intimal thickening that was minimal within the thr oat of the stenosis and reached a plateau at about 1 month. Medial thickening is also mi nimal within the stenosis, with a growth rate that decreases with time. Discussion There have been experimental evidences that relate each individual hemodynamic factor to the remodeling of each wa ll layer of the vein graft, and a set of mathematical equations which define these relationships has been developed in this chapter. Novel in the present approach is the development of a comprehensive PAGE 76 76 A B Axial Position (cm) 02468101214 Intimal Area (cm2) 0.00 0.01 0.02 0.03 0.04 0.05 Axial Position (cm) 02468101214 Medial Area (cm2) 0.00 0.01 0.02 0.03 0.04 0.05 Axial Position (cm) 02468101214 EEL Radius (cm) 0.0 0.2 0.4 0.6 0.8 1.0 Figure 412. Simulation results for human vein graft. A. Vein graft geometry and mesh at implantation. B. Simulation re sults at different time points. mathematical model that includes three layers for vein graft remodeling, in contrast to the less physiological single layer models pr eviously described in the literature (Budu 6 months 3 months 1 month 1 week 6 months 3 months 1 month 1 week 6 months 3 months 1 month 1 week PAGE 77 77 Grajdeanu et al., 2008; Friedman et al., 1986; TranSonTay et al., 2008; Zohdi, 2005). Also unique to the current formulation is the use of a logistic model as the base expression. As a simple equation that descr ibes the Sshaped behavior observed in the experimental data (Figures 42, 43, and 45), this approach has been used by many researchers for modeling cell population gr owth (Britton, 2003; Fujikawa et al., 2004; Goudar et al., 2005), but is novel in its application to soft tissue remodeling. In the formulation of the model, changes in intimal and medial crosssectional area are related to wall shear stress and wall tension, respectively, and EEL radius change is a function of time only. There have been reports that vein graft dilation is dependent on shear stress. Fillinger et al. ( 1994) and Owens et al. (2006) observed a positive correlation between initial shear st ress and percent change in lumen diameter in human vein grafts data. However, the observed shear dependence of lumen size does not depend only on graft dilation but also on wall thickening. Therefore, EEL radius is chosen as a measure of graft dilation to separate the effect of wall thickening. Moreover, Galt et al. (1993), in their r abbit vein graft experiments, observed no significant change in lumen radius of vein gr afts exposed to normal and low flow for 4 weeks. In contrast, intimalm edial thickening was significant ly more pronounced in the low flow case, suggesting a more prominent increase in EEL radius in those grafts exposed to low flow conditions. In the present experimental data, any significant dependence of EEL radius on the s hear stress was not observed. Wall tension has been recognized as one of the hemodynamic factors that has a notable impact on the thickening of the vein gr aft wall (Schwartz et al., 1992; Zwolak et al., 1987). In the present model, wall tension a ffects the medial layer, and the intima is PAGE 78 78 regulated by shear stress. The experimenta l data (Figure 42) show a significant difference in the initial intimal area c hange rates (up to 14 days) between low and high flow vein grafts, even though the experim ents were started under comparable wall tension for both flow groups. The effect of the wall tension on the intima is minimal because the intimal area did not increase much in the high flow vein grafts throughout the 6 month period. The intima model descr ibes an inverse relationship between shear stress and intimal growth rate, which has been observed by many researchers (Jiang el al., 2004; Meyerson et al., 2001; Schwartz et al., 1992). Another char acteristic of the intima model is that it can predict the regression of the intima. Morinaga et al. (1987) reported that the thickened intima in a vein graft exposed to low flow regressed when the vein graft is exposed to higher flow rate in their canine vein graft experiment. In Equation (42), the rate becom es negative when the shear stress is increased such that the ib is less than I A It is understood that hemodynamic forces ar e not the only factors that regulate vein graft remodeling. Incorporation of the bi ologic environment of the arterial system is needed in order to provide a more realistic description of the vein graft behavior. The current model describes vein graft adaptation in the earlier time frame after surgery. Eventually, some vein grafts will become st able after the adaptation and some unstable leading to vein graft failure. The reason wh y some vein grafts are stable and some are unstable is currently unknown. In addition, t he mechanism of vein graft failure seems to be different from that of early adaptation. The ultimate goal of this modeling effort is to develop a model that can predi ct unstable remodeling leading to vein graft failure, so that the model can be used as a pr edictor of vein graft failure. PAGE 79 79 CHAPTER 5 RULEBASED MODEL OF VEI N GRAFT REMODELING In the previous chapters, mathematical models wh ich show the relationship between hemodynamic forces and vein gr aft wall remodeling have been developed. While mathematical models are relatively simple and show the relationship among the variables explicitly, biological mechanisms are often implicitly included in the model. The mathematical models in the previous chapt ers deal with geometrical variables of vein graft remodeling while biol ogical mechanisms are behind the geometrical change. On the other hand, rulebased modeling approach (Hwang et al., 2009) utilizes the biologic knowledge accumulated in the lit erature over the past decades, converts them to rules, and applies the rules to comput ational domain to simulate biological phenomena. In the computational domain, the elem ents can be molecules, cells, organs, or even a group of people. As described in chapter 2, cell level mec hanism of vein graft intimal hyperplasia is mostly smooth muscle cell (SMC) pro liferation and extracellular matrix (ECM) deposition. About 60% to 80% of the volume of intimal hyperplasia is occupied by ECM, and SMC takes up about 20% to 40% of the volume (Kohler et al., 1991; Kraiss et al., 1991; Zwolak et al., 1987; Lemson et al ., 2000). SMC apoptosis and ECM degradation are also involved in the change of the intimal layer (Berceli et al., 2002). The SMC and ECM are reorganized when exposed to arteri al flow, and the reorganization is responsible for the change of graf t diameter (Berceli et al., 2009). The response of vein graft to hemodynam ic forces can be thought of as the sum of the responses of SMCs to the hemody namic forces. SMC has been found to respond to hemodynamic forces (Sumpo et al., 1988a; 1988b). Tissue/organ leve l vein graft wall PAGE 80 80 remodeling results from the cell level acti vities of SMC and ECM, and thus, the rulebased simulation approach is suited for this ty pe of simulation (Boyle et al., 2010). In this chapter, rulebased simulation of hemody namically induced vein graft remodeling is presented. SMC and ECM are two components which occupy the grid elements of rectangular grid system. At each time step, t he probabilities of behaviors are applied to SMC and ECM elements to determine the next states of the elements. Those probabilities were obtained based on experimental data. The simulation is first run on 1dimensional (1D) domain to test the feasibility of t he model. 2dimensional (2D) algorithm is then developed and va lidated against experimental data. The 2D model is applied to focal stenosis simulation and compared with experimental data. Determination of Probabilities of Cellular Activities Probability of Smooth Muscle Cell Division Cell division probability was determi ned based on BrdU (Bromodeoxyuridine) experimental data obtained from the vascular biology laborator y (Dr. Scott A. Berceli). BrdU is a synthetic nucleoside which can subs titute thymidine when DNA is duplicated. In the BrdU experiment, BrdU was injected into vein graft 1 day prior to harvest. It is assumed in this study that BrdU is availa ble for the first 2 hours considering that the halflife of BrdU is on the order of hours (Phuphanich and Levin, 1985). Figure 51 shows a diagram of BrdU injection and cells getting stained. At harvest, the % of BrdU positive cells is det ermined with respect to the total number of cells at harvest. However, during the 1 day from the BrdU injection to vein graft harvest, the total number of cells changes because so me cells divide during the 1 day. Under the assumption that the BrdU is availabl e during the first 2 hours after the BrdU injection, the % of BrdU positiv e cells with respect to the tota l number of cells at the time PAGE 81 81 of BrdU injection would be a more accurate representation of the % of cells dividing. Another thing to note is that BrdU is inco rporated into DNA during Sphase of cell cycle. It is assumed that the cell cycle time is 24 hours, and G1, S, and G2/M phases each takes up onethird of the cell cycle time (Yamamoto et al., 2000) All the factors discussed above put together and based on Figure 51, the space averaged probability of cell division during 1 hour period can be written as follows. ( *, t) = where BU is the % of BrdU positive cells, and is space averaged probability of cell division at 1 day prior to the day when BU is measured, and TB is the time duration (hours) between BrdU injection and tissue har vest, which is 24 hours in the experiment. Figure 51. Schematic of ce lls undergoing cell cycles during 1 day prior to harvest. ( 51 ) PAGE 82 82 Figure 52 shows % of BrdU positive cells measured at harvest. Shear stress is normalized with respect to the preligation shear stress which was measured to be 10 dynes/cm2. The following function is fitt ed to the experimental data. BU = (52) where t is time (day), is normalized shear stress, and the coefficients A and B are found to be 44, 0.16. Figure 52 shows the model curve with the experimental data. Figure 52. % BrdU positive ce lls as a function of normalized shear stress at different time points. Experimental data are from the vascular biology laboratory (Dr. Scott A. Berceli). Figure 53 shows % of BrdU positive ce lls as a function of distance from endothelium. The data are shown for two diffe rent flow conditions. Assuming that the probability of cell division is proportional to the concentration of the molecules released from endothelium and that the c oncentration of the molecule s decreases exponentially, the following form is fitted to the data in Figure 53. (Spatial distribution of cell proliferation) = (53) where x is distanc e from endothelium. When the Equation (53) is fitted to each of the six different data sets (Figure 53) for different flow conditions, coefficient E was found to be 27.39.4 (meanSD) mm1. PAGE 83 83 Day 4Distance from endothelium (mm) 0.00.10.20.30.4BrdU Positive Cells (%) 10 0 10 20 30 40 50 60 High Flow Normal Flow Day 7Distance from Endothelium (mm) 0.00.10.20.30.4BrdU Positive Cells (%) 10 0 10 20 30 40 50 60 High Flow Normal Flow Day 14Distance from Endothelium (mm) 0.00.10.20.30.4BrdU Positive Cells (%) 10 0 10 20 30 40 50 60 High Flow Normal Flow Figure 53. % BrdU positive cells as a f unction of actual distance from endothelium (From Berceli et al. (2002)). PAGE 84 84 Then, the final form of cell division probability can be written as follows. Pdivision ( *, t, x) = (54) where is normalized shear stress, t is time, x is distance from endothelium, and D* is a constant. The constant D* in Equation (54 ) is determined at each time step such that the integration of Equation (54) over t he intimal thickness equals the space averaged probability multiplied by intimal thickness. D* = (55) where IT is intimal thickness (mm). Probability of Smooth Muscle Cell Apoptosis Probability of cell apoptosis is obt ained based on TUNEL (Terminal deoxynucleotidyl transferase dUTP Nick End Labeling) experimental data obtained from the vascular biology laboratory (Dr. Scott A. Berceli). It is assumed that dead cells remain in the tissue for 24 hours before getti ng removed. Then, TUNEL stained cells at harvest died between 1 day prior to harvest and the day of harvest. Using the same correction used in the cell division probability for the total number of cells, the space averaged probability of cell apoptosis for 1 hour can be written as follows. ( *, t) = (56) where TL is the % of TUNEL positive cells, and is space averaged probability of apoptosis at 1 day prior to the day when TL is measured, and TD is time duration (hours) dead cells remain in the tissue before getting removed. TD is assumed in this study to be 24 hours. PAGE 85 85 Figure 54 shows the % of TUNEL positive cells in the same experiment where the BrdU data were obtained. In the Figure 54, the low and high fl ow data are graphed together because the data show no dependence on shear stress. Time (Day) 051015202530 % of TUNEL Positive Cells 0.4 0.2 0.0 0.2 0.4 0.6 0.8 1.0 Experiment Model Figure 54. % of TUNEL positive cells as a function of time. Experimental data are from the vascular biology laborator y (Dr. Scott A. Berceli). The following exponential form is fi tted to the data in Figure 54. TL = (57) where the coefficients F and G are found to be 0.301 and 0.0739. Fi gure 54 shows the model curve with the experimental data. Figure 55 shows the % of TUNEL positive cells as a function of distance from internal elastic lamina (IEL). To determine the spatial distribution of apoptosis, the same procedure used for cell proliferation is used for apoptosis as we ll, and the following exp onential form is fitted to the data (Figure 55). (Spatial distribution of cell apoptosis) = (58) where is distance from IEL. PAGE 86 86 When the Equation (58) is fitted to each of the six different data sets (Figure 55) according to different flow and time condition s, the coefficient I was found to be 6.71.9 (meanSD) mm1. Day 4Distance from IEL (mm) 0.00.10.20.30.40.5TUNEL Positive Cells (%) 10 0 10 20 30 40 50 60 High Flow Normal Flow Day 7Distance from IEL (mm) 0.00.10.20.30.40.5TUNEL Positive Cells (%) 10 0 10 20 30 40 50 60 High Flow Normal Flow Day 14Distance from IEL (mm) 0.00.10.20.30.40.5TUNEL Positive Cells (%) 10 0 10 20 30 40 50 60 High Flow Normal Flow Figure 55. % of TUNEL positive cells as a function of actual distance from IEL (From Berceli et al. (2002)). PAGE 87 87 Then, the final form of cell apoptosis probability can be written as follows. Papoptosis ( *, t, ) = (59) where is normalized shear stress, t is time, is distance from IEL, and H* is a constant. The constant H* in Equation (59) is det ermined at each time step such that the integration of Equation (59 ) over the intimal thickne ss equals the space averaged probability multiplied by intimal thickness. H* = (510) where IT is intimal thickness (mm). Probability of Extracellular Matrix Behaviors It has been reported that the volume % of SM C in vein graft intima ranges from about 20% to 40% regardless of intimal thi ckness (Kohler et al., 1991; Kraiss et al., 1991; Zwolak et al., 1987). In this study, total of 4 ECM elements are produced from each SMC division. Each of t he two daughter cells produce o ne matrix element after 24 hours from the cell division, and one more matrix element after another 24 hours. Sometimes, there is the case where one or both of the daught er cells undergo cell division before producing total of 4 matrix elements. From the onedimensional simulation which will be described later in this chapter, the volume % of SMC at Day 28 resulted in 27%. Since SMC division probability is higher near the lumen compared to the one away from the lumen, there is higher chance of SMC undergoing cell division before producing total of 4 matr ix elements, which would resu lt in higher volume % of PAGE 88 88 SMC near the lumen. When a cell undergoes apoptosis, 4 matrix elements nearest to the cell are removed from the computati onal domain. Table 51 summarizes the assumptions used in this study. Table 51. Assumptions used in the model. Category Assumptions BrdU BrdU is available for the first 2 hours. Cell Cycle Cell cycle time is 24 hours. G1, S, and G2/M phases each takes up onethird of the cell cycle time. Cell Size Cell size is 7m 7 m. Initial volume ratio of cell to matrix is set to be 1:3. Dead Cells Dead cells remain in the tissue for 24 hours before getting removed. Matrix Production Each of the two new daughter cell elements produces two matrix elements with 24 hours time interval. Effects of Monocytes Monocyte Entry Model Based on the data in Jiang et al. (2009b), the same model form used in the BrdU model (Equation (52)) is adopt ed for the model of the net ra te of monocyte entry into vein graft wall. (Net rate of monocyte entry) = (511) where t is time, is normalized shear stress, and A and B are constants. B is set to be the same value (0.16) used in the BrdU model. When A is set to be 175/(mm2hour), the volume occupation of monocyte in the inti ma is simulated to be about 5%, which is consistent with the values f ound in the literature (Lemson et al., 2000). Each monocyte that has entered the computational dom ain is placed at a random position. Monocyte Effects on ECM Synthesis PAGE 89 89 Hoch et al. (1999) showed that macrophage depletion suppresses intimal hyperplasia in vein graft. However, they observed that there was no significant change in the total number of cells, implying decreased synthesis of extracellular matrix. As a first approximation, the effects of macrophages can be incorporated into the model via the amount of matrix created as follows. Nm = [2 ] (512) where Nm is the number of matrix elements generated by each of the 2 new daughter cells, is normalized monocyte entry rate, and [ ] denotes the rounding to the nearest integer. The normalization is with respect to the monocyte entry rate calculated from the model shown above. The two matrix elements from each of the tw o new daughter cells are produced with time interval of 24 hours. Figure 56 shows the simulated intimal area as a function of the normalized monocyte entry rate at shear stress of 1.8 dynes/cm2, where a positive correlation is observed. Average of 104 simulations Normalized Rate of Monocyte Entry 0.00.51.01.52.0 Intimal Area (cm2) 0.000 0.005 0.010 0.015 0.020 0.025 Figure 56. Simulated intimal area with normalized rate of monocyte entry. PAGE 90 90 Algorithm 1Dimensional Algorithm To test the feasibility of t he model including the probabiliti es of cellular activities, the simulation was performed on 1D domain, which is simpler than 2D simulation. 1D domain is a stack of elements in one line wh ere each element is occupied by either SMC or ECM (Figure 57). The size of cell was chosen to be 7m7m. When a cell is divided, the new cell is placed either on the right side or on the left side of the cell with equal probability. When a new cell is placed the new cell pushes all the elements on the luminal side to make r oom for itself. The same procedure is performed when a new ECM element is produced by SMC element. W hen a cell dies, the empty space is filled by the elements on the luminal side by shi fting themselves to f ill the empty space. Figure 58 shows the flow c hart of the general algorithm. Figure 57. Schematic of 1dimensional domain. (C: Cell element, M: Matrix element) 2Dimensional Algorithm In 2D domain, when a new cell or matrix el ement is produced, it is placed at one of the 4 adjacent spaces pushing the existing element in the luminal direction (Figure 59A). The new cell or matrix does not push t he adjacent elements in the direction normal to the luminal direction, because the incr ease of the length of vein graft in the longitudinal direction is not observed experi mentally. Also, the new element does not push the adjacent element in the opposite direction to the luminal direction because PAGE 91 91 Figure 58. Flow chart for the rulebas e simulation of vein graft remodeling. the surrounding tissue is assumed to prev ent such movement. Also not observed experimentally is excessive protrusion of the wall into the lumen. To prevent the excessive protrusion, a redi stribution algorithm is appli ed as shown in Figure 59B. When there is more than 2 elements differe nce between adjacent columns, the last element in the longer column moves to the shorter column. In computational domain for a more arbi trary geometry such as focal stenosis, the redistribution algorithm is a little differ ent (Figure 510). For all the elements facing PAGE 92 92 lumen, the empty spaces ar ound the element is examined, and if there is an empty space closer to the IEL, then the el ement is moved to the empty space. Figure 59. 2D algorithm of rulebased simulation of vein graft remodeling. A. New element placement algorithm. B. Element redistri bution algorithm. Results 1D Model Validation Initial volume ratio of cell to matrix is set to be 1:3. This ratio determines the initial number of cells, which affects the simulation results. Figure 511A shows an example of the change of the number of elements with time in a single simulation, which shows the stochastic nature of the simulation. Figure 511B shows the location of the divided cells in a single simulation. Most of the cells divided near the lumen as expected from Equation (53). PAGE 93 93 Figure 510. Element redi stribution algorithm for stenosis simulation. Figures 511C and 511D show the com parison between experimental data and simulation results for intimal area in radi al crosssection for low (shear stress=1.8 dynes/cm2) and high (shear stress=14 dynes/cm2) flows, respectively. The intimal area in 2D crosssection was calculated from the 1D simulation by multiplying the 1D intimal area by the initial number of cells in 2D domain. The initial number of cells in 2D domain was calculated by dividing the init ial number of all the elements by 4 because 1 out of 4 elements is cell element. The simulation curves are the average of 104 runs, and agree well with experimental data though t he simulation tends to underestimate the intimal area for low flow case. The standard deviation is 2.2 and 0.65 times the average value at Day 28 for low and high flows, respectively. PAGE 94 94 Example of 1 simulationTime (Day) 051015202530 Number of Elements 0 10 20 30 40 50 60 70 Intimal Area (cm2), Low FlowAverage of 104 1D simulations (Average Shear Stress=1.8 dynes/cm2)Time (day) 0102030 0.000 0.005 0.010 0.015 0.020 Experiment Simulation (1D) Intimal Area (cm2), High FlowAverage of 104 Simulations (Average Shear Stress=14 dynes/cm2)Time (day) 0102030 0.000 0.005 0.010 0.015 0.020 Experiment Simulation (1D) Figure 511. 1D simulation results. A. Change of the number of elements with time. B. Location of cell division. C. Intimal area fo r low flow condition. D. Intimal area for high flow condition. Experimental data are from Jiang et al. (2004). 2D Model Validation Figure 512 shows an example of 2D simu lation. The simulation starts with one line of elements composed of SMC (25%) and ECM (75%) elements. The simulation is 0 10 20 30 40 50 60 70 36789111315161723Number of ElementsTime (Day) Location of Divided Cell Cell and Matrix Divided CellLumenA B C D PAGE 95 95 run up to 28 days with time step of 1 hour. Most of the growth of the wall area occurs up to 2 weeks, and the wall area does not change much from 2 weeks to 4 weeks. More monocytes (grey cell) are located away fr om the lumen because m onocyte migration is not taken into account in the simulation. Ac cording to Equation (511), more monocytes enter the wall at earlier time poi nts. Due to the redistribution algorithm (Figure 59B), no two adjacent columns show the difference of the element number greater than 1. Figure 513 shows 2D simulation re sult of intimal area compared with experimental data. Among the simulation para meters, cell size and BrdU availability time are the ones that are uncertain and have significant impact on simulation results. As such, the two parameters were determi ned such that the simulated intimal area matches well with the experimental data. Figure 512. An example of 2D simulation. Empty elements are ECM. The elements with a filled circle at the center ar e SMC. Grey elements are monocytes. Day 1 Day 7 Day 14 Day 28 PAGE 96 96 Figure 513 also compares intimal ar eas from 1D and 2D simulations. 1D simulation predicts less wall area than 2D simu lation because the greater thickness in 1D domain makes the cell proliferation concentrated near the lumen with higher probability which gives cells less time to pr oduce matrix elements. The reason for the concentrated area of higher prob ability near the lumen for 1D simulation is because the same coefficient E in Equation (53) is app lied to both 1D and 2D simulation. Because the space averaged probability of cell divi sion is the same for both 1D and 2D simulation, the probability at the surface of the lumen s hould be higher in 1D case (thicker intima) to maintain t he same coefficient E in Equation (53). The 2D simulation curve in the low flow case is Sshaped becaus e of the initial burst and later decrease of cell division. Intimal Area (cm2), Low FlowAverage of 104 simulations (Average Shear Stress=1.8 dynes/cm2)Time (day) 0102030 0.000 0.005 0.010 0.015 0.020 Experiment Simulation (2D) Simulation (1D) Intimal Area (cm2), High FlowAverage of 104 Simulations (Average Shear Stress=14 dynes/cm2)Time (day) 0102030 0.000 0.005 0.010 0.015 0.020 Experiment Simulation (2D) Simulation (1D) Figure 513. Simulation results and experim ental data of intimal area and comparison of 1D and 2D simulations for low and high flow conditions. Experimental data are from Jiang et al. (2004). PAGE 97 97 Focal Stenosis Simulation Figure 514 shows an example of focal st enosis simulation. The initial curve was reconstructed from experimental data (Figur e 515) using smoothing spline algorithm as described in chapter 3. The initial curv e is placed on a rectangular grid, and the elements through which the curve passes are assigned either cell (25%) of matrix (75%). Shear stress was calculated based on Poiseuille flow assumption. Greater intimal thickening is observed distal to the st enosis where shear stress is relatively low. In the figure, black elements are cells, and grey elements are matrix. Figure 514. An example of focal stenosis si mulation. Black dots are cells, and grey area is extracellular matrix. PAGE 98 98 Figure 515 shows the comparison bet ween the simulation and experimental data. The simulation result overpredicts the in timal thickening distal to the stenosis at Day 28, which is consistent with the overp rediction observed in the 2D validation at Day 28 (Figure 513). This suggests that the experimental data at Day 14 may contain significant amount of error. Figure 515. Comparison between simulation result and experimental data for intimal thickening around focal stenosis. Exper imental data are from the vascular biology laboratory (Dr. Scott A. Berceli). Discussion Rulebased modeling method is appropriate for simulation of biological systems because biological systems are intrinsically multiscale systems in which different scale components such as molecules, cells, ti ssue, and organs comprise the system. The behaviors of greater scale component result from the interactions among the smaller PAGE 99 99 scale components, and the rulebased method enables the observati on of the greater behavior emerging from the interactions am ong the small scale components in the system. Vein graft wall consists of mostly smooth muscle cells and extracellular matrix, and the wall thickening is the result of t he cell proliferation and matrix deposition, and the interactions among them. In this chapter, the wall thi ckening has been simulated by applying rules of behaviors to cell and matrix elements. The rules were determined from the exper imental data when they are available, and those rules that were not available exper imentally were obtaine d either from the literature or by assuming appropriate va lues. Experimental data often contains uncertainties and needs to be analyzed more to be used as quantitative values for simulation. In this study, BrdU experiment al data needed to be analyzed more to extract cell division probability from the data. The mo st uncertain thing about BrdU data is how many hours BrdU is available after it is injected into the body. Phuphanich and Levin (1985) found that the halflife of BrdU is on the or der of hours, but it depended on whether it was administered intravenously or or ally. In this study, the BrdU availability time together with other parameters were dete rmined by fitting the si mulation results to experimental data of wall area (global behav ior). The BrdU availability time determined in this manner is 2 hours which is consis tent with the observation by Phuphanich and Levin (1985). The other parameters determined together with the BrdU availability time were cell size and initial volume ratio of cell to matrix. The initial volume ratio of cell to matrix determines the initial number of cells which signific antly affects the number of cells at later time points. The initial rati o of 1(cell):3(matrix) was tested based on the observation in the literature that the cells oc cupy 20% to 40% of wall volume (Kohler et PAGE 100 100 al., 1991; Kraiss et al., 1991; Zwolak et al., 1987) regardless of wall thickness, and this ratio together with cell size of 7m7m re sulted in simulation results that match experimental data very well (Figure 513). When the rulebased approach was applied to the simulation of wall thickening around focal stenosis, it predicted thicker wall at the regions where wall shear stress is lower. Though this is consistent with exper imental observation, Figure 515 shows that the model overpredicts the wall thickne ss compared to experim ental data where wall shear stress is lower. This overprediction is consistent with that observed in Figure 513 where the intimal area is overpredicted at Day 28 for low flow condition. As shown in the Figure 513, the experimen tal data at Day 14 shows a peak compared with the other data points, but there seems to be no physiological reason that the intimal area should decrease from Day 14 to Day 28. Neverthele ss, the model parameters were determined such that the model prediction fits the exper imental data including t he data point at Day 14. As a result, the model overpredicts t he experimental data at Day 28, and this is also shown in the focal stenosis simulation as shown in Figure 515. Considering that the two graphs (Figure 513 and Figure 515) show the same tendency to overprediction at Day 28, it is very likely t hat the experimental dat a at Day 14 contains significant amount of variation. PAGE 101 101 CHAPTER 6 CONCLUSION Vein graft remodeling is a complex proce ss involving a large number of biological and hemodynamic factors. Thes e factors act together affe cting the outcomes of the other factors. The cells in the vascula r system respond to hemodynamic forces, and collectively generate global pheno mena such as wall thickening. Endothelial cells are in direct contact with blood flow, thus respond to wall shear stress imposed by blood flow, and release signaling molecules which modulate wall remodeling. On the other hand, wall tension acts on the whole wall, and thus smooth muscle ce lls, which are the primary cell type in the vascular system, respond to wall tension. In this dissertation, the effects of hemodynamic forces on vein graft remodeling have been examined through mathematical and computat ional modeling approach. Vein graft experiences higher flow and pr essure environment when implanted in the arterial system compared with the hemodynamic environm ent of the venous system. Vein graft adapts to the arterial hemodynam ic environment by thickening its wall and expanding its diameter. Thickening wall r educes wall tension, and increasing diameter reduces wall shear stress. The size of the lumen is determined by these two remodeling processes. After this adaptation, some ve in grafts proceed to become pathologically occlusive just like the blocked artery that t he vein graft bypasses. The occlusion occurs at focal regions while the neighboring r egions remain patent. The mechanism of occlusive adaptation seems to be different fr om that of early ada ptation, and the focus of this dissertation has been modeli ng the early stage of the adaptation. In the early stage of the ve in graft adaptation, the mo st notable change occurs in the layer of intima. The thickness of the in timal layer become comparable to that of PAGE 102 102 medial layer in the time frame of weeks, while the intima in normal vein is just a single cell layer. The intima is composed of mo stly SMC and ECM. SMC comes from medial layer and proliferate in the intimal layer, and ECM is produced by SMC. A mathematical model of this intimal layer remodeling has been developed in chapter 3. Experimental data show shear stress imposed by blood fl ow has significant impact on the intimal thickening, and the intimal thickening rate is dependent on time. Thus, the mathematical model is a function of wall shear stre ss and time. The mathem atical form has been determined based the behavior of the experiment al data. As simple model form as possible has been selected, while the experimental data has relatively large variations from individual animal to animal. When the model was applied to a focal stenosis simulation, the simulation results agreed well with the focal stenosis experimental data. Certainly, shear stress is not the only fact or influencing the intimal layer, and all the other factors can be considered to be include d in the time term. However, the model successfully describes the change of intimal la yer, where most notable change occurs in the vein graft remodeling, using onl y two variables: shear stress and time. Even though the intimal layer shows the most notable change, medial layer also thickens, and the diameter of EEL layer c hanges as well. The size of lumen can be considered to be determined by the geomet rical relation among these three layers. Mathematical models of these three laye rs have been developed in chapter 4. The mathematical form is based on logistic model and the experimental data at 3 months and 6 months were added to the one used in chapter 3. Wall tension has been added as a hemodynamic factor since experimental dat a showed that it is related to medial layer remodeling. Experimental data show t hat the diameter change of EEL layer is not PAGE 103 103 very related to shear stress. This is so mewhat surprising because it has been believed that vein grafts dilate to reduce wall shear stress. However, at higher shear stress, the intimal layer thickening is reduced, whic h makes the size of the lumen greater. Experimental data of vein gra ft remodeling show Sshaped behavior which is typically observed in cell population growth, and the logist ic model which is one of the simplest form that describes Sshaped behavior succe ssfully modeled the behaviors of each layer of vein graft. When the model was coup led with a flow solver, and applied to the simulation of idealized stenos is and human vein graft, it predicted the remodeling of each layer of vein graft around stenosis regi on depending on the severity of the stenosis at different flow conditions. Even though mathematical model of vein graft remodeling explicitly shows the relationship between the hemodynamic factor s and the geometrical variables of each layer of vein graft, biologic mechanism behi nd the vein graft morphology change is not included in the model equation. To overcome this limitation, rulebased modeling approach has been applied to the simulation of vein graft remodeling in chapter 5. The components in the model are cells and matrix elements, and rule of behaviors were applied to the elements. The rules were selected based on experimental data. The model parameters were determined based on the evidences in the literature, but those that were hard to find in the literature were determined by comparing the simulation results with experimental data. Simulation results show intimal thickening (global behavior) around a focal stenosis, which resulted from the local interactions among the cell and matrix elements. The simulation resu lts also show advanced intimal thickening distal to the throat where wall shear stress is lower, which is ex pected from the model, PAGE 104 104 and agrees with experimental observations. Ho wever, quantitative comparison shows that the model overpredicts the intimal thickness at the low shear region, presumably due to the variations in the ex perimental data the model is bas ed on. It is consistent with 2D model validation where t he simulation results were co mpared with another set of experimental data. The models developed in this dissertat ion are focused on predictive capability rather than explanatory one, though explanatory model usually precedes the predictive one (Engelberg et al., 2011). Despite significant advances in the understanding of the mechanism of vein graft failure and surgical techniques, vein graft failure rate is still significantly high as described in chapter 1, and if it can be predicted whether a vein graft will be successful or not after the surger y, healthcare cost would be significantly reduced. The models in this dissertation are closely linked to experimental data which were obtained within 6 months period, and t he predictive capability beyond 6 months was not confirmed in this study. To make both the mathematical and rulebased models be able to predict longer time period than 6 months without experimental data, more fundamental biologic mechanisms would have to be incorporated into the models, and, because the complete biologic mechanism of ve in graft failure is currently unknown, the investigation of biologic me chanism together wit h theoretical modeling work would be needed for the development of a successful predictor of vein graft failure. PAGE 105 105 LIST OF REFERENCES Alber, M.S., Kiskowski, M.A. Glazier, J. A., Jiang. Y., 2003. On cellular automaton approaches to modeling biological cells. In: Rosenthal, J., Gilliam, D.S. (Eds.), Mathematical Systems Theory in Biol ogy, Communications, Computation, and Finance. SpringerVerlag, New York, 140. Alford, P.W., Humphrey, J.D., Taber, L.A. 2008. Growth and remodeling in a thickwalled artery model: effects of spatial vari ations in wall constituents. Biomechanics and Modeling in Mechanobiology 7, 245262 Alon, R., Kassner, P.D., Carr, M.W., Finger E.B., Hemler, M.E., Springer, T.A., 1995. The integrin VLA4 supports tethering and rolling in flow on VCAM1. Journal of Cell Biology 128, 12431253 Ausk, B.J., Gross, T.S., Srinivasan, S. 2006. An agent based model for realtime signaling induced in osteocytic networks by mechanical stimuli. Journal of Theoretical Biology 39, 26382646. Bailey, A.M., Lawrence, M.B., Shang, H., Ka tz, A.j., Peirce. S. M., 2009. Agentbased model of therapeutic adiposederived st romal cell trafficking during ischemia predicts ability to roll on Pselectin. P LoS Computational Bi ology 5, e1000294. Bailey, A.M., Thorne, B.C., Pe irce, S. M., 2007. Multicell agentbased simulation of the microvasculature to study the dynamics of circulating inflammatory cell trafficking. Annals of Biomedical Engineering 35, 916936. Bartha, K., Rieger, H., 2006. Vascular network remodeling via vessel cooption, regression and growth in tumors. Journal of Theoretical Biology 241, 903918. Berceli, S.A., Davies, M.G., Kenagy, R. D., Clowes, A.W., 2002. Flowinduced neointimal regression in baboon polytetrafluor oethylene grafts is associated with decreased cell prolifer ation and increased apoptosis. J ournal of Vascular Surgery 36, 12481255. Berceli, S.A., TranSonTay, R., Garbey, M., Jiang, Z., 2009. He modynamically driven vein graft remodelig: a systems biology approach. Vascular 17, (Suppl 1)S2S9. Berguer, R., Higgins, R.F., Reddy, D.J., 1980. Intimal hyperplasia. Archives of Surgery 115, 332335. Bhardwaj, S., Roy, H., YlHerttuala, S., 2008. Gene therapy to prevent occlusion of venous bypass grafts. Expert Review of Cardiovascular Therapy 6, 641652. PAGE 106 106 Bonabeau, E., 2002. Agentbased modeling: Me thods and techniques for simulating human systems. Proceedings of the Nati onal Academy of Sciences U.S.A. 99, 72807287. Bongrand, P., Bell, G.I., 1984. Cellce ll adhesion: parameters and possible mechanisms. In: Perelson, A.S., DeLisi, C., Wiegel, F. W. (Eds.), Cell surface dynamics, Concepts and models. Marcel Dekker, New York. Boyle, C.J., Lennon, A.B., Early, M., Kelly D.J., Lally, C., Predergast, P.J., 2010. Computational simulati on methodologies for mec hanobiological modeling: application of a cellcent ered approach to neointima development around stents. Philosophical Transactions. Series A, Mathematical, Physi cal, and Engineering Sciences 368, 29192935. Britton, N.F., 2003. Essential Mat hematical Biology. Springer, London. Bruehl, R.E., Springer, T.A., Bain ton, D.F., 1996. Quant itation of Lselectin distribution on human leukocyte microvilli by immunogol d labeling and electron microscopy. Journal of Histochemistry and Cytochemistry 44, 835844. BuduGrajdeanu, P., Schugart, R.C., Friedman, A., Valentine, C., Agar wal, A.K., Rovin, B.H., 2008. A mathematical model of v enous neointimal hyperplasia formation. Theoretical Biology and Medical Modeling 5, 2. Byrne, H., Drasdo, D., 2009. Individualbased and conti nuum models of growing cell populations: a comparison. Journal of Mathematical Biololgy 58, 657687. Checa, S., Prendergast, P.j., 2009. A mechanobiological model for tissue differentiation that includes angiogenesis: A lattice based modeling approach. Annals of Biomedical Engine ering 37, 129145. Chen, S., Springer, T.A., 2001. Selectin re ceptorligand bonds: Formation limited by shear rate and dissociation governed by the Bell model. Pr oceedings of the National Acadamy of Sci ences U S A 98, 950955. Cheng, G., Youssef, B.B., Ma rkenscoff, P., Zygourakis. K., 2006. Cell population dynamics modulate the rates of tissue growth processes. Biophysical Journal 90, 713724. Chavali, A.K., Gianchandani, E.P., Tung, K.S ., Lawrence, M.B., Peirce, S.M., Papin, J.A., 2008. Characterizing em ergent properties of i mmunological systems with multicellular rulebased computatioanl modeling. Trends in Immunology 29, 589599. Conte, M.S., Bandyk, D.F., Clo wes, A.W., Moneta, G.L., Seel y, L., Lorenz, T.J., Namini, H., Hamdan, A.D., Roddy, S.P., Belkin, M. Berceli, S.A., DeMasi, R.J., Samson, PAGE 107 107 R.H., Berman, S.S., and P. I. Invest igators, 2006. Results of PREVENT III: A multicenter, randomized trial of edifoligide fo r the prevention of vein graft failure in lower extremity bypass surgery. Jour nal of Vascular Surgery 43, 742751. Deutsch, A., Dormann, S., 2005. Cellular Au tomaton Modeling of Biological Pattern Formation. Birkhuser, Boston. Dewitt, S., Hallett, M., 2007. Leukocyte memb rane "expansion": a central mechanism for leukocyte extravasation. Jour nal of Leukocyte Biology 81, 11601164. Dobrin, P.B., Littooy, F.N., Endean, E.D., 1989. Mechanical factors predisposing to intimal hyperplasia and medial thickeni ng in autogenous vein gr afts. Surgery 105, 393400. Dong, C., Lei, X.X., 2000. Biomechanics of cell rolling: shear flow, cellsurface adhesion, and cell deformability. Jour nal of Biomechanics 33, 3543. Dormann, S., Deutsch, A., 2002. Modeling self organized avascular tumor growth with a hybrid cellular automaton. In Silico Biology 2, 0035. Dwir, O., Solomon, A., Mangan, S., Kansas G.S., Schwarz, U.S., Alon, R., 2003. Avidity enhancement of Lselectin bonds by flow: shearpromoted rotation of leukocytes turn labile bonds into functional tethers. Journal of Cell Biology 163, 649659. Engelberg, J.A., Datta, A ., Mostov, K.E., Hunt, C.A., 2011. MDCK Cystogenesis Driven by Cell Stabilization within Computati onal Analogues. PLoS Computational Biology 7(4), e1002030. Engelberg, J.A., Ropella, G. E.P., Hunt, C.A., 2008. Essent ial operating principles for tumor spheroid growth. BMC Systems Biology 2, 110. Ermentrout, G.B., Edelstei nKeshet, L., 1993. Cellular Automata Approaches to Biological Modeling. J ournal of Theoretical Biology 160, 97133. Fernandez, C.M., Goldman, D.R. Jiang, Z., Ozaki, C.K., Tr anSonTay, R., Berceli, S.A., 2004. Impact of shear st ress on early vein graft re modeling: A biomechanical analysis. Annals of Biomedica l Engineeri ng 32, 14841493. Ferreira, Jr., S.C., Martins, M.L., Vilela, M.J., 2002. R eactiondiffusion model for the growth of avascular tumor. Physical Review E 65, 021907. Figueroa, C.A., Baek, S., Taylor, C.A., Humphrey, J.D., 2009. A computational framework for fluidsolidgrowth model ing in cardiovascular simulations. Computational Methods in Applied Me chanics and Engineering 198, 35833602. PAGE 108 108 Fillinger, M.F., Cronenwett, J.L., Besso, S ., Walsh, D.B., Zwolak, R.M., 1994. Vein adaptation to the hemodynamic environment of infrainguin al grafts. Journal of Vascular Surgery 19, 970979. Finger, E.B., Purl, K.D., Alon, R., Lawrence, M.B., Von Andr ian, U.H., Springer, T.A., 1996. Adhesion through Lse lectin requires a threshold hydrodynamic shear. Nature 379, 266269. Fitzgibbon, G.M., Kafka, H.P., Leach, A.J. Keon, W.J., Hooper, G.D., Burton., J.R., 1996. Coronary bypass graft fa te and patient outcome: A ngiographic followup of 5,065 grafts related to survival and reoper ation in 1,388 pati ents during 25 years. Journal of the American Coll ege of Cardiology 28, 616626. Friedman, M.H., Deters, O.J., Bargeron, C. B., Hutchins, G.M., Ma rk, F.F., 1986. Sheardependent thickening of the human arterial intima. Atherosclerosis 60, 161171. Fujikawa, H., Kai, A., Morozumi, S., 2004. A new logistic model for Escherichia coli growth at constant and dynamic temper atures. Food Microbiology 21, 501509. Galle, J., Hoffmann, M., Aust, G., 2009. From single cells to tissue architecturea bottomup approach to modelling the spat iotemporal organisation of complex multicellular systems. Journal of Mathematical Biology 58, 261283. Galle, J., Loeffler, M., Drasdo, D., Modeling the effect of deregulated proliferation and apoptosis on the growth dynamics of epithe lial cell populations in vitro. Biophysical Journal 88, 6275. Galt, S.W., Zwolak, R.M., Wagner R.J., Gilbertson, J.J., 1993. Differential response of arteries and vein grafts to blood flow r eduction. Journal of Vascular Surgery 17, 563570. Galvo, V., Miranda, J.G.V., Ribeirodos Santos, R., 2008. Development of a twodimensional agentbased model for chronic chagasic cardiomyopathy after stem cell transplantation. Bi oinformatics 24, 20512056. Gerlee, P., Anderson, A.R.A ., 2007. An evolutionary hybrid cellular automaton model of solid tumour growth. Journal of Theoretical Biology 246, 583603. Gerlee, P., Anderson, A.R. A., 2008. A hybrid cellular automaton model of clonal evolution in cancer: The emergence of the glycolytic phenotype. Journal of Theoretical Biology 250, 705722. Gervertz, J.L., Torquato, S., 2006. Modeling the effects of va sculature evolution on early brain tumor growth. Journal of Theoretical Biology 243, 517531. Gilbert, N., 2008. AgentBased Models. Sage Publications, Los Angeles. PAGE 109 109 Glagov, S., Weisenbeg, E., Zarins, C.K., St ankunavicius, R., Kolettis, G.J., 1987. Compensatory enlargement of human atherosclerotic coronary arteries. New England Journal of Medicine 316, 13711375. Goetz, D.J., ElSabban, M.E ., Pauli, B.U., Hammer, D.A ., 1994. Dynamics of neutrophil rolling over stimulated endothelium in vitr o. Biophysical Journal 66, 22022209. Goldman, A.J., Cox, R.G., Brenner, H., 1967. Slow viscous motion of a sphere parallel to a plane wall. I. Motion through a quiesc ent fluid. Chemical Engineering Science 22, 637652. Goldman, A.J., Cox, R.G., Brenner, H. 1967. Slow viscous motion of a sphere parallel to a plane wall. II. Couette Flow. Chem ical Engineering Science 22, 653659. Goudar, C.T., Joeris, K., Konstantinov, K. B., Piret, J.M., 2005. Logistic Equations Effectively Model Mammalian Cell Batch and FedBatch Kinetics by Logically Constraining the Fit. Biot echnology Progress 21, 11091118. Grabe, N., Neuber, K., 2005. A multicellular systems biology model predicts epidermal morphology, kinetics and Ca2+ flow. Bioinformatics 21, 35413547. Grant, M.R., Mostov, K.E., Tisty, T.D., Hunt, C.A., 2006. Simulati ng properties of in vitro epithelial cell morphogenesis. PLoS Computational Biology 2, e129. Grondin, C.M., Lepage, G., Castonguay, Y. R., Meere, C., Grondin, P., 1971. Aortocoronary bypass graft: Initial bl ood flow through the graft, and early postoperative patency. Cir culation 44, 815819. Guo, Z., Sloot, P.M.A., Ta y, J.C., 2008. A hybrid agent based approach for modeling microbiological systems. Journal of Theoretical Biology 255, 163175. Haga, M., Yamashita, A., Paszkowiak, J., Su mpio, B.E., Dardik, A., 2003. Oscillatory shear stress increases smooth muscle cell proliferation and Ak t phosphorylation. Journal of Vascular Surgery 37, 12771284. Hammer, D.A., Apte, S.M., 1992. Simulation of cell rolling and adhesion on surfaces in shear flow: general results and analysis of selectinmediated neutrophil adhesion. Biophysical Journal 63, 3557. Hoch, J.R., Stark, V.K., R ooijen, N., Kim, J.L., Nutt M.P., Warner, T.F., 1999. Macrophage depletion alters vein graft in timal hyperplasia. Surgery 126, 428437. Hwang, M., Garbey, M., Berceli, S.A., TranSonTay, R., 2009. RuleBased Simulation of MultiCellular Biological SystemsA Review of Modeling Techniques. Cellular and Molecular Bioengineering 2, 285294. PAGE 110 110 Jacot, J.G., Abdullah, I., Belk in, M., GerhardHerman, M., Gaccione, P., Polak, J.F., Donaldson, M.C., Whittemore, A.D., Cont e, M.S., 2004. Early adaptation of human lower extremity vein grafts: Wall st iffness changes accompany geometric remodeling. Journal of Va scular Surgery 39, 547555. Jadhav, S., Eggleton, C.D., Konstantopoulos K., 2005. A 3D co mputational model predicts that cell deformation affects se lectinmediated leukocyte rolling. Biophysical Journal 88, 96104. Jeremy, J.Y., Gadsdon, P., Shukla, N., Vijay an, V., Wyatt, M., Ne wby, A.C., Angelini, G.D., 2007. On the biology of saphenous vein grafts fitted with external synthetic sheaths and stents. Biom aterials 28, 895908. Jiang, Y., PjesivacGrbovic, J., Cantrell, C ., Freyer, J.P., 2005. A multiscale model for avascular tumor growth. Bioph ysical Journal 89:38843894. Jiang, Z., Tao, M., Omalley, K.A., Wang, D., Ozaki, C. K., Berceli, S.A., 2009. Established neointimal hyperplasia in vein grafts expands via TGFbmediated progressive fibrosis. American Journal of Physiology Heart and Circulatory Physiology 297, H1200H1207. Jiang, Z., Wu, L., Miller, B.L., Goldman, D.R., Fer nandez, C.M., Abouhamze, Z.S., Ozaki, C.K., Berceli, S.A., 2004. A novel ve in graft model: Adapta tion to differential flow environments. American Journal of Physiology Heart and Circulatory Physiology 286, H240H245. Jiang, Z., Yu, P., Tao, M., Fernandez, C., Ifa ntides, C., Moloye, O., Schultz, G.S., Ozaki, C.K., Berceli, S.A., 2007. TGF{bet a}/CTGF mediated fibr oblast recruitment influences early outward vein graft remode ling. American Journal of Physiology Heart and Circulatory Physiology 293, H482H488. Jiang, Z., Yu, P., Tao, M., Ifant ides, C., Ozaki, C.K., Berceli, S.A., 2009. Interplay of CCF2 signaling and local shear force determi nes vein graft neointimal hyperplasia in vivo. FEBS Letters 583(21), 35363540. Jung, U., Bullard, D.C., Tedder, T.F., Ley, K., 1996. Velocity differences between Land Pselectindependent neutrophil rolling in venules of mouse cremaster muscle in vivo. Americal Journal of Physiology Heart and Circulatory Physiology 271, H2740H2747. Kansal, A.R., Torquato, S., Harsh IV, G.R., Chiocca, E.A., Deis boeck, T.S., 2000. Simulated brain tumor growth dynami cs using a threedi mensional cellular automaton. Journal of Theor etical Biology 203, 367382. PAGE 111 111 Kim, S.H.J., Park, S., Mostov, K., Debnat h, J., Hunt, C.A., 2009. Computational Investigation of Epithelial Ce ll Dynamic Phenotype In Vitr o. Theoretical Biology and Medical Modelling 6, 8. Knutton, S., Sumner, M.C.B. Pasternak, C.A., 1975. Role of microvilli in surface changes of synchronized P815Y mastocytom a cells. Journal of Cell Biology 66, 568576. Kohler, T.R., Kirkman, T.R., Kraiss, L.W., Zierler, B.K. Clowes, A.W., 1991. Increased blood flow inhibits neointimal hyperpla sia in endothelialized vascular grafts. Circulation Research 69, 15571565. Korshunov, V. A., Berk, B.C., 2004. Straindependent vascular remode ling: The "Glagov phenomenon" is genetically determi ned. Circulation 110, 220226. Kraiss, L.W., Kirkman, T.R., Kohler, T.R., Zierler, B., Clo wes, A.W., 1991. Shear stress regulates smooth muscle proliferation and neointimal thi ckening in porous polytetrafluoroethylene grafts. Arteri osclerosis and Thrombosis 11, 18441852. Lao, B.J., Kamei, D.T., 2008. Investigation of cellular movement in the prostate epithelium using an agentbased model. J ournal of Theoretical Biology 250, 642654. Lawrence, M.B., Kansas, G.S., Kunkel, E.J. Ley, K., 1997. Threshold levels of fluid shear promote leukocyte adhesion through sele ctin (CD62L, P, E) Journal of Cell Biology 136, 717727. Lawrence, M.B., Springer, T.A., 1991. Leukocytes roll on a selectin at physiologic flow rates: distinction from and prerequisite for adhesion through integrins. Cell 65, 859873. Lee, D., Chiu, J.J., 1992. A numerical simula tion of intimal thickening under shear in arteries. Biorheology 29, 337351. Lemson, M.S., Tordoir, J.H.M., Daemen, M.J. A.P., Kitslaar, P.J.E.H.M., 2000. Intimal Hyperplasia in Vascular Grafts. Eur opena Journal of Vascular and Endovascular Surgery 19, 336350. Mai, J., Sameni, M., Mikkelsen, T., Sloane, B.F., 2002. Degradation of Extracellular Matrix Protein TenascinC by Cathepsin B: An Interaction Involved in the Progression of Gliomas. Bi ological Chemistry 383, 14071413. Mallet, D.G., De Pillis, L.G., 2006. A cellula r automata model of tumorimmune system interactions. Journal of Theor etical Biology 239, 334350. PAGE 112 112 Mann, M.J., 2004. Novel Strategies for the Pr evention of Bypass Graft Failure. BioDrugs 18, 18. Mansury, Y., Deisboeck, T.S., 2003. The impac t of "search precision" in an agentbased tumor model. Journal of Theor etical Biology 224, 325337. Meyerson, S.L., Skelly, C.L., Curi, M.A. Shakur, U.M., Vosicky, J.E., Glagov, S., Schwartz, L.B., 2001. The effects of extremely low shear stress on cellular proliferation and neointimal thickening in the failing bypass graft. Journal of Vascular Surgery 34, 9097. Mills, J.L., Bandyk, D.F ., Gahtan, V., Esses, G.E ., 1995. The origin of infrainguinal vein graft stenosis: A prospective study bas ed on duplex surveillance. Journal of Vascular Surgery 21, 1625. Mitra, A.K., Gangahar, D.M., Agrawal, D,K., 2006. Cellular, molecular and immunological mechanisms in the pat hophysiology of vein graft intimal hyperplasia. Immunology and Cell Biology 84, 115124. Morinaga, K., Eguchi, H., Miyazaki, T. Okadome, K., Sugimachi, K., 1987. Development and regression of intimal thickening of arterially transplanted autologous vein grafts in dogs. Jour nal of Vascular Surgery 5, 719730. N'Dri, N.A., Shyy, W., TranSonTay, R., 2003. Computati onal modeling of cell adhesion and movement using a continuumkineti cs approach. Biophysical Journal 85, 22732286. Neumann, J. v., 1966. Theory of SelfRepro ducing Automata. Edited and completed by Burks, A.W.. University of Illinois Press, Urbana. O'Sullivan, D., 2001. Graphcellular aut omata: a generalised discrete urban and regional model. Environment and Planning B 28, 687705. Owens, C.D., Wake, N., Jacot, J.G., GerhardHerman, M. Gaccione, P., Belkin, M., Creager, M.A., Conte, M.S., 2006. Early biomechanical changes in lower extremity vein grafts Distinct temporal phases of remodeling and wall stiffness. Journal of Vascular Surgery 44, 740746. Palsson, E., 2008. A 3D model used to expl ore how cell adhesion and stiffness affect cell sorting and movement in multicellular systems. Journal of Theoretical Biology 254, 113. Pappu, V., Bagchi, P., 2008. 3D computatio nal modeling and simulation of leukocyte rolling adhesion and deformation. Computers in Biology and Medicine 38, 738753. PAGE 113 113 Park, E.Y.H., Smith, M.J., Stropp, E.S., Snapp, K.R., DiV ietro, J.A., Walker, W.F., Schmidtke, D.W., Diamond, S.L., Lawrenc e, M.B., 2002. Comparison of PSGL1 microbead and neutrophil rolling: Microvillus elongation st abilizes Pselectin bond clusters. Biophysical Journal 82, 18351847. Prez, M.A., Prendergast, P.J., 2007. Randomw alk models of cell dispersal included in mechanobiological simulations of tissue di fferentiation. Journal of Biomechanics 40, 22442253. Peskin, C.S., 1977. Numerical analysis of blood flow in the heart. Journal of Computational Physics 25, 220252. Phuphanich, S., Levin, V.A., 1985. Bioavaila bility of Bromodeoxyuridine in Dogs and Toxicity in Rats. Cancer Research 45, 23872389. Piotrowska, M.J., Angus, S.D., 2009. A quant itative cellular automaton model of in vitro multicellular spheroid tumour growth. J ournal of Theoretical Biology 258, 165178. Qutub, A.A., Popel, A.S., 2009. Elongation, proliferation & migr ation differentiate endothelial cell phenotypes and determine capillary sprouting. BMC Systems Biology 3, 13. Richter, Y., Groothuis, A., Seifert, P., Edelm an, E.R., 2004. Dynami c flow alterations dictate leukocyte adhesion and response to endovascular interventions. Journal of Clinical Investigat ion 113, 16071614. Robertson, S.H., Smith, C.K., Langhans, A.L., McLinden, S.E., Oberhardt, M.A., Jakab, K.R., Dzamba, B., DeSimone, D.W., Papin, J.A., Peirce, S.M., 2007. Multiscale computational analysis of Xenopus laevis morphogenesis reveals key insights of systemslevel behavior. BMC Systems Biology 1, 46. Roger et al., 2011. Heart Disease and Stroke Statistics 2011 Update. Circulation 123, e18e209 Rosamond et al., 2007. Heart Disease and Stro ke Statistics 2007 Update. Circulation 115, e69e171. Schaller, G., MeyerHermann, M., 2007. A modelling approach towards epidermal homoeostasis control. Journal of Theoretical Biology 247, 554573. Schwartz, L.B., O'Donohoe, M.K ., Purut, C.M., Mikat, E.M. Hagen, P.O., McCann, R.L., 1992. Myointimal thickening in experim ental vein grafts is dependent on wall tension. Journal of Va scular Surgery 15, 176186. Schwartz, S.M., 1999. The Intima: A New Soil. Circulation Research 85, 877879. PAGE 114 114 Shao, J.Y., TingBeall, H.P., Hochmuth, R.M., 1998. Static and dynamic lengths of neutrophil microvilli. Proceedi ngs of the National Academy of Sciences U S A 95, 67976802. Shuhaiber, J.H., Evans, A.N., Massad, M.G. Geha, A.S., 2002. Mechanism and future directions for prevention of vein graft fa ilure in coronary bypass surgery. European Journal of Cardiothoracic Surgery 22, 387396. Shyy, W., Francois, M., Udaykumar, H. S., 2001. Moving boundaries in microscale biofluid dynamics. Applied Mechanics Review 5, 405453. Shyy, W., Kan, H.C., Uda ykumar, H.S., 1999. Interact ion between fluid flows and flexible structures. In: Shyy, W., Na rayanan, R. (Eds.), Fluid dynamics at interfaces. Cambridge Universi ty Press, Cambridge, UK. Silverthorn, D.U., Ober, W. C., Garrison, C.W., 2001. Human physiology an integrated approach. Prentice Hall, Upper Saddle River. Simpson, M.J., Merrifield, A., Landman, K. A., Hughes, B.D., 2007. Simulating invasion with cellular automata: Connecting ce llscale and populationscale properties. Physical Review E 76, 021918. Smith, M.J., Berg, E.L., Lawr ence, M.B., 1999. A direct com parison of selectinmediated transient, adhesive events using high tempor al resolution. Biophysical Journal 77, 33713383. Sumpio, B.E., Banes, A.J., 1988. Response of porcine aortic smooth muscle cells to cyclic tensional deformation in culture. J ournal of Surgical Research 44, 696701. Sumpio, B.E., Banes, A.J., Link, W.G., Johnson, Jr., G., 1988. Enhanced collagen production by smooth muscle cells duri ng repetitive mechanical stretching. Archives of Surgery 123, 12331236. Tang, J., Ley, K.F., Hunt, C.A., 2007. Dynami cs of in silico leukocyte rolling, activation, and adhesion. BMC Systems Biology 1, 14. Thorne, B.C., Bailey, A.M., DeSimone, D.W., Peirce, S.M. 2007. Agentbased modeling of multicell morphogenic processes during development. Birth Defects Research C Embryo Today 81, 344353. Thorne, B.C., Bailey, A.M., Peirce, S.M., 2007. Combining exper iments with multicell agentbased modeling to study biologic al tissue patterning. Briefings in Bioinformatics 8, 245257. PAGE 115 115 TranSonTay, R., Hwang, M., Garbey, M., Jiang, Z., Ozaki, C.K., Berceli, S.A., 2008. An ExperimentBased Model of Vein Gra ft Remodeling Induced by Shear Stress. Annals of Biomedical Engineering 36, 10831091. Valentn, A., Humphrey, J.D., 2009. Eval uation of fundamental hypotheses underlying constrained mixture models of arterial growth and remodeling. Philosophical Transacionts of the Roya l Society A 367, 35853606. Varty, K., Allen, K.E., Bell, P.R.F., London, N.J.M., 1993. Infrainguinal vein graft stenosis. British Journal of Surgery 80, 825833. Walker, D.C., Southgate, J., Hill, G., Holcom be, M., Hose, D.R., Wood, S.M., Mac Neil, S., Smallwood, R.H., 2004. T he epitheliome: agentbased modelling of the social behaviour of cells. BioSystems 76, 89100. Walpola, P.L., Gotlieb, A.I., Cybulsky, M.I., Langille, B. L., 1995. Expression of ICAM1 and VCAM1 and monocyte adherence in arteries exposed to altered shear stress. Arteriosclerosis, Thrombosis, and Vascular Biology 15, 210. Wentzel, J.J., Kloet, J., Andhyiswara, I., Oo men, J.A.F., Schuurbier s, J.C.H., de Smet, B.J.G.L., Post, M.J., de Kleijn, D., Pasterka mp, G., Borst, C., Slager, C.J., Krams, R., 2001. Shearstress and wallstress regulation of vascular remodeling after balloon angioplasty : Effect of matrix metalloproteinase inhibition. Circulation 104, 9196. Westerband, A., Crouse, D., Ric hter, L.C., Aguirre, M.L., Wi xon, C.C., James, D.C., Mills, J.L., Hunter, G.C., Heim ark, R.L., 2001. Vein adaptation to arterialization in an experimental model. Journal of Vascular Surgery 33, 561569. Wolfram, S., 2002. A New Kind of Sc ience. Wolfram Media, Champaign. Yamamoto, M., AcevedoDuncan, M., Chalfant, C.E., Patel, N.A., Watson, J.E., Cooper, D.R., 2000. Acute glucoseinduced downr egulation of PKCbeta II accelerates cultured VSMC proliferation. American Journal of Physiology Cell Physiology 279, 587595. Yang, C., Tang, D., Liu, S.Q., 2003. A mult iphysics growth model with fluidstructure interactions for blood flow and restenosis in rat vein grafts A growth model for blood flow and restenosis in grafts. Computers and Struct ures 81, 10411058. Zhang, L., Athale, C.A., Deisboeck, T.S., 2007. Development of a threedimensional multiscale agentbased tumo r model: Simulating geneprotei n interaction profiles, cell phenotypes and multicellular patterns in brain cancer. Journal of Theoretical Biology 244, 96107. PAGE 116 116 Zhang, L., Hagen, P.O., Kisslo, J., P eppel, K., Freedman, N.J., 2002. Neointimal hyperplasia rapidly reaches steady stat e in a novel murine vein graft model. Journal of Vascular Surgery 36, 824832. Zhang, L., Wang, Z., Sagotsky, Z.A., Deisboeck, T.S., 2009. Multiscale agentbased cancer modeling. Journal of Ma thematical Biology 58, 545559. Zohdi, T. I., 2005. A simple model for s hear stress mediated lumen reduction in blood vessels. Biomechanics and Modeling in Mechanobiology 4, 5761. Zou, Y., Dietrich, H., Hu, Y., Metzler, B., Wick, G., Xu, Q., 1998. Mouse Model of Venous Bypass Graft Arteriosclerosis. American Journal of Pathology 153, 13011310. Zwolak, R.M., Adams, M.C., Clowes, A.W., 1987. Kinetics of vein graft hyperplasia: Association with tangential stress. Jour nal of Vascular Surgery 5, 126136. PAGE 117 117 BIOGRAPHICAL SKETCH Minki Hwang received his bachelors and masters degrees in aerospace engineering from the Ko rea Advanced Institute of Scienc e and Technology (KAIST) in 1995 and 1997, respectively. He worked as a research and development engineer for Samsung Techwin Corporation from 1997 to 2002. He received a masters degree in bioengineering from Penn State University in 2004. 