Processing of MRI data for simulation and monitoring of drug delivery

The work in this thesis centers around monitoring and simulation of a novel drug delivery system. The major features are the development of a pipeline for creating realistic tetrahedral geometries from MRI data suitable for finite elements simulation, research into methods for quantifying the impact of the shape of an injected drug depot on drug diffusion, and research into automated methods for segmenting hyperplasia growth from MRI images for monitoring the efficacy of the drugs and delivery methods being tested. The error quantification portion of this work exploits assumptions of smoothness in order to employ the stochastic collocation method to study the effect of drug depot shape variability on the outcome of drug diffusion simulations. Simple parameterizations of depot shape are also explored. We demonstrate that once the underlying stochastic process is characterized, quantification of the introduced variability is quite straightforward and provides an important step in the validation and verification process. The hyperplasia segmentation method described here attempts to exploit minimal user interaction to semi-automatically produce an accurate segmentation. The challenge arises from the need to track a combination of strong and weak features near confounding structures. The 3D segmentation process is broken down into a series of constrained 2D segmentation tasks, which are carried out automatically on preprocessed image slices by an active contour model.

by J. Samuel Preston A thesis submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Master of Science School of Computing The University of Utah August 2009 PROCESSING OF MRI DATA FOR SIMULATION AND MONITORING OF DRUG DELIVERY by t he iJl fC<luircmcnts of i\llastcr COIIIPuting © J. Copyright © Samuel Preston 2009 All Rights Reserved J. Preston and satisfactory. -A a 4- THE UNIVERSITY OF UTAH GRADUATE SCHOOL SUPERVISORY COMMITTEE APPROVAL of a thesis submitted by J . Samuel Preston This thesis has been read by each member of the following supervisory committee and by majority vote has been found to be satisfactory. Chair: Tolga Tazdizen 5/&/01 r ~ 1 Robert M. Kirby .....------.... .;;/&/0'1 ( Guido trig q 0 SCHOOL Utah: J. format) Graduate hlkh') : Chair, Supervisory Committee Department Chair /David S. Chapman THE UNIVERSITY OF UTAH GRADUATE SCHOOL FINAL READING APPROVAL To the Graduate Council of the University of Utah: I have read the thesis of J. Samuel Preston in its final form and have found that (1) its format, citations, and bibliographic style are consistent and acceptable; (2) its illustrative materials including figures, tables, and charts are in place; and (3) the final manuscript is satisfactory to the Supervisory Committee and is ready for submission to The Grad uate School. Date Tolga Tazdizen Approved for the Major Department Martin Berzins Chair/Dean Approved for the Graduate Council D~vid S. Chapman ~ . Dean of The Graduate School The work in this thesis centers around monitoring and simulation of a novel drug delivery system. The major features are the development of a pipeline for creating realistic tetrahedral geometries from MRI data suitable for finite elements simulation, research into methods for quantifying the impact of the shape of an injected drug depot on drug diffusion, and research into automated methods for segmenting hyperplasia growth from MRI images for monitoring the efficacy of the drugs and delivery methods being tested. portion order to employ the stochastic collocation method to study the effect of drug depot shape variability on the outcome of drug diffusion simulations. Simple parameterizations of de­pot shape are also explored. We demonstrate that once the underlying stochastic process is characterized, quantification of the introduced variability is quite straightforward and provides an important step in the validation and verification process. The hyperplasia segmentation method described here attempts to exploit minimal user interaction to semi-automatically produce an accurate segmentation. The challenge arises from the need to track a combination of strong and weak features near confounding structures. The 3D segmentation process is broken down into a series of constrained 2D segmentation tasks, which are carried out automatically on preprocessed image slices by an active contour model. ABSTRACT arc r..'IRJ eicrnents a ll d iffusion, resemch seglnenti ng Mill drngs allel delivt:ry The error quantification port ion of this work exploits assumptions of smoothness in t he si mula tions. de­po~ a rc cxplored. dClllonstrate t hat t be process dIal"acterized , aud process. t.o a utomatically prodllce segment.ation. t.be ueed t.o combi nation stwng st.ruct.ures. Thc bwken int.o arc automat.ically s lices all act ive To my family, for freedom and encouragement To my friends, for inspiration and distraction To my Liz, for love and support To CsHioN 10-2 and ( 1 . for all the ups and downs encouragement inspi ration lll}' CSHlON40 1 C2H;:,OH. CONTENTS iv FIGURES 1. 1 Introduction Segmentation Geometric Discretization 3. 9 3.1 Introduction 9 3.2 Model Equations 9 Numerical Approximation 4. OF ERROR DUE TO DEPOT SHAPE 12 4.1 Introduction 12 4.2 Modeling Uncertainty in Depot Shape 14 4.3 Characterization of Object Variability 4.3.1 Level-set Shape Representation 4.3.2 Parameterization 1: Mean Curvature Flow 21 4.3.3 Parameterization 2: Metamorphosis 22 4.4.1 Domain Resampling 4.4.2 Total Drug Delivery 4.4.3 Nodal Statistics 27 Discussion 34 5.1 Introduction 5.2 Slice Plane Selection and Contour Initialization 36 Data Preprocessing 37 5.4 Active Contours 5.5 Volume Calculation Results CONTENTS ABSTRACT ............................ . .. .. ..... . ...... . .. .. ... . LIST OF FIGURES viii CHAPTERS INTRODUCTION .................. • . •• .................... 2. MODEL GENERATION . . . . . . . .. . .... . . . . . .. .... .. . . . . . . ... . 4 5. 2.1 Introduction..................... . . .. . .. .................. 4 2.2 Image Segmentation. . . . . . . . . . .. ......•.... .. . . . . . . . . . . . . . . . . . 4 2.3 Gcometl'ic Oiscrctiz.ation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 DIFFUSION SIMULATION ................ .. .... .. .. . . ....... . Introduct.ion ................... • .••• •• • .... . . .. . . ...... . . ... r.,'lodcl ................... .... . . 3.3 Nml'lcricai Approximat.ion ..... . QUANTIFICATION OF ERROR DUE TO DEPOT SHAPE ..... . ............................................... . r.. lodcling Unccnaiuty ill Shupe ......... ................ . ....................... . . . sct Reprcscntation .. . ........ ........... . Pan,uueterizatioll I: i\leun ...... .......... . . . l\ lcL<tI liorphosis . .................... • . 4.4 Results ......... .... ................................. . ........... • . ..•.. .............. Tolal ...........• •• • •••• ............. '1.4.3 ~odal ....... .... .......•••••••...... .... ... 4.5 .......................... . .• ... . ...... . . .. ... .. HYPERPLASIA SEGMENTATION ................. ........ . . 5.3 5.6 a nd .....•............... PrcproCt-"Ssing ............ . ............................. . .•............ Volullle Calculat.ion ....................... • • • ••.•............. ............. . 9 12 12 19 19 23 24 24 29 34 39 42 42 6. CONCLUSION ............. ...... ............................. 47 REFERENCES ........•.•• . •.. ...... .. .... ...... . . .. .. . ....... . 49 vii 1.1 Diagram of standard arteriovenous graft 2 1.2 Porcine Model Procedure 3 2.1 Segmentation and discretization pipeline 4.1 Shape statistics using stochastic collocation 18 4.2 MCF gel morph 22 4.3 Metamorphosis gel morph 4.4 Domain resampling results 4.5 Venous tissue drug delivery region 4.6 results, MCF 27 4.7 Total drug delivery results, morph parameterization Total 4.9 ROI MCF results morph statistics parameterization Hyperplasia contour 38 Possible contours for initialization points 5.4 Data 39 43 5.7 Qualitative comparison of segmentations, automatic vs. hand segmentation of anisotropic data 45 5.8 Qualitative comparison of segmentations, automatic vs. hand segmentation of isotropic data 46 LIST OF FIGURES Diagram graft. . . . . . . . . . . . . . • . . . . . . . . . . . . 2 Porcine r-.lodcl Procedure.......... .. . . . . . . .... . ••. •...... . Scgmclltntion d iscretization ..................... . slatistics 1... ICF ............... . MctlUtlorphosis rcsampling ...................... . t issne ..... '1.6 Total drug delivery results. !-.ICF parameterization ............... . Illorph ............ . 4.8 Tolf.II drug delivery ROI ...........................•............. 4.n Total drug delivery, ROJ ~'I C F resl\lts .... 4.10 Total drug delivery, ROI ltIorph results ................ . 4.11 Higher-order moments ... 4.12 Nodal stat is tics results .... ............. . 4 18 22 23 25 26 28 28 29 30 31 32 4.13 Surface area of shapes from shape parameterization. .....••••. ........ 32 5.1 Hyperpl asia growth and imaging ........ . 35 5.2 User contoUl' initialization ...................................... . 38 5.3 rossH,]. contoms fm the same ;n;t;,];,ation po;n" . . . . . . . • • • • . . . . . . . .. 38 5,4 Dnla preprocessing pipeline .............. . 5.5 Contour forces 5.6 The segmentation pipeline ................ . scgmentntions. segmentation 39 43 44 an isotropic data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. Qunlitative nutomatic haud segmentation data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. Over the past two decades, computational techniques have become an indispensable part of medical research. Demand for innovative and improved procedures and phar­maceuticals will always have a particular urgency, as improvements in medicine direct ly lead to longer and higher quality lives. However, procedures used in medical research must be thoroughly vetted to ensure the safety of study participants and the validity of results. When properly applied, computational methods can increase accuracy and reproducibility, and decrease time and cost of studies. This for delivering antiproliferative drugs to dialysis patients with implanted arteriovenous grafts. Repeated or long-term catheterization for dialysis can cause scarring and stenosis of vessels, and therefore, patients undergoing chronic dialysis treatment can quickly render vessels unsuitable for vascular access. longer term solution is to implant an arteriovenous graft that will be used for access (Figure 1.1(a)). As this is a foreign object, the body responds by creating a proliferation of smooth muscle cells, known as neointimal hyperplasia (Figure 1.1(b)). This cellular proliferation is most pronounced at the opening from the graft to the vein (the venous anastomosis). Over time, this tissue causes stenosis of the graft, and is the primary reason for failure of the vascular access [26]. Antiproliferative drugs such as rapamycin are effective at slowing the growth of these cells but treatment is limited to a coating at the time of implantation. The proposed delivery method uses a triblock copolymer gel. liquid at 4° C but a gel at body temperature, as an injectable "drug depot" from which the active drug can naturally diffuse over time. A porcine model with an arteriovenous graft implanted between the common carotid artery and the ipsilateral external jugular vein was used in the study. The procedure is shown in Figure 1.2. Grafts are implanted on both sides of the animal, although CHAPTER 1 INTRODUCTION t.echniques illdispensable innovat ive im proved proc(.'dures und irnprovements directly qualil,y procedttres llsed thorollghly st,udy \Vhen applied , ami T his thesis centers around work done for a study aimed at developing a novel technique anti proliferative arteriovenolls graflH. catheterizat.ion fOI" u nsuitable acceSH. A implaut lISed acc<.'Ss 1.1 (a». AH t his musele ncoint.imal hypcrplasia l.l(b». T his Ovcr t his t issue t he [261. Antiprolifcrative us raprunycin slowi ng of [16], il t ime lnethod lIseS gel, 'drug depot' tillle. ~he arc implalltt.'d 2 Figure 1.1. Diagram of standard arteriovenous graft placement for hemodialysis (a), and diagram of common location for hyperplasia growth (b) only one graft is shown in the figure. Drug and gel are injected at weekly intervals following surgery. The drug is injected outside the vascular structures near the venous anastomosis. Hyperplasia growth is measured at the venous anastomosis to test drug efficacy. MR imaging is performed at the time of gel injection to verify the location of internal structures, including the implanted graft and injected gel. and to monitor the growth of hyperplasia within the graft. The primary contributions of this work to the project were as follows: • The development of a pipeline for segmenting structures of interest from MR data element • Research into novel methods of quantifying error due to drug depot shape variability • Development of methods and software for semi-automatically segmenting hyperpla­sia ( (a) Artery Hyperplasia (b) a). fib'lJl'C. arc Hy perplasia Ivl R Ilt tillle of gel, II ]VIR and producing tetrahedral meshes suitable for use in finite clement simulations eHor hyperpla-sia Graft Hyperplasia Figure 1.2. Porcine Model Procedure: An AV graft is implanted between the common carotid artery and the ipsilateral external jugular vein (top). Gel and drug are injected near the graft site at regular intervals (middle). Hyperplasia growth is monitored at the venous anastomosis to measure drug efficacy (bottom). 3 Artery Vein Graft Artery Vein F igure Porci ne l\lodcl (;urotid externa l jUb'ldar arc Ilcar moni tored anastomos is bottom) . A major goal of the drug delivery study is to develop simulation methods that accurately reflect the in vivo results. Towards this end, this chapter describes a pipeline for quickly generating tetrahedral models of the structures of interest directly from MR images. This includes the vein, artery, and graft walls, as well as the drug depot. Figure 2.1 depicts this processing pipeline for image segmentation and mesh generation of the vascular structures, which is explained below. 2.2 Image Segmentation Lumen and depot geometries were constructed from MR images taken at weekly intervals following surgery. Gel injection was performed at these times, with pre- and postinjection scans acquired. Lumen geometry was obtained via a 'black blood' 3D Turbo Spin Echo pulse sequence. This sequence uses a nonselective inversion pulse followed by a selective reinvcrsion pulse in preparation for imaging in order to suppress signal from blood flowing into the imaging plane (Figure 2.1 (1)). Although the black blood pulse sequence creates a signal void within the lumen, inhomogeneoiis RF penetration Figure 2.1. Segmentation and discretization pipeline for lumen geometry. (1) Black blood MR image. (2) Initial segmentation from user defined seed points. (3) Disconnected segments of lumen (see artery in the left of previous image) are connected using a Minimum Cost Path (MCP) algorithm. (4) Segmentation is refined using a level set-based region growing step. (5) Surface triangulation is generated from refined segmentation. (6) Tetrahedral mesh of vessel walls is generated using an offsetting technique. CHAPTER 2 MODEL GENERATION 2.1 Introduction A major goal of the drng delivery study is to develop simulation methods that accurately reflect the in vivo results. Towards this end, this chapter describes a pipeline for quickly generating tet.rahedral models of t he structures of interest directly from MR images. This includes the vein, artery, and graft walls, as well as the drug depot. Figure 2.1 depicts this processing pi peline for image segmentation and mesh generation of the vascular structures, which is ex plained below, geometri\~s wcre constrllctt."<] ~Ht lakclI injectioll times. postiujection i:lcquired. gcometry 1\u'bo T his uscs select ive reinversion prepa ration imagillg 1». lumeu, inhomogeneolls penetration seg~!~~tion I MCP f.1 Level Set HSUrface MeShH Gen!;!tion I [t1J[EL(]~[E21J~ pipcline .s~:.'(1 Disconnected IUllIcn ().IICP) Segmcntation rcfined gcncratL'(lusing 5 causes a gradient in the signal across the structures of interest. A simple thresholding is therefore not adequate for segmentation. An initial rough segmentation of the lumen was created using user-defined seed points within the lumen as inputs to a, connected component / thresholding algorithm (Figure 2.1 (2)). Due to the extreme narrowing of the vessels near the anastomosis and the common presence of flow and pulsatile motion artifacts in this region, a Minimum Cost Path (MCP) algorithm is necessary to complete the lumen connection through this area (Figure 2.1 (3)). This technique is similar to the stenosis-robust segmentation technique described in [31]. Starting from the initial segmentation, this algorithm defines a graph in which voxels are considered nodes and each voxel has an edge to each immediate neighbor (6-connected). Since lumen should be the darkest region in this black blood scan, the cost associated with each edge is based on the intensity of the neighboring pixel. A search is then performed to find the MCP connecting the disconnected segmentation regions. An active contour-based region growing step (using the edge-based method implemented in the SNAP [36] software package) was then used to refine and 'fill in' the segmentation (Figure 2.1 (4)). Gel images were taken immediately following injection using a specially developed pulse sequence for depot imaging [25]. Due to the improved contrast introduced by this sequence, a simple thresholding algorithm was sufficient for segmentation. The segmented volume agreed well with the known injected volume - taking a maximum and minimum 'reasonable' threshold level changed the segmented volume to approximately ±20% of the known volume. The threshold level used was chosen such that the segmented volume agreed with the known injected volume. 2.3 Geometric Discretization Once segmented, geometric surface models of the structures of interest were created and embedded in a control volume. This multisurface structure is used to define material regions and voids in a surface-conforming tetrahedralization of the control volume. Since the imaged region is smaller than the desired control volume, the vessels are extended to the boundaries of the control volume. The vessels travel in roughly the z-direction within the imaged region, so this extension is accomplished by simply repeating the boundary z slices of the segmentation to the extents of the control volume. A surface mesh of the lumen segmentation was created by slightly blurring the segmen­tation with a Gaussian kernel, then using this gray level image as input to an isosurface Ii !:lignal Ii (}lICP) (t.etrahedralization t.he repeating 6 meshing algorithm. The algorithm used is part of the CGAL library [4], and uses iterative point insertion on the isosurfacc to create a triangle surface mesh with constraints on the maximum size and minimum angle of any member triangle (Figure 2.1 (5)). Each of the vascular structures (vein wall, artery wall, and graft) have a different thickness and different diffusive properties. The portions of the lumen contained in the different structures cannot be clearly distinguished in the available MR images, so the generated surface mesh of the lumen was hand-annotated with this information based on knowledge of the anatomical structures. The resolution of the MRI sequence used to image the lumen was not high enough to directly image the vessel walls or graft. In order to create the geometry of the graft and vessel walls, thickness measurements from histology were used. Using the labelled surface dataset. a tetrahedral model of the vascular structures was generated by duplicating the surface mesh and offsetting each duplicate mesh point in the surface normal direction by a distance equal to the thickness of the vascular structure (See Table 2.1 for thickness values). Connecting each vertex to its offset duplicate produces a mesh of triangular prisms. The maximum triangle size of the surface mesh was chosen with respect to the thickness of the vascular structure in order to produce well-shaped prisms. This mesh is then subdivided such that each prism becomes three tetrahedra and all tetrahedra have coherent faces. Coherent prism splits are computed using the 'rippling' algorithm described in [18]. The tetrahedra generated by this process then have maximum edge lengths roughly equal to the thickness of the structure they represent. 0.3 to 0.7mm. Although the simple offsetting technique used has the possibility of producing intersecting or incorrectly oriented tetrahedra. experimental results show that the curvature of our surface is low enough at all points relative to the offset distance to avoid producing any invalid elements (Figure 2.1 (G)). The surface mesh of a specific drug depot shape is generated from the output of our parameterized depot shape generation algorithm by the same method as used for the lumen surface. algorit.lnn. librar), 41. uscs isosu rfacc mf\..'(imuill t riangle (5)). waH, T he port ions disti uguished t.he hand~al\notatcd all of st ructures. T he of (viRI scquCll(,:e gcornctl'y amI measmements histolo!,,), wcrc labeled dataset, tetraht.-Gral Wlscular wa.':> gellerHted duplicHting offsel.ting smfuce di stance t hickness strllct\ll'C maxinmm of t he sllch tetra hedra hHVC are cornputed using t he ' rippling' algorithm described in [18J. The tetrahedra generated by this process then have maximum edge lengths roughly equal to the thickness of the structure they represent, 0.3 to O.7mm. Although the simple offsetting technique used has the possibili t.y of producing intersecting or incorrect.ly oriented tet rahedra, experimental results show that. the curvature of our surface is low enough at all points relative to the offset dist.ance to avoid producing any invalid clements (Figure 2.1 (6)). The surfRce Hlc.'Sn of a specific drug depot shape is generated from the output of our parameterized depot shape generation algorithm by the same method as used for the lumen surface. 2.1. Gel XVT1 Diffusion (mm2/1 0 - ' 7.80 x 10~8 10-" 10-" 1.0 x 1 0 _ o Struct tire Thickness (nun) 0.6 mm3) 0.0 ( ) . ( ) 1 0 - 9 Table 2 .1 Physical Region Attributes Vein Artery Graft NVT Di ffusion Coefficient (mm:t/ s) 4.09 x 10 7.80 X to- II 2.91 x 10-0 3.89 x 10- ' 1.0 x 10- > Structure mm) 0.3 0.7 n/a n/a Init. Drug Concentration (mol/ mill ) 0.0 0.0 2.73 x 10 0.0 8 The depot surface and the outer (offset) lumen surface are embedded in a cube representing the extent of our test volume. These surfaces are combined and used as the input Piecewise Linear Complex PLC) to tetgen [22]. The cube and the offset lumen surface define the boundaries of the PLC (i.e., the offset lumen surface represents a 'hole' in the volume), and the elements of the depot surface mesh are used as internal constraining faces of the generated Constrained Delaunay Tetrahedralization (CDT), defining the drug depot region. The area outside of the depot in this CDT is considered nonvascular tissue. Elements generated are constrained to have a maximum radius-edge ratio of 2.0. allowing element size to grow away from the structures of interest. Finally, the previously computed tetrahedral mesh of the vascular structures is merged with this tetrahedral volume to create our final test domain containing five regions the graft, arterial wall, venous wall, drug depot, and nonvascular tissue. The tetrahedralization step is constrained to preserve boundary faces along the offset lumen surface, so this merging process requires only that duplicate points be merged. This process yields a discretization containing approximately 250.000 elements, varying somewhat depending on the specific structure geometries. Nonvascular (PLe) 22J. T he und huncn su rface bOlUldarics PLe (i.e. , the offset lumen surface represents a uud elcmcllts arc tL',ed constrai ning Constraiued Delauuay Tetrahedralizat ion COT), d rug t his COT nOllvascular tissue_ Element-s bave radius--ratio of 2.0, allowing element si;-.c to grow away from I.he st.rllctures of interest. Finally, the previously computed tetrahedral mesh of the vascular structures is merged wit.h Lhis tetrahedral volume to create our final test domain containing five regions - the graft, ,-wtcrial wall, venous wall, drug depot, and nonvascular t issue. The tetrahedralizatioll step is constrained to preserve boundary faces along the offset lumen surface, so t.his merging process requires only t.hat. duplicate points be merged, This process yields a discretization conta ining approximately 250,000 clements, varying somewhat depending on the specific structure geometries. ! Nonvn.scular Tissue This chapter describes software developed to solve a standard diffusion equation over the discretized tetrahedral mesh using linear finite elements. Although commercial packages are available for this purpose, development of in-house software allows explicit understanding of the assumptions and methods used in simulating drug diffusion, and allowed this software to be developed as a component in the SCIRun Problem Solving Environment [1] in order to facilitate its use in a meshing / simulation pipeline. For the purposes of the simulation, the software assumes that drug transport is due only to diffusion. Furthermore, the implementation assumes that diffusion rates are isotropic and constant for each material; let be a function that is piecewise-discontinuous and consists of diffusion coefficients D,„ associated with each material m occupying the domain Jim c ft, where Q represents the control volume. The drug concentration c at position x and time t is then governed by the equation: ^ | ^ = V-{D(x)Vc(xU)).xen (3.1) where x) = Dm for € iim, and where appropriate boundary conditions are applied. The assumption that any drug diffusing in'" the interior <>l the vascular structures is immediately removed by bloodflow leads to a known (zero) concentration at the lumen boundary (zero Dirichlet conditions). Drug is disallowed from diffusing across the outer control volume boundary. 3.3 Numerical Approximation Solving Equation 3.1 for any but the simplest geometries requires a numerical approx­imation, typically based on methods such as boundary elements or finite elements. The CHAPTER 3 DIFFUSION SIMULATION 3.1 Introduction This chapter describes software developed to solve a standard diffusion equat.ion over the d iscretized tetrahedral Illesh using linear finite clelllents. Although commercial pac&'lges are available for this purpose, development of in-honse software allows explicit understanding of the assumptions and methods used in simulat ing drug diffw;ion. and allowed this software to be developed as a component in the SCIRull Problem Solving EnvirOlllucut [1] in order to facilitate its LIse in a meshing / simulation pipeline. 3.2 Model Equations for flSStmlCS t.hat arc D piecewise­discontinuous Dm Om C f! , n 8c(axt'. t) ~ V · ( D(x)Vc(x.t)), XE fl (3.1 ) D(x ) Dm x E Om, condit ions an.y into of st ructures irJ'lmediately blood flow concentrat.ion boulldary n .'.< luires Ilumerical ba.<;ed stich 10 software- described uses a standard finite clement formulation with piecewise linear basis functions. For the domain Q representing our control volume, let S1) be the tetrahedralization of the volume as described in Section 2.3. Also, let the set M contain the indices of the vertices of T(il), which are the nodes of our finite element mesh. The set is then decomposed into tun nonintersecting sets. B and I. representing nodal indices thai lie on the boundary and nodal indices of the interior, respectively. B is further decomposed into Bq, the outer boundary of the test volume, and Bc, the lumen boundary. To satisfy our requirements from Section 3.2. we impose zero Neumann boundary conditions at Bq and zero Dirichlct boundary conditions at Bc. 4>i(x) denote the global finite element interpolating basis functions, which have the property that <pi(Xj) = <5jj where Xj denotes a node of the mesh for j £ Af. Solutions are then of the form: u(x) = $^«f c<MaO (h A" = 5^wf e ( M x ) + $ 3 (3.3) k&t k€Bc where U = X\jB,,. which represents the unknown values of the problem. The first term of Equation 3.3 then denotes the sum over the degrees of freedom (the values at the unknown vertices weighted by the basis functions) and the second term denotes the same sum for the (known) Dirichlet boundary conditions of the solution. Substituting the expansion (Equation 3.3) into the differential equation (3.1). multiplying by a function from the test space {4>j(x);j EU}. taking inner products, integrating by parts, and applying boundary conditions yields a linear system of the form: M | f + KtZ = / (3.4) where il denotes a vector containing the solution of the system (i.e., concentration values at the nodal positions in U) and M. K and / denote the mass matrix, stiffness matrix, and right- hand-side function, respectively, given by the following expressions: M = Mjk = (<)>s>4>k) (3.5) K = Kjk = ( V 0 j . D • V d / , ) . (the expressions above, j,k € U, and (•, •) denotes the inner product taken over the entire spatial domain ft. For the described conditions, all iii for G Bc are zero and the flux on the exterior boundary was taken to be zero, so / , 0 Wj. software Il fmitc element Iincar loaf n tet T(O) tctmhedraliz<ltion N T (f2}, arc t he uodes N decol'nposed into two lIonintersecting sets, and I , representing: nodal indices that lie on the boundary and nodal indices of the interior, respectively_ is further decomposed into Sq , the outer boundary of the test volullIe, and Be. the lumen boundary. To S<:'\tisfy Ollr requirements from Section 3.2, we impose zero Koumann boundary conditions at Bq and zero Dirichlet boundary conditions at Be> Let ¢i(X) denote t be global finite elemellL interpolating basis functions, whicb have propcl'ty ¢i(Xj) oij X j denott:s tuesh E N. arc then of t he form: . (x ) ~ L u,~, (x) k EN (3.2) ~ 2:>,¢,(x) + L ",¢,·(x) , (3 .3) kEU k € Bc wherc U = I u B!!. which represents the unknown valucs of the problem. The first term of Equntiolt 3.3 thcn denot<.:S the stlln over the degrees of freedom (Lhc values at thc unknown vertices weighted by t he basis fnnctions) and the st'Con(\ tenn denotes the same sum for the (k uowll) Dirichlet boundary conditions of the solution. Substituting the expansion (Equation 3.3) into I,he difl:"erential equation (3.1 ), multiplying by a fUllction from the test space {¢j(x):.i E U}, taking iuner products, integrating by parts, and applying boundary conditions yields a linear system of the form: MNail + K u -- I­( 3 .4) 'ii vt.'Ctor cOlltailling (i .e .. conc~ntrat i on Ht M . j haud-side function. M ~ Mj' ~ (¢j,¢k) f{] ~. Vtbj, · Vtbd. (3.5) (3.6) In j ,k E U , (-, .) t he n. 'Uj ·i E Be zero. Ij = Vj. The time-dependent equation (3.4) can be discretized using the Crank-Nicolson scheme. Arranging the unknown (future timestep) values on the left-hand side yields the timestep-ping equation: (M - fK) u"+l = (M + fK) u" (3.7) At timestep, represents timestep n. this scheme require an initial condition x, t = 0), which in this case represents the time of injection with a known, constant drug concentration within the depot and /cn> concentration elsewhere. 11 dependent discrctizcd Crank-)Jicolson timcstcp) hand s ide tilllcstcp­ping (3.7) where 6.t is the timcstcp, and u" t'Cprescnts the solution at timcstep Solutions using c(x. 0). t his t he ;\ known. zero 4.1 Introduction In any numerical simulation, assumptions arc made that do not precisely reflect the physical reality of the problem. These simplifications yield tractable problems, but also introduce error. As part of the traditional simulation validation and verification (V&V) process, it is important to understand the impact of a given assumption on the final result. Assuming that a given modeling decision represents a specific choice from a set of possible decisions, it is possible to quantify the impact of that choice on the simulation. Knowing the variability introduced by the range of possible choices allows confidence in the final result to be adjusted accordingly. Probabilistic characterization of the effect of modeling decisions on simulation results is a powerful tool for validation, but is generally prohibitively costly due to the number of simulations that must be run to obtain a reasonable confidence level for the statistics. To quantify the impact of a given modeling decision using a standard Monte Carlo approach, individual samples (each representing a specific modeling choice) are randomly taken from the population of possible choices ac­cording to the underlying probability distribution. However, this method converges quite slowly (as \/s/N [19]). often requiring thousands of samples (and therefore thousands of simulations), which generally proves a prohibitive requirement. If instead we can model this population as arising from an underlying stochastic process, it is possible to apply the generalized Polynomial Chaos - Stochastic Collocation (gPC-SC) method [33. 32] to efficiently sample the probability space, yielding accurate statistics from a nearly minimal number of samples [23]. The Stochastic Collocation method can be seen as a "sampling' extension to generalized Polynomial Chaos, which represents the stochastic process as a linear combination of orthogonal polynomials of random variables. This gPC representation has been used to solve stochastic differential equations in a number of fields [34, 35, 17, 11]. Stochastic collocation has been successfully used to quantify the CHAPTER 4 QUANTIFICATION OF ERROR DUE TO DEPOT SHAPE [n are reffect simplificatioll8 givell it poStiible resu lt pO\verful ~ullples l j .jJii 19]), If ns it pO$.<;ible the generalized Polynomial Chaos - Stochastic Collocation (gPC-SC) method [33, 32] to efficiently sample the probability space, yielding accurate statistics from a nearly minimal number of samples [23]. The Stochastic Collocation method can be seen as a 'sampling' extension to generalized Polynomial Chaos, which represellts the stochastic process a<; it linear combination of orthogonal polynomials of random variable-<;. This gPC representation has been used to solve stochastic differential equations in a number of fields [34, 35, 17, 11J. Stochastic collocation has been successfully used to quantify the 13 effects of heart position in ECG forward simulation [12], and in modeling uncertainty in diffusion simulation due to microstructure variability [10]. Since stochastic collocation builds statistics based on deterministic solutions for sampled stochastic parameter values, it requires only a standard discrete solver for the problem of interest. This makes for easy implementation and allows its use on problems with complicated governing equations for which a nonsampling gPC formulation would be difficult or impossible. The specific experiment we are simulating in this wrork attempts to test the efficacy of the antiproliferative drug rapamycin delivered via an injected triblock copolymer gel in limiting neointimal hyperplasia growth at the venous anastomosis of an implanted arteriovenous expanded polytetrafruoroethylene (ePTFE) graft. This experiment is part of a study analyzing the efficacy of such a procedure for possible future use in hemodialysis patients with implanted ePTFE grafts. The current study uses a porcine model with grafts implanted between the common carotid artery and ipsilateral external jugular vein [15]. We are specifically interested in drug transport via diffusion from the injected gel deposit containing drug (the 'drug depot') to the venous anastomosis (the 'target site'). The geometry of the vein, artery, graft, and depot are modeled from MR images and diffusion coefficients of each are experimentally determined [6] (see Table 2.1). All other structures in the test region are combined and labeled as homogeneous 'nonvascular tissue.' Additionally, the vein, artery, and graft wall thicknesses can not be measured reliably from MR images. Therefore, we use measurements obtained from histological studies. One of the modeling decisions involved in this simulation is the geometric shape of the drug depot. Although the volume of gel and amount of drug is known and constant from experiment to experiment, each injection is a nonrepeatable process that yields a different depot shape based on the properties of the local tissue and variability in the injection technique. For an individual simulation, a single depot shape must be chosen from the range of possible shapes. This choice introduces variability to the experimental outcome; hence, it is important to study. In this chapter, we propose a method for studying the effect of depot shape on the outcome of the experiment. In order to determine the variability in drug concentrations at some simulation time due to the choice of depot shape, we develop a shape parameterization to represent the underlying population of shapes. We then use the gPC-SC method to sample this population at well-chosen and in modeling uncertainty in impossible. work anti proliferative polytetrafluoroethylene [15]. \Vc :VIR sec inclividual \Ve well-14 chosen simulation time. Although the specific shape models developed in this chapter are characterized by a single parameter and are too simple to capture the entire variability observed from all injections, the contribution of this chapter is the development of a method for studying the effect of shape variability on simulation outcomes. More complex shape models with multiple parameters can be easily inserted into the proposed framework in future studies. To begin our study, we must first understand and clearly articulate the mathematical problem we are attempting to solve. We accomplish this by first describing our "exper­iment" (in the sense that a probabilist would), then by articulating a probability space that captures the salient features of the experiment about which we want to reason, and finally by describing how the stochastic collocation method can be employed to allow computation of statistics of interest. Our "experiment" of interest is as follows. We assume that we have been given a test volume containing the static structures of interest - in our case the vein, artery, and graft. These objects are modeled as remaining fixed in both position and in material properties (see Table 2.1) throughout the course of all experiments. We assume that we have also been given a parameterization of drug depot shapes such that a specific parameter value corresponds to a specific geometric depot model. Having both the static test volume and a specific depot model, diffusion simulations can bc run yielding drug concentrations for positions within the test volume at any future time, respecting assumptions regarding our simulation method. For our probabilistic experiment, we limit the set of admissible depot shapes to those represented by a univariate parameterization. In Section 4.3, we will formulate two specific parameterizations: one which generate shapes with varying regularity in surface curvature and another which generates shapes by morphing between two realistic depot shapes as segmented from a Magnetic Resonance image. Initial concentration of drug in the depot is assumed to be uniform, and the volume of all depot geometries generated by this parameterization is constant. The volume, position, and gross shape of the depot were chosen to represent a 'standard' injection result as seen by postinjection MRI. Variation in depot shape can now be expressed by assuming that any specific shape is drawn from a set consisting of shapes resulting from our parameterization. We assume that there is an parameter values and derive statistics of the effect of shape on drug concentration at our 4.2 Modeling Uncertainty in Depot Shape mathematical \Ve - arc \Ve be \Ve 15 equal probability of obtaining any specific depot shape from our parameterization; that is, the probability distribution is uniform. Using this parameterization, we aim to compute statistics of drug concentration at some time tj induced by the variations in depot shape expressed by this simple (but quantifiable) probabilistic experiment. To formalize this process, let (x, A. fx) be a complete continuous probability space that expresses variation in the depot shape where \ ' s the event space consisting outcomes corresponding to depot shapes, A C 2* the cr-algebra used to define measurable events, and fi the probability-measure expressing the distribution from which outcomes are drawn: in our case we choose a uniform distribution. We can now express the depot shape as a function of the uniform random variable £ where the length of the distribution interval represents the range of admissible depot shapes. We note that although the uniform distribution is used throughout this study, the stochastic: collocation approach can be used for any compactly-supported second-order random process. The random field of interest (and in particular, its statistical characterization) in this study is the final drug concentration at a time tj. Given that the depot shape in our experiment can be completely expressed in terms of £ (given a single-parameter diffeomorphisni of admissible shapes), and given that the drug concentration follows as a direct consequence of a depot shape, the drug concentration can be expressed as a function of £. Let the distribution of drug concentration due to £ be denoted by f(x. tf,£) where denotes points within our test volume and tj denotes the final time of interest. Mathemat­ically, the quantities we are interested in computing are given by the following expressions: mean(f)(x,tf) = E[f(x.tf,0] (4-1) var(f)(x.tf) E[(f(x.tf.O-mean(f))2} (4.2) (f(x,t/,^) - mean{f))2dfi. IV or more generally for the p"' moment: E[(f(x,if,0-mean(f))*]= (4.3) j^nx.tf.O-meanif))",!,, where mean(f) and var(f) (and other high-order moments) are themselves fields defined over our test volume giving the mean and variance of drug concentration at any x due to obtailling a ny is, probab i !i~y t his paramcterizat,ion. conccntration t f the vnriations in depot shape X, A, /-I.) be a complete continuous probability space that expn.'sses variation ill X is spru.-e c.:onsisting shape!;, c 2); t he t7 - algcbra dcfiuc Jl probability measure arc drawn; ill om we c~m of ~ t.he represents \Ve ulliform distribution stochastic collocatioll a pproach us(.'{\ any T he interes t, in particula r, its statistical characterization) in this study is the final drug concentration at tf . t.hat ~ diffeomorphism dl"llg concentration can be expressed as a function of~. ~ deuoted 1(x , tf ' ~) x \vithin a nd t f fi nal j\'lathemat­ically, t.he intcre;ted fol lowing mean([)(x ,/.,) ~ E[J(x, t"OJ ~ I f(x,t,.()d~ V",([)( x , t,) ~ 1E[([ (x .t,.(j - meau([))2J ~ lrr ( [(x , 'f, () - mean(j))2d~, or more generally for the pilI moment: n<:[([(x , t, ,(j - mean([))"J ~ r (f( x , tf' ~) - m.ea.n(f))l'dll. lr (4. 1) (4.2) (4.3) meanU) and V(lTU) (and other high-order moments) arc themselves fields defined x due 16 depot shape variability and T is the domain over which the random variable varies. Notice also that mean(f) and var(f) are functions of tf. Based upon these quantities, we can make steps toward quantifying the impact of depot shape variability on drug transport in our test volume. We note that the process of moving from the uniform shape distribution to the distribution of f(x,tf,£) (the drug concentration distribution) is highly nonlinear; the resulting distribution is unlikely to be uniform. We present means and variances (i.e.. the first two moments) as the lowest order statistical information one would examine, but higher order moments may (and in general, will) be necessary to fully characterize the statistics of the resulting process. An example analysis showing higher order moments is presented in Section 4.4. The generalized polynomial chaos/stochastic collocation approach takes advantage of quadrature rules to integrate the stochastic process of interest over the stochastic domain, thus allowing the efficient computation of quantities such as means and variance. Under assumptions of smoothness, which in this case equate to the assumption that drug concentration varies smoothly with changes to our parameterized depot shape model, we can gain exponential convergence in the accuracy of our computed statistics as a function of the number of diffusion simulations wc are required to run. In the discussion of our experiment, we will refer to the simplified diagram in Figure 4.1. In the stochastic collocation approach, a collection of points at which the random field is to be sampled is specified and a set of corresponding weights are computed that account for the probability density function characteristics underlying the set from which the points (or outcomes) are drawn. In the context of our study, each collocation point £ represents a particular gel shape selected from our outcome set. Figure 4.1 (1) shows an example of a univariate smoothing parameterization sampled at three collocation points, and Figure 4.1 (2) shows the resulting depot shapes. Since the embedding process is deterministic given a depot shape, the entire volume (shown in Figure 4.1 (3)) is a consequence of the chosen collocation point. At each collocation point, wc can then compute the drug concentration f(x,tf,£) through traditional finite element diffusion simulation techniques, as shown in Figure 4.1 (4). Now, unlike traditional Monte Carlo in which large numbers of collocation points are used in an attempt to form accurate statistics, only a limited number of samples are used. The gPC-SC approach exploits smoothness assumptions on the random process to orchestrate the selection of points and weights such that accurate depot shape variability and r is the domain over which the random variable varies. Notice also that mean(j) and var(f) are functions of tf. Based upon thcse quantities, we can make steps toward quantifying the impact of depot shape variability on drug transport in our test volume. We note that the process of moving from the uniform shape distribution to the distribution of f(x, tf'~) (the drug concentration distribution) is highly nonlinear; the resulting distribution is unlikely to be uniform. We present means and variances (i. e., the first two moments) as the lowest order statistical information one would examine, but higher order moments may (and in general, will) be necessary to fully characterize the statistics of the resulting process. An example analysis showing higher order moments is presented in Section 4.4. The generalized polynomial chaos/stochastic collocation approach takes advantage of quadrature rules to integrate the stochastic process of interest over the stochastic domain, thus allowing the efficient computation of quantities such as means and variance. Under assumptions of smoothness, which in this case equate to the assumption that drug concentration varies smoothly with changes to our parameterized depot shape model, we can gain exponential convergence in the accuracy of our computed statistics as a function of the number of diffusion simulations we are required to run. In the discussion of our experiment, we will refer to the simplified diagram in Figure 4.1. In the stochastic collocation approach, a collection of points at which the random field is to be sampled is specified and a set of corresponding weights are computed that account for the probability density function characteristics underlying the set from which the points (or outcomes) are drawn. In the context of our study, each collocation point ~ represents a particular gel shape selected from our outcome set. Figure 4.1 (1) shows an example of a univariate smoothing parameterization sampled at three collocation points, and Figure 4.1 (2) shows the resulting depot shapes. Since the embedding process is deterministic given a depot shape, the entire volume (shown in Figure 4.1 (3)) is a consequence of the chosen collocation point. At each collocation point, we can then compute the drug concentration f (x, t f'~) through traditional finite element diffusion simulation techniques, as shown in Figure 4.1 (4). Now, unlike traditional Monte Carlo in which large numbers of collocation points arc used in an attempt to form accurate statistics, only a limited number of samples are used. The gPC-SC approach exploits smoothness assumptions on the random process to orchestrate the selection of points and weights such that accurate 17 results can be obtained. We use third-order Smolyak points that require (/ = 9 points for a single random dimension, as this level of integration is sufficiently accurate for the needs of this study (note that Figure 4.1 shows only the three first-order Smolyak sample points). In contrast, the number of realizations required for equivalent solution accuracy via the Monte Carlo method is significantly larger. Whereas in the one-dimensional case, the gPC-SC approach reduces to Clenshaw-Curtis quadrature [27], we will continue to formulate our method as stochastic collocation iir order to motivate the extension to multidimensional parameterizations. Once; diffusion computations have been done for each depot, shape dictated by the collocation sampling, the statistics of drug concentration as a consequence of depot shape is given by the following expressions: mean{f)(xAf) = E[f(x,tf,0] (4.4) 9 var(f)(x,tf) = E[(f(x,tf,0-mean(f))2} (4.5) 9 i=\ where q is the number of points used and Wj denotes the weights. As discussed further in Section 4.4.1, calculation of these statistics over the distinct diffusion solutions requires points between Figure 4.1 The 6))- Similarly, higher-order centralized moments can be approximated by: E[(f(x.tf.O-mean(f))»}* (4.6) 9 Y,wj(f(x,tf,Zj)-mean(f))'> j=\ Discussion of the formulation and computation of further statistics based upon the stochastic collocation method can be found in [32]. Of note is that increasing the order of the moments under examination may lead to an increase in the number of collocation points required. \Vc lise t hi rd-order q t his noLe order III t he Ilutnber re<luired solut ion lurger. 27J, ollr ill Ilmlt.idimcnsional pal'aInctcrizations. Ollce each depot shape dictated the collocation sampling, the statistics of drug concentration a<; t he meaaU)(x ,t,) ~ E[f(x ,t,,{)] , " L wjf(x , t" {j) j=i " ,,-(f)(x ,t,) ~ E[U(x ,t"O - meanU))'] q :::::: L 'Wj(f( x , tf : {j) - menn(f))2 j=i ( 4.5) t he a nd 'Wj 4.4.1 , calcu lat.ion stat.istics sampling the solutions at a set of point.s uniform het.w(.'Cn simulation solutions, as shown in Figurc 4. 1 (5). Thc means and variances can then be computed at these points (Figure 4.1 (0)). eNl E[U(x , t" {) - me.aU))"] " <, L Wj(f( x , t'f , ~j) - me(tn(f»" j=1 (40) furt her 32}. t hat uoder examinat ion increaBC 18 Shape Parameterization Shape Generation 7 t V (3) Shape Embedding / Mesh Generation l i t ^5 $ i Diffusion Simulation (4) (5) Compute Statistics t (6) Figure for shape parameterization Smolyak points collocation. (1) The univariate shape parameterization and three collocation points. (2) The three gel shapes corresponding to the collocation points. (3) Cutaway of tetrahedral mesh showing embedding of gel shapes. (4) Cutaway of tetrahedral mesh colored by concentration solution at day 14 of simulation. (5) Sampling of concentration values on a regular grid spacing over a region of interest. (6) Mean concentration at each sample point is calculated as a weighted siun of concentration values from (5). Para meterization smoother • (1 ) • • • .. .. .. (2) ; I .. .. .. Simulation .. Sampling .. .. Mean F ig ure 4.1. Shape statistics using stochastic collocation: calculation of the mean drug concentration fol' a sha pe paramcteri:-:ation using 1st order SmoJyak POillt8 for collocntion. collocat.ion Illcsh concent ration spacing: 11 ~'I ean cOllcellLration weightt.>ci sum concellLration (19 In order to run our diffusion simulation, we must first have a discretized representation of the tissue in which we are interested. For our statistical experiment, the model will be a combination of 'static" components not changed from simulation to simulation, and a model of the injected drug depot generated from our parameterization that changes based on the parameter value chosen. The depot shape is assumed to be unchanged within a single diffusion simulation. Although this ignores degradation of the gel over time, changes to the shape or diffusive properties of the depot would be difficult to experimentally determine and greatly increase the complexity of the simulation. Both the models of the static components and the depot shape parameterization rely on segmentations of MRI volumes. Once specific structures are segmented, geometric models of these objects are generated. This process was described in Chapter 2. These models are used to create a finite elements diffusion simulation for each parameter value (/.(..each gel shape) as described in Chapter 3. For the rapamycin study, we retain the solutions to the diffusion simulation at a set of timepoints T at which we will calculate our statistics. Our simulation is carried out over a period of 80 days with a timestep of 2 hours, which has shown sufficiently accurate results for our purposes. As mentioned previously, to apply the stochastic collocation method, we require a parameterization of depot shape - the variable component of the geometric model. To physically model the injection process would require a complex, high-dimensional parameterization based on properties of surgically damaged tissue, gel properties, and injection techniques, and would result in a more complex simulation than the diffusion problem we are attempting to solve. Instead, our parameterization is based on the results we have - namely segmentations of the injected depot from post injection MRI. In this way. we can develop a simple, low-dimensional parameterization that still respects the available data. We explore two pararneterizations, both based on the evolution of existing depot segmentations using level-set methods [21]. Since level-set methods operate directly on a volume of data, they are an obvious choice for working with images acquired from MRI, and have the advantages of easily allowing deformations based on differential properties of a surface and automatically handling changes in surface topology. For both 4.3 Characterization of Object Variability rUll discreti'l.ed representation combinatioll static' t his grently t.he parametcri>':ation Mill volumes. Oncc specific structures are segmented, geometric models of these objects are generated. This process was described in Chapter 2. These models are used to create a finite clements diffusion simulation for each parameter value (i. e.,each gel shape) as described in Chapter 3. For the rapamycin study, we retain the solutions to the diffusion simulation at a set of timepoints T at which we will calculate our statistics. Our simulat ion is carried out over a period of 80 days with a timestep of 2 hours, which has shown sufficiently accurate results for our purposes. parameteri>':ation 1.\ tt."Chniques, paramctcri>':ation - postilljection MRl t.his way, pl1rameteri>':l1tion 4.3.1 Level-set Shape Representation \Ve parameterizations, segmentfltions all voltune acquimd l\'lill, deforlIlations 011 differential ' . 20 parameterizations we require that the volume of the depot be constant for any parameter value in order to keep the amount of drug uniform for all simulations while still using the known value for drug concentration in the injected gel. Our parameterization will be formed by deforming an initial depot surface over time. Note that 'time' in this section does not refer to the diffusion simulation time, but an arbitrary unit describing the amount of evolution of our level set surface has undergone. In order to set up our level set methods, assume that we start with some initial surface <So at time t = 0. We will find a scalar function such that <So.a- is an isosurface of 0 with value A-; that is, So,* = {x\e{x,to) = k}. We can then modify this surface changing the function 0 due to properties of <S,. More generally, any point in a volume dataset at time can be thought of as a point on the isosurface St,o(x)- The function 0 can be modified at all points based on properties of the isosurface passing through that point (the local isosurface"). The most obvious choice for our embedding 0 is a distance transform to the isosurface St- In this case, the surface we are interested in is St.o, which we will now refer to as simply <Sj. For surface evolution, we move the isosurface St in the surface normal direction n: - IvSii <4-8> by modifying the values of 0 over a series of tiniesteps: St+i(x)=St(x) + V(x)-n (4.9) where x) is the 'speed function." a function of local and global properties of the local isosurface and of independent forcing functions. As described in Section 4.2. this parameterization be sampled at specified collocation points (' times' in our current formulation). In order to ensure stability, we use a variable timestep for our level set evolution based on the fastest moving point on our surface. We do not attempt to sample at the exact collocation point specified, but use a 'nearest neighbor' approach, ensuring that our surface is sampled within one timestep of the prescribed sample time. Since the purpose of the variable timestep is to limit the movement of any point on the surface! to one paramcLerization illitial SectiOIl evolutiou lllldergone. So O. 0 SO,A- () k; s", .. ~ (xlO(x , to) ~ k). (4.7) by fUllction () St. ~'ilore x tillie t isosurfacc S I,O( 3:)' flLnction () 'isosurface' ). () transforlll St. t his o , St. surfacc din:"ction (4.8) by tllodifying the valucs of () over a series of timesteps: SH'(X) ~ S, (x) + V(x ) · n (4.9) V(x ) function ,' 4.2, will times ' enslire \ve lise timcstep but, ensuring t he t he t imestep thc movemcnt surface 21 voxel length in a single timestep, the maximum error introduced by this approximation is one grid spacing, and for most of the surface, the actual error is much less. Since we are only interested in a single isosurface of 8, it is unnecessary to update the entire volume at each timestep. Only voxels near St at timestep t are needed to calculate St+i- We therefore use a Sparse-Field method [29] for calculating updates. This greatly reduces the computational burden while giving sub-voxel accuracy for surface position. 4.3.2 Parameterization The first parameterization takes an existing gel shape chosen as an example of a particularly irregular surface and smooths it using a volume-preserving mean curvature flow transformation. Mean curvature flow smooths a surface at each timestep by moving each point on the surface in the surface normal direction by an amount proportional to the mean curvature. In its original formulation (lacking volume preservation), mean curvature flow is the gradient descent for the surface area. The mean curvature H of the local isosurface is proportional to the divergence of the local isosurface normal: H = \(V-n) (4.10) This technique, however, reduces the interior volume of the surface, generally resulting in a singularities at 1^. In order to enforce volume preservation, the movement of the surface is adjusted by the average mean curvature so that the total offset of the surface-in the normal direction is zero: x) = x) H) where »J&JB* sdx The effect of volume preserving mean curvature flow is to create progressively smoother shapes that approach a sphere, as shown in Figure 4.2. For a more rigorous treatment, see [8]. Using this parameterization in conjunction with stochastic collocation, we can study the effect of depot shape smoothness on experimental outcomes of drug diffusion. Equivalently. for a fixed depot volume, this can be seen as studying the effect of surface area. ill maxillltUll Jess. arc s ingle 0, llllllccessary vaxels S t SH1- We therefore usc a Sparse-Field method [29] for calculating updates. This greatly wbile 4 .3 .2 P a rameterization 1: Mean Curvature Flow paramctcriJ';ution existillg part icularly aud SlIIooths timv poiut 011 t.he ill lUcan CUlvdture isosurfacc 1 ~ -('J. n ) 2 reduccs geuerally ill at. too. III prescrvat.ioll, t he t he t hat surface in lIormal V(x ) ~ (H(x ) - H ) H ~ j , H(x )dx fsdx (4.11) (4.12) T he ri b'OI"OUS sec conjuIlction t hc dcpot all ex perimental Figure 4.2. MCF gel morph: gel shape smoothing using volume preserving mean curvature flow Parameterization Metamorphosis The second parameterization performs a volume-preserving metamorphosis between two depot segmentations chosen to represent the extremes of surface regularity for the existing segmentations. The less-regular shape was taken as the 'source' surface and the more-regular shape as the ' target' surface for the morph. In order to allow volume preservation, rcsampled have same interior volume as the source surface. The morphing technique used is based on Brecn and Whitaker's work on level set shape metamorphosis [5]. where transitioning between the two surface models occurs via a level-set formulation with the speed function given by: V(x) y(x) where ^{x) is the distance transform of the target image. In this way. points on the source surface will move towards the target surface whether they are inside or outside the target surface. This shape metamorphosis, however, is not volume preserving. Although we know that the interior volumes at the beginning and end of the transformation are equal, the volume of intermediate1 surfaces is not constrained. As Brecn and Whitaker noted, the volume inside an intermediate surface $2/ can be divided into two regions - a region of voxels inside the target surface that will grow to match the target, and Q". a region of voxels outside the target surface that will shrink to zero volume. For any x € the two regions can be distinguished by the sign of ~/(x). In order to ensure constant volume throughout the transformation, we introduce two scaling factors, pi and p0, which scale the speed function of voxels in ilj and Qa, respectively. The values of the scaling factors rely on the total movement of the surface of the two regions, Sl t = {x e St\j(x) > 0} and S° = {x € 5(|7(x) 0}. Let V{ = j s , V(x)dx and v0 = / s „ V(x)dx. We then constrain the change in the two regions to be equal at each timestep by updating p\ and p0: 22 Fig u re gcl 4.3.3 P a rameterization 2: Me tamorphosis T he t.wo cbosen rcglllari~y cxistiug mor<.'-fcgular ta rget' morpho preservation. the target surface was scaled and rcsamplcd to hove the sallie iutcrior volume slII'fa.cc. mOl'phing ust."(i 011 Breen 5], betwccn forrnulatioll sp("ed fHflctioil Vex) = >(x ) (4.13) 'Y(x ) disl<.U!ce trausform way, ou Ulove This intermediate lIot Breen IIOt.ed, lin int.enn(..uillte stnface nl call - nL 11 illside I.hat nr, illly E nt, c;m 'Y(x ). Pi Po, fllll etion "oxels ni no, 011 St E Sd'Y (x ) > O} /llld Sf = {x E Sd'Y{x ) < OJ. Let Vi = Is' , V{x)dx and Vo = Is", V(x )dx . We then constrain hvo t imestcp Pi Po: Pi = { 1.0 otherwise % Xvi>v0 { 1.0 if Vi > v, otherwis'eo (4.14) where the level set formulation now becomes: S(+i(x) = Si(x) + ptV(x) • n , l e f lj St(x) + PoV{x) n. x G Q 0 Figure 4.3 shows the intermediate shapes obtained by the metamorphosis method between two segmented gel shapes. One of the gel shapes is less regular than the other, providing a second way of studying the effect of surface area for fixed volume depot shapes. The data presented in this section is computed from simulation results for input data representing the discretization of the simulation domain for depot shapes obtained by sampling our shape parameterizations at the nine Smolyak collocation points specified for third-order convergence. Note that although these results give accurate characterization of the variability based on our modeling and simulation assumptions, validation based on in vivo experimental studies has not yet been performed. Three methods will be described below for obtaining statistics based on these sim­ulation results. Section 4.4.1 deals with obtaining statistics over an area for which the location and number of nodes of the solution u(x) differ for different shape parameters. Section 4.4.2 shows the mean and standard deviation of total drug delivery to a region of interest over time, as well as an example analysis of higher order statistics on this data. Section 4.4.3 shows the simpler case of computing statistics over a static region of the discretized domain. Results Figure 4.3. Metamorphosis gel morph: volume preserving gel shape metamorphosis { '" = /0 if 'Ui > Vo t,he fa nnulation 5 ,+,( x) -_ { SSt,((Xx) ) ++ Ppo;V\f ((xx )) '· nn., Uj >vo otherwise 23 (4.15) hy rnetamorphosis Gue t he depot. 4.4 R esults da ta pn.> . scnkxi sectioll froUl tiirnulation reprcscnt ing discret ization simulat ion depot. obtaint>d Ollr lIille t hird-order convergence. t.ha t alt hough \~ariabiliLY obtnining sim­ula tion rcsults. dea ls wi th ovcr t he loca.t.ion 'x ) different. panuneters. t he me~U1 s tandard ~wal ys is s tatis tics t. his compu t.ing s tatis tics F igure i\'letamorphos is metamorphos is 24 AAA Domain Resampling Since a different finite clement mesh is generated for each depot shape (as described in Section 2.3). the solution vectors u from different discretizations do not represent weights for the same basis functions, or even the same number of basis functions. However, the solutions u(x) are defined everywhere over the same domain O. Samples can therefore be taken at a constant set of points in $2, and compute statistics at these points. Since the control volume is specifically chosen to be larger than the area where consequential levels of drug are expected to diffuse during our simulation period, there is no reason to sample the entire control volume. For simplicity, a minimum level of drug concentration in which we are interested is set, and the axis-aligned rectangular region encompassing all areas with at least this drug concentration in any of the simulation solutions is found. This region is then sampled on a regular grid, and the mean and standard deviation calculated at these points. This process is shown in Figure 4.1 (5) and (6). Figure 4.4 shows an 11mm x 11mm 2D region of the venous anastomosis for the solution at simulation day (tf) 14 sampled at a resolution of 512 x 512. The location of the slice is shown in (a), with the structures containing each sample point shown in (b). Note that the gel region changes for each choice of parameter £, and the specific depot region shown in yellow in (a) and (b) is for reference only. The mean and standard deviation at each sample point is shown in (c) and (d), respectively, for the shape parameterization based on mean curvature How (MCF). The mean and standard deviation are similarly shown in (e) and (f) for the metamorphosis-based shape parameterization. Since we are interested in the amount of drug delivered to the venous anastomosis over time, the mean and standard deviation of total drug in the vein wall is computed at each timepoint e T. Each element of our tetrahedralization s € T(it) represents one of the materials in our domain. In particular, we define £ v . the set of elements representing the tissue of the vein wall (See Figure 4.5). From this subset of elements, the total amount of drug in the vein wall Vj(t) = f£u(x)dx is calculated for each solution timepoint t e T in each simulation arising from a parameter value £j. The mean and standard deviation at each timepoint is then given by: 4.4.1 Domain R esampling Since a different finite element mesh is gCllcratl.'CI for each depot shape (as described in Section 2.3), the solution vectors u from different discretizat ions do not represent weights for t he same basis fu nctions, or even the samc munber of basis functions. However, the solutions u(x } are defined everywhere over the same domain n. Samples can therefore be taken at a constant set of points in n, and compute statistics at these points. Since the control volume i1) specifically chosen to be larger than the area where consequentiallevcls of drug are expected to diffuse during our simulation period, there is no reason to sample the entire cont rol volume. For simplicity, a minimum level of drug concentration in which we arc interested is set, and t he axis-aligned rectangular region encompassing a ll areas with at least t his drug concent ration in any of the simulation solutions is found. This region is then sampled 011 a regular grid, and the mean and standard deviation calculated at t hese points. This process is shown in Figure 4.1 (5) and (6). Figure 4.4 shows an llmm x llmm 20 region of the venous anastomosis for the solution at simulation day (tl) 14 sampled at a resolut ion of 512 x 512. The location of the slice is shown in (a), with the structures containing each sample point shown in (b). Note that the gel region changes for each choi<.:e of parameter ~, and the specific depot region shown in yellow in (a) and (b) is for reference only. The mean and standard deviat ion at each sample point is shown in (e) and (d), respectively, for the shape parameterizat ion based on mean curvature flow (MCF). The mean and standard deviation are similarly shown in (e) and (f) for the metamorphosis-based shape parameterization. 4.4.2 Total Drug Delivery t he t imepoint t E T . E E O) £'H clements t he t his \tj(t) = J£~ u(x )dx is calculated for each solut ion timepoint t E T ~j. t imepoillt (f) resampling plane at day 14. (a) Location of collocation slice; (b) Material locations in slice plane, colored to match the regions as shown in (a): (c) MCF diffusion mean values; (d) MCF diffusion std. dev. values; (e) Morph diffusion mean values; (f) Morph diffusion std. dev. values 9 mean(V(t)) = ]Tu/,V}(f) (4-16) var(V(t)) = f^wjiVjW-meaniVit)))2. (4.17) 25 o (e) (f) Figure 4.4. Domain I"csampling results: Slice in the Z-planc of sampled collocation data n) LocatiOIl ~Il ateriallocaiiolls a); rvlCF MeF s td . val ues: i'vlorph diH'usion , ntea.n(V(t» ~ I: Wj Vj(t) (4. 16) j=1 , VM(V(t)) ~ I: Wj (Vj(t) - "wan(V (t )))'. (4 .17) j=1 2(i Figure 4.5. Venous tissue drug delivery region. The opaque geometry shows the target are to Equation 4.17. 4.6 time simulations as well as the mean and standard deviation as calculated above for our two shape parameterizations. To look more specifically at the amount of drug reaching the target site, we define a spherical region of interest surrounding the venous anastomosis (Figure 4.8). We then define the set of elements £ s as those falling completely wit Inn I his sphere (Figure 4.8(b)). Of interest to us arc the elements c", By the same equations (( 1.16) and (4.17)) we can compute similar statistics on this reduced region. The results are shown in Figures and 4.10. As noted in Section 4.2. the resulting drug concentration distribution due to our shape parameterization is not completely specified by a mean and variance. Further statistics can be calculated, as shown in Figure 4.11. This shows the skewness and kurtosis of the distributions of total drug delivered to the venous tissue at our sampled timepoints under the MCF shape parameterization. The skewness and kurtosis are calculated from the normalized third and fourth moments of the distribution, as given in [24]. The PDF (whose domain is normalized to [0.1]) of the distribution at day 10 is also shown. 26 F igure 4 .5. T he region (venous tissue) over which total drug delivery is calculated. while as previously noted, higher order moments arc easily computed in a similar manner Figures 4.G and 4.7 show the amount of drug over tillle for individual shape parameter paramctcri:r.ations. spccifically sile, surroundi ng £.~ t hose within t his sphere (Figure 4.8( b) ). £81) = £!j n £1/. ~he « 4.16) 4.17» t.his rcsults showll 4.9 Clnd 4.2, t.he concent.ration d istriblLtion slwpe cOInplctely statistics t.he \·enous ti;;sue tlmepoints ~'ICF T he arc T he wbose 0, 1]) t he Total Drug in Vein, MCF Parameterization 30 40 50 time (days) Figure 4.6. Total drug delivery results, MCF parameterization. Results for the nine individual collocation points (depot shapes) are shown in blue to green, with the mean and variance shown in red. The means and variances for this simulation are shown in Figures 4.6 and 4.7. Com­paring these figures, we see that the skewness and kurtosis correspond to the increased 'lopsidedness' of the distributions from the mean, as we would expect. The PDF is also as we would expect from the given timepoint's distribution, which has several samples grouped at the lower end of the concentration range - note that the mean is near the lowest concentration at the timepoint. Whereas in general, the discretization Q) is different for different parameter values of f, using the discretization method described in Section 2.3, the elements comprising the static geometry are always the same. In particular, the elements Sv comprising the tissue we are interested in are the same for any value of £. Since the concentration field within a tetrahedral linear finite element is given by a simple linear weighting of the values at the element's nodes, the mean and standard deviation fields within £ v can be completely described by using computing the means and standard deviations at the x lO-e 1.6 1.. 1.2 \i S 2' 0.8 0 ~ 0.6 0.' 0.2 0 0 '" 20 '" lime' ("da ys) 50 collocation point 1 collocation point 9 --mean value 60 70 27 80 F igure J'vICF v/\fiances arc sec lopsidedness ' grouped at the lower end of the concentration range - note that the mean is near the lowest concentration at the t imcpoint. 4.4.3 Nodal Statistics Whercas gcneral , T(O) {, comprisi ng particular. t:1i arc saIne ally {. t:.., dcscrib(.x\ comput ing Figure 4.7. Total drug delivery results, morph parameterization. Results for the nine individual collocation [joints (depot shapes) are shown in blue to green, with the mean and variance shown in red. tissue (shown in red) contained within the spherical region shown around the anastomosis. x 10-' 1.6 1.' 1.2 'ii S ~ 0,8 ~ "li 0.6 OA 0.2 0 0 " Total Drug in Vein, Morph Parameterization time" (d ays) 50 collocation point 1 collocation point 9 --mean value 60 70 28 80 Illorph points ill (a) (b) Figure 4.8. Total drug delivery ROI. Total drug delivery is computed within the venous tissUe showlI x 10~9 T o , a l D r u 9 i n V e ' n N e a r Anastomosis, MCF Parameterization collocation point 1 time (days) Figure delivery. ROI MCF the collocation points (depot shapes) are shown in blue to green, with the mean and variance shown in red. points £ v and using these as the weights of the linear basis functions. Since sampling as described in Section 4.4.1 can be computationally expensive at high sampling resolutions or fine discretizations of the test volume due to the search for containing cell of each sample point, this method can improve performance when applicable. Figure 4.12 shows the mean and standard deviation on the surface of £ a v as computed by this method. 4.5 Discussion The results of Section 4.4 explore both the variability introduced by changes in depot shape and the effect of a specific parameterization of gel shape on this variability. From the results of Section 4.4.2 we see that the total amount of drug delivered to the venous tissue is quite similar for both shape parameterizations. We also see that the variance due to shape variation within each parameterization is fairly small, and that the amount of drug delivered decreases for 'smoother' depot shapes (See Figures 4.6 and 4.7). This agrees with our notion that increased smoothness corresponds to decreased surface area, which should lead to lower total diffusion rates. Figure 4.13 shows how depot surface 5 X 10-9 Total Drug in Vein Near Anastomosis, MCF Parameterization 30 40 50 1 collocation point 9 --mean value 60 70 29 80 F ig ure 4.9. Total drug delivery, ROT lv[CF results. Results for t he nine individual bllle ill Ev fUllctions. ill call computatiollully fille discretizatiollll poillt, t his performHnce all El1v t his t his S(."Ction S(.'C varia tion 30 time (days) Figure 4.10. Total drug delivery, ROI morph results. Results for the nine individual collocation points (depot shapes) are shown in blue to green, with the mean and variance shown in red. area changes over our two parameterizations. Whereas surface area decreases for both MCF see of the two parameterizations. This may be due to changes in the initial location and distribution of drug relative to the tissue in which we are interested. This effect is more pronounced if we restrict our interest to a more specific target site. Figures 4.9 and 4.10 show that if we look at only drug delivered to venous tissue near the anastomosis, the two parameterizations are not equivalent. The MCF parameterization now has higher peak delivery as well as higher variance compared to the morph parameterization. A side effect of the 'smoothing' caused by mean curvature flow is to reduce the convex hull of the depot, concentrating more drug near the center of mass. We believe that this explains the increased delivery seen here, as more drug is located near the target site. As can be seen in Figure 4.3, the depot retains a larger positional 'spread' during our metamorphosis parameterization, and therefore, this effect is not as apparent there. The x 10~9 T o , a l D r u 9 i n V e ' n N e a r Anastomosis. Morph Parameterization collocation point 1 collocation point 9 5 mean value 4 2 ~ 10-9 Total Drug in Vein Near Anastomosis, Morph Parameterization collocation point 1 colloca tion point 9 --°OU."-"'0C---~~--C~C-_-"~"---~~--O~O --C7~O--"80 aJ'C showli parameterizations, the j'vlCF parameterization has a smaller final surface area. From this, we would expect a lower diffusion rate; however, as we can sec in Figures 4.6 and 4.7, the curves for total delivery over time are nearly identical for the later collocation points duc locat ion T his parameteri;-;ations arc not. bas a<; depot , d rug \"'e ncar ta rget sccn sprcad' 31 Skewness and Kurtosis time (days) Figure 4.11. Higher-order moments: for sampled timepoints, skewness and kurtosis are given for concentration distributions of total drug in venous tissue due to shape variability. The PDF of the distribution at day 10 is shown inset. results from Section 4.4.1 show that the variance at any given point is closely related the local mean concentration. This also agrees with expectations about a smooth process such as diffusion varying by a smooth process such as our shape parameterization. These results are also quite helpful in determining how drug concentration varies across the target site and how the ' front' of diffusing drug moves through the test volume. Although the parameterization of properties such as shape can be cumbersome and still requires judgment in determining that the appropriate range of variation is expressed, we have demonstrated that through methodologies such as the stochastic collocation method the problem of studying the impact of shape variability on a particular forward simulation is tractable. The basic requirements are that: 1. The variation can be expressed by a parameterized model, and 2. The resulting stochastic field varies smoothly with changes to model parameters. Once the underlying stochastic process is characterized, quantification of the introduced timepoillts) a re resu lts variancc ;:It ally givcn va rying requircs rauge t he plll'tieular parameterizL-ci !itochastic proCCS!i 7.8e-9 5.9e-9 3.9e-9 2.0e-9 •o.o Figure 4.12. Nodal statistics results: Drug concentration mean and standard deviation at day 14 given by nodal statistics for (a) MCF parameterization and (b) morph parameterization Surface Area 4000 3 4 5 6 7 shape parameterization Figure 4.13. Surface area of shapes from shape parameterization 32 7.8e·g (a) (b) mcrul aL pal'amctcri,.ation "E 2500 g •~ 2000 i ii! 1500 ' 000 soo Surlace Area o L---7---7---~---7---7--~--~~~ 123456789 shape parameteri.talion sample point Surfacc parmnctcrizatioll 4.6e·g 33 variability is quite straightforward and provides an important step in the validation and verification process, without which modeling (geometric, parameter, and mathematical) and simulation error bounds on the final solution must be viewed with some skepticism. Although our univariate shape parameterization is simple, in many cases, such a simple model will suffice. As noted earlier, development of this model require judgment on the part of the modeler and affect the accuracy of the computed statistics. At the very least, however, a formal statement of confidence contingent on experimental realities respecting assumptions of the parameterization can be made. It is worth noting that we are not limited to a univariate shape parameterization. such a simple parameterization is not adequate, a multidimensional shape representation can be created. In this case, we will have a multidimensional random variable £ , and will use multidimensional sample points as described in [32]. The primary restriction on this method is that the entries of £ must be independent. If this cannot be shown directly, the covariance matrix of £ can be generated in a similar manner to the statistics previ­ously described, and independence generated through a Gramm-Schmidt-like procedure. Although requiring more sample points than the univariate case, this extension scales well to a limited number of dimensions and is quite straightforward as the simulation process is identical regardless of the dimensionality of the underlying parameterization. all vulidation a lld OUI' univariat.e asimple will uud affect. fo rmal If s lLch paramet.erization ~ llIult idimensional ~ usc ill 32). ~ ~ shaWl] ~ ~ to previ-ously t hrough H."<jlJiring marc poi nts t he ident.ical lUldcrlying SEGMENTATION As the overall goal ol' the drug delivery project is to prevent the buildup of hyperplasia and therefore the stenosis of the implanted graft, monitoring the growth of hyperplasia is a central concern. Automated and semiautomated segmentation methods not only reduce the time necessary for an expert to segment the structures of interest, but can also reduce artifacts of independently tracing the 2D contours of a 3D structure, and can reduce intra- and interoperator variability [28]. MR angiography (without a contrast-enhancing agent) is the imaging modality used for this study. Although computed tomography (CT) provides higher resolution images with fewer artifacts, patients are subjected to a nontrivial amount of ionizing radiation. It also requires the administration of a contrast agent, making it a somewhat invasive procedure;. Although the invasiveness is not an overriding concern for the purposes of this study, the exploration of noninvasive imaging techniques is an important secondary goal of this project. Furthermore, although CT angiography provides clear images of the lumen, soft tissues are not directly imaged, making it unsuitable for measuring hyperplasia volume. However, as noted in Section 2.2, the limited resolution, imaging artifacts, and nonuniform RF penetration of MRI can cause problems even in the simple task of segmenting the lumen. A 2D MR image of hyperplasia growth at the anastomosis, along with the location of the slice, is shown in Figure 5.1. Hyperplasia segmentation attempts to measure the volume of smooth muscle cells (SMCs) growing at the intima (inner wall) of the vessels. This growth occurs in response to the foreign object (the graft) implanted in the body. The majority of this growth occurs at the venous anastomosis, as shown in Figure 5.1(b). As a large portion of the vessel wall is made of SMCs. there is very little contrast between the hyperplasia and the vessel walls. Additionally, the SMCs have gyromagnetic properties similar to the surrounding muscular tissue, making this CHAPTER 5 HYPERPLASIA SEGMENTATION 5.1 Introduction A'S the overall goal of the drug delivery project is to prevent the buildup of hyperplasia and therefore the stenosis of the implanted graft, monitoring the growth of hyperplasia is a central concern. Automated and seOliautoma.ted segmentation methods not only reduce the time necessary for an expert to segment the structures of interest, but can also reduce artifacts of independently tracing the 2D contours of a 3D structure, and can reduce intra- and interoperator variability [28]. MR angiography (without a contrast-enhancing agent) is the imaging modality used for this study. Although computed tomography (CT) provides higher resolution images with fewer artifacts, patients are subjected to a nontrivial amount of ionizing radiation. It also requires the administration of a contrast agent, making it a somewhat invasive procedure. Although the invasiveness is not an overriding concern for the pmposcs of this study, the exploration of noninvasive imaging techniques is an important secondary goal of t his project. Fmthermore, although CT angiography provides clear images of the lumen, soft tissues are not directly imaged, making it unsuitable for measuring hyperplasia volume. However, as noted in Section 2.2, the limited resolution, imaging artifacts, and nonuniform HE penetration of MRI can cause problems even in the simple task of segmenting the lumen. A 20 MR image of hyperplasia growth at the anastomosis, along with the location of the slice, is shown in Figure 5.1. Hyperplasia segmentation attempts to measure the volume of smooth muscle cdls (SMCs) growing at the intima (inner wall) of the vessels. This growth occurs in response to the foreign object (the graft) implanted in the body. The majority of this growth occms at the venous anastomosis, as shown in Figure 5.1(b). As a large portion of the vessel wall is made of SMCs, there is very little contrast between the hyperplasia and the vessel walls. Additionally, the SMCs have gyrornagnetic properties similar to the surrounding muscular t issue, making this distinction problematic as well. Fortunately, the graft is made of MR-nonresponsive polytetrafluoroethylene. MR However, the limited resolution of MRI compared to the thickness of the graft means that the void is not always a distinct contour. The graft thickness (0.5mm) is roughly a standard in-plane voxel size (0.5mm2). In order to make the graft boundary more distinct, a pulse sequence with higher in-plane resolution (0.25nun2) was used at the expense of slice thickness (2mm). This leads to a clear signal in cases where the slice direction lies in (or nearly in) the tangent plane of the graft wall at a given point. Where the slice direction is closer to the normal of the tangent plane, the graft cuts through only a part of the voxel region and partial-voluming occurs, causing a slightly gray region instead of a strong signal void (see Figure 5.1) Although the lumen can be segmented in a fairly straightforward manner by relying on voxel intensity differences as described in Section 2.2. intensity information alone cannot bc used to segment hyperplasia. Although histology tells us that some hyperplasia growth occurs in the native vessels outside the graft, the imaging sequences used were tumble to provide any contrast between normal and hyperplastic tissue. For this reason, hyperplasia Figure 5.1. Hyperplasia growth and imaging. Figure (a) shows the MR image from the plane shown in (b). 35 t he ~"'I R-nonres ponsivc polytctrafluoroet uylcne. It therefore shows up as a void in the reconstructed MIl images. i\IRJ dist,inct contou!". O.Smm) vaxel si:-:c O .5mm2). d isti nct, O.2511un2 ) expensc 21ll111). ill t<Ulgcnt partifll-voluming causi ng of sec 5.1 ) cau Oil voxcl canllot be segnlent t he unable norn111.1 tiss\le. t his reason. Vein Graft / t I---\~;'-_--;;-- Hyperplasia \ Slice Plane (a) (b) I,he ~ m. quantification was limited to growth within the graft. Even in this limited problem, region growing algorithms such as level sets are ineffective due to 'bleeding out' in areas where there is little contrast at the graft border. In order to effectively segment the growth, more constraints must be added. Since the tissue of interest is by definition the tissue between the lumen and the graft wall, the segmentation technique used was to segment the graft wall, subtract the already-segmented lumen from this volume, and count the remaining volume as hyperplasia. Since several properties of the graft are known, such as the expected circumference and smoothly-varying surface, these can be used to further constrain the segmentation. The method used for segmenting the graft surface consists of manually or automati­cally finding the orientation of the graft lumen, taking a number of regularly-spaced 2D slices perpendicular to the vessel axis, and finding the graft contour in each 2D slice. The use of 2D slices along the structure of interest avoids the problematic and computationally expensive 3D surface fitting problem, while allowing smoothness constraints along the vessel by penalizing large change from one contour to the next. Although the cross section of an unimplanted graft would bc a near-perfect circle, the actual cross-section is changed by contact with the internal structures of the pig and the suturing of the graft to the vessel, as well as extraction of 2D planes that are not perfectly oriented with respect to the vessel. Of these, suturing causes the most deviation from a circular cross section. Since the hyperplasia growth is concentrated near the anastomosis, the warping of the graft is most severe exactly where accurate tracking is most needed. Attempts at automatic contour initialization based on the lumen centerline were attempted, but often failed near the anastomosis where both hyperplasia development and graft distortion are most severe. order to provide a reliable initialization, a semi-automated technique is used in which a user selects a slice plane approximately containing the graft axis (Figure 5.2(a)). The user then draws two lines showing the intersection of the graft with the slice plane in the region of hyperplasia growth. These lines are smoothed and moved towards a local area of minimum image intensity, using an open-ended version of the basic active contour model described below in Section 5.4 (Figure 5.2(b)). A set number of regularly-spaced points are generated along the length of each user-drawn line. For each pair of corresponding points on the two lines, a plane 36 already­segmented snrface, automati­cally 5.2 Slice Plane Selection and Contour Initialization be In is created passing through these points and perpendicular to the current slice plane. An initial circular contour is then generated on this plane, assuming the two initial points used to create the plane define the diameter of the circle. It is worth noting that the user may not have selected a slice plane exactly containing the axis of the graft. Assuming that the graft cross-section is circular, in general, the two contour points define a chord. The assumed case in which this chord is the diameter leads to the circle of minimum circumference, meaning that the circumference of the initial contour generated gives a rough lower bound for the length of the desired contour (See Figure 5.3). 5.3 Data Preprocessing Tracking of the graft wall is complicated by the fact that the boundary is a strong edge between tissue and lumen where there is no hyperplasia growth, but is a weak linear structure near a strong edge where hyperplasia occurs. Since the lumen boundary is also smooth from slice to slice, the contour-finding algorithm would generally prefer this stronger edge to the weaker structure showing the graft. The data is therefore prepro-cessed to enhance the boundaries of interest, as shown in Figure 5.4. The contour created by the graft material in the presence of hyperplasia creates a curvilinear structure. A common strategy for detecting elongated structures uses the shape information contained in the hessian matrix calculated at each pixel ( [14], [20], [9]). Our method uses the ratio of large eigenvalue over small eigenvalue of the hessian matrix to detect linear structures (Figure 5.4(3)). However this detector still produces a response to edges, and in the case of strong edges the response is greater than to the curvilinear structures of interest. The filter is therefore modified to reduce the response in areas of high gradient (Figure 5.4(4)). Finally, noting that the filter should be picking up structures showing a signal void (dark area) in an area of relatively high signal, the data is further processed to remove structures brighter than their surroundings (by checking the sign of the hessian's large eigenvalue when computing the linear structure response, Figure 5.4(5)), and by removing structures whose neighboring pixels show too low a signal to be hyperplasia or surrounding tissue (Figure 5.4(7)). Although this procedure adequately enhances the graft boundary in areas of hyperpla­sia growth, it removes the more obvious boundary between surrounding tissue and lumen. Since this is the graft contour in areas without hyperplasia growth, this boundary must be included in the final preprocessed data. A simple Canny edge filter creates a binary map 37 prepro­cessed t.he whose neighboring pixels show too low a signal to be hyperplasia or surrounding tissue (Figure 5.4(7)). Figure 5.3. Possible contours for the same initialization points of these strong edges quite nicely (Figure 5.4(2)). In order to combine these volumes in a tunable way. the curvilinear detector response is thresholded and inverted, creating uniform dark contours along the detected structures. The Canny response (already binary) is scaled down and inverted, creating less-dark contours. The two volumes are then combined by voxel-wise multiplication of values. In this way. extraneous features in the volume are removed, noise is eliminated, and both edge and line features within volume are turned into line features (Figure 5.4(8)). This will facilitate contour detection in the next sections, as only one type of feature must be detected. A standard active contour ('snakes') model 13 was used for finding the graft contour in each 2D slice. This method performs well at finding a strong but incomplete boundary. It also takes natural advantage of the smooth nature of the graft contour by penalizing curvature, and since the initial contour is guaranteed to be 'close to' the desired final contour, the problems of capture range and local minima is at least partially avoided. 38 (u) (b) Figure 5 .2. User contour init.ializat.ion: selection of graft axis slice plane (a) and lIser drawn initialization contour (b). ini t ializat.ion F igure 5.1l(2» . ill way, t.he responsc thrcsholdcd creating unifonn detccted I'CSpOtl!:iC arc voxeJ-th is way, features arc eliminated. fcaturf'J> arc 8)) . fac ilitate t.ype detect.ed. act ive snake:;') [13) list-xi T his rnethod but. nat.ural guarantccd contolll", 39 (6) (7) if detection. (3) Vesselness filter output (4) Truncation and scaling of the vesselness filter, including division by the gradient magnitude. (5) Removal of ' bright' feature lines. (6) Thresholding. (7) Removal of features detected in 'dark' areas. (8) Weighted combination of linear structure detector and edge detector. 5.4 Active Contours The standard active contour model as proposed by Kass et al. [13] represents the boundary of the object to be detected as a parameterized closed curve C(t), t G [0,1]. The contour is iteratively updated to find a configuration minimizing the internal and external energy of the contour. The internal energy is the weighed sum of the first and second derivatives, integrated over the length of the curve: /'a-C'(t)2 + f3-C"(t)2dt Jo (5.1) This attempts to create a smoothly parameterized contour with minimal curvature. In general, a and 0 could be parameterized as well, allowing corners or discontinuities to occur by setting their value to zero at specific points. Since such features are not expected in the contour representing the graft cross-section, constant values are used for these parameters. The external force represents the image informal ion. In the case of line detection (dark line on light background), the external force can simply be image value. The direction of energy minimization is then -V/: Figure 5.4. Data preprocessing pipeline. (1) The initial slice. (2) Result of Canny edge vessel ness filter. induding hy of linear structure detector and edge detector. contolll' paramet.erized t), E O, IJ. iuternai external euergy of Lhe contour. The internal energy is the weighed sum of the first and :;ccond derivatives. integrated over the length of the curve: 5 .1) curvature. 0 {J parametcrized comcrs Jlel'o featuf(.'S arc arc extel'llal fo rce information. background). - \1 /: Eexl = f I(C{t))dt 5.2) Jo VEexl = f - VI(C(i))dt (5.3) Jo The energy minimization should force the contour to be smooth (low Ei„t), and be located on lines in the image (low Eext). The total energy of the contour is therefore given by: E = Einl + Eexl (5.4) The contour is discretized in order to minimize this energy. It is represented as a finite number of ordered points. Although there are several methods proposed for minimizing this energy term [3]). this method employs the original method by Kass as it is simple, fast, and robust, and allows subpixel accuracy. With these energy terms alone, however, it was not possible to accurately locate the contour in each slice. Balancing the internal and external forces (through tuning the parameters a and 0) proved problematic, and it was common for contours to either expand beyond the region of the graft, attracted by distant features, or stay circular, not deforming enough to accurately represent the graft boundary. Further constraints representing known properties of the graft were therefore added. The first of constraints concerns the contour length. The graft is known to have a constant circumference, and therefore, the contours discovered should have a specific, known length L. A soft length constraint is therefore imposed on the contour as a force: Flenit) = 11. -(L-L)- n{t) (5.5) where 7/, is a tunable weight, L is the current contour length, and n{t) is the normal of the contour at *. This gives a force that can be calculated at each point on the contour, and moves the contour in the normal direction at each point (a "balloon force' [7]) in an attempt to either increase or decrease the boundary length. Although the value of L is theoretically known and fixed, in practice we see that at the anastomosis, the graft is not cut perpendicular to its axis, and therefore, after being sutured to the vessel, the graft cross-sections can have an apparent nonuniform circumference. The length of the initial contour still defines a rough lower bound for the actual contour length, however, so a value slightly larger than this initial radius is E", ~ 10' I (C(t))dt 'VE", ~ 10' 'VI(C(t))dt 40 (5 .2) E;nt), liucs ill {Ee:l:d . (5.4) Ilumber arc ([30], Kas..., et al. sim ple, fast , robust. 'vVith it 0' (3) reprcscnt represclltiug c;onstraiuts kuown L. lengt.h (5.5) "fl. Cllnent 1"i(at. t. 'attempt. Ot· L fixed , t he however. vallie 41 clioosen as the desired length of each contour, e.. L, = 1.2 • D,. where D, is the initial diameter of the contour. The second constraint comes from the fact that two points on the boundary are already known, namely the points used to determine the initial contour diameter. Since the used on the contour to this fixed point: where P is a fixed point and C* is the attached point on the contour. The point on the initial contour closest to P is selected as the corresponding C,. This method of constraining a single point on the curve to stay "near' P is not ideal, as it assumes that the initial diameter is very close to a bisector of the final contour. Trials show that this assumption is not severely violated, but a more general constraint that 'some' point on the enhancement. The final constraint enforces smoothness from contour to contour along the graft. the vessel center and whose x and y axes correspond to the initial contour's diameter vector and the vector perpendicular to the graft axis plane on which the user drew the initialization boundary, respectively. Were the graft cross-section perfectly circular, and the center of the vessel accurately specified, we would expect all contours to be coincident in the 2D plane. Even though these assumptions do not hold in practice, the contours should vary gradually over the course of the contours. This concept of 'smoothness' is enforced by penalizing the distance between a point on a contour and the corresponding point on neighboring contours. A spring energy is used as in the initialization point constraints above: contour point, excluding" the first and last contours, which only have one neighboring contour. These forces are therefore added to the external energy gradient, giving: Ept = lPt • | |P - Q choosen ill) desi red i.e., L; · Di , D; i,ll, t he t,he arc known. poilltS llSed t he contour should pass through these points, a spring energy is tilled to attach a point (5.6) Ci attachcd t he 011 Ci. T his 11 ' near ' idea l, illisumes init.ial t he curve should be near the initialization point would be preferable and a possible future The slice extraction along t he ves..,el axis creates a 2D image whose origin is at the vessel lhe vect.or t he initialization ilt contou rs contOHrs. T his betwccn a.'S constraints E ~ ",oot" = I~ ",ooth . IICi' - CiH I li' (5.7) where ct is the 'i tl ! point Oll the J.P' contour. T his leads to two spring forces on each excludi ng neighboriug coutour. a re VEcxt(t) = - V / ( C ( t ) ) + Flen(t) + VEpl{t) + VEsmoolh(t)dt (5.8) Figure 5.5 shows these forces after 100 iterations on a chosen active contour. It is worth noting that the total energy is minimized when the sum of the internal and external forces is zero as the contour reaches equilibrium, then, the internal and external force should be approximately equal in magnitude and opposite in orientation. 5.5 Volume Calculation Figure 5.6 shows the steps in the segmentation process. Once the contours have reached a stable state (in practice the contours iterate for a fixed number of steps empirically determined to reach a stable state), they are projected hack into the 3D space of the imaged volume, shown in Figure 5.6(5). Since each of the contours have an equal number of points, and point corresponences between contours are known, it is simple to build a closed polygonal model by building a triangle strip connecting each pair of neighboring contours and capping the ends somehow this algorithm inserts a new vertex at the center of the first and last contours, and builds a triangle fan connecting this central point to each pair of adjacent points on the corresponding contour (Figure 5.6(6)). As long as care is taken in orienting the triangles correctly, a simple inside/ outside test with respect to this mesh can be run on the centroid of each voxel contained in the bounding box of the mesh. This gives a simple binary segmentation of all voxels contained in the graft structure within the region of interest (Figure 5.6(7)). This volume contains the lumen as well, but over the small region in question, the lumen can be segmented by a simple thresholding operation (Figure 5.6(8)). After finding the intersection of the lumen segmentation and the segmentation of the interior of the graft, these voxels can be removed from our graft segmentation yielding our hyperplasia segmentation (Figure 5.6(9) and ( 10)). The volume of this segmentation can then be trivially calculated. 5.6 Results showing hyperplasia growth, and the lack of expert segmentations of hyperplasia with which to compare our results. Figure 5.7 shows a comparison of an automatic segmen­tation (yellow) and a hand-segmentation performed on the original 0.25x0.25x2.0mm voxels (orange). We can see that the segmentations do not line up exactly, as the limited 42 \7 E",(t) ~ \7/(C(t)) + F',n(t) + \7 EpI(t) + \7E",.oo/h(t)dt (- iu practi(.'C back corrcsponcnces betwl:cn - t.he segmcntation 8) ). yieldillg OUl' 10)) . Quantitative analysis of this method has been limited by the availability of datasets Figme II O.25xO.25x2.0mm voxeis \\le sec t hat Figure 5.5. Contour forces: Selected active contour after 100 iterations. Current forces acting on the contour are shown, as well as closeups of two regions. The 'External Force' is the sum of the gradient, spring, and length forces. The spring force add smoothness constraints and binds the contour to the initial user-defined points ('Constraint Points'). The large spring force shown in closeup (a) is a result of this binding. - I Internal - Forces External Gradient - Length - x Constraint Point Spring - 43 Seicctl.,{\ nre hvo stun l.'Ontour uscr-(,Points'). fOfce t his Figure 5.6. The segmentation pipeline. ( 1) The initial slice plane showing the region of hyperplasia. (2) Initial points along user-defined graft boundary. (3) Initial circular contours generated from points. (4) refinement using active contours. (5) Final contours after refinement. (6) Closed mesh generated from contours. (7) Segmentation of voxels inside mesh geometry. (8) Graft segmentation shown with lumen obtained by simple thresholding. (9) Hyperplasia segmentation generated as the difference of the graft and lumen segmentations. (10) Closeup view of the resulting segmented hyperplasia. 44 '" (OJ (0) (10) T he ci rcular gCllcmtcd Contour lIsing flfter Illesh silllple threshold ing. STart Closcup 45 Figure 5.7. Qualitative comparison of segmentations, automatic vs. hand segmen­tation of anisotropic data using Seg3D [2]. The hand segmentation on the original 0.25x0.25x2.Omm voxels is shown in orange, and the automatic segmentation on the isotropically resampled data is shown in yellow. Volume of the hand segmentation is 120m7?i2, volume of the automatic segmentation is 132m?n2 resolution of the anisotropic dataset tends to over-segment the hyperplasia by extending the segmentation in the slice plane direction. However, we do see some over-segmentation using the automated method in the axial view, where the lower left portion of the contour was attracted by i he nearby dark linear structure. Figure 5.8 shows a comparison of t he automatic segmentation (red) and a hand segmentation of the isotropically resampled data (green). In this case, the hand segmentation seems to be more conservative in labeling voxels as hyperplasia as opposed to lumen, as the coronal and axial view shows. Although adjusting the thresholding level for designating a voxel as lumen may bring our method closer to the hand segmentation, it is a subjective decision without ground truth about which voxels actually represent hyperplasia. With our automated method and guidelines for thresholding, we would at least expect a more consistent segmentation. SegJD T he humi t he O. 25xO.25x2.0mm voxc!s resam p[ed yeHow. Vol\lme l20mm2, 132n17u2 I'esolution hypel'plasia sHce sec the nearby dark linear structure. Figure 5.8 shows a comparison of the alltomatic sq,rrflelltation isotropicaUy r€Sampled t hresholding OUI' doser grouud t.ruth whicb and b'1.tideHnes for thresholding, we would at lea:~t ex pect a more consistent segmentation. Figure 5.8. Qualitative comparison of segmentations, automatic vs. hand segmentation of isotropic data using Seg3D [2]. The hand segmentation on the isotropically resampled voxels is shown in green, and the automatic segmentation on the resampled data is shown in red. Volume of the hand segmentation is 96mm2 , volume of the automatic segmentation is 46 Qllalilative Scg30 rcsarnpled t,hc 96m'm2. vohunc 3ulomatic 132mm2 The work in this thesis features elements from various segments of the drug delivery experiment; work for both simulation and experimental workflows, and work concerning both theoretical questions and practical engineering solutions. The finite element model generation and simulation pipeline from Chapters 2 and 3 outline a practical solution used to quickly generate finite element models from MRI data. Chapter 4 proposes a method it is shown that accurate results can bc calculated given an appropriate simulation method and shape parameterization, the examples presented are intended only as a framework and stop short of claiming specific results for this experiment due to the unknown accuracy of the current shape parameterizations and simulation methods. Chapter 5 introduces a semi-automated method for segmenting hyperplasia from MRI data, allowing the monitoring of hyperplasia volume over time. Preliminary results are presented, and seem promising, but a true 'ground truth' segmentation is needed. Although this work solved many of the problems initially identified, there are many open questions left as future work. Although the segmentation and simulation pipeline generates realistic finite element meshes, it has not been determined whether this added fidelity significantly impacts the final drug delivery solutions. Also, although results are given for the variability of drug delivery due to depot shape, we know that the proposed shape parameterizations are inadequate, and a multidimensional shape parameterization that better captures the underlying population of shapes is needed. Finally, a better ground-truth is needed for quantitative analysis of the hyperplasia segmentation results, and refinement of the method is necessary to better deal with the uncertainty due to the anisotropic nature of the data. Although many of the techniques developed in this work are specific to the experiment being discussed, the problems in formulating accurate simulations and using computa­tional techniques to expedite human-performed tasks are similar to those encountered in CHAPTER 6 CONCLUSION t hesis clements experi tllcnt, experiment.al workflow!;, solutiolls. T he aud nlld clement i\'IRJ for quantifying error in simulation results due to drug depot shape differences. Although be siulUlation sha pe experimellt current. parameterizatiom; si lllulat.ioll ::vlRI thc truth ' solvt.>d Blany rcalistic clement thc Also. a lthough arc proposed shapc arc inadequate. parameterization ca pl.tII'cs Heeded. l.l 1lCl.·ded t he hypcrpla.sia scgmclltation anbotropic t.he fonnull.lting simulat ions perfonned t.hose enco\lntered 48 integrating computer-assisted simulation and monitoring in a wide range of research. This thesis provides a case study of problems and solutions encountered in one such experiment. Further, the error quantification outlined in Chapter 4 can easily be modified to provide systematic sensitivity analysis for a wide range of problems. I hope that this and other work in this area demonstrates the ease of computing rigorous error bounds, even with a deterministic solver and limited computation time. "8 iHtegratillg assisted ami a nd solut.ions olle slich Further. t.he quant.ification ana lysis fol' dcmons t,rates ti l'l1e. [1] SCIRun: A Scientific Computing Problem Solving Environment, Scientific Comput­ing and Imaging Institute (SCI). [2] Seg3D: Volumetric Image Segmentation and Visualization. Scientific Computing and Imaging Insi itutc ( SCI). [3] A M I N I , A . , TBHRANi, S., A N D W E Y M O U T H . T. Using dynamic programming for minimizing the energy of active contours in the presence of hard constraints, pp. 95 99. [4] B O I S S O N N A T . J.-D.. A N D O U D O T , S. Provably good sampling and meshing of surfaces. Graphical Models 67, 5 (September 2005), 405 451. [5] B R E E N , D. E., A N D WmiAKER, A level-set approach for the metamorphosis of solid models. IEEE Transactions on Visualization and Computer Graphics 07, 2 (2001). 173-192. [6] C H R I S T O P H E R S O N , R. .J., T E R R Y . C. M., K I R B Y . R. M.. C H E U N G . A. K.. A N D SlHU', Y.-T. Finite element modeling of the continuum pharmacokinetics of perivascular drug delivery at the anastomoses of an arterio-venous graft. In preparation. 7] COHEN, L. D. On active contour models and balloons. CVGIP: Image Underst. 53, 2 ( 1991). 211 218. 8] E S C H E R , J., A N D S I M O N E T T , The volume preserving mean curvature flow near spheres. Proceedings of the American Mathematical Society 126. 9 ( 1998), 2789- 2796. [9] F R A N G I . A. F.. N I E S S E N . W. J., V I N C K E N , K. L., A N D V I E R G E V E R , M. A. Muliscale vessel enhancement filtering. In MICCAI '98: Proceedings of the First International Conference on Medical Image Computing and Computer-Assisted In- h ret iihtiii (London. UK. 1998), Springer-V'erlag, pp. 130 137. [10] G A N A P A T 1 1 Y s u I i R A M A N l A N . B . . A N D Z A B A R A S . N. Modeling diffusion in random heterogeneous media: Data-driven models, stochastic collocation and the variational multiscale method. Journal of Computational Physics 226. 1 (2007). 326 353. G E N E S E R . S. E.. K I R B Y . M., A N D M A C L E O D . S. Application of stochastic finite element methods to study the sensitivity of ecg forward modeling to organ conductivity. IEEE Transactions on Biomedical Engineering 55, 1 (January 2008). 31-40. REFERENCES 11 SCIRuu: SciClltific hnaging Visuali:r.ation, Institute SCI) . AMI :-II , A. , TEIIRAN I , 5. , AND WEY:'IOUTI I, T . t he constraints. BOISSON"''''!'. J .-0 .. AND Q UDOT, P rovably 1llt.'Shiug Gmphical 405- 5J BREEN , O. E .. AND WIHTAKER, R. T. A metnmorphosis Tmnsactions OI~ anti 07. 2001 ).173- 192. [6J C IiHISTOPII EIlSON , R . J. , T ERRY, C . lVI. , Kml3Y, R.. lVI. , C HEUNG , A . K., AND SIIIU, Y. -T. Finitc elclllcnt continumn pharmacokinetics artcrio-preparfl tioll. [7J CO Ill~ :\" Imaye Undo 'st. 1991 ). 211- 21 8. [8J ESCII EH, J ., AND SI:vIONETT, G. votlUne l1ear ATltc1-ican 126, 199S), [91 FRANC I, A. F. , N I ESS I~ N , W . J. , VJ :\'CKE1'. 1<' L. , .. \ ND VIERCEVEH, M. A. i'.lutiscale Fi1'st lrltemational ConfeH!1!ce Compu.ter·Assi,stcd In­tC11Jcntion (London, UK, 1998), Springer-Verlag, pp. 130- 137. [10] G ANAPATHYSUI3HAMAN IA i\ , 13 . , AND Z ABAHAS, i'.'lodeli ng hetcl'Ogellcotis driven modeb, lllut t isca te Comp'ut(ltion(ll 326- [11] G I~ N ES ER , E. , KIRBY. R. rd. , AND MACLEOD, R. S. clemcnt (.'Cg o rgan cOllcluctivity. Tmnsactions Enyinee1i.ny 2008), [12] G E N E S E R , S. E., S T I N S T R A , J. G.. K I R B Y , R. M.. A N D M A C L E O D . R. S. Cardiac position sensitivity study in the electrocardiographic forward problem using stochastic collocation and hem. IEEE Transactions on Biomedical Engineering. Submitted. [13] K A S S , M., W I T K I N , A., A N D T E R Z O P O U L O S , D . Snakes: Active contour models. International Journal of Computer Vision VI, (January 1988). [14] K O L L E R , T., G E R I G . G . . S Z E K E L Y , G.. A N D D E T T W I L E R , D . Multiscale detection of curvilinear structures in 2-d and 3-d image data. pp. 864- 869. [15] K U J I , T., M A S A K I . T.. G O T E T I , K.. LI. L.. Z I I U P L A T O V , S., T E R R Y , C. M., Z H U , W . , L E Y P O L D T . J. K.. R A T I I I . R.. B L U M E N T H A L , D. K.. K E R N . S. E., A N D C H E U N G , A. K. Efficacy of local dipyridamole therapy in a porcine model of arteriovenous graft stenosis. Kidney International 69. 12 (2006), 2179 2185. 16] M A R X , ( ) . . A N D M A R K S . A. Bench to Bedside: The Development of Rapanrycin and Its Application to Stent Restenosis. Circulation 104- 8 ( 2001). [17] N A R A Y A N A N , V. A. B . , A N D Z A B A R A S , N. Stochastic inverse heat conduction using a spectral approach. International Journal for Numerical Methods in Engineering 60, 9 (2004). 1569 1593. 18] P O R U M B E S C U , D., B U D G E , B., F E N G , L., A N D J O Y , K. Shell maps. SIGGRAPH '05: ACM SIGGRAPH 2005 Papers (New York, NY. USA. 2005). ACM, pp. 626-633. 1!) R O B E R T . C. P.. A N D C A S K I . I . A . G. Monti Carlo statistical methods. Springer. 2004. [20] S A T O . Y., N A K A J I M A , S., S H I R A G A , N., A T S U M I , H., Y O S H I D A , S., K O L L E R, T., G E R I G , G., A N D K I I < I N I S . R. Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images. Medical Image Analysis 2. 2 ( 1998), 143 - 168. [21] S E T H I A N , J. A. Level Set Methods and Fast Marching Methods. Level Set Methods and Fast Marching Methods, by A. Sethian. pp. ISBN Cam­bridge, UK: Cambridge University Press. June 1999.. June 1999. [22] Si, H., A N D G A R T N E R , K. Meshing piecewise linear complexes by constrained dclaunay tetrahedralizations. In Proceedings. 14th International Meshing Roundtable (2005). Springer-Verlag, pp. 147- 164. S M O L Y A K . S. A. Quadrature and interpolation formulas for tensor products of certain classes of functions. Akad. Nauk SSSR. (1963), S N E D E C O R . G. W . . A N D C O C H R A N , W . G. Statistical Methods. Blackwell Publishing, 1989. T E R R Y . C. M.. LI. H.. K H O L M O V S K I , E., K I M , S.. A N D Z H U P L A T O V , I. Serial magnetic resonance imaging (MRI) to monitor the precise location of a biodegradable thermal gelling copolymer for drug delivery after ultrasound-guided injection in a porcine model of arteriovenous (AV) hemodialysis graft stenosis, in preparation. 50 [12J GENESEH , S. E. , STINSTH!\ , J. G. , KlHl3Y, R. ivl. , AND MACLEOD, R. S. sensit ivity forWilrd lIsing Biomc(lical Engineering. 13J i< ASS, 1',,1. , 'WITKIN, A. , AND T EHZOPOULOS, D, coutour lntcnwtional COfllJHLler 4 1988), 321- 331. KOLLER, T .. GEHIG , G. , SZE I'~~LY , G .. AND DETTW1LEIt, D. l'vlultiscale curvil inear [15J K UJI , T ., i'viASA 1<I, T .. GOTETl , K., LI, L., ZII UPLATOV, S. , TERRY, C. NI. , ZHU, \V. , L EYPOLOT .• J. 1(. , RAHII , R. , BLU!I'IENTHAL, O. K. , K ERN , S. E. , :\1"0 CIIEUNG , }(idney 2179- [16J ,MARX, S. 0. , A:-ID i'vIAHKS. R. Rapalllycin Restcnosis. 104, 2001 ), 852- 855. NARAYANAN , B., AND ZABARAS, cOlHluction Jounutl Engineering 60,9 2004 ), [I 8] PORUl\-IBESC U, S. D. , BUDGE, B., F ENG , L. , Ai'."!) JOY , 1(. I. In SIGGRA PfI Pape1"S York. USA , AC:'vl, 626 [19] ROBEIlT, P. , M'm eASEl-LA , G. Monte [20] SATO. Y. , NAKAJl~IA , 5. , SIIIRAGA , 1'\. , ATSUl\II , H., YOSlIlDA, S. , KOLLER, T. , GERIG , G. , AND « IKINIS. mlllti~scale segmcntaLion alld visllulizatioll Imllye 2, 1998). SETIIlAN , MlU"Chirtg MethOlls. )' h"ll"ching J\'lethods. J. Sethian, 400. 0521645573. Cum~ bridge. Press, Junc 1999 .. J une 22J SI, H., AND C"'lTNER, 1<. ).'Ieshing piecew ise construillcd delaunay tctrahcdralizations. Intemational Roun{itable 2005), Springe r~Vcrlag, 164. [23] SMOLYAK. of SSSR, 4 240 243. [24] SNEDECOH, G. W .. A:'m COCHlb\N , W. G. [25] TERHY, t..,!. , LI. H., KllOLl\JOVSKI , E. , KIM, S. , ANI) ZHUPLI\TOV, Mill) ultl"asollnd~guided stenosis. preparation. 51 26] T O R D O I K . .1. H., V A N D K R SANDS, . M., A N D D E H A A N , M. W . Current topics on vascular access for hemodialysis. Minerva Urol Nefrol 56 (Sep 2 0 0 4 ) . 223 235. [27] T R E F E T H E N . L. N . Spectral SIAM. 2000. [28] WARFIELD, S-. Z O U , K . , A N D WELLS, W. Simultaneous truth and performance level estimation (staple): an algorithm for the validation of image segmentation. Medical Imaging. IEEE Transactions on 23. (July 2004). 9 0 3 - 9 2 1. [29] WlHTAKER, R. T . A level-set approach to 3 D reconstruction from range data. Int. Comput. Vision 29, 3 ( 1 9 9 8 ) , 2 0 3 - 2 3 1. WILLIAMS, D . J., A N D S H A H , M. A fast algorithm for active contours and curvature estimation. CVGIP: Image Underst. 55. 1 ( 1 9 9 2 ) . 14-26. [31] W I N K . ( ) . , F R A N G I , A. F.. V E R D O N C K . B.. V I E R G E V E R , M. A., A N D N I E S S E N , W. J. 3 D MR A coronary axis determination using a minimum cost path approach. Magnetic Resonance in Medicine 47. 6 2 0 0 2 ) . 1169-1175. [32] Xiu, D . Efficient collocational approach for parametric uncertainty analysis. Comm. Comput. Phys. ( 2 0 0 7 ) , 293-309. [33] XlU, D . . A N D HESTHAVEN. S. High-order collocation methods for differential equations with random inputs. SIAM .1. Comput. 27, 3 ( 2 0 0 5 ) . 1118 1139. [34] Xiu, D.. A N D K A R N I A D A K I S , G . E. Modeling tmcertainty in steady state diffusion problems via generalized polynomial chaos. Computer Methods in Applied Mechanics and Engineering 191, 4 3 2 0 0 2 ) , 4 9 2 7 - 4 9 4 8 . [35] XlU, D . . A N D K A R N I A D A K I S , G. E. Modeling uncertainty in flow simulations via generalized polynomial chaos. Journal of Computational Physics 187. 1 ( 2 0 0 3 ) . 137- [36] Y U S H K E V I C H , P. A.. P I V E N , J.. C O D Y H A Z L E T T , H., G I M P E L S M I T H , R., HO, S.. G E E . C , A N D G E R I G , G. User-guided 3 D active contour segmentation of anatomical structures: Significantly improved efficiency and reliability. Neuroimage 81, 3 ( 2 0 0 6 ) . 1116 1128. [26J TOlwom, ,/. 1-1 ., VAN DEH SMWE, F . \1.. AND DE I-I AI\:\', ;\1. W. 011 56' SCI' 2004), 223 235, [271 TREFETHEN. N. Spec:tHll Methods in MATLAB. SIA:'>.1. 2000. \OVAHFIELD, S .. Zou, 1< •• AND WELLS. Silllu!tnllcous pCrfOrlllllIlCe t.hc vlJlidatioll segmentation. Medical Imaging. IEEE Tmnsactions on 7 2004 ), 903- 921. WIHTA I{~~ It, T. 3D ['(!construction runge I nt. J. COHlput. 29. 3 (1998). 203- 231. [30] \VI LLIA;"'S, D. J .. ,\NO 511AII. \ 1. acti\'c CVG'IP: Underllt. 1 (1992), 14 26. [:31] WINK , 0., FRA:>iGI. A. P. , V ~: RI)O:-.lCK, B.: VI~;HGT.;vEH, M. A. , ,\1\1) :\IESSEX, \V .. 1. 3D l\IRA coronflry fl.:ds determinatioll 1.1 lle$o71unce -in 47, (j ( 2002). tW9 1175. X It' . D. ullcerlfdnty Gomut. Gomp"Ul. Pllys. 2. 2 (2007). 293- 309. XIU. .. ANI) H ESTIIAVE1\ .. J. colloc<ltioll differential J. Sci. C01nTJlLt. 27.3 (2005), IllS 1139. [34] XlV, D., A1\1) I<AHN[AI)AK[S, G. }. Iodding uncertainty Computc1" Mechan.ics 191,43 ( 2002). 4927 ·19-1S. [35] Xu..:. D .. ,\NI) J(A II.NIAD,\ KIS. C. l\lodeling ullcenainty ill Joumal I (2003) . 137- 167. [36] YUSHKEVICII , P. A .. PIVE>:. J .. COI)Y H AZ L I~TT. rI .: GJ:\IPEL SI"o.[[T[I, R ., H o. 5 .. GEE. J. C., ANI) G~;IUC. G. 3D actiyc Neumimuye 31. J (2006). Ill6 1128.