<%BANNER%>

A Computational Study of Flow and Mass Transport in Intrathecal Drug Delivery

Permanent Link: http://ufdc.ufl.edu/UFE0042230/00001

Material Information

Title: A Computational Study of Flow and Mass Transport in Intrathecal Drug Delivery
Physical Description: 1 online resource (54 p.)
Language: english
Creator: Fletcher, Adam
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2010

Subjects

Subjects / Keywords: cerebrospinal, computational, concentration, delivery, drug, dynamics, element, finite, fluid, human, intrathecal, mass, morphine, nerve, project, roots, sas, space, subarachnoid, transfer, transport, visible
Biomedical Engineering -- Dissertations, Academic -- UF
Genre: Biomedical Engineering thesis, M.S.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation

Notes

Abstract: A COMPUTATIONAL STUDY OF FLOW AND MASS TRANSPORT IN INTRATHECAL DRUG DELIVERY The detailed geometry of the spinal cord and its surrounding tissues make flow of the cerebrospinal fluid highly variable. Intrathecal drug delivery requires that drugs be injected into the subarachnoid space, which surrounds the spinal cord. In the present work we simulated the flow of cerebrospinal fluid and the mass transfer of morphine in a 29 mm section of the spine, by employing computational fluid dynamics. The Visible Human Project dataset is used to generate a mesh of the subarachnoid space and MRI data is used to apply a physiological pulsatile velocity condition at the inlet of the model. Morphine injection is modeled as a saturated point source with a concentration of 1.0. Results from the study show velocities peaking at 20 cm/s with phase shifts up to 5 degrees. The average Womersley number was calculated to be 3.5 with a deviation of about 1. Spatial velocity profiles show varying flow development in different regions, with a tendency to be blunted during systole. Reynolds numbers were calculated between 200 and 500 during peak systole, and agree well with the literature. After 1.65 seconds, the average relative morphine concentration was 4.4E-4. Dispersion was highly complex, involving circulation about the spinal cord, redirection around nerve roots, and a compression of the subarachnoid space. The Peclet and Schmidt numbers reveal that diffusion was dominated by inertial and viscous effects.
General Note: In the series University of Florida Digital Collections.
General Note: Includes vita.
Bibliography: Includes bibliographical references.
Source of Description: Description based on online resource; title from PDF title page.
Source of Description: This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility: by Adam Fletcher.
Thesis: Thesis (M.S.)--University of Florida, 2010.
Local: Adviser: Tran-Son-Tay, Roger.

Record Information

Source Institution: UFRGP
Rights Management: Applicable rights reserved.
Classification: lcc - LD1780 2010
System ID: UFE0042230:00001

Permanent Link: http://ufdc.ufl.edu/UFE0042230/00001

Material Information

Title: A Computational Study of Flow and Mass Transport in Intrathecal Drug Delivery
Physical Description: 1 online resource (54 p.)
Language: english
Creator: Fletcher, Adam
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2010

Subjects

Subjects / Keywords: cerebrospinal, computational, concentration, delivery, drug, dynamics, element, finite, fluid, human, intrathecal, mass, morphine, nerve, project, roots, sas, space, subarachnoid, transfer, transport, visible
Biomedical Engineering -- Dissertations, Academic -- UF
Genre: Biomedical Engineering thesis, M.S.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation

Notes

Abstract: A COMPUTATIONAL STUDY OF FLOW AND MASS TRANSPORT IN INTRATHECAL DRUG DELIVERY The detailed geometry of the spinal cord and its surrounding tissues make flow of the cerebrospinal fluid highly variable. Intrathecal drug delivery requires that drugs be injected into the subarachnoid space, which surrounds the spinal cord. In the present work we simulated the flow of cerebrospinal fluid and the mass transfer of morphine in a 29 mm section of the spine, by employing computational fluid dynamics. The Visible Human Project dataset is used to generate a mesh of the subarachnoid space and MRI data is used to apply a physiological pulsatile velocity condition at the inlet of the model. Morphine injection is modeled as a saturated point source with a concentration of 1.0. Results from the study show velocities peaking at 20 cm/s with phase shifts up to 5 degrees. The average Womersley number was calculated to be 3.5 with a deviation of about 1. Spatial velocity profiles show varying flow development in different regions, with a tendency to be blunted during systole. Reynolds numbers were calculated between 200 and 500 during peak systole, and agree well with the literature. After 1.65 seconds, the average relative morphine concentration was 4.4E-4. Dispersion was highly complex, involving circulation about the spinal cord, redirection around nerve roots, and a compression of the subarachnoid space. The Peclet and Schmidt numbers reveal that diffusion was dominated by inertial and viscous effects.
General Note: In the series University of Florida Digital Collections.
General Note: Includes vita.
Bibliography: Includes bibliographical references.
Source of Description: Description based on online resource; title from PDF title page.
Source of Description: This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility: by Adam Fletcher.
Thesis: Thesis (M.S.)--University of Florida, 2010.
Local: Adviser: Tran-Son-Tay, Roger.

Record Information

Source Institution: UFRGP
Rights Management: Applicable rights reserved.
Classification: lcc - LD1780 2010
System ID: UFE0042230:00001


This item has the following downloads:


Full Text





A COMPUTATIONAL STUDY OF FLOW AND MASS TRANSPORT IN INTRATHECAL
DRUG DELIVERY




















By

ADAM D. FLETCHER


A THESIS PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
MASTER OF SCIENCE

UNIVERSITY OF FLORIDA

2010

































2010 Adam D. Fletcher

































I dedicate this work to all of my family and friends.









ACKNOWLEDGEMENTS

I would like to thank Dr. Malisa Samtinoranont, Dr. Roger Tran-Son-Tay, and Dr. Paul

Carney for serving as my supervisory committee. They have provided ample support, counsel,

and encouragement throughout my studies. Without them, I could not have accomplished this

work.











TABLE OF CONTENTS

A CK N O W LED G EM EN T S ......................................................................... ........................ 4

L IS T O F T A B L E S ...................................... ..................................................... .. 7

LIST OF FIGURES ........................ ....... .. .. .... .... ...... ...........

L IST O F A B B R E V IA TIO N S............................................................................... .. .............10

CHAPTER

1 IN TR O D U C TIO N ............... .............................. ..................... ........ .. 13

B ack g rou n d ................... ...................1...................3..........
L iteratu re R ev iew .............................................................................14

2 M E T H O D S .......................................................................................................19

M odel Generation ................................................................ ..... ..... ......... 19
Im age Preparation ............................. ......................................... ........ 19
Stereolithography (STL) Surface Generation ............................................................... 20
Non-Uniform Rational B-Spline (NURBS) Surface Generation ..................................20
U nits and Coordinate System ............................................... ............... 21
N um erical Im plem entation ................................................................................................ 2 1
F lu id M o d el ..............................................................................2 1
M a ss M o d el ..............................................................................2 3
M esh G generation ....................................................... 24

3 R E SU L T S .............. ... ................................................................29

D efin itio n s .............. ... ................................................................2 9
F lu id C h aracteristic s ............................................................................................................... 2 9
T em poral V velocity P profiles ........................................................................................ 29
Spatial V elocity Profiles.................................................. 30
D im en sionless P aram eters .......................................................................................... 30
V olu m etric F low R ate ................................................................................................ 3 1
M a ss C h aracteristic s ............................................................................................................... 3 2
C o n c e n tratio n s ............................................................................................................ 3 2
Dim ensionless Param eters ................................. .......................... ............ 32
D ispersion ............................................. 33

4 D IS C U S S IO N ........................................................................................................4 5

Anatom y .................. ........................................ ........ 45
Subarachnoid Space (SA S) Com pression ....................................................... 45
N e rv e R o o ts ........................................................................................................4 5









Flow Characteristics ................................. .. ... .... ................... 46
T em poral Phenom ena ..............................................................................................46
Flow Development ................................... ..... .. ...... ............... 46
M odel V alidation ......... ............................................................................ ...... ......4 7
M ass Characteristics .................. ............ ...................... ..........48
D isp version ............... ...................................................... ......... ....... 4 9
D iffusive Properties ................ ................ ........................... ..... ....... 49

5 C O N C L U SIO N S ................................................................52

L IST O F R EFER EN CE S .......................................................... ...................... ............... 53

B IO G R A PH IC A L SK E T C H ................................................ .............................. .....................54









































6









LIST OF TABLES

Table page

3-1 Values for relative morphine concentrations and their average, at each region
sam p led (t= 1.64 7 s)..................................................................... .. 34

3-2 Maximum and average Peclet numbers over the cardiac cycle for each sampling
location. The Schmidt number is included as well, and is a constant based on
m material properties .................. ........... .................... ........... 34

4-1 Temporal parameters relating CSF and cardiac systoles, measured using phase-
contrast M RI. (Source: Henry-Feugeas et al. [13])................................ ............... 51









LIST OF FIGURES


Fire pae

1-1 Pump and catheter implantation for intrathecal drug delivery (L: posterior view; R:
saggital view. Sources: M edtronic and M ayfield Clinic) ............................................. 16

1-2 Saggital close-up of catheter tip and surrounding anatomy. (Source: Mayfield Clinic) ...16

1-3 Spatial velocity profiles for the circular annulus at different times during the cardiac
cycle. A: external and internal radius are 10 and 3 mm, respectively, a=16. B:
external and internal radius are 10 and 9 mm, respectively, a=2.3. (Adapted from
[6 ])............. ....... ......................... ................................................1 7

1-4 Comparison of analytical velocity profiles along the minor axis of the annular cross
section with CFD over a complete oscillatory cycle with period T. [7]..........................17

1-5 Stockman's 4 models used to show the difference in CSF flow by varying fine
structures. Only A and D were used in his dispersion study. (Adapted from [8])............18

1-6 Dispersion differences after the time shown in the right-hand column. Models D and
A refer to those in Fig. 1-5. (Adapted from [9]). .................................... .................18

2-1 The region modeled for this study is represented as the section between the two
horizontal lines that intersect the T10 and T11 vertebrae ............................................. 26

2-2 Patches and grids used to construct a NURBS surface. Grid density lowered for
better visualization ......................................................................... ..... ..... .. 27

2-3 Cerebrospinal fluid velocity over one cardiac cycle at the T10 vertebral level ...............27

2-4 Spread of analgesia after intrathecal administration in the lumbar cistern. Values
shown in parentheses are octanol:water partition coefficients; a larger value indicates
that a drug is m ore lipid soluble [12] ........................................... .......................... 28

2-5 Finite element mesh of SAS model, using 4-node tetrahedral FCBI elements. The
velocity and mass loadings are shown in pink .................................................28

3-1 Sampling locations for velocity and mass values denoted by colored elements ...............35

3-2 Maximum velocity nodes depicted as triangles in these band plots................................35

3-3 Y-Velocity plots at the locations shown in Figure 3-1. ............................................. 36

3-4 Spatial velocity profile sampling locations. The caudal-left reference does not lie on
an edge as depicted. A body is absent from this view, and the true shape of the
section can be seen in Figure 3-1 ..................................................................... 37









3-5 Spatial velocity profiles in each region. All series are representative of the y-
velocity at one tenth of a cardiac cycle (T/10, 2T/10, etc.) ............................................ 38

3-7 Flow entrance length at the dorsal aspect of the model inlet and exit.............................39

3-8 Reynolds plots at each sam pling location..................................... ......................... 40

3-9 Womersley number at the cephalic, center, and caudal cross sections............................41

3-10 Error between the model's inlet and exit volume flow rates ..........................................41

3-11 Comparison of the model's volume flow rate with Loth's MRI measured flow rate........42

3-12 Relative concentrations of morphine at each of the 12 sampling regions. Values are
over the course of 1.4 cardiac cycles. ........................................ .......................... 43

3-13 Relative concentration of morphine at time 1.647 seconds. This was the final mass
transfer time step, and represents 1.4 cardiac cycles. .............................. ......... ...... .44

3-14 Cephalic (left) and caudal (right) cross sections showing the velocity field at peak
systole. Note the counterclockwise and clockwise trends. .............................................44

4-1 VHP images showing compression of the SAS at the T10-T11 vertebral disc. Also
note the presence of nerve roots, which obstruct the flow. LEFT =model inlet,
RIGHT = 12 mm inferior. (Source: Visible Human Project [5])...................................51









LIST OF ABBREVIATIONS


Acs cross sectional area

C cervical level

CSF cerebrospinal fluid

cm centimeters

D characteristic length scale (annular gap)

Dh hydraulic diameter

Di diffusivity tensor of ith species

f body force on fluid

f temporal frequency

FCBI flow-condition-based-interpolation

ITDD intrathecal drug delivery

L lumbar level

le entrance length

m meters

mi mass of ith species

mi'" mass creation rate of ith species

MRI magnetic resonance imaging

mm millimeters

ms milliseconds

n number of species

NURBS non-uniform rational b-spline

p fluid pressure









Pe Peclet number

Pwet wetted perimeter

Re Reynolds number

Rh hydraulic radius

s seconds

SAS subarachnoid space

Sc Schmidt number

STL stereolithography (format)

t time

T thoracic level

u axial component of fluid velocity

v fluid velocity

vi velocity of ith species

VHP visible human project

a Womersely number

0 phase angle

It fluid viscosity

p fluid density

(Di mass ratio of ith species

co angular frequency









Abstract of Thesis Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Master of Science

A COMPUTATIONAL STUDY OF FLOW AND MASS TRANSPORT IN INTRATHECAL
DRUG DELIVERY

By

Adam D. Fletcher

August 2010

Chair: Roger Tran-Son-Tay
Major: Biomedical Engineering

The detailed geometry of the spinal cord and its surrounding tissues make flow of the

cerebrospinal fluid highly variable. Intrathecal drug delivery requires that drugs be injected into

the subarachnoid space, which surrounds the spinal cord. In the present work we simulated the

flow of cerebrospinal fluid and the mass transfer of morphine in a 29 mm section of the spine, by

employing computational fluid dynamics. The Visible Human Project dataset is used to generate

a mesh of the subarachnoid space and MRI data is used to apply a physiological pulsatile

velocity condition at the inlet of the model. Morphine injection is modeled as a saturated point

source with a concentration of 1.0.

Results from the study show velocities peaking at 20 cm/s with phase shifts up to 5 degrees.

The average Womersley number was calculated to be 3.5 with a deviation of about 1. Spatial

velocity profiles show varying flow development in different regions, with a tendency to be

blunted during systole. Reynolds numbers were calculated between 200 and 500 during peak

systole, and agree well with the literature. After 1.65 seconds, the average relative morphine

concentration was 4.4E-4. Dispersion was highly complex, involving circulation about the

spinal cord, redirection around nerve roots, and a compression of the subarachnoid space. The

Peclet and Schmidt numbers reveal that diffusion was dominated by inertial and viscous effects.









CHAPTER 1
INTRODUCTION

Background

Implantable intrathecal drug delivery (ITDD) systems are becoming increasingly

common [1]. They work by administering drugs directly to the intrathecal space surrounding the

spinal cord. This eliminates some side effects that are present with intravenous and oral

administration. In addition, smaller doses can be used which is cheaper and healthier for the

patient [2]. Using implantable drug pumps also allows doctors to precisely program the size and

timing of doses which increases the quality of life for the patient.

Two common types of drugs that are delivered intrathecally are muscle relaxants and

analgesics. Morphine is a common analgesic that is used in treatment of a variety of pain

conditions. Sometimes when morphine is injected in a nonhomogeneous fashion, the result can

be a high dose to a localized part of the spinal cord, which may result in injury [3]. Baclofen is

the commonly prescribed muscle relaxant for ITDD. It is used to treat muscle spasticity in

conditions such as cerebral palsy. We investigated morphine in this study because it is

prescribed more often and has a larger incident of side effects.

A typical procedure for implanting the pump and catheter involves insertion of an

introducer needle between the third and four lumbar levels (L3-L4) and delicately guiding the

catheter through the needle and up to the eleventh vertebral level (T11). Figure 1-1 shows a

posterior and saggital view of the implant. Once inserted, the catheter is anchored to the

thoracolumbar fascia and any excess catheter is trimmed away. A pocket for the drug pump is

created above the ilium. The catheter is connected to the pump, and the pump is then sutured

into the pocket. An additional illustration of the catheter tip and surrounding anatomy is

provided in Figure 1-2.









Literature Review

One of the first numerical investigations into the distribution of anesthesia in the intrathecal

space was performed by Myers [4]. His technique involved using the finite element (FE) method

to simulate the dispersion of anesthetics during spinal anesthesia. Spinal anesthesia is performed

between the third and fourth lumbar vertebral levels, and has the potential to result in cauda

equina syndrome. This occurs when high concentrations of the analgesic invade the cauda

equina. His results showed that the size of the needle used for injection was a critical parameter

in dispersion of the anesthesia. Larger needles resulted in higher dispersion.

Francis Loth has been a strong contributor to the study of ITDD. For one of his works,

his group used the visible human project (VHP) [5] magnetic resonance imaging (MRI) datasets

to measure the geometry of the intrathecal space [6]. From these measurements they made

concentric and eccentric elliptical cross sections as well as concentric circular cross sections for

their numerical computations. To define the CSF flow waveform, dynamic phase-contrast MRI

was performed on a healthy adult at the second cervical vertebral level (C2). This provided the

CSF volumetric flow rate over the cardiac cycle. They found the flow in their models to be

laminar and inertially dominated with Reynolds numbers on the order of 102. Simulations of

larger annular gaps between the cord and dura mater resulted in velocity profiles that were

blunted at the center and peaked at the walls. This type of flow is typical of high Womersley

numbers (a) and probably due to flow reversal before the CSF can accelerate to a full developed

state. Smaller annular gaps resulted in lower Womersley numbers and a more typical profile.

Figure 1-3 displays these phenomena. Loth's results for the concentric elliptical annulus were

verified by Gupta et al. [7], who found the analytical solution for pulsatile flow in an elliptical

annulus. Figure 1-4 shows his comparison to Loth's work.









Previous studies by Stockman have shown that fine anatomic structures in the

subarachnoid space (SAS), such as dorsal and ventral nerve roots, do not significantly alter the

velocity profiles of CSF [8]. However, they are responsible for a more rapid dispersion of

injected solutes in an oscillatory environment [9]. Figure 1-5 displays the geometry used in these

experiments. Stockman used MRI data from Loth [6] to define his cross sectional geometry. For

his flow simulations the length of the model was 1 cm. His dispersion simulations were done

with a model length of 4-8 cm. Figure 1-6 provides the reader with a visual of the difference in

dispersion between the models containing and absent of fine structures. It was concluded that

fine structures enhance dispersion by 5 to 10 times over the range of Schmidt numbers (Sc) that

were used.

Our mission was to simulate an ITDD system using the high resolution VHP dataset.

This allowed us to more accurately model the intrathecal space with fine structures. The high

resolution dataset has not been used in previous literature and should provide a more realistic

response due to better accuracy of the model geometry. A pulsatile flow was used because this

best emulates the oscillatory motion of CSF. This periodic motion occurs due to expansion of

the cerebrovascular bed during cardiac systole. To achieve the proper physiological oscillation,

Loth's MRI-generated waveform for CSF flow at the C2 level was used for our velocity input.

Morphine injection was simulated by a mass concentration point source on the dorsal aspect of

the intrathecal space. With this combination of techniques we attempt to take a new step in

ITDD research.


























Figure 1-1. Pump and catheter implantation for intrathecal drug delivery (L: posterior view; R:
saggital view. Sources: Medtronic and Mayfield Clinic)





spinal
.:--cord
SI sarsrachnoid

catheter

g hw mater

+; n Ilr ttfepkdanl

SMayfed Cliime

Figure 1-2. Saggital close-up of catheter tip and surrounding anatomy. (Source: Mayfield Clinic)











S0.01

0.02
0.015 --- r-----.......----...

0.015



-0.005

-0.01
0.003 0-004 000 0.006 0.007 0.008 0.009 001
r Lnm


B o0.



CJ.I




-0.05

0.
-0-1
0-0


Figure 1-3.


0
s

d


i9


O.0092 0.0094 0.0096 0.0098 0.01
r [n]


Spatial velocity profiles for the circular annulus at different times during the cardiac
cycle. A: external and internal radius are 10 and 3 mm, respectively, a=16. B:
external and internal radius are 10 and 9 mm, respectively, a=2.3. (Adapted from
[6])


t= 0
t=3T/5 -6
0 0 f .. : ---- ,--
SI t=2T/5


j -0.01

\t=T/5


-0.02
4 6 8 10
y (m) x 103

Figure 1-4. Comparison of analytical velocity profiles along the minor axis of the annular cross
section with CFD over a complete oscillatory cycle with period T. [7]


--------------- ----------------













tx A

1r
Nerve
Bundle
zf B

Denticulate
Ligament

x C




zD


Figure 1-5. Stockman's 4 models used to show the difference in CSF flow by varying fine
structures. Only A and D were used in his dispersion study. (Adapted from [8]).







0.24 s
11.04 s




Model D Model A
Figure 1-6. Dispersion differences after the time shown in the right-hand column. Models D and
A refer to those in Fig. 1-5. (Adapted from [9]).










18









CHAPTER 2
METHODS

Model Generation

The first phase of this study was to construct a surface for the finite element model. The

surface was based upon high resolution male cadaveric images from the VHP. The region spans

from the bottom one-fourth of T10 to half of T 1, as seen in Figure 2-1.

Image Preparation

The high resolution VHP dataset consists of images taken from a cryogenically frozen male

cadaver. All of the images are axial cross sections that are spaced 1 mm apart. Each slice was

photographed with 70 mm film and subsequently scanned into digital RAW images. The pixel

size for each image is 0.144 x 0.144 x 1.000 mm and the color quality is 24 bit. 30 images were

used to create a model 29 mm long.

In order to convert the images from RAW to bitmap, a batch process was performed

using Photoshop (Adobe Systems Inc., San Jose, CA). The images contained fiducial markers

from a rod situated in the frozen gel that surrounded the cadaver. However, the rod had a bow to

it and was not perfectly true, rendering the fiducial markers invalid. A custom program was

written to register the images by sampling the edge of the frozen gel and performing a least

squares regression. It then aligned the images rotationally and translated them so the edge of the

gel was flush with the left side of the image. Minor translational corrections were made

manually to finalize the registration. These were done using the cord as a fiducial; therefore our

model does not have any curvature to it. However, in this section of the spine, the cord is linear

and it was determined that the centroid of the cord only moved by 14 pixels over the length of

the model, which is about 4.6 mm, and assumed to be negligible.









After registering the dataset, a small MATLAB (The MathWorks Inc., Natick, MA) script

was written to filter out the excess anatomy that did not define the fluid space. This was done by

reading through the pixels and converting to 0 those that did not fall into predefined ranges of

red, green, and blue. MATLAB's imtool was used to assess the color values of anatomy that

defined the fluid space. This provided the range of values for the filter to exclude. After the

filters were run, manual corrections were performed on pixels that did not agree with the filtered

values but were indeed fluid space.

Stereolithography (STL) Surface Generation

A three dimensional model was created from the 30 two dimensional images using Amira

(Visage Imaging, Inc., San Diego, CA) software. Each image was segmented in order to define

the fluid space to be used in the model. Much of this was already completed in the filtering

process, but a few modifications were necessary because the software did not handle thicknesses

of 2 or less pixels very well. Sharp corners were also removed with the segmentation editor so

that it would be more practical to generate a NURBS surface after exporting the STL. Next, the

surface was generated, consisting of 98740 triangles. Unconstrained smoothing was used to

soften the surface, and the ends of the annulus were closed off. The surface was exported in

little-endian STL format.

Non-Uniform Rational B-Spline (NURBS) Surface Generation

In order to apply boundary conditions to the model, we must define areas known as

"patches" or "faces" to refer to. A NURBS surface is a collection of patches that consist of B-

splines. The patches are populated with grids that precisely represent the surface to be built.

Geomagic (Geomagic, Inc., Durham, NC) was used to construct the patches on our model.

Please see Figure 2-1 for a visualization of what this looks like on our model. Prior to adding the

patches the model had to be refined. This was done by increasing the number of triangles that









make up the surface to 278212. Several smoothing techniques were then utilized to eliminate the

striations caused by the Imm gap between each image. The ends of the annulus were cut to

create flat surfaces that would be easily identifiable as slip boundary conditions, and also to

identify loading conditions. After these steps were completed, the model was patched and grids

were placed to create the final NURBS surface.

Units and Coordinate System

Using Geomagic, the length units of the surface were established to be SI, and a

Cartesian coordinate system was defined for the model. The origin was placed in the center of

the cord at the mid point of its length. The spinal cord was aligned to the y-axis with positive

being in the caudal direction, and the positive z-axis aligned to the ventral direction. Positive x-

axis is directed to the cadaver's anatomical left. The surface was exported in IGES file format.

Numerical Implementation

The model was solved numerically using the FE package ADINA-F (ADINA R&D Inc,

Watertown, MA). The FE equations were discretized using the Galerkin method. Due to the

transient nature of the problem, the Euler method was used to perform time integration. Solving

the discretized FE equations was done using the Newton-Raphson method.

Fluid Model

To solve for the flow field within the reference and anatomic annuli, the Navier-Stokes mass

continuity and momentum equations were used. The mass continuity equation is shown below

(p density, t=time, v=velocity)

+ V (pv) = 0 (2-1)
at

The fluid was approximated as water since, according to Stockman, viscosity data on CSF is

clustered toward being more like water [9], and the density of water is within 0.4% of values









measured by Schiffer et al [10]. Being that water is a liquid, and liquids at this low of a speed (<

ap
8 cm/s) can be assumed incompressible, the density is constant and = 0. This reduces the
at

continuity equation to

V v = 0 (2-2)

Due to the low velocities encountered in this problem, the flow field was assumed to be laminar

and the conservative form of the momentum equation was used (p=pressure, u=axial velocity,

f=body force)


p -+v* Vv -Vp +uV2v+f (2-3)
(at

To simplify the problem, the patient was assumed to be lying in the supine position so that body

forces could be neglected. Due to the pulsatile nature of CSF flow, the viscous nature of the

fluid, and varying geometry of the intrathecal space, all terms were kept except body force,

resulting in


p v- + Vv -Vp + uV2v (2-4)


At the walls of all the anatomic structures, including the spinal cord, nerve roots, and arachnoid

mater, a no-slip condition (v=0) was assumed. The energy equation was not solved due to the

assumption of incompressible flow without heat generation.

To simulate motion of the CSF, a velocity loading condition was applied at the cephalic end

of the model. The velocities used were adapted from the volume flow rates measured at the C2

vertebral level by Loth et al. using phase-contrast MRI [6]. Two programs were written to

satisfy this adaptation. The first was to digitize the CSF waveform that was measured via MRI.

This allowed us to convert the image from his publication into discrete values of time and flow









rate. The second program analyzed the VHP image that represented the cephalic end of the

model. It summed all of the pixels that defined the fluid space and, with knowledge of the pixel

size, computed the cross sectional area of the inlet. From here, the digitized waveform was

divided by the inlet area to produce the velocity waveform used for loading. Figure 2-2 shows

the converted CSF velocity waveform. The period (T) is equal to 1.176 seconds. The velocity

was prescribed uniformly at the cephalic inlet of the model, in the axial (caudal) direction. The

simulation was set to run for 5 cardiac cycles, totaling 250 time steps at intervals At = 0.0235 s.



Mass Model

In order to solve for the mass transport of morphine, the mass transport equation had to be

coupled with the flow equations. The mass transport equation is derived by combining Fick's

Law and the species conservation equation. The latter is shown below


O +V .(pv,)=m" (2-5)
at

Where p, is the density of the ith species and m,' is the mass creation rate. Fick's Law represents

the diffusion of a species in a fluid

~ + V. (pv, D, V )= m, (2-6)
aSt

Where Di is the diffusion tensor and (pi is the mass ratio for the ih species. For this study we

assume the diffusivity to be equal in all directions. Clark et al. measured the diffusivity of

morphine in CSF to be in the region of 3 to 5 x 10-6 cm s-1 [11]. We took the average value and

set D = 4 x 10-10 m2s-1 to keep with the SI standard. The mass ratio is a ratio of the density of the

ith species to the total, or "bulk" density of all species combined. The definition is shown below.










n= A (2-7)
ZP,
1=0

Since we are only concerned with one species, morphine, (pi = 1.0 by default. This value also

represents the concentration of the species at the site where the condition is applied. In this case

it is 1.0, or 100% concentrated. Finally the conservation equation and Fick's Law combine to

form the mass transfer equation

p 'V +V. (pv, -pD, V)= m" (2-8)
at

To best approximate the conditions of ITDD, the species was input on the dorsal side of the

model, directly on the arachnoid mater. This was done because there were no geometric

references in the center of the annulus to use. In reality the injected morphine is in the form of a

solution; however, this was neglected to simplify the problem and dismiss any fluid-fluid

interactions. The spinal cord, nerve roots, and arachnoid mater were assumed to be inert (non-

porous) to morphine because according to Rathmell et al. [12], morphine has an octanol-water

partition coefficient of one, meaning that it is not attracted to lipids. Figure 2-3 shows the

comparison of morphine distribution with other common intrathecal drugs. Note that morphine

diffuses throughout the whole spine while the others do not, implying that it does not readily

react with the tissue. Due to restricted computation time, the mass model was only solved for 1.4

cardiac cycles, but with the same step size as the fluid model.

Mesh Generation

The model was sectioned into 17 bodies in order to align nodes along model lines in regions

of interest. This was done by placing an axial section cut through the center of the model, two

additional axial cuts at 10 mm fore and aft of center, a saggital cut along the center of the model,

and a coronal cut 1.5mm from center in the ventral direction. Each body was subdivided into 4-









node tetrahedral Flow-Condition-Based-Interpolation (FCBI) elements, with a characteristic

length 0.5 mm. This made for a total of 129542 elements. A Delaunay meshing algorithm was

used on all elements and boundaries. In one case Advancing Front had to be used for boundary

meshing. The total number of nodes generated was 27316. Figure 2-4 shows the meshed model

including velocity and species loading.

















It rlrer.I kw


Figure 2-1. The region modeled for this study is represented as the section between the two
horizontal lines that intersect the T10 and T11 vertebrae.




































Patc- 2.42
Cuent Tn- ,e 278.212


Figure 2-2. Patches and grids used to construct a NURBS surface. Grid density lowered for
better visualization


10


8


6


4
,a

2
S2
C)


Time (s)


Figure 2-3. Cerebrospinal fluid velocity over one cardiac cycle at the T10 vertebral level































Figure 2-4.


Spread of analgesia after intrathecal administration in the lumbar cistern. Values
shown in parentheses are octanol:water partition coefficients; a larger value
indicates that a drug is more lipid soluble [12].


Figure 2-5. Finite element mesh of SAS model, using 4-node tetrahedral FCBI elements. The
velocity and mass loadings are shown in pink.









CHAPTER 3
RESULTS

Definitions

Much of the flow and mass transfer data was collected from three axial cross sections of the

model and used as reference planes for extracting data. These coincide with the axial cuts

described in the section "Mesh Generation", and the center cut coincides with the morphine

source at the dorsal wall. For the remainder of the paper, the data related to the center cut will be

prefixed with "center", and the two cuts fore and aft of center will be prefixed with cephalicc"

and "caudal", respectively. Figure 3-1 shows these three cuts and the locations were the velocity

was sampled. The element colors correspond to the colors used in the figures that follow. The

terms "left" and "right" refer to the cadaver's anatomical left and right sides. All data series with

a "max" suffix refer to the node within the cross section that had the maximum velocity during

peak systole. Their locations are shown in Figure 3-2 and are represented as triangles on the

velocity band plot. This is a plot of the axial velocity at peak systole. By comparing Figure 3-1

and 3-2, the reader can also identify where the sampling locations are in respect to the velocity

field.

Fluid Characteristics

The results of the fluid portion of the simulation are presented here. The solution required 98

hours to process 250 time steps, equaling 5 cardiac cycles. Temporal and spatial velocity

profiles from 12 sites along the length of the model were constructed over one cardiac cycle.

Dimensionless numbers are given, as well as model validation by the volumetric flow rate.

Temporal Velocity Profiles

Figure 3-3 presents graphs of the temporal velocity profiles. The velocity prescribed at the

inlet is included as a reference. One remarkable feature of these plots is that the waveforms are









out of phase with one another. Due to the low temporal resolution of this simulation (23.5 ms),

cubic spline interpolation was used to determine the times of maximum axial velocity at each

location. The phase shift (0) for each location was computed relative to the prescribed velocity.

0 = 360fAt (3-1)

Wherefis the frequency of the cardiac cycle and At is the time shift at peak velocity. The results

are shown in Figure 3-4.

Spatial Velocity Profiles

Spatial velocity profiles were constructed from the locations highlighted in Figure 3-5. These

planes are the same as those shown in Figure 3-1. For each location, ten profiles of the axial

velocity were plotted over the cardiac cycle. Figure 3-6 displays these plots. Dorsal/ventral and

left/right locations are overlaid to conserve space. The x-axis of all plots begins at zero and

increases in the model's positive x- or z-direction, which are the cadaver's left and ventral

directions, respectively.

The entrance length (le) was computed as a supplement to help assess the flow

development. The following equation for entrance length in pipe flow was used.

le = 0.06ReD (3-2)

The prescribed velocity was used to calculate the Reynolds number (Re) and the annular gap at

the dorsal aspect, between the cord and arachnoid mater, was used as the characteristic length

scale (D=2.0-2.1 mm). Figure 3-7 plots the entrance length against the input velocity. The

prescribed velocity was integrated and plotted to provide the reader with an idea of how far the

fluid moves, neglecting viscous effects.

Dimensionless Parameters

Reynolds number. This number describes the ratio of inertial to viscous forces acting on a

flow. It was computed to help validate the model against the literature.










Re= pVD (3-3)
/1

The hydraulic diameter (Dh) was used as the length scale for consistency with Loth et. al [6] and

Stockman [8]. The formula for the hydraulic diameter was borrowed from Loth. It is a function

of the annulus cross sectional area (Acs) and the wetted perimeter of the cross section (Pwet).

4A
Dh = (3-4)
Pwet

The wetted perimeter and cross sectional area were calculated using Solidworks 2008 (Dassault

Systemes, Velizy, France). Figure 3-8 gives the varying Reynolds number over one cardiac

cycle.

Womersley number. This value (a) represents the relationship between an oscillatory fluid's

frequency and its viscosity.


Rh P Y2 (3-5)


Rh is the hydraulic radius (Rh = Dh/2 ), p=density, p=dynamic viscosity, and co is the frequency

of the flow oscillation in rad/s. For our simulation o=5.34 rad/s. The Womersley number can

indicate whether a pulsatile flow will reach a fully developed state before it reverses, and is a

common component ofbiofluid analysis. The values for this study were computed at the

cephalic, center, and caudal cross sections referred to previously. Figure 3-9 displays the three

values.

Volumetric Flow Rate

As a means of validating the model, volumetric flow rates were extracted from the solution.

The flow rates were computed 1 mm from the inlet and exit of the annulus. Figure 3-10 shows

the percent error between the two values, plotted against the prescribed velocity. As another









check, the simulation's volumetric flow rate was compared to Loth's MRI measured flow rate.

The result is shown in Figure 3-11.

Mass Characteristics

The mass transport model gives the reader an idea of how morphine will disperse in the

intrathecal space of the cadaver, and possibly other patients with similar anatomy. All of the

mass solutions reported in this text are relative concentrations. That is, they represent the

concentration of morphine relative to the fully saturated point source, which is p=1.0.

Concentrations

The transient nature of the morphine concentrations were plotted over the solution time of

1.647 seconds and can be found in Figure 3-12. For reference, the first systole occurs at 0.35

seconds and the second systole occurs at 1.53 seconds. Table 3-1 summarizes the mass

distribution at the 70h (final) time step for each region. The dorsal and right regions are more

steady, while the ventral and left regions vary more.

Dimensionless Parameters

Schmidt number. To assist in the analysis of the mass transport results, two dimensionless

parameters were calculated. The first is the Schmidt number


Sc = (3-6)
pD

Which describes the ratio of viscosity to diffusivity in a flow. In our case, the diffusive tensor is

a constant for all directions, and CSF was modeled as a Newtonian, incompressible fluid, so t

and p are both constants as well. Therefore the Schmidt number for this problem is a constant

value itself, Sc = 2450.

Peclet number. The second parameter calculated was the Peclet number









VD
Pe = ReSc = h (3-7)
D

This value represents the ratio of inertial to diffusive forces in a flow. The maximum and

average Peclet numbers over one cardiac cycle were tabulated and can be found in Table 3-2.

Graphs were not constructed because in our case, the Peclet number is just the Reynolds number

multiplied by a constant, and temporally it would look the same as the Reynolds and velocity

plots.

Dispersion

As an additional tool in understanding how the morphine dispersed, a band plot was

produced for each cross section and can be found in Figure 3-13. This plot is representative of

the last time step in the series (t=1.647 s). The morphine disperses in a complex manner due to

the anatomy in the SAS. Noting that the dorsal and right regions consistently have higher

concentrations than the ventral and left regions, a vector plot was constructed to help identify

why the morphine dispersed in this way. Figure 3-14 presents the velocity fields at the cephalic

and caudal cross sections during systole.










Table 3-1. Values for relative morphine concentrations and their average, at each region
sampled (t=-.647 s).

cross section \ region dorsal right ventral left average
cephalic 6.4E-04 2.3E-07 2.6E-09 3.1E-08 1.6E-04
center 7.7E-04 2.3E-07 9.2E-10 1.8E-08 1.9E-04
caudal 1.8E-03 3.8E-06 7.5E-08 3.6E-10 4.4E-04



Table 3-2. Maximum and average Peclet numbers over the cardiac cycle for each sampling
location. The Schmidt number is included as well, and is a constant based on
material properties


Peclet Number


dorsal


Region Sampled


ventral


right


max


Schmidt Number 2450 constant


9.2E+05 6.1E+05 8.8E+05 8.6E+05 9.4E+05 max
cephalic
3.1E+05 1.8E+05 3.2E+05 2.9E+05 3.0E+05 avg
Cross
Cross 9.0E+05 1.0E+06 1.0E+06 8.4E+05 1.2E+06 max
Section center
Location 3.0E+05 3.2E+05 3.2E+05 2.6E+05 3.9E+05 avg
caudal 1.1E+06 4.9E+05 1.0E+06 5.3E+05 1.1E+06 max
caudalE+05 1.3E+05 3.4E+05 1.8E+05 4.1E+05 avg
4.0E+05 1.3E+05 3.4E+05 1.8E+05 4.1E+05 avg


2450


constant


Schmidt Number


























Figure 3-1. Sampling locations for velocity and mass values denoted by colored elements


Figure 3-2. Maximum velocity nodes depicted as triangles in these band plots.





















- cephalic-dorsal
- cephalic-ventral
- cephalic-right
cephalic-left
- cephalic-max
- input-velocity


Time (s)


- center-dorsal
- center-ventral
- center-right
center-left
- center-max
- input-velocity


Time (s)


- caudal-dorsal
- caudal-ventral
- caudal-right
caudal-left
- caudal-max
- input-velocity


Time (s)


Figure 3-3. Y-Velocity plots at the locations shown in Figure 3-1.


25


20


15


10


5


0


-5


-10


































Figure 3-4. Velocity phase shifts relative to the prescribed velocity at the inlet


Figure 3-5.


Spatial velocity profile sampling locations. The caudal-left reference does not lie on
an edge as depicted. A body is absent from this view, and the true shape of the
section can be seen in Figure 3-1.


60

40

?2 2 mcephalic
( m center




-4


Cross Section Location


















10 10



5 5



0 0
>0 5

-5



-10 -10
Distance (mrro cephalic-dorsal Distance (nrn) cephalic-right
cephalic-ventral cephalic-left

20 25




15
*; 10



5

0


-5


-10 -10
Distance (inm) center-dorsal Distance (wn) center-right
center-ventra center-lef

14 12

12 10

10
8
8
6



2





-4 -4

-6 -6
Distance (nm) caudal-dorsal Distance (mn) caudal-right
-caudal-ventral caudal-left



Figure 3-6. Spatial velocity profiles in each region. All series are representative of the y-

velocity at one tenth of a cardiac cycle (T/10, 2T/10, etc.)





















-Particle Distance
-Prescribed Velocity
Entrance Length


Time (s)


Figure 3-7. Flow entrance length at the dorsal aspect of the model inlet and exit.


35

30

E|E 25
aE
20

> 2 15
1-
10







-10
U 5

a-



-5

-10






















- cephalic-dorsal
- cephalic-ventral
- cephalic-right
cephalic-left
- cephalic-max


400


z 300
N


S200


100


0
0.0


- center-dorsal
-center-ventral
- center-right
center-left
- center-max


0.0 0.2 0.4 0.6 0.8 1.0 1.2
Time (s)


- caudal-dorsal
- caudal-ventral
- caudal-right
caudal-left
-caudal-max


0.2 0.4 0.6 0.8
Time (s)


1.0 1.2 1.4


Figure 3-8. Reynolds plots at each sampling location


0.2 0.4 0.6 0.8 1.0 1.2
Time (s)


z 300

C,
200


100


500

450

400

350

300

250

200

150

100



50
0.0





























Cephalic
3.3


Center


Caudal
4.4


Womersley number at the cephalic, center, and caudal cross sections


Figure 3-10. Error between the model's inlet and exit volume flow rates


, n .


Figure 3-9.


0.45% 10

040% 8

0.35%
6
0 30%

0.25% 4
L 4
W 0 20% 2

0.15% >
0
0 10%

0.05% -2

0.00% -4
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4
Time (s) Flow Rate Error
Prescnbed Velocity


I


I













7

6

5

4

3

2


E

O
-1

-2

-3
Tirme (s) MRI Measured
-Model Solution


Figure 3-11. Comparison of the model's volume flow rate with Loth's MRI measured flow rate















1 00E+01
1 00E-01
1 00E-03
1 00E-05
c 1.00E-07
1 DE-09
1 00E-09 cephalic-dorsal

8 1 00E-11 cephalc-ventral
S 1 00E-13 cephalic-right
1.00 E-15 cephallc-left
1.00E-17
1.00E-19
1.00E-21
1.00E-23
1.00E-25
Time (s)



1 00E+01
1 00E-01 8
1 00E-03
1 00E-05

1 .00E-07
S1.00E-09.
1 DD-D -center-dorsal
1. 00E-11 -c enter-ventral
1.00E-13 center-right
1.00E-15 center-left
S1.00E-17
1.00E-19
1 00E-21
1.00E-23
1.00E-25
Time (s)



1 00E+01
1 00E-01 8
1 00E-03
1 00E-05
c 10E-07
ODE09 cauda-dorsal

1.00E-11 --- caudalvetral
1 00E-13 cauda-right
1 00E-15 caudakleft
S100E-17
1 00E-19
1 00E-21
1 ODE-23

1 00E-25
Tim e (s)


Figure 3-12. Relative concentrations of morphine at each of the 12 sampling regions. Values are

over the course of 1.4 cardiac cycles.
































Figure 3-13. Relative concentration of morphine at time 1.647 seconds. This was the final mass
transfer time step, and represents 1.4 cardiac cycles.


Figure 3-14. Cephalic (left) and caudal (right) cross sections showing the velocity field at peak
systole. Note the counterclockwise and clockwise trends.









CHAPTER 4
DISCUSSION

Anatomy

To fully understand the analysis of this model, it is important to note the anatomical features

of the cadaver. Those that are particularly important are the nerve roots found within the SAS,

and equally as important, the compression of the SAS at the T10-T11 vertebral disc.

Subarachnoid Space (SAS) Compression

Taking a look at Figure 3-1, the reader should note that the ventral side of the center cross

section has an annular gap much smaller than the other sections. For better visualization of why

the cross section narrows so much at the center, two of the original VHP images used to generate

this model are displayed in Figure 4-1. The left-hand image represents the model inlet, and the

other is 12 mm from the inlet. In the first image there is pink trabecular bone ventral (bottom) to

the SAS, and in the second image the nucleus pulposus of the T10-T11 disc can be seen on the

ventral side. The presence of the disc at the location of compression may indicate a geometry

typical of a slightly herniated disc.

Nerve Roots

Also apparent in Figures 3-1 and 4-1 is the large nerve root on the left side of the SAS. It

entirely blocks off flow between the dorsal and left regions in this section. Therefore, the CSF

that is trying to enter the left side chamber is redirected to faster moving, lower pressure fields.

Also, the thin clearance at the ventral aspect further isolates the CSF in the left chamber. The

fluid in the left chamber, therefore, is more stagnant and plays less of a role in mixing. The large

nerve root and the cord compression will be the dominant factors in determining the morphine

distribution in this segment of the cord.









Flow Characteristics


Temporal Phenomena

One of the remarkable features of the temporal velocity profiles is their phase delay. Henry-

Feugeas et al. measured similar types of phase delay in her MRI phase-contrast study [13].

Table 4-1 gives the temporal parameters of their work. These values represent the relationship

between the patient's CSF systole and cardiac systole at the thoracolumbar junction. Their

measurements are not intended to agree with this simulation, since we do not know what the

geometries consisted of. Also, one must consider that the standard deviations on Table 4-1 are

somewhat large, which of itself implies that this is a highly variable characteristic of the flow.

These phase shifts are likely the product of flow obstruction due to the cadaver anatomy. This

hypothesis is based on the spatial and temporal velocity data that showed the dorsal and right

regions of the SAS to have higher velocities most of the time. Also, the phase angles calculated

in Figure 3-4 show the dorsal and right regions to have positive shifts, while the ventral and left

aspects have negative shifts.

Flow Development

The spatial velocity profiles in Fig. 3-6 demonstrate the degree of flow development in each

region that was sampled. When viewing these plots note that entrance effects take place at both

the inlet and exit of the model, due to the oscillatory nature of the flow. Thus at the beginning of

the cycle, the velocity is negative and the exit has entrance effects imposed upon it, but after

systole begins the flow becomes positive (caudal) and the inlet is subjected to entrance effects.

Due to the irregular varying geometry, the entrance length was not able to be determined.

The entrance length is a function of both the size and shape of the cross section that the fluid

passes through. Since the geometry varies continuously, so does the entrance length, making it









difficult to estimate. For the sake of a reference, using the solution for pipe flow, the entrance

length was computed and is shown in Figure 3-7.

Profiles with a width greater than 1.5 mm showed signs of being more developed during

diastole. This is due the lower Reynolds numbers during this phase of the cardiac cycle,

effectively reducing the entrance length. During systole the profiles are more blunted, and the

entrance length may be equal to or greater than the length of the model. However, the plots do

agree favorably with those calculated by Loth in Fig. 1-3b. Similar velocities are reported and

the Womersley numbers are close to Loth's, with his a=2.3, and our a =2.8-4.4 (Fig. 3-9). In

this study the annular gap was quite small compared to Stockman [8], and as such his

simulations had a higher Womersley number (a=16) and the profiles were under-developed such

as that in Fig. 1-3a.

Model Validation

Reynolds Number. The Reynolds numbers plotted in Fig. 3-8 were validated against

those in Loth's study. Along the length of the spine, the maximum Reynolds numbers that he

reported were on the range of 175-450. This is consistent with the maximum values in the local

regions of this model, which vary between 200 and 500.

Volumetric Flow Rate. To check the simulation for continuity, the volumetric flow rate

was calculated at the inlet and exit of the annulus and compared. Under the assumption of

incompressible flow, the flow rate entering the annulus should be equal to that exiting the

annulus. Fig. 3-10 displays the error between these two values. The error spikes to 0.26%

during systole when the flow accelerates, and again to 0.41% during diastole when the velocity

waivers around zero. The former can be explained by a low temporal resolution, relative to the

abrupt increase in velocity. During systole, the axial velocity increases by 2-5 cm/s in one time

step, depending on the location sampled. The 0.41% error may be explained by a combination of









the velocity oscillating around zero, and the intrinsic nature of the regions being out of phase. In

this period the prescribed flow reverses twice over four time steps, and discrepancies about the

direction of flow between individual regions can clearly be seen in Fig 3-3. Indeed, when

comparing regions, there are large differences in time for when the flow finally reverses during

diastole.

An additional check on the model was done by comparing the calculated volumetric flow

rate with that measured by Loth. This is shown in Fig. 3-11. The MRI measured flow rate is the

net flow rate after all losses are considered. In this study those losses were not accounted for,

and as such the prescribed velocity was calculated with inherent losses, only to undergo more

loss after entering the annulus. This explains the difference in the flow rates.

A final note about the volumetric flow rate is that over one period it has a non-zero mean.

This implies that there is a net movement of fluid, and in this case in the caudal direction. This is

more easily visualized by the particle distance plotted in Fig. 3-7. By continuity there must be a

sink that accommodates this excess fluid, or a secondary cephalic-directed flow that begins at the

bottom of the cord. The MRI measured the mean velocity over an average slice thickness of 0.5

cm [8]. It may be that slow, secondary flows, that are cephalically directed, do not produce

enough radio frequency (RF) signal to be detected by MRI.

Mass Characteristics

In this simulation, we were only able to solve for a time of 1.65 seconds due to a lack of

computational power and time. As such, we cannot determine what the concentration of

morphine will be at a more useful times of 1-10 minutes after injection. However, we can get a

glimpse at how anatomical structures effect morphine dispersion, and given more computational

time, a follow-up study will report on the concentrations and distribution over longer periods.









Dispersion

The concentration of morphine tended to disperse in a net clockwise fashion around the SAS,

due primarily to the large nerve root on the left side. This is most clearly seen in Figure 3-13.

However, emphasis should be placed on "net" since there were counterclockwise flows in the

cephalic region, which is shown in Figure 3-14, and more subtly in Figure 3-13. Also, because

of the SAS compression, flow in the ventral region became choked, prohibiting circulation

around the cord in this area, and thereby reducing the morphine concentration in this region.

Another remarkable feature about the morphine distribution are the shifts in

concentration over time. In Figure 3-12 there is a time delay between when each region's

concentration began to reach their more "steady" orders of magnitude. The ventral aspect

consistently took the longest to catch up with the rest. However, in the cephalic and center cross

sections, its concentration surpasses the left side at the onset of the second cardiac cycle, and

falls below it after systole of the second cycle. This could be an artifact of the concentrations

having not yet settled, and increasing the solution time will help clear this up.

Another point worth mentioning is how the morphine does not disperse symmetrically

about the source, along the longitudinal axis. On average, the caudal cross section had higher

concentrations than the cephalic or center sections, which were very close in their averages.

Table 3-1 best summarizes this. The asymmetric distribution is attributed to being a product of

the net volume flow rate being in the caudal direction.

Diffusive Properties

Morphine has a very low coefficient of diffusion (D = 4.0e-10 m2/s) relative to the size of our

model's geometry. To put this in perspective, morphine diffuses 0.04% of a square millimeter in

one second. In order to assess the contribution that diffusion has on the dispersion of Morphine,

the Schmidt and Peclet numbers were employed. These values can be found in Table 3-2. For









all of the locations, the Peclet number averages on the order of 105, which means the inertial

effects of the fluid are 105 greater than the diffusive effects. Thus, the distribution of morphine

is governed almost entirely by the motion of the fluid.

When the CSF goes through diastole, there are times when the Peclet number drops to 0(0),

since it is a function of fluid velocity. When this happens, diffusion can contribute to the process

of dispersion; however, the Schmidt number reveals a different story. The Schmidt number is

not dependent on velocity, and is constant at Sc = 2450 throughout the course of the solution.

This indicates that viscous effects are 2450 times greater than diffusive effects. So, despite the

lack of CSF advection during diastole, the viscous forces inherent to the fluid are still

overpowering the diffusive properties of morphine.










Table 4-1. Temporal parameters relating CSF and cardiac systoles, measured using phase-
contrast MRI. (Source: Henry-Feugeas et al. [13])

Key temporal parameters of the CSF systoles
Locations Onset Peak End Duration
Thoracolumbar junction
Dorsal SAS 22.9 45 363 3.9 53.9 10.2 311 t B.1
Lateral SAS 24.1 + 5.1 37.0 5.0 603 10.3 36.1 7.4
Ventral SAS 20.9 7.4 37.6 6.3 61.3 15.6 40.2 13.1
Values are mean times + SD in percentage of the cardiac cycle. They are obtained from measurements inthe intracranial cisterns and spinal subarachnoid
spaces (SAS).


&L


VHP images showing compression of the SAS at the T10-T11 vertebral disc. Also
note the presence of nerve roots, which obstruct the flow. LEFT =model inlet,
RIGHT = 12 mm inferior. (Source: Visible Human Project [5])


Figure 4-1.









CHAPTER 5
CONCLUSIONS

This study was intended to be a preliminary evaluation of using the finite element method for

computing the distribution of morphine in an anatomically correct geometry. Clinically, this

type of work is not practical, since cadaveric data had to be used in order to generate the model.

However, it does provide an understanding about the role that fine structures and possibly

pathologies in the case of the SAS compression, play in the distribution of morphine. If imaging

modalities such as MRI progress to a high enough image resolution, it is possible that this

technique could be clinically applicable in the future.

Overall this work confirms the notion that CSF flow around the spinal cord is highly variable,

and can actually promote mixing under some anatomical conditions. It was found that morphine

dispersion was inertially and viscously dominated. It distributed in an irregular pattern, but what

appeared appropriate considering the SAS geometry. Velocities in different regions varied

greatly due to flow obstructions, promoting flow circulation and phase shifts in the flow.

Reynolds and Womersley numbers were on target with those found in literature.

Future work will include simulations of longer duration and longer lengths of the SAS.

Different sections of the spine will also be assayed. Models absent of nerve roots and other SAS

anatomy will be constructed to act as controls. The input velocity waveform will be re-evaluated

so that the volumetric flow rate agrees better with the published data. Baclofen and other species

such as antibiotics may also be investigated as the mass transport emphasis.









LIST OF REFERENCES


[1] Murphy, P. M., and Cousins, M. J., "Neural Blockade and Neuromodulation in Persistent
Pain Management," The Paths of Pain: 1975-2005, H. Merskey, and et al., eds., IASP
Press, Seattle, pp. 447-468.

[2] Simpson, K., 2008, "Interventional techniques for pain management in palliative care,"
Medicine, 36(2), pp. 72-74.

[3] Follett, K. A., 2003, "Intrathecal analgesia and catheter-tip inflammatory masses,"
Anesthesiology, 99(1), pp. 5-6.

[4] Myers, M. R., 1996, "A numerical investigation into factors affecting anesthetic
distribution during spinal anesthesia," J Biomech, 29(2), pp. 139-149.

[5] National Library of Medicine, The Visible Human Project, 1997,
http://www.nlm.nih.gov/research/visible/visible_human.html

[6] Loth, F., Yardimci, M. A., and Alperin, N., 2001, "Hydrodynamic modeling of
cerebrospinal fluid motion within the spinal cavity," J Biomech Eng, 123(1), pp. 71-79.

[7] Gupta, S., Poulikakos, D., and Kurtcuoglu, V., 2008, "Analytical solution for pulsatile
viscous flow in a straight elliptic annulus and application to the motion of the cerebrospinal
fluid," Phys Fluids, 20(9), pp. -.

[8] Stockman, H. W., 2006, "Effect of anatomical fine structure on the flow of cerebrospinal
fluid in the spinal subarachnoid space," Journal of Biomechanical Engineering, 128(1), pp.
106-114.

[9] Stockman, H. W., 2007, "Effect of anatomical fine structure on the dispersion of solutes in
the spinal subarachnoid space," Journal of Biomechanical Engineering, 129(5), pp. 666-
675.

[10] Schiffer, E., Van Gessel, E., and Gamulin, Z., 1999, "Influence of sex on cerebrospinal
fluid density in adults," Brit J Anaesth, 83(6), pp. 943-944.

[11] Clark, S. L., Edeson, R. O., and Ryall, R. W., 1983, "The relative significance of spinal and
supraspinal actions in the antinociceptive effect of morphine in the dorsal horn: an
evaluation of the microinjection technique," British Journal of Pharmacology, 79(3), pp.
807-818.

[12] Rathmell, J. P., Lair, T. R., and Nauman, B., 2005, "The role of intrathecal drugs in the
treatment of acute pain," Anesthesia and Analgesia, 101(5), pp. S30-S43.

[13] Henry-Feugeas, M. C., Idy-Peretti, I., Baledent, O., 2000, "Origin of subarachnoid
cerebrospinal fluid pulsations: a phase-contrast MR analysis," Magn Reson Imaging, 18(4),
pp. 387-395.









BIOGRAPHICAL SKETCH

Adam Daniel Fletcher was born on April 15, 1984 in Gainesville, Florida. He was raised in

Gainesville and graduated in 2002 from the International Baccalaureate program at Eastside

High School. He attended the University of Florida that fall and graduated in spring 2008 with

Bachelor of Science degrees in aerospace engineering and mechanical engineering, as well as a

minor in biomechanics. He continued his education at the University of Florida that fall of 2008

and completed his Master of Science in biomedical engineering in the summer of 2010. During

his study at the University of Florida, Adam worked for UF Computing and Network Services as

a technical consultant, as well as the Department of Astronomy as a mechanical engineer.





PAGE 1

1 A COMPUTATIONAL STUDY OF FLOW AND MASS TRANSPORT IN INTRATHECAL DRUG DELIVERY By ADAM D. FLETCHER A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 2010

PAGE 2

2 2010 Adam D. Fletcher

PAGE 3

3 I dedicate this work to all of my family and friends

PAGE 4

4 ACKNOWLEDGEMENTS I would like to thank Dr. Malisa Sarntinoranont, Dr. Roger Tran SonTay, and Dr. Paul Carney for serving as my supervisory committee They have provided ample support, c ounsel and encouragement throughout my studies Without them, I could not have accomplished this work.

PAGE 5

5 TABLE OF CONTENTS ACKNOWLEDGEMENTS .............................................................................................................4 LIST OF TABLES ...........................................................................................................................7 LIST OF FIGURES .........................................................................................................................8 LIST OF ABBREVIATI ONS ........................................................................................................10 CHAPTER 1 INTRODUCTION ..................................................................................................................13 Background .............................................................................................................................13 Literature Review ...................................................................................................................14 2 METHODS .............................................................................................................................19 Model Generation ...................................................................................................................19 Image Preparation ............................................................................................................19 Stereolithography (STL) Surface Generation ..................................................................20 Non Uniform Rational B Spline (NURBS) Surface Generation ....................................20 Units and Coordinate System ..........................................................................................21 Numerical Implementation .....................................................................................................21 Fluid Model .....................................................................................................................21 Mass Model .....................................................................................................................23 Mesh Generation .............................................................................................................24 3 RESULTS ...............................................................................................................................29 Definitions ..............................................................................................................................29 Fluid Characteristics ...............................................................................................................29 Temporal V elocity Profiles .............................................................................................29 Spatial Velocity Profiles ..................................................................................................30 Dimensionless Parameters ...............................................................................................30 Volumetr ic Flow Rate .....................................................................................................31 Mass Characteristics ...............................................................................................................32 Concentrations .................................................................................................................32 Dimensionless Para meters ...............................................................................................32 Dispersion ........................................................................................................................33 4 DISCUSSION .........................................................................................................................45 Anatomy .................................................................................................................................45 Subarachnoid Space (SAS) Compression .......................................................................45 Nerve Roots .....................................................................................................................45

PAGE 6

6 Flow Characteristics ...............................................................................................................46 Temporal Phenomena ......................................................................................................46 Flow Development ..........................................................................................................46 Model Validation .............................................................................................................47 Mass Characteristics ...............................................................................................................48 Dispersion ........................................................................................................................49 Diffusive Properties .........................................................................................................49 5 CONCLUSIONS ....................................................................................................................52 LIST OF REFERENCES ...............................................................................................................53 BIOGRAPHICAL S KETCH .........................................................................................................54

PAGE 7

7 LIST OF TABLES Table page 31 Values for relative morphine concentrations and their average, at each region sampled ( t =1.647 s). ...........................................................................................................34 32 Maximum and average Pclet numbers over the cardiac cycle for each sampling location. The Schmidt number is included as well, and is a constant based on ma terial properties .............................................................................................................34 41 Temporal parameters relating CSF and cardiac systoles, measured using phasecontrast MRI. (Source: Henry Feugeas et al. [13]) ...........................................................51

PAGE 8

8 LIST OF FIGURES Figure page 11 Pump and catheter implantation for intrathecal drug delivery (L: posterior view; R: saggital view. Sources: Medtronic and Mayfield Clinic) ................................................16 12 Saggital close up of catheter tip and surrounding anatomy. (Source: Mayfield Clinic) ...16 13 Spatia l velocity profiles for the circular annulus at different times during the cardiac [6]) ......................................................................................................................................17 14 Comparison of analytical velocity profiles along the minor axis of the annular cross section with CFD over a complete oscillatory cycle with period T. [7] ............................17 15 Stockmans 4 models used to show the difference in CSF flow by varying fine structures. Only A and D were used in his dispersion study. (Adapted from [8]). ...........18 16 Dispersion differences after the time shown in the right hand column. Models D and A refer to those in Fig. 15. (Adapted from [9]). ...............................................................18 21 The region modeled for this study is represented as the section between the two horizontal lines that intersect the T10 and T11 vertebrae. .................................................26 22 Patches and grids used to construct a NURBS surface. Grid density lowered for better visualization .............................................................................................................27 23 Cerebrospinal fluid velocity over one cardiac cycle at the T10 vertebral level .................27 24 Spread of analgesia after intrathecal administration in the lumbar cistern. Values shown in parentheses are octanol:water partition coefficients; a larger value indicates that a drug is more lipid soluble [12]. ................................................................................28 25 Finite element mesh of SAS model, using 4node tetrahedral FCBI elements. The velocity and mass loadings are shown in pink. ..................................................................28 31 Sampling locations for velocity and mass values denoted by colored elements ...............35 32 Maximum velocity nodes depicted as triangles in these band plots. .................................35 33 Y Velocity plots at the locations shown in Figure 31. .....................................................36 34 Spatial velocity profile sampling locations. The caudal left reference does not lie on an edge as depicted. A body is absent from this view, and the true shape of the section can be seen in Figure 3 1. ......................................................................................37

PAGE 9

9 35 Spatial velocity profiles in each region. All series are representative of the y velocity at one tenth of a cardiac cycle (T/10, 2T/10, etc.) ...............................................38 37 Flow entrance length at the dorsal aspect of the model inlet and exit. ..............................39 38 Reynolds plots at each sampling location ..........................................................................40 39 Womersley number at the cephalic, center, and caud al cross sections ..............................41 310 Error between the models inlet and exit volume flow rates .............................................41 311 Comparison of the models volume flow rate with Loths MRI measured flow rate ........42 312 Relative concentrations of morphine at each of the 12 sampling regions. Values are over the course of 1.4 cardiac cycles. ................................................................................43 313 Relative concentration of morphine at time 1.647 seconds. This was the final mass transfer time step, and represents 1.4 cardiac cycles. ........................................................44 314 Cephalic (left) and caudal (right) cross sections showing the velocity field at peak systole. Note the counterclockwise an d clockwise trends. ...............................................44 41 VHP images showing compression of the SAS at the T10T11 vertebral disc. Also note the presence of nerv e roots, which obstruct the flow. LEFT =model inlet, RIGHT = 12 mm inferior. (Source: Visible Human Project [5]) ......................................51

PAGE 10

10 LIST OF ABBREVIATION S Acs cross sectional area C cervical level CSF cerebrospinal fluid cm centimeters D characteristic length scale (annular gap) Dh hydraulic diameter Di diffusivity tensor of ith species f body force on fluid f temporal frequency FCBI flow c onditionbased interpolation ITDD intrathecal drug delivery L lumbar level le entrance length m meters mi mass of ith species mi mass creation rate of ith species MRI magnetic resonance imaging mm millimeters ms millisecond s n number of species NURBS nonuniform rational bspline p fluid pressure

PAGE 11

11 Pe Peclet number Pwet wetted perimeter Re Reynolds number Rh hydraulic radius s seconds SAS subarachnoid space Sc Schmidt number STL stereo lithography (format) t time T thoracic level u axial component of fluid velocity v fluid velocity vi velocity of ith species VHP visible human project Womersely number phase angle fluid viscosity fluid density i mass ratio of ith species angular frequency

PAGE 12

12 Abstract of Thesis Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Master of Science A COMPUTATIONAL STUDY OF FLOW AND MASS TRANSPORT IN INTRATHECAL DRUG DELIVERY By Adam D. Fletcher August 2010 Chair: Roger Tran SonTay Major: Biomedical Engineering The detailed geometry of the spinal cord and its surrounding tissues make flow of the cerebrospinal fluid highly variable. Intrathecal drug delivery requires that drugs be injected into the subarachnoid space, which surrounds the spinal cord. In the present work we simulated the flow of cerebrospinal fluid and the mass transfer of morphine in a 29 mm section of the spine, by employing computational fluid dynamics. The Visible Human Project dataset is used to generate a mesh of the subarachnoid space and MRI data is used to apply a physiological pulsatile velocity condition at the inlet of the model. Morphine inject ion is modeled as a saturated point source with a concentration of 1.0. Results from the study show velocities peaking at 20 cm/s with phase shifts up to 5 degrees. The average Womersley number was calculated to be 3.5 with a deviation of about 1. Spati al velocity profiles show varying flow development in different regions, with a tendency to be blunted during systole. Reynolds numbers were calculated between 200 and 500 during peak systole, and agree well with the literature. After 1.65 seconds, the a verage relative morphine concentration was 4.4E 4. Dispersion was highly complex, involving circulation about the spinal cord, redirection around nerve roots, and a compression of the subarachnoid space. The Peclet and Schmidt numbers reveal that diffusi on was dominated by inertial and viscous effects.

PAGE 13

13 CHAPTER 1 INTRODUCTION Background Implantable intrathecal drug delivery (ITDD) systems are becoming increasingly common [1] They work by administering drugs directly to the intrathecal space surrounding the spinal cord. This eliminat es some side effects that are present with intravenous and oral administration. In addition, smaller doses can be used which is cheaper and healthier for the patient [2] Using implantable drug pumps also allows doctors to precisely program the size and timing of doses which increases the quality of life for the patient. Two common types of drugs that are delivered intrathecally are muscle relaxants and analgesics. Morphine is a common analgesic that is used in treatment of a variety of pain conditions. Sometimes when morphine is injected in a nonhomogeneous fashion, the result can be a high dose t o a localized part of the spinal cord, which may result in injury [3] Baclofen is the commonly prescribed muscle relaxant for ITDD. It is used to trea t muscle spasticity in conditions such as cerebral palsy. We investigated morphine in this study because it is prescribed more often and has a larger incident of side effects. A typical procedure for implanting the pump and catheter involves insertion of an introducer needle between the third and four lumbar levels (L3 L4) and delicately guiding the catheter through the needle and up to the eleventh vertebral level (T11). Figure 11 shows a posterior and saggital view of the implant. Once inserted, the c atheter is anchored to the thoracolumbar fascia and any excess catheter is trimmed away. A pocket for the drug pump is created above the ilium. The catheter is connected to the pump, and the pump is then sutured into the pocket. An additional illustrati on of the catheter tip and surrounding anatomy is provided in Figure 1 2.

PAGE 14

14 Literature Review One of the first numerical investigations into the distribution of anesthesia in the intrathecal space was performed by Myers [4] His technique involved using the finite element (FE) method to simulate the dispersion of anesthetics during spinal anesthesi a. Spinal anesthesia is performed between the third and fourth lumbar vertebral levels, and has the potential to result in cauda equina syndrome. This occurs when high concentrations of the analgesic invade the cauda equina. His results showed that the size of the needle used for injection was a critical parameter in dispersion of the anesthesia. Larger needles resulted in higher dispersion. Francis Loth has been a strong contributor to the study of ITDD. For one of his works, his group used the visi ble human project (VHP) [5] magnetic resonance imaging (MRI) datasets to measure the geometry of the intrathecal space [6] From these measurements they made concentric and eccentric elliptical cross sections as well as concentric circular cross sections for their numerical computations. To define the CSF flow waveform, dynamic phase contrast MRI was performed on a healthy adult at the second cervical vertebral level (C2). This provided the CSF volumetric flow rate over the cardiac cycle. They found the flow in their models to be laminar and inertially dominated with Reynolds numbers on the order of 102. Simulations of larger annular gaps between the cord and dura mater resulted in velocity profiles that were blunted at the center and peaked at the walls. This type of flow is typical of high Womersley nd probably due to flow reversal before the CSF can accelerate to a full developed state. Smaller annular gaps resulted in lower Womersley numbers and a more typical profile. Figure 1 3 displays these phenomena. Loths results for the concentric ellipti cal annulus were verified by Gupta et al. [7] who found the analytical solution for pulsatile flow in an elliptical annulus. Figure 1 4 shows his comparison to Loths work.

PAGE 15

15 Previous studies by Stockman have shown that fine anatomic structures in the subarachnoid space (SAS), such as dor sal and ventral nerve roots, do not significantly alter the velocity profiles of CSF [8] However, they are responsible for a more rapid dispersion of injected solutes in an oscillatory environment [9] Figure 1 5 displays the geometry used in these experiments. Stockman used MRI data from Loth [6] to define his cross sectional geometry. For his flow simulations the length of the model was 1 cm. His dispersion simulations were done with a model length of 4 8 cm. Figure 1 6 provides the reader with a visual of the difference in dispersion between the models containing and absent of fine structures. It was concluded that fine structures enhance dispersion by 5 to 10 tim es over the range of Schmidt numbers (Sc) that were used. Our mission was to simulate an ITDD system using the high resolution VHP dataset. This allowed us to more accurately model the intrathecal space with fine structures. The high resolu tion dataset h as not been used in previous literature and should provide a more real istic response due to better accuracy of the model geometry. A pulsatile flow was used because this best emulates the oscillatory motion of CSF. This periodic motion occurs due to expa nsion of the cerebrovascular bed during cardiac systole. To achieve the proper physiological oscillation, Loths MRI generated waveform for CSF flow at the C2 level was used for our velocity input. Morphine injection was simulated by a mass concentration point source on the dorsal aspect of the intrathecal space. With this combination of techniques we attempt to take a new step in ITDD research.

PAGE 16

16 Figure 1 1. Pump and catheter implantation for intrathecal drug delivery (L: posterior view; R : saggital view. Sources: Medtronic and Mayfield Clinic) Figure 1 2. Saggital close up of catheter tip and surrounding anatomy. (Source: Mayfield Clinic)

PAGE 17

17 Figure 1 3. Spatial velocity profiles for the circular annulus at different times during the cardiac external and internal radius are 10 and 9 mm, respectively (Adapted from [6] ) Figure 1 4. Com parison of analytical velocity profiles along the minor axis of the annular cross section with CFD over a complete oscillatory cycle with period T. [7]

PAGE 18

18 Figure 1 5. Stockmans 4 models used to show the difference in CSF flow by varying fine structures. Only A and D were used in his dispersion study. (Adapted from [8 ] ). Figure 1 6. Dispersion differences after the time shown in the right hand column. Models D and A refer to those in Fig. 15. (Adapted from [9] ).

PAGE 19

19 CHAPTER 2 METHODS Model Generation The first phase of this study was to construct a surface for the finite element model. The surface was based upon high r esolution male cadaveric images from the VHP. The region spans from the bottom one fourth of T10 to half of T11, as seen in Figure 21. Image Preparation The high resolution VHP dataset consists of images taken from a cryogenically frozen male cadaver. All of the images are axial cross sections that are spaced 1 mm apart. Each slice was photographed with 70 mm film and subsequently scanned into digital RAW images. The pixel size for each image is 0.144 x 0.144 x 1.000 mm and the color quality is 24 bi t. 30 images were used to create a model 29 mm long. In order to convert the images from RAW to bitmap, a batch process was performed using Photoshop (Adobe Systems Inc ., San Jose, CA). The images contained fiducial markers from a rod situated in the f rozen gel that surrounded the cadaver. However, the rod had a bow to it and was not perfectly true rendering the fiducial markers invalid. A custom program was written to register the images by sampling the edge of the frozen gel and performing a least squares regression. It then aligned the images rotationally and translated them so the edge of the gel was flush with the left side of the image. Minor translational corrections were made manually to finalize the registration. These were do ne using the cord as a fiducial; therefore our model does not have any curvature to it. However in this section of the spine, the cord is linear and it was determined that the centroid of the cord only moved by 14 pixels over the length of the model, which is about 4.6 mm, and assumed to be negligible.

PAGE 20

20 After registering the dataset, a small MATLAB (The MathWorks Inc., Natick, MA) script was written to filter out the excess anatomy that did not define the fluid space. This was done by reading through the pixels and converting to 0 those that did not fall into predefined ranges of red, green, and blue. MATLABs imtool was used to assess the color values of anatomy that defined the fluid space. This provided the range of values for the filter to exclude. After the f ilters were run, manual corrections were performed on pixels that did not agree with the filtered values but were indeed fluid space. Stereolithography (STL) Surface Generation A three dimensional model was created from the 30 two dimensional images using Amira (Visage Imaging, Inc., San Diego, CA) software. Each image was segmented in order to define the fluid space to be used in the model. Much of this was already completed in the filtering process, but a few modifications were necessary because the software did not handle thicknesses of 2 or less pixels very well. Sharp corners were also removed with the segmentation editor so that it would be more practical to generate a NURBS surface after exporting the STL. Next, the surface was generated, cons isting of 98740 triangles. Unconstrained smoothing was used to soften the surface, and the ends of the annulus were closed off. The surface was exported in little endian STL format. Non Uniform Rational B Spline (NURBS) Surface Generation In order to apply boundary conditions to the model, we must define areas known as patches or faces to refer to. A NURBS surface is a collection of patches that consist of B splines. The patches are populated with grids that precisely represent the surface to be bui lt. Geomagic (Geomagic, Inc., Durham, NC) was used to construct the patches on our model. Please see Figure 2 1 for a visualization of what this looks like on our model. Prior to adding the patches the model had to be refined. This was done by increasi ng the number of triangles that

PAGE 21

21 make up the surface to 278212. Several smoothing techniques were then utilized to eliminate the striations caused by the 1mm gap between each image. The ends of the annulus were cut to create flat surfaces that would be easily identifiable as slip boundary conditions, and also to identify loading conditions. After these steps were completed, the model was patched and grids were placed to create the final NURBS surface. Units and Coordinate System Using Geomagic, the leng th units of the surface were established to be SI, and a Cartesian coordinate system was defined for the model. The origin was placed in the center of the cord at the mid point of its length. The spinal cord was aligned to the y axis with positive being in the caudal direction, and the positive z axis aligned to the ventral direction. Positive x axis is directed to the cadavers anatomical left. The surface was exported in IGES file format. Numerical Implementation The model was solved numerically using the FE package ADINA F (ADINA R&D Inc, Watertown, MA) The FE equations were discretized using the Galerkin method. Due to the transient nature of the problem, the Euler method was used to perform time integration. Solving the discretized FE equations was done using the NewtonRaphson method. Fluid Model To solve for the flow field within the reference and anatomic annuli, the Navier Stokes mass continuity and momentum equations were used. The mass continuity equation is shown below ( density, time, v=velocity) 0 v t (2 1) The fluid was approximated as water since, according to Stockman, viscosity data on CSF is clustered toward being more like water [9] and the density of water is within 0.4% of values

PAGE 22

22 measured by Schif fer et al [10] Being that water is a liquid, and liquids at this low of a speed (< 8 cm/s ) can be assumed incompressible the density is constant and 0 t This reduces the continuity equation to 0 v (2 2) Due to the low velocities encountered in this problem, the flow field was assumed to be laminar and the conservative form of the momentum equation was used ( p =pressure, u=axial velocity, f =body force) f v v v v 2u p t (2 3) To simplify the problem, the patient was assumed to be lying in the supine position so that body forces could be neglected. Due to the pulsatile nature of CSF flow, the viscous nature of the fluid, and varying geometry of the intrathecal space, all terms were kept except body force, resulting in v v v v2 u p t (2 4) At the walls of all the anatomic structures, including the spinal cord, nerve roots, and arachnoid mater, a no slip condition ( v=0) was assumed. The energy equation was not solved due to the assumption of incompressible flow without heat generation. To simulate motion of the CSF, a velocity loading condition was applied at the cephalic end of the model. The velocities used were adapted from the volume flow rates measured at t he C2 vertebral level by Loth et al. using phase contrast MRI [6] Two programs were written to satisfy this adaptation. The first was to digitize the CSF waveform that was measured via MRI. This allowed us to convert the image from his publication into discrete values of time and flow

PAGE 23

23 rate. The second program analyzed the VHP image that represented the cephalic end of the model. It summed all of the pixels that define d the fluid space and, with knowledge of the pixel size, computed the cross sectional area of the inlet. From here, the digitized waveform was divided by the inlet area to produce the velocity waveform used for loading. Figure 22 shows the converted CSF velocity waveform. The period (T) is equal to 1.176 seconds. The velocity was prescribed uniformly at the cephalic inlet of the model, in the axial (caudal) direction. The simulation was set to run for 5 cardiac cycles, totaling 250 time steps at inter Mass Model In order to solve for the mass transport of morphine, the mass transport equation had to be coupled with the flow equations. The mass transport equation is derived by combining Ficks Law and the species conservation equ ation The latter is shown below i i i im t v (2 5) Where i is the density of the ith species and mi is the mas s creation rate. Ficks Law represents the diffusion of a species in a fluid i i i i im t D v (2 6) Where Di i is the mass ratio for the ith species. For this study we assume the diffusivity to be equal in all directions. Clark et al. measured the diffusivity of morphine in CSF to be in the region of 3 to 5 x 106 cm2s1 [11] We took the average value and set D = 4 x 1010 m2s1 to keep with the SI standard. The mass ratio is a ratio of the density of the ith species to the total, or bulk density of all species combined. The definition is shown below.

PAGE 24

24 n i i i i 0 (2 7) Since we are only concerned with i = 1.0 by default. This value also represents the concentration of the species at the site where the condition is applied. In this case it is 1.0, or 100% concentrated. Finally the conservation equation and Ficks Law combine t o form the mass transfer equation i i i i im t D v (2 8) To best approximate the conditions of ITDD, the species was input on the dorsal side of the model, directly on the arachnoid mater. This was done because there were no ge ometric references in the center of the annulus to use. In reality the injected morphine is in the form of a solution; however, this was neglected to simplify the problem and dismiss any fluidfluid interactions. The spinal cord, nerve roots, and arachnoi d mater were assumed to be inert (nonporous) to morphine because according to Rathmell et al. [12] morphine has an octanol water partition coefficient of one, meaning that it is not attracted to lipids. Figure 2 3 shows the comparison of morphine distribution with other common intrathecal drugs. Note that morphine diffuses throughout the whole spine whil e the others do not, implying that it does not readily react with the tissue. Due to restricted computation time, the mass model was only solved for 1.4 cardiac cycles, but with the same step size as the f luid model. Mesh Generation The model was sectioned into 17 bodies in order to align nodes along model lines in regions of interest. This was done by placing an axial section cut through the center of the model, two additional axial cuts at 10 mm fore an d aft of center, a saggital cut along the center of the model, and a coronal cut 1.5mm from center in the ventral direction. Each body was subdivided into 4-

PAGE 25

25 node tetrahedral Flow ConditionBased Interpolation (FCBI) elements, with a characteristic length 0.5 mm. This made for a total of 129542 elements. A Delaunay meshing algorithm was used on all elements and boundaries. In one case Advancing Front had to be used for boundary meshing. The total number of nodes generated was 27316. Figure 2 4 shows the meshed model including velocity and species loading.

PAGE 26

26 Figure 2 1. The region modeled for this study is represented as the section between the two horizontal lines that intersect the T10 and T11 vertebrae.

PAGE 27

27 Figure 2 2. Patches and grids used to construct a NURBS surface. Grid density lowered for better visualization Figure 2 3. Cerebrospinal fluid velocity over one cardiac cycle at the T10 vertebral level

PAGE 28

28 Figure 2 4. Spread of analgesia after intrathecal administrati on in the lumbar cistern. Values shown in parentheses are octanol:water partition coefficients; a larger value indicates that a drug is more lipid soluble [12] Figure 2 5. Finite element mesh of SAS model, using 4node tetrahedral FCBI elements. The velocity and mass loadings are shown in pink.

PAGE 29

29 CHAPTER 3 RESULTS Definitions Much of the flow and mass transfer data was collected from three axial cross sections of the model and used as reference planes for extracting data. These coincide with the axial cuts described in the section Mesh Generation and the center cut coincides with the morphine source at the dorsal wall. For the remainder of the paper, the data related to the center cut will be prefixed with center, and the two cuts fore and aft of center will be prefixed with cephalic and caudal, respectively. Figure 3 1 shows these three cuts and the locations were the velocity was sampled. The element colors correspond to the colors used in the figures that follow. The terms left and right refer to the cadavers anatomical left and right sides. All data series with a max suffix refer to the node within the cross section that had the maximum velocity during peak systole. Their locations are shown in Figure 32 and are represented as triangles on the velocity band plot. This is a plot of the axial velocity at peak systole. By comparing Figure 3 1 and 32, the reader can also identify where the sampling locations are in respect to the velocity field. Fluid Characteristics The results of the fluid portion of the simulation are presented here. The solution requir ed 98 hours to process 250 time steps, equaling 5 cardiac cycles. Temporal and spatial velocity profiles from 12 sites along the length of the model were constructed over one cardiac cycle. Dimensionless numbers are given, as well as model validation by the volumetric flow rate. Temporal Velocity Profiles Figure 3 3 presents graphs of the temporal velocity profiles. The velocity prescribed at the inlet is included as a reference. One remarkable feature of these plots is that the waveforms are

PAGE 30

30 out of phase with one another. Due to the low temporal resolution of this simulation (23.5 ms), cubic spline interpolation was used to determine the times of maximum axial velocity at each location. The phase shift ( ) for each location was computed relative to the prescribed velocity. t f 360 (3 1) Where f t is the time shift at peak velocity. The results are shown in Figure 34. Spatial Velocity Profiles Spatial velocity profiles were constructed from the locations highlighted in Figure 3 5. These planes are the same as those shown in Figure 3 1. For each location, ten profiles of the axial velocity were plotted over the cardiac cycle. Figure 3 6 displa ys these plots. Dorsal/ventral and left/right locations are overlaid to conserve space. The x axis of all plots begins at zero and increases in the models positive x or z direction, which are the cadavers left and ventral directions, respectively. T he entrance length ( le ) was computed as a supplement to help assess the flow development. The following equation for entrance length in pipe flow was used. D le Re 06 0 (3 2) The prescribed velocity was used to calculate the Reynolds number (Re) and the annular gap at the dorsal aspect, between the cord and arachnoid mater, was used as the characteristic length scale (D=2.0 2.1 mm). Figure 3 7 plots the entrance length against the input velocity. The prescribed velocity was integrated and plotted to provide the reader with an idea of how far the fluid moves, neglecting viscous effects. Dimensionless Parameters Reynolds number This number describes the ratio of inertial to viscous forces acting on a flow. It was computed to help va lidate the model against the literature.

PAGE 31

31 hVD Re (3 3) The hydraulic diameter (Dh) was used as the length scale for consistency with Loth et. al [6] and Stockman [8] The formula for the hydraulic diameter was borrowed from Loth. It is a function of the annulus cross sectional area (Acs) and the wetted perimeter of the cross section (Pwet). wet cs hP A D 4 (3 4) The wetted perimeter and cross sectional area were calculated using Solidworks 2008 (Dassault Systmes, Vlizy, France). Figure 3 8 gives the varying Reynolds number over one cardiac cycle. Womersley number. frequency and its viscosity. 2 1 hR (3 5) Rh is the hydraulic radius ( 2h h D R ), =density, =dynamic viscosity, and is the frequency of the flow oscillation in rad/s. For our simulation =5.34 rad/s. The Womersley number can indicate whether a pulsatile flow will reach a fully developed state before it reverses, and is a common component of biofluid analysis. The va lues for this study were computed at the cephalic, center, and caudal cross sections referred to previously. Figure 3 9 displays the three values. Volumetric Flow Rate As a means of validating the model, volumetric flow rates were extracted from the solu tion. The flow rates were computed 1 mm from the inlet and exit of the annulus. Figure 310 shows the percent error between the two values, plotted against the prescribed velocity. As another

PAGE 32

32 check, the simulations volumetric flow rate was compared to Loths MRI measured flow rate. The result is shown in Figure 311. Mass Characteristics The mass transport model gives the re ader an idea of how morphine will disperse in the intrathecal space of the cadaver, and possibly other patients with similar anat omy. All of the mass solutions reported in this text are relative concentrations. That is, they represent the Concentrations The transient nature of the morphine concentrations were plotted over the solution time of 1.647 seconds and can be found in Figure 3 12. For reference, the first systole occurs at 0.35 seconds and the second systole occurs at 1.53 seconds. Table 31 summarizes the mass distribution at the 70t h (final) time step for each region. The dorsal and right regions are more steady, while the ventral and left regions vary more. Dimensionless Parameters Schmidt number To assist in the analysis of the mass transport results, two dimensionless paramet ers were calculated. The first is the Schmidt number D Sc (3 6) Which describes the ratio of viscosity to diffusivity in a flow. In our case, the diffusive tensor is a constant for all directions, and CSF was modeled value itself, Sc = 2450. Pclet number. The second parameter calculated was the Peclet number

PAGE 33

33 DhVD Sc Pe Re (3 7) This value represents the ratio of inertial to diffusive forces in a flow. The maximum and average Peclet numbers over one cardiac cycle were tabulated and can be found in Table 32. Graphs were not constructed because in our cas e, the Peclet number is just the Reynolds number multiplied by a constant, and temporally it would look the same as the Reynolds and velocity plots. Dispersion As an additional tool in understanding how the morphine dispersed, a band plot was produced for each cross section and can be found in Figure 313. This plot is representative of the last time step in the series ( t =1.647 s). The morphine disperses in a complex manner due to the anatomy in the SAS. Noting that the dorsal and right regions consistently have higher concentrations than the ventral and left regions, a vector plot was constructed to help identify why the morphine dispersed in this way. Figure 3 14 presents the velocity fields at the cephalic and caudal cross sections during systole.

PAGE 34

34 Table 3 1. Values for relative morphine concentrations and their average, at each region sampled ( t =1.647 s). cross section \ region dorsal right ventral left average cephalic 6.4E 04 2.3E 07 2.6E 09 3.1E 08 1.6E 04 center 7.7E 04 2.3E 07 9.2E 10 1.8E 08 1.9E 04 caudal 1.8E 03 3.8E 06 7.5E 08 3.6E 10 4.4E 04 Table 3 2. Maximum and average Pclet numbers over the cardiac cycle for each sampling location. The Schmidt number is included as well, and is a constant based on material properties Region Sampled Peclet Number dorsal ventral right left max Cross Section Location cephalic 9.2E+05 6.1E+05 8.8E+05 8.6E+05 9.4E+05 max 3.1E+05 1.8E+05 3.2E+05 2.9E+05 3.0E+05 avg center 9.0E+05 1.0E+06 1.0E+06 8.4E+05 1.2E+06 max 3.0E+05 3.2E+05 3.2E+05 2.6E+05 3.9E+05 avg caudal 1.1E+06 4.9E+05 1.0E+06 5.3E+05 1.1E+06 max 4.0E+05 1.3E+05 3.4E+05 1.8E+05 4.1E+05 avg Schmidt Number 2450 constant

PAGE 35

35 Figure 3 1. Sampling locations for velocity and mass values denoted by colored elements Figure 3 2. Maximum velocity nodes depicted as triangles in these band plots.

PAGE 36

36 Figure 3 3. Y Velocity plots at the locations shown in Figure 3 1.

PAGE 37

37 Figure 3 4. Velocity phase shifts relative to the prescribed velocity at the inlet Figure 3 5. Spatial velocity profile sampling locations. The caudal left reference does not lie on an edge as depicted. A body is absent from this view, and the true shape of the section can be seen in Figure 3 1.

PAGE 38

38 Figure 3 6. Spatial velocity profiles in each region. All series are representative of the y velocity at one tenth of a cardiac cycle (T/10, 2T/10, etc.)

PAGE 39

39 Figure 3 7. Flow entrance length at the dorsal aspect of the model inlet and exit.

PAGE 40

40 Figure 3 8. Reynolds plots at each sampling location

PAGE 41

41 Figure 3 9. Womersley number at the cephalic, center, and caudal cross sections Figure 3 10. Error between the models inlet and exit volume flow rates 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 Womersley Number 3.3 2.8 4.4 Cephalic Center Caudal

PAGE 42

42 Figure 3 11. Comparison of the models volume flow rate with Loths MRI measured flow rate

PAGE 43

43 Figure 3 12. Relative concentrations of morphine at each of the 12 sampling regions. Values are over the course of 1.4 cardiac cycles.

PAGE 44

44 Figure 3 13. Relative concentration of morphine at time 1.647 seconds. This was the final mass transfer time step, and represents 1.4 cardiac cycles. Figure 3 14. Cephalic (left) and caudal (right) cross sections showing the velocity field at peak systole. Note the counterclockwise and clockwise trends.

PAGE 45

45 CHAPTER 4 DISCUSSION Anatomy To fully understand the a nalysis of this model, it is important to note the anatomical features of the cadaver. Those that are particularly important are the nerve roots found within the SAS, and equally as important, the compression of the SAS at the T10T11 vertebral disc. Subarachnoid Space (SAS) Compression Taking a look at Figure 3 1, the reader should note that the ventral side of the center cross section has an annular gap much smaller than the other sections. For better visualization of why the cross section narrows so much at the center, two of the original VHP images used to generate this model are displayed in Figure 4 1. The left hand image represents the model inlet, and the other is 12 mm from the inlet. In the first image there is pink trabecular bone ventral ( bottom) to the SAS, and in the second image the nucleus pulposus of the T10T11 disc can be seen on the ventral side. The presence of the disc at the location of comp ression may indicate a geometry typical of a slightly herniated disc. Nerve Roots Also apparent in Figures 3 1 and 41 is the large nerve root on the left side of the SAS. It entirely blocks off flow between the dorsal and left regions in this section. Therefore, the CSF that is trying to enter the left side chamber is redirected to faster moving, lower pressure fields. Also, the thin clearance at the ventral aspect further isolates the CSF in the left chamber. The fluid in the left chamber, therefore, is more stagnant and play s less of a role in mixing. The large nerve root and the cord compression will be the dominant factors in determ ining the morphine distribution in this segment of the cord.

PAGE 46

46 Flow Characteristics Temporal Phenomena One of the remarkable features of the temporal velocity profiles is their phase delay. Henry Feugeas et al. measured similar types of phase delay in her MRI phasecontrast study [13] Table 4 1 gives the temporal parameters of their work. These values rep resent the relationship between the patients CSF systole and cardiac systole at the thoracolumbar junction. Their measurements are not intended to agree with this simulation, since we do not know what the geometries consisted of. Also, one must consider that the standard deviations on Table 41 are somewhat large, which of itself implies that this is a highly variable characteristic of the flow. These phase shifts are likely the product of flow obstruction due to the cadaver anatomy. This hypothesis is based on the spatial and temporal velocity data that showed the dorsal and right regions of the SAS to have higher velocities most of the time. Also, the phase angles calculated in Figure 3 4 show the dorsal and right regions to have positive shifts, whi le the ventral and left aspects have negative shifts. Flow Development The spatial velocity profiles in Fig. 3 6 demonstrate the degree of flow development in each region that was sampled. When viewing these plots note that entrance effects take place at both the inlet and exit of the model, due to the oscillatory nature of the flow. Thus at the beginning of the cycle, the velocity is negative and the exit has entrance effects imposed upon it, but after systole begins the flow becomes positive (caudal) and the inlet is subjected to entrance effects. Due to the irregular varying geometry, the entrance length was not able to be determined. The entrance length is a function of both the size and shape of the cross section that the fluid passes through. Since the geometry varies continuously, so does the entrance length, making it

PAGE 47

47 difficult to estimate. For the sake of a reference, using the solution for pipe flow, the entrance length was computed and is shown i n Figure 37. Profiles with a width greater than 1.5 mm showed signs of being more developed during diastole. This is due the lower R eynolds numbers during this phase of the cardiac cycle effectively reducing the entrance length. During systole the profiles are more blunted, and the entrance length may be equal to or greater than the length of the model. However, the plots do agree favorably with those calculated by Loth in Fig. 13b. Similar velocities are reported and =2.8 4.4 (Fig. 39). In this study the annular gap was quite small compared to Stockman [8] and as such his simulations had a the profiles were under developed such as that in Fig. 1 3a. Model Validation Reynolds Number. The Reynolds numbers plotted in Fig. 38 were validated against those in Loths study. Along the length of the spine, the maximum Reynolds numbers that he reported were on the range of 175 450. This is consistent with the maximum values in the local regions of this model, which vary between 200 and 500. Volumetric Flow Rate. To check the simulation for continuity, the volumetric flow rate was calculated at the inlet and exit of the annulus and compared. Under the assumption of incompressible flow the flow rate entering the annulus should be equal to that exiting the annulus. Fig. 310 displays the error between these two values. The error spikes to 0.26% during systole when the flow accelerates, and again to 0.41% during diastole when the veloc ity waivers around zero. The former can be explained by a low temporal reso lution, relative to the abrupt increase in velocity. During systole, the axial velocity increases by 2 5 cm/s in one time step, depending on the location sampled. The 0.41% error may be explained by a combination of

PAGE 48

48 the velocity oscillating around zero, and the intrinsic nature of the regions being out of phase. In this period the prescribed flow reverses twice over four time steps, and discrepancies about the direction of flow between individual regions can clearly be seen in Fig 3 3. Indeed, when comparing regions, there are large differences in time for when the flow finally reverses during diastole. An additional check on the model was done by comparing the calculated volumet ric flow rate with that measured by Loth. This is shown in Fig. 3 11. The MRI measured flow rate is the net flow rate after all losses are considered. In this study those losses were not accounted for, and as such the prescribed velocity was calculated with inherent losses, only to undergo more loss after entering the annulus. This explains the difference in the flow rates. A final note about the volumetric flow rate is that over one period it has a nonzero mean. This implies that there is a net move ment of fluid, and in this case in the caudal direction. This is more easily visualized by the particle distance plotted in Fig. 3 7. By continuity there must be a sink that accommodates this excess fluid, or a secondary cephalicdirected flow that begi ns at the bottom of the cord. The MRI measured the mean velocity over an average slice thickness of 0.5 cm [8] It may be that slow, secondary flows, that are cephalically directed, do not produce enough radio frequency (RF) signal to be detected by MRI Mass Characteristics In this simulation, we were only able to solve for a time of 1.65 seconds due to a lack of computational power and time. As such, we cannot determine what the concentration of morphine wi ll be at a more useful times of 110 minutes after injection. However, we can get a glimpse at how anato mical structures effect morphine dispersion, and given more computational time, a follow up study will report on the concentrations and distribution over longer periods.

PAGE 49

49 Dispersion The concentra tion of morphine tended to disperse in a net clockwise fashion around the SAS, due primarily to the large nerve root on the left side. This is most clearly seen in Figure 3 13. However, emphasis should be placed on net since there were counterclockwise flows in the cephalic region, which is shown in Figure 3 14, and more subtly in Figure 313. Also, because of the SAS compression, flow in the ventral region became choked, prohibiting circulation around the cord in this area, and thereby reducing the mor phine concentration in this region. Another remarkable feature about the morphine distribution are the shifts in concentration over time. In Figure 3 12 there is a time delay between when each regions concentration began to reach their more steady ord ers of magnitude. The ventral aspect consistently took the longest to catch up with the rest. However, in the cephalic and center cross sections, its concentration surpasses the left side at the onset of the second cardiac cycle, and falls below it after systole of the second cycle. This could be an artifact of the concentrations having not yet settled, and increasing the solution time will help clear this up. Another point worth mentioning is how the morphine does not disperse symmetrically about the source, along the longitudinal axis. On average, the caudal cross section had higher concentrations than the cephalic or center sections, which were very close in their averages. Table 3 1 best summarizes this. The asymmetric distribution is attributed to be ing a product of the net volume flow rate being in the caudal direction. Diffusive Properties Morphine has a very low coefficient of diffusion ( D = 4.0e 10 m2/s) relative to the size of our models geometry. To put this in perspective, morphine diff uses 0.04% of a square millimeter in one second. In order to assess the contribution that diffusion has on the dispersion of Morphine, the Schmidt and Pclet numbers were employed. These values can be found in Table 32. For

PAGE 50

50 all of the locations, the Pe clet number averages on the order of 105, which means the inertial effects of the fluid are 105 greater than the diffusive effects. Thus, the distribution of morphine is governed almost entirely by the motion of the fluid. When the CSF goes through diastole, there are times when the Peclet number drops to O (0), since it is a function of fluid velocity. When this happens, diffusion can contribute to the process of dispersion; however, the Schmidt number reveal s a different s tory. The Schmidt number is not dependent on velocity, and is constant at S c = 2450 throughout the course of the solution. This indicates that viscous effects are 2450 times greater than diffusive effects. So, despite the lack of CSF advection during di astole, the viscous forces inherent to the fluid are still overpowering the diffusive properties of morphine.

PAGE 51

51 Table 4 1. Temporal parameters relating CSF and cardiac systoles, measured using phasecontrast MRI. (Source: Henry Feugeas et al. [13] ) Figure 4 1. VHP images showing compression of the SAS at the T10T11 vertebral disc. Also note the presence of nerve roots, which obstruc t the flow. LEFT =model inlet, RIGHT = 12 mm inferior. (Source: Visible Human Project [5] )

PAGE 52

52 CHAPTER 5 CONCLUSIONS This study was intended to be a preliminary evaluation of using the finite element method for computing the distribution of morphine in an anatomically correct geometry. Clinically, this type of work is not practical, since cadaveric data had to be used i n order to generate the model. However, it does provide an understanding about the role that fine structures and possibly pathologies in the case of the SAS compression, play in the distribution of morphine. If imaging modalities such as MRI progress to a high enough image resolution, it is possible that this technique could be clinically applicable in the future. Overall this work confirms the notion that CSF flow around the spinal cord is highly variable, and can actually promote mixing under some anatomical conditions. It was found that morphine dispersion was inertially and viscously dominated. It distributed in an irregular pattern, but what appeared appropriate considering the SAS geometry. Velocities in different regions varied greatly due to flow obstructions, promoting flow circulation and phase shifts in the flow. Reynolds and Womersley numbers were on target with those found in literature. Future work will include simulations of longer duration and longer lengths of the SAS. Different sections of the spine will also be assayed. Models absent of nerve roots and other SAS anatomy will be constructed to act as controls. The input velocity waveform will be re evaluated so that the volumetric flow rate agrees better with the published data. Baclofen and other species such as antibiotics may also be investigated as the mass transport emphasis.

PAGE 53

53 LIST OF REFERENCES [1] M urphy, P. M., and Cousins, M. J., "Neural Blockade and Neuromodulation in Persistent Pain Management," The Paths of Pain: 19752005, H. Merskey, and et al., eds., IASP Press, Seattle, pp. 447468. [2] Simpson, K., 2008, "Interventional techniques for pain management in palliative care," Medicine, 36(2), pp. 72 74. [3] F ollett, K. A., 2003, "Intrathecal analgesia and catheter tip inflammatory masses," Anesthesiology, 99(1), pp. 56. [4] Myers, M. R., 1996, "A numerical investigation into factors affecting anesthetic distribution during spinal anesthesia," J Biomech, 29(2), pp. 139149. [5] National Library of Medicine, The Visible Human Project 1997, http://www.nlm.nih.gov/research/visible/visible_human.html [6] Loth, F., Yardimci M. A., and Alperin, N., 2001, "Hydrodynamic modeling of cerebrospinal fluid motion within the spinal cavity," J Biomech Eng, 123(1), pp. 7179. [7] Gupta, S., Poulikakos, D., and Kurtcuoglu, V., 2008, "Analytical solution for pulsatile viscous flow in a straight elliptic annulus and application to the motion of the cerebrospinal fluid," Phys Fluids, 20(9), pp. [8] Stockman, H. W., 2006, "Effect of anatomical fine structure on the flow of cerebrospinal fluid in the spinal subarachnoid space," Journa l of Biomechanical Engineering, 128(1), pp. 106114. [9] Stockman, H. W., 2007, "Effect of anatomical fine structure on the dispersion of solutes in the spinal subarachnoid space," Journal of Biomechanical Engineering, 129(5), pp. 666675. [10] Schiffe r, E., Van Gessel, E., and Gamulin, Z., 1999, "Influence of sex on cerebrospinal fluid density in adults," Brit J Anaesth, 83(6), pp. 943944. [11] Clark, S. L., Edeson, R. O., and Ryall, R. W., 1983, "The relative significance of spinal and supraspinal actions in the antinociceptive effect of morphine in the dorsal horn: an evaluation of the microinjection technique," British Journal of Pharmacology, 79(3), pp. 807818. [12] Rathmell, J. P., Lair, T. R., and Nauman, B., 2005, "The role of intrathecal drugs in the treatment of acute pain," Anesthesia and Analgesia, 101(5), pp. S30S43. [13] Henry Feugeas, M. C., Idy Peretti, I., Baledent, O., 2000, "Origin of subarachnoid cerebrospinal fluid pulsations: a phase contrast MR analysis," Magn Reson Imaging 18(4), pp. 387395.

PAGE 54

54 BIOGRAPHICAL SKETCH Adam Daniel Fletcher was born on April 15, 1984 in Gainesville, Florida. He was raised in Gainesville and graduated in 2002 from the International Baccalaureate program at Eastside High School. He attended the University of Florida that fall and graduated in spring 2008 with Bachelor of Science degrees in aerospace engineering and mechanical engineering, as well as a minor in biomechanics. He continued his education at the University of Florida that fall of 2008 and completed his Master of Science in biomedical engineering in the summer of 2010. During his study at the University of Florida, Adam worked for UF Computing and Network Services as a technical consultant, as well as the Department of Astronomy as a mechanical engineer.