USFDC Home  USF Electronic Theses and Dissertations   RSS 
Material Information
Subjects
Notes
Record Information

Full Text 
xml version 1.0 encoding UTF8 standalone no
record xmlns http:www.loc.govMARC21slim xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.loc.govstandardsmarcxmlschemaMARC21slim.xsd leader nam Ka controlfield tag 001 001917304 003 fts 005 20071119155135.0 006 med 007 cr mnuuuuuu 008 071119s2007 flu sbm 000 0 eng d datafield ind1 8 ind2 024 subfield code a E14SFE0002035 035 (OCoLC)181592668 040 FHM c FHM 049 FHMM 090 QC21.2 (ONLINE) 1 100 Hersh, Lawrence T. 0 245 Mathematical techniques for the estimation of the diffusion coefficient and elimination constant of agents in subcutaneous tissue h [electronic resource] / by Lawrence T. Hersh. 260 [Tampa, Fla.] : b University of South Florida, 2007. 520 ABSTRACT: The purpose of this work was to develop methods to estimate the diffusion coefficient and elimination constant for dexamethasone in subcutaneous tissue. Solutions to the diffusion equation were found for different conditions relevant to implantation and injection. These solutions were then used as models for measured autoradiography data where the unknown model parameters were the diffusion coefficient and the elimination constant. The diffusion coefficient and elimination constant were then estimated by curve fitting the measured data to these models. Having these estimates would be of practical importance since inflammation surrounding implantable glucose sensors may be controlled through local release of dexamethasone at the site of implantation. Derivation of the appropriate model, how the model was used to estimate D and k, and various specific profile examples were investigated in detail.Osmotic pumps containing [3H] dexamethasone were implanted into the subcutaneous tissue of rats. Digital autoradiography was used to measure the distribution of the [3H]dexamethasone within the subcutaneous tissue at 6, 24, and 60 hours after implantation. Measured concentration profiles, near the catheter tip through which the agent was released, were compared to solutions of the diffusion equation in order to characterize drug diffusion coefficients and elimination constants. There was good agreement between the experimental data and the mathematical model used for estimation. The diffusion coefficient for dexamethasone in subcutaneous tissue was found to be D = 4.11+1.77x10E10 m2/s, and the elimination rate constant was found to be k = 3.65+2.24x10E5/s. Additionally, [3H]dexamethasone was injected into the subcutaneous tissue of rats.Digital autoradiography was again used to measure the distribution of the [3H] dexamethasone within the subcutaneous tissue at 2.5 and 20 minutes after injection. Measured concentration profiles were again compared to a mathematical model of drug diffusion for injection. There was good agreement between the experimental data and the mathematical model. The diffusion coefficient found using this simple injection method was 4.01+2.01x10E10 m2/s. The simple method given here for the determination of the diffusion coefficient is general enough to be applied to other substances and tissues as well. 502 Thesis (M.S.)University of South Florida, 2007. 504 Includes bibliographical references. 516 Text (Electronic thesis) in PDF format. 538 System requirements: World Wide Web browser and PDF reader. Mode of access: World Wide Web. 500 Title from PDF of title page. Document formatted into pages; contains 73 pages. 590 Advisor: Yvonne Moussy, Ph.D. 653 Mathematical models. LevenbergMarquardt algorithm. Osmotic pump. Injection. Radiolabeling. 690 Dissertations, Academic z USF x Physics Masters. 773 t USF Electronic Theses and Dissertations. 4 856 u http://digital.lib.usf.edu/?e14.2035 PAGE 1 Mathematical Techniques for the Estimation of the Diffusion Coefficient and Elimination Constant of Agents in Subcutaneous Tissue by Lawrence T. Hersh A thesis submitted in partial fulfillment of the requirements for the degree of Master of Science Department of Physics College of Arts and Sciences University of South Florida CoMajor Professor: Yvonne Moussy, Ph.D. CoMajor Professor: Wei Chen, Ph.D. Francis Moussy, Ph.D. Date of Approval: February 27, 2007 Keywords: Mathematical Models, Levenbe rgMarquardt Algorithm, Osmotic Pump, Injection, Radiolabeling Copyright 2007, Lawrence T. Hersh PAGE 2 Acknowledgments It is with much appreciation and respect that I thank Dr. Yvonne Moussy for her guidance and insight in doing th is research. She was an excellent advisor and leader. Being a part of her research group for the last two years has been an enjoyable and enlightening experience. Secondly, I would like to thank Paul Dungel. He was crucial in obtaining and initially processing the autora diographic data. Without his diligence and skill, this work would have been significantl y more difficult to do. Thirdly, I would like to thank Dr. Wei Chen whose cooperation and encouragement are much valued. It was his openness and help that made the comple tion of my endeavors in the USF Physics Department possible. Finally, I would like to de dicate this thesis to my past teachers and particularly my grandmother Roselle Karrer Hersh, whose love of science and reasoning has always been a great motivating factor. This work was supported by NIH Grant R01 EB0164003 and by the General Electric Em ployee Tuition Assistance Program. PAGE 3 i Table of Contents List of Tables ii List of Figures iii Abstract iv Introduction 1 Experimental Methods 4 Materials 4 Subcutaneous Implantation 4 Preparation of Osmotic Pumps 5 Subcutaneous Injection 5 The 20 Minute Experiment 5 The 2.5 Minute Experiment 6 Autoradiographic Imaging and Analysis 6 Mathematical Formulation 10 Mathematical Formulation, Implantation 11 Mathematical Formulation, Injection 15 SteadyState Solution With Elimination, Implantation 21 General Solution, Implantation 25 General Solution, Injection 28 Marquardt Curve Fitting 33 Results 35 Discussion 50 References 58 Appendices 61 Appendix A: The LevenbergMarquardt Algorithm 62 Appendix B: Example MatLab Marquardt Program 67 PAGE 4 ii List of Tables Table 1. Resultant Estimated D and k Values Using Implantation 45 Table 2. Penetration Distance of Radioac tivity from the Tip of the Catheter 46 Table 3. Diffusion Coefficients for Dexamethasone in Various Media 47 Table 4. Elimination Constants of Various Agents in Subcutaneous Tissue 48 Table 5. Estimated D from Injection Experiments 49 PAGE 5 iii List of Figures Figure 1. Example Autoradiographic Image, Implantation 8 Figure 2. Example Autoradiographic Image, Injection 9 Figure 3. Theoretical Concentration vers us Distance Profiles for Implantation 14 Figure 4. Theoretical Concentration ve rsus Distance Profiles for Injection 19 Figure 5. Example Measured Injection Profiles 20 Figure 6. Implantation Profile with Cu rve Fit, 24 Hours, 45 Degrees.. 38 Figure 7. Implantation Profile with Curve Fit, 24 Hours, 0 Degrees. 39 Figure 8. Implantation Profile with Curve Fit, 60 Hours, 195 Degrees. 40 Figure 9. Implantation Profile with Curve Fit, 6 Hours, 60 Degrees. 41 Figure 10. Implantation Profile with Curve Fit, 6 Hours, 75 Degrees. 42 Figure 11. Implantation Profile with Curve Fit, 6 Hours, 30 Degrees. 43 Figure 12. Injection Prof ile with Curve Fit. 44 Figure 13. Flow Chart Showi ng the Marquardt Algorithm 66 PAGE 6 iv Mathematical Techniques for the Estimation of the Diffusion Coefficient and Elimination Constant of Agents in Subcutaneous Tissue Lawrence T. Hersh ABSTRACT The purpose of this work was to develop methods to estimate the diffusion coefficient and elimination constant for dexame thasone in subcutaneous tissue. Solutions to the diffusion equation were found for differe nt conditions relevant to implantation and injection. These solutions were then used as models for measured autoradiography data where the unknown model parameters were th e diffusion coefficient and the elimination constant. The diffusion coefficient and elim ination constant were then estimated by curve fitting the measured data to these m odels. Having these estimates would be of practical importance since inflammation surr ounding implantable glucose sensors may be controlled through local re lease of dexamethasone at the si te of implantation. Derivation of the appropriate model, how the model was used to estimate D and k, and various specific profile examples were investigated in detail. Osmotic pumps containing [3H]dexamethasone were implanted into the subcutaneous tissue of rats. Digital autora diography was used to measure the distribution of the [3H]dexamethasone within the subcutan eous tissue at 6, 24, and 60 hours after implantation. Measured concentration profile s, near the catheter tip through which the agent was released, were compared to solu tions of the diffusion equation in order to characterize drug diffusion coefficients a nd elimination constants. There was good PAGE 7 v agreement between the experimental data and the mathematical model used for estimation. The diffusion coefficient for de xamethasone in subcutaneous tissue was found to be D = 4.11 1.77 1010 m2/s, and the elimination rate constant was found to be k = 3.65 2.24 105 s1. Additionally, [3H]dexamethasone was injected in to the subcutaneous tissue of rats. Digital autoradiography was again used to measure the distribution of the [3H]dexamethasone within the subcutaneous tis sue at 2.5 and 20 minutes after injection. Measured concentration profiles were again compared to a mathematical model of drug diffusion for injection. There was good agreemen t between the experimental data and the mathematical model. The diffusion coeffici ent found using this simple injection method was 4.01 2.011010 m2/s. The simple method given here for the determination of the diffusion coefficient is general enough to be applied to other substa nces and tissues as well. PAGE 8 1 Introduction Several recent reports sugge st that controlled local re lease of dexamethasone may be useful for preventing inflammation around an implantable glucose sensor ( 1 3 ) This decrease in inflammation is expected to increase glucose sensor function and lifetime. Local drug delivery may be achieved usi ng biodegradable polymer implants ( 4 ), hydrogels ( 5 ), and osmotic pumps ( 6 ) Local delivery of dexamethasone would permit high interstitial drug concentrations at the s ite of glucose sensor implantation without producing high systemic drug levels. For su ccessful local treatment, dexamethasone must be released and penetrate through the tissue surrounding the implanted glucose sensor. Additionally, the concentration of dexamethasone in the subcutaneous tissue surrounding the implanted glucose sensor must be high enough to prevent inflammation caused by an implant. In a previous study using dexamethasone to suppress inflammation due to an implant, local distribu tion of the drug in s ubcutaneous tissue was not determined ( 1 ) Although dexamethasone is a comm only used antiinflammatory agent, its local concentration, diffusion coeffici ent, and rate of elimination have not been reported following subcutaneous release. The ability of dexamethasone to penetrate subcutaneous tissue can be measured, quantif ied, and compared to mathematical models ( 4 ). Therefore, the controlle d delivery of dexamethasone in normal rat subcutaneous tissue was used in order to help develop a fundamental understanding of how the drug is transported in the subcutaneous tissue. Becau se the efficacy of controlled interstitial delivery depends on the distance the drug can penetrate into the tissue surrounding the PAGE 9 2 implantable glucose sensor, [3H]dexamethasone was delivered from osmotic pumps implanted into the subcutaneous tissue of rats. Digital autoradiographic imaging was used to quantify the spatial distribution of ra dioactivity in the subc utaneous tissue at 6, 24, and 60 hours after subcutaneous implantation. Both the extent of penetration of dexamethasone and the effectiveness of simp le transport models for quantification of penetration were investigated. Additionally, many transport experiments are based on th e injection of a finite volume of substance into the tissue of interest which then diffuses away. Some examples of injectionbased diffusion experiments are th e determination of the diffusion coefficient of small molecules in the brain ( 7 ), the determination of th e diffusion coefficient of growth factors in the brain ( 8 ), and the determination of th e diffusion coefficient of drugs in tumors( 9 10 ). Knowledge of the diffusion of a s ubstance of interest in the tissue of interest is important for treatment efficacy. Therefore, it is also important to develop a method in which the diffusion coefficient of an injected substan ce in tissue can be determined in a relatively simple manner. Such a technique is investigated in this work by finding the diffusion coefficient of [3H]dexamethasone in rat subcutaneous slices after an injection. The purpose of this project was to devel op a technique to estimate the diffusion coefficient and the rate of elimination of dexamethasone in subcutaneous tissue. The technique includes finding the so lutions to relevant forms of the diffusion equation, and then using these solutions as models with unknown parameters: the diffusion coefficient and the elimination constant. Curve fitting measured autoradiography data to the models was then undertaken to estimate the diffusion co efficient and the elimination constant. In PAGE 10 3 this way, this newly developed technique aids in providing an understanding of the diffusion of dexamethasone in subcutaneous tissue. All animal experiments were performed unde r the approval of the University of South Florida Animal Care and Use Program. PAGE 11 4 Experimental Methods Materials [3H]Dexamethasone (392.46 MW), specifically [1,2,4,6,73H]dexamethasone, was obtained from Amersham Biosciences Corp (Piscataway, NJ). The specific activity was 88.0 Ci/mmol. Alzet osmotic pumps ( 1003D model) were obtained from Durect Corp. (Cupertino, CA). Subcutaneous Implantation Six male Sprague Dawley rats (Harlan, Indianapolis, IN, 375399 g) were used. The rats were initially anesthetized by pl acing each rat in an induction chamber filled with a 5% mixture of isoflurane in oxygen. During surgery anesthetization of the rats was maintained using a 2.5% mixture of is oflurane in oxygen. Two pumps containing radiolabeled dexamethasone were implanted subcutaneously on either side of the shoulders of the rat. A 34 cm incisi on was made between the shoulder blades. A hemostat was inserted into the incision on th e lateral aspect. By opening and closing the jaws of the hemostat, a pocket in the subcut aneous tissue just larg e enough for the pump was created. A tunnel to insert the tu bing was made using a blunt probe. Excess bleeding was removed with sterile cotton ga uze. The osmotic pump was implanted tubing end first. The wound was closed with 46 surgical staples. Two rats were sacrificed at 6, 24, and 60 hours after implanta tion. The rats were euthanized using CO2. The tissue around the tip of th e catheter was removed, quickly frozen on dry ice, and PAGE 12 5 stored at 80 C to immobilize the tracers within the tissue sample. The frozen tissue samples were mounted on a cryostat chuck and cut in 10 m sections. Sections taken at every 200 m were used for autoradiographic imaging. Preparation of Osmotic Pumps A solution of [3H]dexamethasone and sterile 0.9% (w/v) saline was loaded into the osmotic pumps (total volume 114 L ) using the protocol provided by the manufacturer. Each pump cont ained a total activity of 127 Ci. The pumps provided a controlled delivery at a rate of 1.0 L/hour. To prevent the pump from causing a tissue reaction at the site of drug delivery, a 4 cm length of polyethylene tubing was connected to the body of the pump. Subcutaneous Injection Three male Sprague Dawley rats (Harlan, Indianapolis, IN; 375399 g) were used. Again, the rats were euthanized using CO2 prior to the experiment. It took an average of 86 seconds for the tissue sections to freeze and is included in the total time of 2.5 minutes or 20 minutes. All tissue samples were then stored at 80 C to immobilize the tracers within the tissue sample. As before, the frozen tissue samples were mounted on a cryostat chuck and cut in 10 m thick sections. Sec tions taken at every 200 m were used for autoradiographic imaging. The 20 Minute Experiment Three 0.04 mL solutions of [3H]dexamethasone in sterile 0.9% (w/v) saline were used. Injections using insulin syringes were made into subcutaneous tissue on the backs of the rats. The approximate duration of the injection was 1 second. Each solution contained a total activity of 0.65 Ci. The tissue around the injec tion site (1 cm x 1 cm x PAGE 13 6 0.5 cm) was removed and frozen on dry ice to immobilize the [3H]dexamethasone in the tissue. The average time from injection to when the tissue froze, as measured using a surface thermometer (Mannix Testing & Measurement, Lynbrook, NY), was approximately 20 minutes after injection. The 2.5 Minute Experiment Three subcutaneous sections were harv ested (1 cm x 1 cm x 0.5 cm) from the backs of the rats. Each section was injected with a 0.04 mL solution of [3H]dexamethasone in sterile 0.9% (w/v) saline us ing an insulin syringe and then frozen on dry ice. The injection duration was appr oximately 1 second. The average time from injection to when the tissue froze was a pproximately 2.5 minutes after injection. Autoradiographic Imaging and Analysis Autoradiographic images of the tissue s ections were obtained using a recently developed realtime digital ra dioactivitydetection system, the MicroImager (Biospace Mesures, Paris, France) ( 11 12 ) With the MicroImager, acquisition of events can be visualized in realtime on a monitor screen. Each event is individually analyzed by the computer. An event is a radioactivity decay event ( 11 ). The acquisition of events needs only to proceed for as long as is necessary to obtain a good image. In the osmotic pump implantation case, autoradiographic imag es with between 380,715 to 686,390 events were acquired over 2445 hours to obtain good images. An optical image of the same tissue sample using the MicroImager was al so obtained. The spatial variation in drug concentration from the osmotic pumps was qua ntified in the following way. The areas of subcutaneous tissue were identified on the optical image and then superimposed onto the corresponding autoradiographic image. The co ncentration profiles in the subcutaneous PAGE 14 7 tissue surrounding the catheter tip were dete rmined directly from the autoradiographic images using Beta Vision+ software (Biosp ace Mesures, Paris, France). A line profile tool (1 mm wide) was used from the center of the catheter tip to the periphery of the subcutaneous tissue to obtain a number of events versus distance profile. The background number of events was subtracted from the number of events acquired. A number of events versus distance profiles were performed at 15 increments around the catheter tip on each section select ed for analysis. The number of events at the catheter tip opening was calibrated to the known concentra tion of the agent in the pump to obtain concentration versus distance profiles at 6, 24, and 60 hours after implantation. An example implantation autoradiographic image is shown in Figure 1. Similar processing was done for the injection experiments providing concentration versus distance profiles at 2.5 minutes and 20 minutes. For the injection experiments, autoradiographic images with between 1,593,815 and 1,918,869 events were acqui red over 71 hours 32 minutes to 72 hours 43 minutes to obtain good images. Exam ple injection autoradiographic images are shown in Figure 2. PAGE 15 8 Figure 1. Example Autoradiographic Image, Implantation. This figure shows an autoradiographic image from rat subcutan eous tissue obtained using the MicroImager after implantation of an osmotic pump containing [3H]dexamethasone for 6 hours. The arrow shows the location and di rection of the catheter tip. Each dot represents a radioactivity decay event. Li ghter shades indicate higher activity. The bar represents a distance of 1 mm. PAGE 16 9 Figure 2. Example Autoradiographic Image, Injection. Panels a and b show a typical autoradiographic image obtained using the MicroImager 2.5 minutes and 20 minutes after injection of [3H]dexamethasone. Lighter areas indicate higher activity. The bar represents a distance of 2 mm. PAGE 17 10 Mathematical Formulation In the diffusion equation, the diffusion coefficient (D) and the elimination constant (k) are parameters. Using the con centration versus distance profile information, estimated values for these parameters can be found by using a nonlinear curve fitting technique. Essentially, the measured concentr ation profiles can be f it to the appropriate solution to the diffusion equation to optim ally estimate D and k. The MarquardtLevenberg algorithm was chosen for the curve fitting due to its simplicity of programming, stability in searching, and rapid convergence to a good solution. The development of the required relevant solutions to the diffusion equation for implantation and injection results in three major forms, which are detailed below. First, the steady state solution for implantation is need ed. This solution essential requires that the time derivative be set to zero and that th ere is a constant con centration surface. Elimination is included for the steadystate case. Additionally, th e steadystate solution will be useful for expected and theoretical ex amination of concentration profiles. Also, the steadystate solution is essential as an in termediate step in deriving the full transient solution for the implantation experimental m odel. The steadystate solution was found using a standard series technique. Secondly, the full solution to the diffusion equation with elimination is needed for modeling the implantation experiments. This solution is quite involved and requires Laplace transform t echniques, which are also detailed below. Third, the full transient solution with no elimin ation but with a given initial concentration in a confined volume is needed for modeli ng the injection experi ments. Armed with PAGE 18 11 these solutions, curve fitting can be undertaken to estimate the diffusion coefficients and elimination constants for the various experiments in this work. Mathematical Formulation, Implantation The concentration profiles of [3H]dexamethasone obtai ned using the MicroImager were compared to mathematical m odels of drug diffusion and elimination. The model assumed constant drug concentration at the catheter tip/tissue interface, firstorder elimination of drug; homogeneous and isotr opic diffusion transport of drug through the subcutaneous tissue, negligib le fluid convection, and spheri cal symmetry. The governing equation for the diffusion and elimination of a drug in subcutaneous tissue is: C k r C r r C D t C 22 2. (1) C is the concentration of the dr ug in the subcutaneous tissue; D is the diffusion coefficient of the drug in subcutaneous tissue; r is the radial distance from the center of the catheter tip; k is the firstorder elimination constant for the drug from the subcutaneous tissue; and t is the time after implanta tion. The boundary and initial conditions are: , 0 0 C , 0 C C , 0 0 C0 r t a r t a r t (2) where a is the radius of the catheter and C0 is the concentration at the catheter tip. The solution ( 13 ) of equation 1 using the boundary and in itial conditions of equations 2 is: t k t D a r erfc D k a r r aC t r C 2 ) ( exp 2 ,0 PAGE 19 12 t k t D a r erfc D k a r 2 ) ( exp. (3). The methods for solution of the diffusion equation are provided below. Assuming steadystate and applying the boundary conditio ns from equations 2, equation 1 can be solved using a series technique as is show n below. Alternatively, the steadystate solution can also be found from equa tion 3 by letting time go to infinity: 1 exp0a r D k a r a C r C (4) However, this would require knowing the full so lution to begin with. The series solution technique for the steadystate case and the methods to find the general solution of the diffusion equation are complex and are not easily found in the literature. They are provided here to show their corre ctness and detailed structure. The Brownian diffusion coefficient for dexamethasone in water was estimated from the StokesEinstein equation: A sN r T R D 6 (5) rs = 0.657 M1/3 [x 1010 m] and is the equivalent spherical solute radius; M is the molecular weight of dexamethasone (392.46 MW); R is the ideal gas constant 8.314 J K1 mol1; T is temperature; is the dynamic viscosity; and NA is Avogadros number. The calculated diffusion constant of de xamethasone in water at 37 C is D = 6.82 x1010 m2/s. The StokesEinstein equation under pred icts the actual diffusion coefficient for small solutes of molecular weight less than several hundred, and over predicts it for large solutes of molecular weight gr eater than several thousand ( 14 ). Measured concentration profiles for [3H]dexamethasone at t = 6 hours, t = 24 hours, and t = 60 hours were PAGE 20 13 compared to the transient equation 3. Typical concentration profiles predicted by equation 3 for two values of for various times are shown in Figure 2, a and b. The dimensionless parameter where = D k a, is analogous to the Thiele modulus obtained in analysis of heterogeneous catalysis ( 4 ) and is a predictor of the extent of drug penetration from the catheter tip. The radius of the catheter was appr oximated at 0.6 mm. Values for D and k were found in the following manner. First, initial estimates for D and k were found. For the initial estimate of D the diffusion constant for dexamethasone in water was used, D= 6.82x1010 m2/s. This D was used in the steadystate solution of the diffusion equation (equation 4) to find an initial estimate for k which was found by using the MarquardtLevenberg technique ( 15 ) with two independent variables, k and C0, to minimize the residual of the sumsquarederror between the predicted and experiment al concentrations. Note that the steadystate solutions are genera lly expected to be close to the profiles for the 6, 24, and 60 hour data. (Figure 3.) Second, these initial estimates for the k and D values were then used as the starting point s for the MarquardtLevenberg algorithm using the transient equation (equation 3) with the k D and C0 being the three independent variables over which the resi dual of the sum squared erro r between the predicted and experimental concentrations was to be minimized. The initial value for C0 was always the maximum concentration in the measured data set. The MarquardtLevenberg algorithm efficiently searched over the k D and C0 space to find the point which best fits the data ( 16 ) This technique was repeated to find k and D for 6, 24, and 60 hours. For each of these times, the calculations were repeated for the autoradiographic scans at various angles. The MarquardtLevenberg algorithm was written in MATLAB. PAGE 21 14 Figure 3. Theoretical Concentration ve rsus Distance Profiles for Implantation. Concentration versus distance profiles can be obtained by solving the transient diffusion and elimination equation (Equation 3) for various times until steadystate is reached. Panels a and b demonstrate the depe ndence of the penetration depth with the modulus D k a / Panel a (modulus 0.2) has a larger penetration depth than panel b (modulus1). PAGE 22 15 Mathematical Formulation, Injection In order to investigate th e diffusion coefficient for the injection experiments, the model required four specific assumptions. The model assumed (1) that the diffusing substance is deposited with in a sphere of radius a at t =0, (2) isotropic and homogeneous diffusion transport of drug through the subc utaneous tissue, (3 ) negligible fluid convection, and (4) negligible elimination. A ssuming that the elimination is negligible is justified as tissue samples were obtained from a sacrificed rat. Note that the absence of blood flow eliminates most clearance mechanisms normally present in vivo ( 8 ). Hence, the governing equation for diffusion of a drug in the subcutaneous tissue is r C r r C D t C 22 2. (6) The following initial conditions al so come from the assumptions: 0 0 C , 0 C C , 0 0 C0 r t a r t a r t (7) The analytic solution for equation 1 using the above initial and boundary conditions is ( 13,17 ): 4 / ) ( exp 4 / ) ( exp 2 2 2 2 ,2 2 2 / 1 0 t D a r t D a r r t D t D a r erf t D a r erf C t r C (8) where a is the radius of the sphere. If r is very much greater than a then expression 8 becomes ( 13 ): PAGE 23 16 t D a t D r t D r t D m C 40 6 1 4 exp 82 2 2 2 / 3 (9) where m = VC o = 4/3 a3C o and V is the injected volume. If the radius of the sphere tends to zero, a0, with m remaining constant ( 13 ), t D r t D m C4 exp 82 2 / 3, (10) or t D r t D a C t r C4 exp 6 ,2 2 / 3 2 / 1 3 0. (11) Equation 11 has the same solution as that for the instantaneous point source in 3 D ( 18 ). However, Nicholson ( 19 ) suggests that, at measurement locations sufficiently far from the source, equation 10 (or equation 11) will provide a useful approximation. Moreover, Thorne et al. ( 8 ) suggest that when the in jection time is very brief compared to the time of the subsequent diffusion measurements the concentration can be described by equation 11 (or equation 10). Typica l concentration profiles for [3H]dexamethasone predicted by equation 11 at t =2.5 minutes, t =5 minutes, t =10 minutes, and t =20 minutes are shown in Figure 4. The radius of the injected spherical volume was 2.1 mm. For equation 11 to be a useful approximation, data away from the source were used ( 19 ). For the 20 minute experiment, a portion of the c oncentration profile from the tail end was used in the mathematical model. This porti on ranged from the ta ilend to a position 3 mm toward the source from the firstzero event value. The tailend was defined as the furthest position on the curve from the source. The firstzero event value was defined as the position where the events first decreas e to a zero value. For the 2.5 minute PAGE 24 17 experiment, a firstzero event value could not be used as a reference point. Instead, a location on the profile where th e profile bends from a steep curve to a plateau region was used as a reference. This bend was defi ned to occur at a position where the number of events was 100. The number of events at the bend is dependent on the detection time in the MicroImager. If the detection time is increased, more events are detected. However, the detection time was the same for all profiles in the 2.5 min experiment. Therefore, the portion of the concentration pr ofile used in the mathematical model was from the tailend to a distance 0.7 mm toward the source after the bend. The reason that the firstzero event value coul d not be used as a referen ce point is discussed further below. The value for D was found in an iterative manner. First, an initial estimate for D was needed. The diffusion constant fo r dexamethasone in water was used, D =6.82x1010 m2/s. This initial estimate for D was then used as the star ting points for the MarquardtLevenberg algorithm ( 15 ) with the D and C0 being the two indepe ndent variables over which the residual of the sumsquarederro r between the predicted and experimental concentrations was to be minimized. The initial value for C0 was always the maximum concentration in the measured data set. The initial values for D and C0 should be physiologically reasonable to avoid final estim ates that are based on meaningless local minima in the sumsquarederror function. While sensitivity to the initial guess was not specifically investigated, careful thought was given as to wher e to start the searching in the D and C0 space. Quick convergence, along with physiologically reasonable results using the Marquardt algorithm, led us to be lieve that our strategy was generally sound. The MarquardtLevenberg algorithm efficiently searched over the D and C0 space to find the point which best fits the data ( 16 ). This technique was repeated to find D at 2.5 and PAGE 25 18 20 minutes for each scan. The MarquardtLe venberg algorithm was again written in MATLAB. PAGE 26 19 Figure 4. Theoretical Concentration ve rsus Distance Profiles for Injection. Concentration profiles for diffusion when a c oncentrated bolus of solute is deposited within a small region are shown. The curv es shown result from equation 5 with D = 4.11x1010 m2/s, a =2.1 mm, and t =2.5 minutes, 5 minutes, 10 minutes, and 20 minutes. PAGE 27 20 Figure 5. Example Measured Injection Profiles. Typical number of events versus distance profiles obtained usi ng the MicroImager at (a) 2.5 and (b) 20 minutes after injection are shown. Data from only one s can are shown. The ordinate represents the location of the center of injection. For (a) all data shown was used for diffusion coefficient estimation. PAGE 28 21 SteadyState Solution With Elimination, Implantation The diffusion partial differential equation with spherical symmetry that is needed as the starting expression for this work is: kC r C r r C D / 2 / 02 2. (12) Now divide through by D and set D k / to 1k. C k r C r r C1 2 2/ 2 / 0 (13) Next, for clarity we shift to the following notational form which is the same equation but closer to classical notation used for this category of problem. Let C be replaced by y and the independent variable r by x. We then get: 0 2 '1 2 2 y k x y x y x (14) This differential equation can be solved using a power series. By taking successive derivatives of the general series fo rm assumed as the solution we get: n r n n r n r n nx a x a x a y 1 0 0, (15) 1 0) ( n r n nx a n r y 1 1 1 0) ( n r n n rx a n r x ra, (16) 2 0) 1 )( ( ' n r n nx a n r n r y 2 1 2 0) 1 )( ( ) 1 ( n r n n rx a n r n r x a r r. (17) PAGE 29 22 Now substitute these derivatives into the differential equation. n r n n rx a n r n r x a r r 1 0) 1 )( ( ) 1 ( (18) Combining terms and simplifying, a sequence of straightforward steps can be taken: rx a r r r02 ) 1 ( n r n nx a n r n r n r 1) ( 2 ) 1 )( ( (19) rx a r r0) 1 ( n r n nx a n r n r 1) 1 )( ( (20) rx a r r0) 1 ( n r n nx a n r n r 1) 1 )( ( (21) 1 1 0) 2 )( 1 ( ) 1 ( r rx a r r x a r r (22) 0 ) 1 )( (2 2 1 n r n n nx a k a n r n r n r n n rx a n r x a r 1 0) ( 2 2 02 1 1 1 2 0 n r n n rx a k k x a. 02 0 1 n r n nx a k. 02 0 1 n r n nx a k. 02 2 1 n r n nx a k PAGE 30 23 At this point there are many different path s that can be taken to provide the two independent solutions needed for our second or der situation. For our case the best and most general choice will be 1 r and0 ,1 0 a a We now need to develop a recurrence relation for 2 n. Therefore, if both 0a and 1aare taken as nonzero values, for 2 n the recurrence relationship that must always be true is: ) 1 )( (2 1 n n a k an n (23) Specifically, for example, this would mean: ) 1 )( 2 (0 1 2a k a (24) ) 2 )( 3 (1 1 3a k a (25) 40 2 1 4a k a (26) 51 2 1 5a k a (27) 60 3 1 6a k a (28) 71 3 1 7a k a (29) 80 4 1 8a k a (30) and 91 4 1 9a k a (31) PAGE 31 24 For the even series, in general for ,... 3 2 1 n we get )! 2 (0 1 2n a k an n. (32) For the odd series, in general for ,... 3 2 1 n we get )! 1 2 (1 1 1 2 n a k an n. (33) Going back to the original series form with 1 r we get: ...5 6 4 5 3 4 2 3 1 2 1 0 0 x a x a x a x a x a a x a x a yn r n n or (34) ... 6 5 4 3 25 0 3 1 4 1 2 1 3 0 1 2 2 1 1 0 1 1 0 x a k x a k x a k x a k x a k a x a y (35) If we now let 10 a and 1 1k a this general solution becomes: ... 6 5 4 3 2 15 3 1 4 1 2 1 3 1 2 2 1 1 1 1 x k x k k x k x k k x k k x y (36) This can now obviously be written as: exp1x x k y (37) Similarly, if we chose 10 a and 1 1k a the solution becomes: exp1x x k y (38) Returning to our original not ation, these 2 independent so lutions can be combined to form the full solution: exp exp r r D k B r r D k A r C (39) PAGE 32 25 Next, the boundary conditions 0C a C and 0 C must e met. Therefore, the steadystate solution is: 1 exp 1 exp0 0 a r r a C a r D k a r a C r C (40) This is the first important equation that we will need to model the diffusion behavior in some of our experimental situations. This equation will be particularly useful in finding a first estimate of k, as described earlier. This equation will also be useful in deriving the more sophisticated transient cases as will now be shown. General Solution, Implantation We will now take recourse to the Laplace method to derive a solution for the full diffusion differential equation. The ultimate solution will involve the elimination term, but first we must go back to the solution fo r the diffusion equation without elimination. Laplace transform methods will be used to de rive the full solution. However, if the solution for the diffusion equation without elim ination is understood, it will make finding the full solution easier. Taking the Laplace tr ansform of Equation 1 with a zero initial condition can quickly identify the Laplace transform for the diffusion equation without elimination. It must simply be the st eadystate solution with k replaced by s: a r D se r a C r s C 0) (. (41) When returning to the time domain this formula will handle the 0 0 r C initial condition. Next consider some of the si mple properties of the Laplace Transform: PAGE 33 26 t f at a s F exp, where .t f s F However, the differential equation we really want to solve is: C k r C r r C D t C 22 2. (42) By going through a similar argument to that ju st given for the no elimination case we can make progress toward the desired full solu tion. Taking the Laplace transform of both these equations with the zero initial conditions (whi ch is our case) we get: 22 2C k r C r r C D C s (43) The bar over the concentration term is the notation for the Laplace transform. Now rearranging we can say: 2 02 2C s k r C r r C D (44) This means the transform we want to work with is: ) (0 a r D k se r a C r s C (45) Next, the boundary condition 0) ( C t a C for t>0 must be taken into account. This can be done by multiplying by a s/ 1 factor. The need for doing this is clear since there can only be a solution for time greater than zero. This requirement comes from the basic properties of the Laplace transform as can be verified by examining any Laplace transform table. Therefore, the transform we want to invert is: PAGE 34 27 ) (0a r D k se s r a C r s C (46) Note that we can expect to get the final a nd full solution to the diffusion equation because all the initial and boundary conditions will be met. However, an algebraic manipulation is needed before we will find the proper expressions in Laplace transform tables. This algebraic manipulation is: a r D k s a r D k se k s k s s r a C e s r a C r s C 2 2 1 ) (0 0 a r D k se k s s k s r a C 1 2 20 (47) a r D k se k s k k s k k s r a C 1 1 1 20 (48) 1 1 1 1 20 a r D k s a r D k se k s k k s e k s k k s r a C (49) The two terms in the previous expression are commonly available in Laplace transform tables (20) This will allow us to return to the time domain. The needed transform pairs are: t c t c erfc c c e c s c c sc s c2 exp 1 11 2 2 1 2 2 22 1, (50) t c t c erfc c c e c s c c sc s c2 exp 1 11 2 2 1 2 2 22 1. (51) PAGE 35 28 Remember that because of the exponential te rms, we are making a shift in the complex plane. That is, when coming back to th e time domain there must be a factor of ) exp(2t c if the s is shifted by 2c This is why the ) exp(2t c term drops out of the transform pair found in the tables. 2 exp exp 1 11 2 2 2 1 21 t c t c erfc t c c c e s c ss c (52) The t c2exp term drops out due to the 2c s terms. To get back to the time domain just let D a r c 1 and k c 2. Note that we have assumed k>0 in this development. So finally we get: t k t D a r erfc D k a r r aC t r C 2 ) ( exp 2 ,0 t k t D a r erfc D k a r 2 ) ( exp (53) General Solution, Injection Using the assumptions described earlier for the injection st udies, the governing equation for diffusion of a drug in the subcutaneous tissue is: r C r r C D t C 22 2. (54) The following initial conditions al so come from the assumptions: 0 0 C , 0 C C , 0 0 C0 r t a r t a r t (55) PAGE 36 29 A transformation of the diffusion equation can be used to simplify the problem by making a simple variable change. 2 2 / ) ( ) (2 3 2 2 2r u r u r u r C r u r u r C r r t u r t C (56) By substitution, this reduces the problem to the following simpler case: 2 2r u D t u (57) The initial conditions change to: 0 0 u , 0 C u , 0 0 u 0 r t a r t r a r t (58) The point source solution to equation 54 is easily found and verified by substitution. This solution is: ) 4 /( 2 / 12) ( 2 1 ) (t D r pe t D r t u (59) If equation 59 is a solution, then the follo wing equation must also be true for some general function r f, which can also be quickly verified by substitution: ) ( ) ( 2 1 ) () 4 /( ) ( 2 / 12dr e r f t D r t ut D r r (60) In our particular case the f unction we want to chose is r f = C0 r C0 will be ignored for now. No w, let C(t,r) = C0 v(t,r) = C0 u(t,r)/r. Using the in itial conditions (58), we get: PAGE 37 30 ' ) ( 2 1 ) () 4 /( ) ( 2 / 12dr e r t D r t ut D r r a a (61) Next, it is just a matter of simplifying the ma thematic expressions and getting them into a more tractable form. First, just split the integral: ' ' ') 4 /( ) ( 0 ) 4 /( ) ( 0 ) 4 /( ) (2 2 2dr e r dr e r dr e rt D r r a t D r r a t D r r a a (62) Now just substitute r for r in the second integral, changing the integral and its limits appropriately one gets: " ') 4 /( ) ( 0 ) 4 /( ) ( 02 2dr e r dr e rt D r r a t D r r a (63) Reversing the limits of th e second integral, we get: " ') 4 /( ) ( 0 ) 4 /( ) ( 02 2dr e r dr e rt D r r a t D r r a (64) Then reducing these integr als into one, we get: ') 4 /( ) ( ) 4 /( ) ( 02 2dr e e rt D r r t D r r a (65) Returning to the full expression for the solution, but still disregarding C0, we get: ' ) ( 2 1 ) () 4 /( ) ( ) 4 /( ) ( 0 2 / 12 2dr e e r t D r r t vt D r r t D r r a (66) Doing the indicated integration is straightforwar d. By simple change of variables the first part of the integral can be easily evaluated. Letting r r r where r is taken as a PAGE 38 31 constant, changing integration limits appropriately, and splitting the integral into two parts, gives the following sequence of equalities: ) ( ') 4 /( ) ( ) 4 /( ) ( 02 2dr e r r dr e rt D r a r r t D r r a (67) ) () 4 /( ) (2dr e r rt D r r a r (68) " ") 4 /( ) ( ) 4 /( ) (2 2dr e r dr e rt D r r a r t D r r a r (69) Next, the precise definitions of the error func tion needs to be understood and adjusted for these integrals. The error function is: d e x erfx022 ) (. (70) Now making the variable change: t D r 4 ', (71) we get: Dt dr e x erfDt x Dt r4 2 ) (4 0 4 '2 (72) With further simplification the following equality results. 1 40 4 '2dr e t D Dt r erfr Dt r (73) Using equation 69, further progress ca n now be made with equation 66. PAGE 39 32 " ") 4 /( ) ( ) 4 /( ) (2 2dr e r dr e rt D r r a r t D r r a r " ") 4 /( ) ( ) 4 /( ) (2 2dr e r dr e rt D r r a r r a r t D r (74) Now, just doing the integration and using the sp ecified limits the resulting expression is: Dt a r erf t D r Dt r erf t D r 4 4 Dt a r Dt re t D e t D4 ) ( 42 22 4 2 4 (75) By symmetry a further integral re sult can be immediately written: ') 4 /( ) ( 02dr e rt D r r a Dt a r erf t D r Dt r erf t D r 4 4 Dt a r Dt re t D e t D4 ) ( 42 22 4 2 4 (76) Equations 74 and 76, when put into equati on 66, provide the solution needed for the injection studies. This final result is: 4 / ) ( exp 4 / ) ( exp 2 2 2 2 ,2 2 2 / 1 0 t D a r t D a r r t D t D a r erf t D a r erf C t r v (77) PAGE 40 33 where C0 is the concentration of the drug in subcutaneous tissue, D is the diffusion coefficient of the drug in subcutaneous tissue, r is the radial distance from the center of the injection, and t is time. The initial concentration is C o in the sphere 0 r < a and zero for r > a Additionally, a boundary condition is C (+ t ) = 0. The analytic solution for equation 54 using the above ini tial and boundary conditions is ( 13,17 ): 4 / ) ( exp 4 / ) ( exp 2 2 2 2 ,2 2 2 / 1 0 t D a r t D a r r t D t D a r erf t D a r erf C t r C(78) where a is the radius of the sphere. This completes the needed derivations and mathematical background for this research. Marquardt Curve Fitting Equipped with the solutions for the st eadystate and general solutions to the diffusion equation, curve fitting techniques can be employed to do the estimation of the D and k values. The Marquardt technique offers a way to do the needed curve fitting (13, 27) The Marquardt method is based on two prin ciples. It combines a gradient search with a GaussNewton technique. This method ba lances these two prin ciples to provide a stable yet efficient search in the soluti on space for a mulivariable nonlinear modeling formula. In this particular cas e, the desire is to find estimate s of the diffusion coefficient, the elimination constant, and the maximum count value. In essence, the desire is to find the best values for these parameters to f it the measured data. The Marquardt method minimizes a sumsquarederror expression to get the best estimate for the desired parameters. A brief description of the form ulation of the Marquard t Method is given in PAGE 41 34 Appendix I. An example of the actual Ma rquardt program in MATLAB is provided in Appendix II. To start the Ma rquardt algorithm, a beginning point in the solution space must be specified. The choice of the star ting point for the least squares optimization process is crucial for stability and quick c onvergence. How the ini tial point is found for the various cases related to implantation a nd injection have already been described. PAGE 42 35 Results Radiolabeled dexamethasone spread through the subcutaneous tissue after implantation of the osmotic pump (Figure 1). The local concentration of drug within the tissue was quantified from the autoradiographi c images using the Beta Vision+ software. The Beta Vision+ software was used to cons truct the number of events as a function of distance profiles. An event is a radioactivity decay event ( 12 ) The number of events was greatest at the tip of the catheter. A high nu mber of events on the autoradiographic image represent a high drug concentration. The number of events at the tip of the catheter can be calibrated to the known con centration in the pump. Hence, the local concentration of the drug in the subcutaneous tissue surr ounding the catheter can be estimated by comparing the local number of events to the number of events at the catheter tip. In general, at distances more than a few millim eters from the catheter tip, the radioactivity was not significantly different from bac kground. Figure 1 is repr esentative of the autoradiographic images obtained using the MicroImager after implantation of the osmotic pumps for 6, 24, or 60 hours. C oncentration profiles obtained from the autoradiographic images of the subcutaneous tissue surrounding the catheter tip were examined and compared to the mathematical model of diffusion and firstorder elimination to find the best estimates for D and k Figures 6 to 11 are representative of the measured profiles with their final curv e fit. The best estimates obtained for D and k are given in Table 1. A si nglefactor analysis of va riance (ANOVA) indicated (26) that there was no significant difference betw een the 6 hour and 24 hour data for k ( p > 0.05) PAGE 43 36 or for D ( p > 0.05). There was not enough data at 60 hours for comparison. The average, based on the 6, 24, and 60 hour data, for the diffusion coefficient is D = 4.11 1.77x1010 m2/s and for the elimination constant is k = 3.65 2.24 x 105 s1. To quantify differences in drug penetration with time after releas e from the osmotic pump, the best fit concentration profiles were used to find the distance where the local concentration drops to 10% of its maximum value. For the 6 hour cas e, the majority of the drug was confined to a region within 2.22 0.42 mm from the tip of the catheter. For the 24 and 60 hour cases, the majority of the drug wa s confined to a region within 2.70 0.38 mm and 1.80 mm from the tip of the catheter, respectivel y. (Table 2.) The penetration distance of [3H]dexamethasone increased from 6 to 24 hours but decreased from 24 to 60 hours. Radiolabeled dexamethasone spread through the subcutaneous ti ssue after injection. Figure 2 is representative of the autora diographic images obtai ned using the MicroImager after injection of a radiolabeled dr ug. The local concentration of drug within the tissue was again quantified from the autora diographic images using the Beta Vision+ software. Figure 5 is representative of th e number of events versus distance profiles obtained from the autoradiographic images for the injection experiments. The number of events was greatest at the center of the injec tion. The number of events can be calibrated to concentration to obtain concentration vers us distance profiles. For the 20 min case, the concentration profile from the tailend to 3 mm toward the source from the firstzero event value was compared to the mathemati cal model of diffusion to find the best estimate for D For the 2.5 min case, the concentra tion profile from 0.7 mm toward the source from the bend to the tailend was co mpared to the mathematical model of PAGE 44 37 diffusion to find the best estimate for D The best estimates obtained for D are given in Table 5. PAGE 45 38 Figure 6. Implantation Profile with Curve Fit, 24 Hours, 45 Degrees. A concentration versus distance example profile obtaine d by solving the transient diffusion and elimination equation (equation 3) is shown. This curve fit is for the data obtained from the 24 hour case at an angl e of 45 degrees. Note that th ere were points that went into the curve fit in producing this graph that fell below the zero line due to background adjustments that are not shown. PAGE 46 39 Figure 7. Implantation Profile with Curve Fit, 24 Hours, 0 Degrees. A concentration versus distance example profile obtaine d by solving the transient diffusion and elimination equation (equation 3) is shown. This curve fit is for the data obtained from the 24 hour case at an angl e of 0 degrees. Note that th ere were points that went into the curve fit in producing this graph that fell below the zero line due to background adjustments that are not shown. PAGE 47 40 Figure 8. Implantation Profile with Curve Fit, 60 Hours, 195 Degrees. A concentration versus distance example pr ofile obtained by so lving the transient diffusion and elimination equation (equation 3) is shown. This curve fit is for the data obtained from the 60 hour case at an angle of 195 degrees. Note that there were points that went into the curv e fit in producing this graph that fell below the zero line due to background adjustments that are not shown. PAGE 48 41 Figure 9. Implantation Profile with Curve Fit, 6 Hours, 60 Degrees. A concentration versus distance example profile obtaine d by solving the transient diffusion and elimination equation (equation 3) is shown. This curve fit is for the data obtained from the 6 hour case at an angle of 60 degrees. Note that there were points that went into the curve fit in producing this graph that fell below the zero line due to background adjustments that are not shown. PAGE 49 42 Figure 10. Implantation Profile with Curve F it, 6 Hours, 75 Degrees. A concentration versus distance example profile obtaine d by solving the transient diffusion and elimination equation (equation 3) is shown. This curve fit is for the data obtained from the 6 hour case at an angle of 75 degrees. Note that there were points that went into the curve fit in producing this graph that fell below the zero line due to background adjustments that are not shown. PAGE 50 43 Figure 11. Implantation Profile with Curve Fit, 6 Hours, 30 Degrees. A concentration versus distance example profile obtaine d by solving the transient diffusion and elimination equation (equation 3) is shown. This curve fit is for the data obtained from the 6 hour case at an angle of 30 degrees. Note that there were points that went into the curve fit in producing this graph that fell below the zero line due to background adjustments that are not shown. PAGE 51 44 Figure 12. Injection Profile with Curve Fit. Concentration profiles near the tail end at (a) 2.5 and (b) 20 min after injection are shown. Data from only one scan per time period are shown. Combining data from all other scans would make the figure unreadable. The solid lines show the diffusion model in which D and C0 were varied to minimize the residual of the sums quarederror between the predicted and experimental values. PAGE 52 45 Time after implantation [hours] k [1/s] x 105 D [m2/s] x 1010 6 4.80 2.56 3.63 1.06 24 2.52 1.65 4.92 1.97 60 4.70 1.73 Table 1. Resultant Estimated D and k Values Using Implantation. The diffusion coefficient and elimination constant were determined by fitting a model of diffusion and elimination to the concentration profile s measured near the tip of a catheter attached to an osmotic pump. For the 60 hour data, only was very thin, allowing measurement without boundary effects only in one case. The outliers and obviously questionable numeric values were excluded to arrive at these statistical results. PAGE 53 46 Time after implantation [hours] Penetration distance [mm] D k a 6 2.22 0.42 0.22 24 2.70 0.38 0.14 60 1.80 0.31 Table 2. Penetration Distance of Radioac tivity from the Tip of the Catheter. The penetration distance is the distance where th e local concentration drops to 10% of the concentration at the cathete r tip. This radial distance was found using the best fit curve through the data and corre sponds to the location where C / C o = 0.1. The dimensionless parameter, D k a, determines the extent of drug penetration and was found using the corresponding k and D values in Table 1. PAGE 54 47 Medium D [m2/s] Reference Water 6.82 x1010 StokesEinstein equation Subcutaneous tissue 4.11 1.77x1010 This study Brain 2.0 x1010 Saltzman and Radomsky, 1991 ( 21 ) Cellulose acetate membrane 3.15 x1010 Barry and Brace, 1977 ( 22 ) Table 3. Diffusion Coefficients for Dexamethasone in Various Media. The diffusion coefficient of dexamethasone in subcutan eous tissue was compared to the diffusion coefficient for dexamethasone in other media from the literature. PAGE 55 48 Agent k [1/s] Reference RSAa 6.42 1.19 x 105 Kim and Burgess, 2002) ( 24 ) Dexamethasone 3.65 2.24 x 105 This study VEGFb 3.50 1.03 x 105 Kim and Burgess, 2002) ( 24 ) Table 4. Elimination Constants of Vari ous Agents in Subcutaneous Tissue The elimination constant of dexamethasone in subcutaneous tissue was compared to the elimination constant of other agents in subcutaneous tissue from the literature. a rat serum albumin. b vascular endothelial growth factor. PAGE 56 49 Time after injection (min) D x 1010 (m2/s) 2.5 2.69 1.08 20 4.01 2.01 Table 5. Estimated D from Injection Experiments The diffusion coefficient was determined by fitting a model of diffusion to the concentration profiles from the tailend of the profiles for injection analysis. PAGE 57 50 Discussion The purpose of this research was to i nvestigate the capabi lities of using the solutions to the diffusion equation with s pherical symmetry along w ith a nonlinear curve fitting technique to estimate the diffusion co efficient and the elimination constant for dexamethasone in subcutaneous tissue. While the estimates were in the expected range and the general technique did work, there is much than can be done to improve the overall estimation calculations. One of the shortcomings was that some data had to be considered outliers and eliminated for doing the final estimation (Table 1). Specifically, the data at some angles was such that the tails were clearly not representative of a diffusion phenomenon. One explanation for this discrepancy is that, in this analysis, convection was not taken into account. Du e to channeling by anatomical structures convection could have been significant in so me directions. Future research could investigate the importance of convection to get a better result in the modeling of the subcutaneous mass transfer. A second weakness in what has been done here is that only the tails of the concentration profiles were used to get the estimates of the diffusion coefficients and elimination constants. Th is was appropriate because of the focus on diffusion and elimination. However, an improve ment in this technique could be made if different criteria were used to decide which data represente d the diffusion process. Here the maximum looking back from the tail was us ed as the cutoff point for curve fitting. Other criteria might be established for achievi ng better results. For example, a curvature criterion could be used to establish which data on the profile tails shou ld be used in curve PAGE 58 51 fitting. Clearly, much research could be done to investigate how to better estimate the diffusion and elimination properties. Another area that could be investigated is the control of the Marquardt algorithm. While the Marquardt algorithm worked quite well, there are some controlling parameters that could be investigated to improve the re sults. Specifically, th e initial lambda setting could be investigated in detail to make sure it is being set appr opriately. There were instances of curve fitting for which the initial lambda was changed so that the initial first step would not be so large as to put th e calculation into a mode which provided nonphysiological results. While, in this work, the Marquardt algorithm generally provided physiological results, it may be possible to discover strategi es for picking lambda which are more methodical than what was done here Similarly, better starting points for the Marquardt search may be able to be found w ith more work. The strategies used for picking the initial lambda and starting point were carefully considered, but no investigation was actually done to see if other strategies might be more appropriate as this issue was beyond the scope of this research. The diffusion coefficient, D of dexamethasone in subcutaneous tissue at 6 and 24 hours after implantation was 3.63 1.06 x1010 m2/s and 4.92 1.97 x1010 m2/s, respectively. The 60 hour data suggests a D of 1.73 x1010 m2/s. There was no significant difference between the 6 and 24 hour data for D ( p > 0.05). A comparison with the 60 hour data was not made as the sample size was too small. Even though the concentration profile at 6 hours has not yet reached stea dystate (Figure 3a), the value found for D should not be different from that found for the 24 hour case, which is very close to steadystate. (Note that at 6 hours = 0.22 (Table 2), and Figure 3a shows concentration PAGE 59 52 profiles for = 0.2 for various times.) The concentr ation profile at 60 hours has reached steadystate (Figure 1a). D should have similar values at 6, 24, and 60 hours because the best estimate for D and k for all cases was achieved using the transient diffusion and elimination equation (equation 3). Since the tr ansient equation takes time into account, be it for a short time period or for a long time period, the D and k values for the same agent in the same tissue should be the same. D and k are assumed to be constants. As time becomes large, the transient equation (equati on 3) reduces to the steadystate equation (equation 4). Hence, the average diffusion coefficient D = 4.11 1.77 x1010 m2/s, based on the 6, 24, and 60 hours data, results in a reasonable value for dexamethasone in subcutaneous tissue. The diffusion coefficient of dexamethasone in subcutaneous tissue is slightly less than in water but slightly greater than in brain tissue (Table 3). Our diffusion coefficient for dexamethasone in rat subcutan eous tissue is slightly greater than the diffusion coefficient of sodium fluorescei n (MW 376) in rat subcutaneous tissue D = 2.35 0.24 x1010 m2/s (23). Sodium fluorescein has a molecula r weight similar to that of dexamethasone (MW 392). The elimination constant, k at 6 and 24 hours was 4.80 2.56 x 105 s1and 2.52 1.65 x 105 s1, respectively. The 60 h data suggests a k of 4.70 x 105 s1. There was no significant difference between the 6 and 24 hour data for k ( p > 0.05). A comparison with the 60 hour data wa s not made, as the sample size was too small. The average, based on the 6, 24, and 60 hour data, for the elimination constant is k = 3.65 2.24 x 105 s1. This value is quite reasonable de spite the fact that the 6 hour case has not yet reached steadystate for the reasons given in the paragraph above. Table 4 shows values for k of other agents in subcutaneous tis sue. Our elimination constant for dexamethasone in rat subcutaneous tissue is sl ightly greater than th at of dexamethasone PAGE 60 53 in rat brain k = 1.19 x 105 s1 (25). Although only two rats were used for each time point, variation was not observed between the two rats as they were the same age, sex, size, and strain and were all from the same ve ndor. A detailed study would be useful to demonstrate that the age, sex, size, strain and vendor have no si gnificant effect on the values of D and k When a substance is injected into tissue in a period that is effectively instantaneous, it may exhibit two distinct behavi ors: (1) form a fluidfilled cavity or (2) infiltrate the extracellu lar space of the tissue ( 19 ). The subsequent diffusion from each case can be described by its own set of expressions ( 19 ). In this study, we have assumed that the substance does not form a cavity but infiltrates the extracellular space and then diffuses away. Hence, the appropriate soluti ons and approximations have been used for this case. The approximations to the case wh ere substance infiltrates the extracellular space lead to equation 11. The two criter ia for equation 11 to provide a useful approximation are that the measurement locations be sufficiently far from the source ( 19 ) and that the injection time is very brief comp ared to the time of the subsequent diffusion measurements ( 8 ). To comply with these criteria, the data near the tailend of the concentration profiles were used as describe d below. The measurement distance was kept as small as possible while large enough to provide meaningful data. To investigate criterion 2, two diffusion times were chosen, t =2.5 minutes and t= 20 minutes. For this study, radiolabeled dexamethasone was intr oduced into the subcutaneous tissue by injection. The highest concentrations of the agent were assumed to be at the location of the center of the injection. This assumption is supported by our theoretical curves (Figure 4). The local distribution of the ag ent in the subcutane ous tissue surrounding the PAGE 61 54 center of injection was measured (Figure 5). Th e local distribution of the agent at the tailend of the distribution was compared to th e mathematical model of diffusion. For the 20 minute case, the mathematical model was compar ed to the local distribution from the tailend to a distance 3 mm toward the center of th e injection from the firstzero event value. For the 2.5 minute case, the mathematical mode l was compared to the local distribution from the tailend to a distance 0.7 mm toward the center of the inje ction from a bend. The concentration profile bends from a steep cu rve to a plateau regi on. The bend was defined to be the position where the number of events had a value of 100 (Figure 5a). The plateau region was defined as having a relatively flat profile where the events values were between 0 and 100. The position of the first ze ro event value could not be used as a reference as the plateau region varied greatly in length. Hence, it woul d not be possible to set a specified measurement distance from th e first zero event value. A plateau region was not seen with the 20 minute data. Th e distribution of the agent within the subcutaneous tissue near the ta ilend of the concentration prof ile was consistent with the mathematical model of diffusion (Figure 12). The mathematical model was compared to the experimental data in order to obta in values for the diffusion coefficient, D at 2.5 and 20 minutes after injection. The diffusion coefficient, D of dexamethasone in subcutaneous tissue slices at 2.5 a nd 20 minutes after injection was 2.69 1.08 x1010 and 4.01 2.01 x1010 m2/s, respectively. As mentioned a bove, there were two criteria for equation 5 to provide a useful approximation. Also, to comply with the criteria, the data near the tailend of the concentration prof iles were used. However for a few of the concentration profiles for the t =2.5 minute case, using data 0.7 mm toward the source from the bend meant using all the data as the profile was very steep (Figure 12a). Hence, PAGE 62 55 at t =2.5 minutes the criteria could not be co mplied with. For the 20 minute case, there was an offset ranging from 0.77 to 2.25 mm from the center of the inje ction (Figure 12b). Although this offset is not large, it may be sufficient enough to comply with criterion 1. In addition, data in Nicholson ( 19 ) showed that the accur acy of equation 5, at measurement distances near the source, increases with time. Further, the criteria required that the injection time be very brief compar ed to the time of the subsequent diffusion measurements. The approximate duration of the injection was 1 second. Hence, the injection of 0.04 mL of substance was very brief. The two diffusion times were t =2.5 minutes and t =20 minutes. The t =2.5 minute concentration prof ile had a plateau region that was not seen in the t =20 minute concentration profile (F igure 12). It could be that, at t =2.5 min, the injected substance both formed a fluid filled cavity and infiltrated the extracellular space to some degree, producing the plateau region. If this were the case, then equation 5 would not be the appropriate expression. This is a phenomenon that needs to be investigated further and is beyond the scope of this present study. The concentration profile at t =20 minute is similar in shape to the theoretical curves realized by using equation 11 (Figure 4), wher eas the concentration profile for t =2.5 minutes is not similar due to the plateau region. Hence, we assume that for t =2.5 minutes that the diffusion time was not long enough and that equation 11 does not provide a useful approximation. A diffusion time of t =20 minutes probably provides an adequate diffusion time, and, hence, equa tion 11 does provide a useful a pproximation in this case. Therefore, the best estimate for the diffusion coefficient of dexamethasone in subcutaneous tissue slices based on the t= 20 minute data is D= 4.01 2.01 x1010 m2/s. This compares with D= 4.11 1.77x1010 m2/s for the implantation case. These values PAGE 63 56 for D are very similar, suggesti ng that equation 11 provides an adequate approximation as long as the two criteria are met. Our math ematical model assumed that the diffusing substance was deposited within a sphere at t =0. Figure 2a shows th at the shape of the injection at 2.5 minutes is rela tively spherical. Hence, th e assumption that the injected volume at t= 0 was spherical is acceptable. The mathematical model also assumed isotropic diffusional transport of drug thr ough the subcutaneous ti ssue. Figure 2b shows the diffusion of the drug at t =20 minutes, and diffusion is re latively spherical away from the site of injection. Hen ce, our assumption of isotr opic diffusional transport is reasonable. However, the actual shape is not perfectly spherical, and this resulted in a relatively large standard deviation. The elimination constant was assumed to be negligible as the injections were made in ei ther harvested subcutaneous tissue or in a sacrificed rat so that the nor mal clearance processes that depend on circulation of blood were eliminated ( 8 ). It is recommended that for fu ture diffusion experiments with other substances in other tissue that (1) the injection volume be sm all, i.e., 0.04 mL or less; and (2) higher radioactivity be used, i.e., 0.65 Ci or higher. The experimental duration time can be estimated a priori by (1) finding a di ffusion coefficient for a substance similar to the one of interest in the ti ssue, (2) making theoretical curv es like these shown in figure 4, and (3) choosing a curve that shows a larg e diffusion distance away from the source (e.g., the 20 minute curve in Figure 5 shows a large diffusion distance away from the source). In conclusion, equation 11 pr ovides an adequate approximation for measurement locations sufficiently far from the source and for diffusion times much longer than the injection time ( 8 19 ). The main advantages of this injection technique to determine an approximation for the diffusion coe fficient are that it is a relatively simple PAGE 64 57 technique and can be ap plied to any radiolabeled substanc e of interest injected into any tissue of interest. PAGE 65 58 References (1) Hickey, T.; Kreutzer, D.; Burgess, D. J.; Moussy, F. In vivo evaluation of a dexamethasone/PLGA microsphere system desi gned to suppress the inflammatory tissue response to implantable medical devices. J. Biomed. Mater. Res. 2002, 61 180187. (2) Wisniewski, N.; Moussy, F.; Reichert W. M. Characterization of implantable biosensor membrane biofouling. Fresenius J. Anal. Chem 2000, 366 611621. (3) Norton, L. W.; Tegnell, E.; Topore k, S. S.; Reichert, W. M. In vitro characterization of vascular endothelial gr owth factor and dexamethasone releasing hydrogels for implantable probe coatings. Biomaterials 2005, 26 32853297. (4) Strasser, J. F.; Fung, L. K.; Eller, S.; Grossman, S. A.; Saltzman, W. M. Distribution of 1,3bis(2chloro ethyl)1nitrsourea and tracer s in the rabbit brain after interstitial delivery by biode gradable polymer implants. J. Pharmacol. Exp. Ther. 1995, 275 16471655. (5) Kurisawa, M.; Chung, J. E.; Yang, Y. Y.; Gao, S. J.; Uyama, H. Injectable biodegradable hydrogels com posed of hyaluronic acidtyramine conjugates for drug delivery and tissue engineering. Chem. Commun. ( Cambridge ) 2005, 34 43124314. (6) Huang, T. S. Concomitant infusion of ovine corticotrophinreleasing hormone does not prevent suppression of the h ypothalamuspituitaryadrenal axis by dexamethasone in male rats. J. Endocrinol. In V est 1997, 20 ( 7 ), 393396. (7) Xin, Q.; Wightman, M. R. Transpor t of choline in rat brain slices. Brain Res. 1997, 776 126132. (8) Thorne, R. G.; Hrabetova, S.; Nic holson, C. Diffusion of epidermal growth factor in rat brain extrace llular space measured by in tegrative optical imaging. J. Neurophysiol. 2004, 92 34713481. (9) Zhou, X.; Pogue, B. W.; Chen, B.; Hasan, T. Analysis of effective molecular diffusion rates for Verteporfin in subcutan eous versus orthotopic Dunning prostate tumors. Photochem. Photobiol. 2004, 79 (4), 323331. (10) Lankelma, J.; Luque, R. F.; Dekker: H.; Schinkel, W.; Pinedo, H. M. A mathematical model of drug tran sport in human breast cancer. Micro V asc. Res. 2000, 59 149161. PAGE 66 59 (11) Rapkin, E. Realtime radioima ging of biological tissue sections. Am. Biotechnol. Lab 2001, 19 3436. (12) Laniece, P.; Charon, Y.; Car dona, A.; Pinot, L.; Maitrejean, S.; Mastrippolito, R.; Sandkamp, B.; Valentin, L. A new highresolution radioimager for the quantitative analysis of radiolab eled molecules in tissue section. J. Neurosci. Methods 1998, 86 15. (13) Carslaw, H. S.; Jaeger, J. C. Conduction of Heat in Solids 2nd ed.; Oxford University Press: Oxford, 1959. (14) Colton, C. K. Analysis of me mbrane processes fo r blood purification. Blood Purif. 1987, 5 202251. (15) Press: W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes in C ; Cambridge University Press: Cambridge, 1988. (16) Hersh, L.; Friedman, B.; Medero, R. Method for oscillometric blood pressure determination employing curve fitting U.S. Patent Number 5704362, 1998. (17) Crank, J. The Mathematics of Diffusion 2nd ed.; Clarendon Press: Oxford, U.K., 1975; pp 2829. (18) Deen, W. M. Analysis of Transport Phenomena ; Oxford University Press: New York, 1998; pp 187188. (19) Nicholson, C. Diffusion from an in jected volume of a substance in the brain tissue with arbitrary volu me fraction and tortuosity. Brain Res. 1985, 333 325329. (20) CRC Standard Mathematical Tables, 25th Edition, CRC Press, 1978. (21) Saltzman, W. M.; Radomsky, M. L. Drugs released from polymers: diffusion and elimination in brain tissue. Chem. Eng. Sci 1991, 46 ( 10 ), 24292444. (22) Barry, B. W.; Brace, A. R. Permeati on of oestrone, oestradiol, oestriol and dexamethasone across cellulose acetate membrane. J. Pharm. Pharmacol 1977, 29 ( 7 ), 397400. (23) Sharkawy, A. A.; Klitzman, B.; Trus key, G. A.; Reichert, W. M. Engineering the tissue which encapsulates subcutaneous implants. I. Diffusional properties. J. Biomed. Mater. Res. 1997, 37 401412. (24) Kim, T.K.; Burgess, D. J. Pharmacokinetics characterization of 14Cvascular endothelial growth factor contro lled release microspheres using rat model. J. Pharm. Pharmacol. 2002, 54 897905. PAGE 67 60 (25) Reinhard, C. S.; Radomsky, M. L. ; Saltzman, W. M.; Hilton, J.; Brem, H. Polymeric controlled release of dexamethasone in normal rat brain. J. Controlled Release 1991, 16 331340. (26) Moussy, Y.; Hersh, L.; Dungel, P. Distribution of [3H]dexamethasone in rat subcutaneous tissue after delivery from osmotic pumps. Biotechnol. Prog. 2006, 22 819824. (27) SAS Institute Statistical Analysis System. SAS Users Guide 1979 Edition. SAS Institute Inc. Cary, North Carolina, 1979. PAGE 68 61 Appendices PAGE 69 62 Appendix A: The LevenbergMarquardt Algorithm The LevenbergMarquardt method of optimization is a minimum sum squared error numerical technique for fitting a nonlin ear expression to a set of data. The optimization can be easily extended to several va riables. This makes it an ideal technique for the problem addressed in this research (27) The solution of the diffusion equation contains a set of parameters, C0, D, and k, over which one wishes to search for the best possible values so that the experimental measurements are adequately summarized. More generally, we have a function F which relates a single i ndependent variable to a single dependent variable. ) (i ix F y The independent and depe ndent variables form an ordered pair ( xi yi ). The measured data that we work with is a set, m, of these ordered pairs. F contains a set of parameters,i that must be found to best represent the data with the function. The goal becomes finding the best set {i } that minimizes the sum squared error between the function and the raw data: 2) ( . .i i iy y E S S (79) where i i iy y i2 1and ) , ( 2 1 0i ix F y The first strategy that the Marquardt technique uses is to find the gradient of the sumsquarederror with respect to PAGE 70 63 the function parameters: E X F X Y X ) ( 2 1 (80) Note that X is an m by n matrix containing the partial deri vatives of F with respect to the parameters, F X and E is an n by 1 matrix containing the error at each data point in the set. The gradient is used to help determine the direction to move in the i space to achieve the smallest sumsquarederror: E X ki i 1. (81) The variable k controls how fa r to move in the direction opposite to the gradient in updating the parameter values. The disadvantag e of the gradient method is that while it can tell the best direc tion to move for finding a solution, it does not specify how far to move. The Marquardt algorithm rectifie s this problem by using a GaussNewton technique. The Gauss technique assumes that the f unction of interest can be expanded in a Taylor series around the present location in the i space. Only the linear terms are kept in this approach: 0 0 X F F (82) Now, by assuming there is an exact i location that will make the above equation exact or, in other words, when the error is at zero, the following algebra leads to an update formula. Y X X F X 0 0 (83) PAGE 71 64 0 0 F X Y X X X (84) So that finally, the Gauss method gives an update formula: E X X Xi i 1 1 (85) Note that this update formula te lls us everything we need to ge t to a better solution in the i space. However, as opposed to the st eepest decent gradient method it requires a matrix inverse calculation which can be unstable under certain conditions. The Marquardt algorithm combines the philosophy of these two approaches to provide a stable yet complete way to move to a better position in the solution space. Marquardt uses the following update formula: E X I X Xi i 1 1 (86) This formula balances the gradientsteep estdecent and the Gauss approaches by essentially forcing the two together by use of the scaling parameter. Note that if goes to zero the formula become s the Gauss expression and if becomes large the formula goes to the steepest decent expr ession with very small movement. So by controlling we can keep the inverse calculation stable and yet still move in a good direction even when the inverse is not stable. The key then becomes picking a good st arting point in the solution space. Additionally, developing a good strategy for updating is important. The Marquardt algorithm handles this issue by increasing the if the error gets larg er on that particular iteration of the updating loop and decreasing if the error got better. In this way, the distance moved gets more aggressive as the e rror gets smaller but us es the Gauss concept to know precisely how far to move. This t echnique is then able to provide a stable, efficient, easily programmable, easily understo od, and mathematically elegant approach PAGE 72 65 for finding the optimal solution for a nonlinear function containing unknown parameters. (Figure 13.) PAGE 73 66 Figure 13. Flow Chart Show ing the Marquardt Algorithm. This diagram is a detailed flowchart for the Marduardt Algorit hm used for the curve fitting. Start Set the initial point at which to start search. Set the error limit (epsilon) and the initial lambda. Calculate the initial sum(SQE) squared error. Calculate the gradient. Calcualte the least squares matrix. Is RSQE PAGE 74 67 Appendix B: The Marquardt Program % % This program will fit a curve to autoradiogrphic data to determine % diffusion and elimination constants. % % The initializations. % clear all; format short e; format compact; warning off all % % Get or generate the data needed. % The number of data points should be specified as n. % This data should have the following form for each sample row: % 1. First Column; distance from source sample. % 2. Second Column; time of the sample. % 3. Third Column; concentration value of the sample. % % R is the distance from the source. % % % k is the elimination constant. % % % D is the diffusion constant. % % % C is concentration level. C0 is the source concentration. % % % Thin line 45 degrees background already subtracted. % a=0.0006; % a is the radius of the source. % InData=[ 1.47336 40.692426 1.4944 32.927726 1.51545 35.045326 1.5365 35.045326 1.55755 23.751226 1.5786 29.398226 1.59964 29.398226 1.62069 42.104126 1.64174 26.574726 1.66279 26.574726 1.68384 23.751226 1.70488 23.751226 1.72593 21.633526 1.74698 21.633526 1.76803 23.045326 1.78908 28.692426 1.81012 28.692426 1.83117 26.574726 PAGE 75 68 1.85222 20.927726 1.87327 23.751226 1.89432 23.751226 1.91536 14.574726 1.93641 17.398226 1.95746 17.398226 1.97851 14.574726 1.99955 23.045326 2.0206 17.398226 2.04165 17.398226 2.0627 23.045326 2.08375 13.868826 2.10479 13.868826 2.12584 6.104126 2.14689 12.457126 2.16794 12.457126 2.18899 10.339426 2.21003 8.927726 2.23108 19.515926 2.25213 19.515926 2.27318 18.810026 2.29423 24.457126 2.31527 24.457126 2.33632 28.692426 2.35737 12.457126 2.37842 21.633526 2.39947 21.633526 2.42051 19.515926 2.44156 20.221826 2.46261 20.221826 2.48366 13.868826 2.50471 11.751226 2.52575 3.280626 2.5468 3.280626 2.56785 15.986526 2.5889 7.515926 2.60995 7.515926 2.63099 21.633526 2.65204 13.868826 2.67309 13.868826 2.69414 11.045326 2.71518 12.457126 2.73623 5.398226 2.75728 5.398226 2.77833 11.045326 2.79938 6.810026 2.82042 6.810026 2.84147 12.457126 2.86252 10.339426 2.88357 8.927726 2.90462 8.927726 2.92566 16.692426 2.94671 11.045326 2.96776 11.045326 2.98881 3.280626 3.00986 1.163026 3.0309 9.633526 3.05195 9.633526 3.073 6.104126 3.09405 3.986526 3.1151 3.986526 3.13614 3.280626 PAGE 76 69 3.15719 6.810026 3.17824 6.810026 3.19929 1.868826 3.22034 1.163026 3.24138 0.457126 3.26243 0.457126 3.28348 10.339426 3.30453 6.810026 3.32558 6.810026 3.34662 10.339426 3.36767 10.339426 3.38872 0.457126 3.40977 0.457126 3.43082 8.927726 3.45186 3.778174 3.47291 3.778174 3.49396 0.248774 3.51501 5.398226 3.53605 5.189974 3.5571 5.189974 3.57815 11.751226 3.5992 13.163026 3.62025 13.163026 3.64129 5.398226 3.66234 8.927726 3.68339 8.927726 3.70444 6.104126 3.72549 8.013474 3.74653 0.248774 3.76758 0.248774 3.78863 3.280626 3.80968 4.692426 3.83073 4.692426 3.85177 3.986526 3.87282 3.280626 3.89387 0.954674 3.91492 0.954674 3.93597 1.868826 3.95701 5.398226 3.97806 5.398226 3.99911 4.692426 4.02016 1.163026 4.04121 6.601774 4.06225 6.601774 4.0833 1.868826 4.10435 5.398226 4.1254 5.398226 4.14645 2.574726 4.16749 0.248774 4.18854 0.248774 4.20959 0.248774 4.23064 1.868826 4.25168 8.221826 4.27273 8.221826 4.29378 0.954674 4.31483 6.810026 4.33588 6.810026 4.35692 5.398226 4.37797 1.163026 4.39902 5.398226 4.42007 5.398226 4.44112 8.221826 PAGE 77 70 4.46216 4.692426 4.48321 4.692426 4.50426 6.810026 4.52531 1.868826 4.54636 2.574726 4.5674 2.574726 4.58845 4.692426 4.6095 3.072274 4.63055 3.072274 4.6516 11.045326 4.67264 13.163026 4.69369 13.163026 4.71474 3.986526 4.73579 8.221826 4.75684 6.601774 4.77788 6.601774 4.79893 4.484074 4.81998 5.895874 4.84103 5.895874 4.86208 3.280626 4.88312 0.457126 4.90417 0.457126 4.92522 0.457126 4.94627 0.248774 4.96732 0.954674 4.98836 0.954674 5.00941 0.457126 5.03046 5.398226 5.05151 0.248774 5.07255 0.248774 5.0936 0.954674 5.11465 0.954674 5.1357 0.954674 5.15675 3.778174 5.17779 3.072274 5.19884 3.072274 5.21989 0.457126 5.24094 8.013474 5.26199 6.810026 5.28303 6.810026 5.30408 3.986526 5.32513 1.163026 ]; [n,m]=size(InData); R=zeros(1,n); C=zeros(1,n); for i=1:n T(i)=24*60*60; % T is the time. R(i)=(InData(i,1)InData(1,1))*0.001+a; C(i)=InData(i,2); % if C(i)<0.0 % C(i)=0.0; % end end % % The following creates the initial estimates for the optimization starting % point. % estC0=max(C); PAGE 78 71 C0=estC0; k=4.5825e5; D=6.82e10; Dold=1.0e10; % % figure(1); plot(Ra,C,'r+'); grid on; title('Concentration vs. Distance from Source'); xlabel('Position (meters)'); ylabel('Concentration'); % % Creation of the necessary matrices for the Marquardt algorithm. % X=zeros(n,3); XT=zeros(3,n); XTX=zeros(3,3); INV_XTX=zeros(3,3); DELTA=zeros(n,1); ADJUST=zeros(3,1); INIT=zeros(3,1); NEW=zeros(3,1); Cest=zeros(1,n); LAMBDA=eye(3); % % The following are the initial estim ates with other ne eded settings. % epsilon=0.0001; % Small numbers needed for the curve fit. lambda=100000000000000; oldSSQE=0.0; % Storage for the sum squared error and its update. newSSQE=0.0; INIT(1,1)=C0; INIT(2,1)=k; INIT(3,1)=D; % % Calculate the sum squared error of the concentrtion with the latest % concentration model. % iteration=0; newSSQE=0.0; for i=1:n Cest(i)=C0*a/(2*R(i))*(exp((R(i)a)*sqrt(k /D))*erfc((R(i)a)/(2*sqrt(D*T(i)))sqrt(k*T(i))) ... +exp((R(i)a)*sqrt(k/D))*erfc((R(i)a)/(2*sqrt(D*T(i)))+sqrt(k*T(i)))); newSSQE = newSSQE + (Cest(i)C(i))^2; end % % Iterate until convergence, but loop at least 5 times. while (abs(newSSQEoldSSQE)/oldSSQE>epsilon) &(iteration<100)&(abs(DDold)/D>epsilon) newSSQE oldSSQE=newSSQE; Dold=D; iteration=iteration+1 % % Fill the X matrix and do the matix calculations. % The derivatives have been split up for efficiency and clarity. % for i=1:n fparm1=(R(i)a)/(2*sqrt(D*T(i)))sqrt(k*T(i)); fparm2=(R(i)a)/(2*sqrt(D*T(i)))+sqrt(k*T(i)); fparm3=(R(i)a)*sqrt(k/D); fparm4=0.5*sqrt(T(i)/k); PAGE 79 72 fparm5=0.5*sqrt(k/(D^3)); fparm6=(R(i)a)*0.5/sqrt(k*D); fparm7=a/(2*R(i)); X(i,1)=fparm7*(exp(fparm3)*erfc(fparm1) ... +exp(fparm3)*erfc(fparm2)); X(i,2)=C0*fparm7*((fparm6)*exp(fparm3)*erfc(fparm1) ... +(2.0/sqrt(pi)*exp((fparm1^2)))*(fparm4)*exp(fparm3) ... +(fparm6)*exp(fparm3)*erfc(fparm2) ... +(2.0/sqrt(pi)*exp((fparm2^2)))*fparm4*exp(fparm3)); X(i,3)=C0*fparm7*(((R(i)a )*fparm5)*exp(fparm3)*erfc(fparm1) ... +(2.0/sqrt(pi)*exp((fparm1^2)))*((R(i)a)/(4*sqrt(D^3*T(i))))*exp(fparm3) ... +((R(i)a)*fparm5)*exp(fparm3)*erfc(fparm2) ... +(2.0/sqrt(pi)*exp ((fparm2^2)))*((R(i)a)/(4*sqrt(D^3*T(i))))*exp(fparm3)); DELTA(i,1)=C(i)C0*fparm7*(exp(fparm3)*erfc(fparm1) ... +exp(fparm3)*erfc(fparm2)); end % % The essence of the Marquardt algorithm. % XT=X'; % LAMBDA(1,1)=lambda; LAMBDA(2,2)=lambda; LAMBDA(3,3)=lambda; % XTX = XT*X + LAMBDA; INV_XTX=inv(XTX); % ADJUST=INV_XTX*XT*DELTA % NEW=ADJUST+INIT % INIT=NEW; % C0=INIT(1,1); k=INIT(2,1); D=INIT(3,1); % newSSQE=0.0; for i=1:n Cest(i)=C0*a/(2*R(i))*(exp((R(i)a)*sqrt (k/D))*erfc((R(i)a)/(2*sqrt(D*T(i)))sqrt(k*T(i))) ... +exp((R(i)a)*sqrt(k/D))*erfc((R(i)a)/(2*sqrt(D*T(i)))+sqrt(k*T(i)))); newSSQE = newSSQE + (Cest(i)C(i))^2; end % if oldSSQE>newSSQE lambda=lambda/2; else lambda=lambda*2; end lambda end % % Create the found curve for plotting with the data. % Rsg=zeros(1,300); Cestg=zeros(1,300); Restg=zeros(1,300); kn=4.5825e5 Dn=6.82e10 k D PAGE 80 73 % for i=1:300 Rsg(i)=(i1)*0.000015+a; Cestg(i)=a/(2*Rsg(i))*(exp((Rsg(i)a)*sqrt(kn/ Dn))*erfc((Rsg(i)a)/(2*sqrt(Dn*T(1)))sqrt(kn*T(1))) ... +exp((Rsg(i)a)*sqrt(kn/Dn))*erfc((Rsg(i )a)/(2*sqrt(Dn*T(1)))+sqrt(kn*T(1)))); Restg(i)=a/(2*Rsg(i))*(exp((Rsg(i)a)*sqrt(k /D))*erfc((Rsg(i)a)/(2*sqrt(D*T(1)))sqrt(k*T(1))) ... +exp((Rsg(i)a)*sqrt(k/D))*erfc((Rsg(i)a)/(2*sqrt(D*T(1)))+sqrt(k*T(1)))); end % for i=1:300 if Restg(i)<0.1 r10=Rsg(i)a break; end end C=C/C0; % figure(2); %plot((Ra)*1000,C,'b*',(Rsga)*1000,Cestg,(Rsga)*1000,Restg); plot((Ra)*1000,C,'b*',(Rsga)*1000,Restg); grid on; title('Concentration vs. Distance from Source, Thinline 24Hr. 45 Degrees'); xlabel('Position (mm)'); ylabel('Concentration/Counts'); text(3.0,0.6,... ['D = ',num2s tr(D)],'FontWeight','bold'); text(3.0,0.8,... ['k = ',num2str(k)],'FontWeight','bold'); % text(2.8e03,C(1),... % ['Dn = ',num2str(Dn)],'FontWeight','bold'); % text(2.8e03,C(25),... % ['kn = ',num2str(kn)],'FontWeight','bold'); % text(4.8e03,C(30),... % ['r10 = ',num2str(r10)],'FontWeight','bold'); % k/D 