Ind. Eng. Chem. Res. 2007, 46, 6157-6163
6157
Measurements of the Effective Diffusion Coefficient of Oxygen in Pancreatic Islets Efstathios S. Avgoustiniatos,†,‡ Keith E. Dionne,†,§ David F. Wilson,# Martin L. Yarmush,†,⊥ and Clark K. Colton*,† Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, and Department of Biochemistry and Biophysics, UniVersity of PennsylVania, Philadelphia, PennsylVania 19104
The effective diffusion coefficient of oxygen in pancreatic islets is an important parameter that is needed for the development of models of viability and function of pancreatic islets that are transplanted to cure diabetes, especially when encapsulated in devices. We have developed a method for estimating the maximum oxygen consumption rate (Vmax) and the oxygen permeability in tissue (RtDt) in a suspension of spherical aggregates of different sizes from measurements of the bulk oxygen partial pressure in batch oxygen consumption experiments. The method made use of a model of oxygen consumption and diffusion in the spheres and a parameter estimation scheme. Measurements made with three different batches of rat islets yielded batch averages of Vmax ) 3.40 ( 0.77 × 10-8 mol/(cm3 s) and RtDt ) 1.34 ( 0.36 × 10-14 mol/(cm mm Hg s). Using a literature estimate for the solubility (Rt) yields Dt ) 1.31 ( 0.36 × 10-5 cm2/s at 37 °C. These oxygen permeability and diffusion coefficient values are ∼38% and ∼47% of their values in water, respectively. Introduction Islets of Langerhans are spheroidal aggregates of cells that are located in the pancreas and secrete hormones that are involved in glucose metabolism, most notably, insulin. Transplantation of isolated islets is a promising treatment for some forms of Type 1 diabetes.1-3 Islets removed from the pancreas lose their internal vascularization and are dependent on the diffusion of oxygen from the external environment and through the oxygen-consuming islet tissue to satisfy the metabolic requirements of all the cells. The supply of oxygen is a limiting factor for the viability and functionality of islets that are (a) encapsulated or placed in devices for implantation,4-6 (b) cultured, especially at high density,7 and (c) generally used in situations where the islets are exposed to low oxygen conditions.8,9 The development of theoretical models of oxygen diffusion and consumption to describe these phenomena requires estimates of the oxygen consumption rate, solubility, and effective diffusion coefficient in islet tissue. No measurements of the oxygen diffusion coefficient have been reported for islet tissue, and estimates for other tissues vary widely.10-18 In this study, we have developed a method for estimating the oxygen diffusion coefficient in spheroidal cellular aggregates using experimental data from batch oxygen uptake experiments and applied the method to rat islet preparations. In these experiments, the bulk oxygen concentration in a well-stirred microreactor containing a suspension of islets was measured by monitoring the quenching of the phosphorescence of Pdcoproporphyrin by oxygen.19 We used a model that incorporated oxygen consumption and diffusion in islets and allowed us to * To whom correspondence should be addressed. Tel.: 617-2534585. Fax: 617-252-1651. E-mail address:
[email protected]. † Department of Chemical Engineering, Massachusetts Institute of Technology. ‡ Currently with Diabetes Institute for Immunology and Transplantation, Department of Surgery, University of Minnesota, Minneapolis, MN. § Currently with Alantos Pharmaceuticals, Cambridge, MA. # Department of Biochemistry and Biophysics, University of Pennsylvania. ⊥ Currently with Department of Biomedical Engineering, Rutgers University, Piscataway, NJ.
extract kinetic and transport parameters by fitting theoretical prediction to experimental data. The islets were present in a distribution of sizes that complicated the analysis because intraislet oxygen profiles and the mass transfer coefficient between the bulk and islet surface (evaluated from an empirical correlation) varied with islet size. Thus, it was necessary to simultaneously solve for oxygen profiles in each group of islets of about the same size in an iterative optimization scheme to find the values of two adjustable parameters, the maximum oxygen consumption rate and the effective oxygen diffusion coefficient. Experimental Procedures Islet Isolation, Culture, and Sizing. Islets were isolated from male Sprague-Dawley rats, using a modified collagenase digestion/ficoll purification technique,8,20-22 handpicked under a dissecting microscope, and cultured for 1 day in non-attaching polystyrene Petri dishes (Catalog No. 8-757-12, Fisher Scientific Co., Pittsburgh, PA) in Dulbecco’s modified Eagles medium (DMEM, D5523, Sigma Chemical Co.), which contained 5.6 mM (100 mg/dL) glucose, 50 U/mL penicillin, 50 µg/mL streptomycin, and 10% (v/v) newborn calf serum (NBCS, 2006010, GIBCO). The material was placed in an incubator at a temperature of 37 °C with a mixture of 95% air and 5% CO2. On the day of oxygen uptake measurements, aliquots that contained ∼50-100 islets were photographed under a dissecting microscope for later sizing using a calibrated reticule. The size distribution and total tissue volume of the islets were estimated by individual measurements of the external dimensions of each islet. Deviations from a spherical shape were taken into consideration by assuming that each islet could be represented by an ellipsoid, measuring both the largest islet diameter (a) and the diameter normal to it (b), and estimating the third dimension (c) from the ratio c/b ) 0.87 (found in a parallel study22). The islet volume was calculated from V ) πabc/6. The equivalent spherical diameter (d) was then given by d ) (abc)1/3. Oxygen Consumption. A known number of islets were placed in a tube that contained ∼45 mL of culture medium and ∼5 mL of air. The tubes were intermittently rotated to prevent
10.1021/ie070662y CCC: $37.00 © 2007 American Chemical Society Published on Web 08/18/2007
6158
Ind. Eng. Chem. Res., Vol. 46, No. 19, 2007
settling and aggregation of the islets, and to enhance oxygen transfer; they were carried by hand and transferred by commercial air flight from MIT to the University of Pennsylvania (total transit time of ∼4 h), where they were incubated at a temperature of ∼35 °C for 2 h before testing. The oxygen uptake chamber was equipped for measuring the oxygen-dependent lifetime of Pd-coproporphyrin phosphorescence to provide rapid and accurate measurements of oxygen concentration down to values as low as 0.05 µM.19,23 The chamber was a glass cuvette (0.275 mL volume, square inside cross-section, ∼0.5 cm wide on each side) that contained a small Teflon-coated magnetic stirring bar (0.4 cm in length, 0.1 cm in diameter) rotated at a speed of 20 rev/s. An aliquot that contained ∼1500 islets was loaded in the cuvette, which was then completely filled with phosphate-buffered saline (pH 7.4) containing 0.35 g/L Hepes buffer, 0.5 g/L bovine serum albumin, and 300 mg/dL glucose supplemented with 0.01 µM palladium coproporphyrin and 1-5 U/mL catalase. The cuvette was capped with a ground glass stopper to eliminate the gas phase. The stopper had a 2-cm-long, 1-mm-diameter hole that was bored through its center; this hole functioned as an unstirred oxygen diffusion barrier. To repeat an experiment after the oxygen was depleted, H2O2 was immediately injected into the cuvette through the hole without removing the cap to produce dissolved oxygen and water and thus reoxygenate the medium to between 80 and 100 mm Hg. The temperature of the cuvette was controlled by wrapping the aluminum block of the cuvette holder with heating tape that was connected to a variable voltage source. The temperature was monitored continuously with a thermistor and was maintained at ∼37 °C throughout the oxygen uptake experiment. In two of the runs, the oxygen uptake was uncoupled from ATP generation via the addition of carbonylcyanide-p-trifluoro-methoxyphenylhydrazone (FCCP, E.I. DuPont de Nemours Co., Wilmington, DE) to a final concentration of 15 µM and inhibited by the addition of amytal (amobarbital sodium, Eli Lilly, Indianapolis, IN). Additional details are available in ref 22. Extraction of Kinetic and Diffusion Parameters Oxygen Reaction and Diffusion Model. We assume that the islet preparation is a suspension of tissue spheres that can be divided into m groups; each sphere in group i has the same equivalent radius Ri, which differs from group to group. We assume that the tissue is uniform with constant physical properties, independent of position. Because oxygen is poorly soluble in aqueous media, the species conservation equation for transient oxygen diffusion and reaction in spherical coordinates with azimuthal symmetry is
() ( )
∂c 1 ∂ 2 ∂c ) Dt 2 r -V ∂t ∂r r ∂r
(1)
where c is the oxygen concentration in the tissue (in units of mol/cm3), t the time (in seconds), Dt the effective diffusion coefficient of oxygen in the tissue (in units of cm2/s), r the radial distance from the center of the sphere (in centimeters), and V the intrinsic oxygen consumption rate per unit volume of tissue (in terms of mol/(cm3 s)). A simplified but reasonable representation of the dependence of the oxygen consumption rate on the oxygen concentration is given by Michaelis-Menten kinetics:
V)
Vmaxc K′m + c
(2)
where Vmax is the maximum oxygen consumption rate and K′m is the Michaelis-Menten constant (in terms of mol/cm3). Because partial pressure is continuous across an interface, it is convenient to express eq 1 in terms of the oxygen partial pressure (p, in units of mm Hg), using
c ) Rtp
(3)
where Rt is the Bunsen solubility coefficient of oxygen in the tissue (in units of mol/(cm3 mm Hg)). Substituting eqs 2 and 3 into eq 1 leads to
Rt
() ( )
Vmaxp ∂p 1 ∂ 2 ∂p r ) R t Dt 2 ∂t ∂r ∂r K r m+p
(4)
where the product RtDt (in terms of mol/(cm mm Hg s)) is the oxygen permeability in the tissue and Km is expressed in units of mm Hg. The initial condition is
p ) p(r,0)
(for t ) 0)
(5)
The boundary conditions are symmetry about the center of the sphere,
∂p )0 ∂r
(at r ) 0)
(6)
and the requirement at the surface that the oxygen diffusive flux within the sphere be equal to the flux Ji for convective transport across the mass transfer boundary layer surrounding each sphere:
R t Dt
∂p ) kiRm[pm - p(Ri)] ) Ji ∂r
(at r ) Ri)
(7)
where p(Ri) is the oxygen partial pressure at the surface of a sphere of radius Ri, ki the oxygen mass transfer coefficient between the bulk medium and the surface of the sphere (in units of cm/s), and Rm the oxygen solubility in the medium. Finally, an expression is needed to relate the total rate of oxygen transfer (N) from the bulk medium to all of the spheres:
N ) -VmRm
dpm dt
m
)
Ji(4πRi2)ns fi ∑ i)1
(8)
where Vm is the volume of the medium (in units of cm3), ns the total number of spheres, and fi the number fraction of spheres in group i. The associated initial condition is
pm ) pm(t0)
(for t ) t0)
(9)
Known Parameters in the Model. The Bunsen solubility coefficient in the tissue at 37 °C was estimated by averaging values for five different tissues to yield Rt ) 1.02 × 10-9 mol/ (cm3 mm Hg); the coefficient for the medium at 37 °C is Rm ) 1.27 × 10-9 mol/(cm3 mm Hg).16,22,24 Measurements of oxygen consumption rates with isolated rat liver mitochondria at 25 °C yielded an estimate of K′m ) 0.7 µM,23 which corresponds to Km ) 0.44 mm Hg from eq 3 with Rm ) 1.60 × 10-9 mol/(cm3 mm Hg) at 25 °C.25 We assumed that Km has little temperature dependence and used this value at 37 °C. The mass transfer coefficient (k) for spherical particles in an agitated tank was estimated from correlations of the Sherwood number Sh with the agitation power per unit fluid mass.26 For particles with diameter di in the range of 100-300
Ind. Eng. Chem. Res., Vol. 46, No. 19, 2007 6159
µm (that of islets used in this study) and for little difference in the densities of the fluid medium and the suspended particles,26,27
() ( )
ki di dI ) Shi ) 2 + Sc1/3 Dm dT
di
f
(10)
ν3
k i di ) Shi ) 2 + 101di0.65 Dm
(11)
where the diameter of the sphere di is measured in centimeters and the coefficient 101 has units of cm-0.65. The first term on the right-hand side of eq 11 is the value of Sh in an unstirred medium, whereas the second is the increase in Sh due to stirring. The value of Sh ranged from ∼7 to ∼13 for islet diameters in the range of di ) 100-300 µm. The oxygen partial pressure drop across the mass transfer boundary layer may be estimated from a balance equating the oxygen consumption rate within an individual sphere of radius Ri ) di/2 to the rate of convective mass transfer of oxygen at the surface of the sphere. This is a good approximation when the transient term in eq 4 is small, as will be discussed later. For sufficiently high pm values, such that p(0,t) . Km and V ) Vmax throughout the sphere,
kiRm[pm - p(Ri)](4πRi2) ) Vmax
(4π3 R ) 3
i
(12)
and
Vmax Ri kiRm
VmRm
4
0.15
where Dm is the oxygen diffusion coefficient in isotonic saline at 37 °C (Dm ) 2.78 × 10-5 cm2/s);18,22 ν is the kinematic viscosity of the medium at 37 °C (ν ) 6.89 × 10-3 cm2/s);22 Sc is the Schmidt number (Sc ) ν/Dm ) 248); dI and dT are the impeller and tank (cuvette) diameters, respectively; is the power input per unit fluid mass; and f is a function derived from experimental data.27 The power input was estimated from correlations of the impeller power number, impeller Reynolds number, and geometrical parameters for systems similar to the cuvette and stirring bar used here.28 Details of these calculations are provided elsewhere.7 From eq 10, a simplified approximate result for the conditions of the oxygen uptake experiments in this study is
pm - p(Ri) )
eq 8 can be replaced by
(13)
Using a typical value of Vmax ) 3 × 10-8 mol/(cm3 s) and a ki value estimated from eq 11, the partial pressure drop calculated from eq 13 is ∼10 mm Hg for d ) 300 µm and ∼2 mm Hg for d ) 100 µm. Model Solution and Estimation of Parameters. Order-ofmagnitude scaling analysis can identify the requirements for dropping the transient term on the left side of eq 4. By equating the transient term to the diffusion term, the time scale for relaxation of diffusion gradients within the sphere is tD ≈ R2/Dt. Using Dt ) 1.3 × 10-5 cm2/s (determined in this study) and the largest sphere radius (R ) 1.5 × 10-2 cm), tD ≈ 17 s, which is small compared to the experimental time scale (tE) of ∼200 s for essentially complete depletion of oxygen in our experiments. The parameter tE can be estimated in a more general way by noting that the rate of change of pm is determined by the oxygen consumption reaction in the spheres. Assuming the transient term in eq 4 is small, and under conditions where V ≈ Vmax throughout all spheres,
dpm ≈ VtVmax dt
(14)
where Vt is the total volume of tissue. The time scale of the experiment is then given by
tE ≈
VmRmpm(0) VtVmax
(15)
The transient term in eq 4 can be neglected if tD , tE, which leads to
Vt ,1 DtRmpm(0) Vm R2Vmax
(16)
In this study, the order of magnitude of R2Vmax/DtRmpm(0) is 1 and that of Vt/Vm is 10-2, and eq 16 is satisfied. Consequently, the transient term in eq 4 and the initial condition eq 5 may be dropped, quasi-steady-state conditions prevail in the spheres, and all partial derivatives become ordinary derivatives. Because of the nonlinearity of Michaelis-Menten kinetics, eqs 4-8 must be solved numerically. The solution occurs within two nested iteration loops, which are contained within an optimization procedure. The first step is to establish a value for the bulk oxygen partial pressure pm(t0) at the selected initial time t0. Because the parameter estimation procedure was sensitive to scatter in the experimental data, the data (∼50 measurements) in the vicinity of t0 were smoothed by fitting them with a fourth- or higher-order polynomial, and pm(t0) was evaluated from that polynomial. For the ith size group, a value for p(Ri) was selected, eq 4 was solved by Euler’s method, dp/ dr at Ri was evaluated, and eq 7 used to establish a new estimate of p(Ri). This iterative procedure was repeated until the solution for p(Ri) converged. The procedure was conducted with each of the m size groups to establish the quasi-steady-state partial pressure field in the spheres in each group at t0. Each flux Ji was calculated from eq 7 and used in eq 8 to estimate the total oxygen transfer rate N at t0. The value of pm at the end of the first time step t1 is estimated from the integrated form of eq 8, which is expressed generally as
pm(tj) - pm(tj-1) ) )-
∫t t
j
j-1
N dt VmRm
N h (t - tj-1) VmRm j
(17)
Using N evaluated at t0, pm(t1) at the end of the first time step is estimated. Using this estimated value of pm(t1), the same calculations performed at t0 are repeated at t1, and the oxygen transfer rate, N(t1) from the bulk to the spheres at the end of the first time step t1 is estimated. The average transfer rate during the first time step, N h ) [N(t0) + N(t1)]/2 is used to re-evaluate pm(t1) from eq 17. This second iteration loop is repeated until successive estimates of pm(t1) converge. The calculation procedure previously described is a double iteration loop, with the inner loop involving the calculation of p(Ri) for all islet sizes and the outer loop involving the calculation of pm as a function of time t. By successively repeating this procedure, the bulk oxygen partial pressure time profile for the entire duration of the experiment was predicted as it decreased from its initial value to almost zero. Because the values calculated from the theoretical model were spaced
6160
Ind. Eng. Chem. Res., Vol. 46, No. 19, 2007
equally in time (usually 1 s apart), the theoretical predictions at the exact times at which experimental measurements were taken were obtained by interpolation. Although the scaling analysis indicated that a quasi-steadystate approximation was valid, we investigated whether the very small error that was incurred at each time step amassed to a larger error by solving the complete transient problem using the same double iteration scheme applied to the quasi-steadystate model. The only difference was that the quasi-steady-state solution was used as an initial estimate of the initial condition p(r,0) at t0, and the solution at t1 was obtained by solving eq 4 with a finite difference method. Values of the flux Ji at t0 and t1 were averaged to produce the next estimate of p(Ri), and the procedure was updated until the average Ji converged. To determine the set of parameters that resulted in the closest match between theoretical prediction and experimental data, we defined an objective function s to be minimized: nd
s)
wj(pm,j - p*m, j)2 ∑ j)1
(18)
where nd is the number of experimental data points, pm,j the measured value of pm at time tj, p*m,j the (interpolated) predicted value, and wj the weighting coefficient:
wj )
Figure 1. Frequency distribution of equivalent islet diameter for batch 3; the upper panel shows number fraction data, and the lower panel shows volume fraction data.
|∆pm| nd
(19)
|∆pm| ∑ j)1 where ∆pm ) pm,j - pm,j-1. This definition of wj was used to give greater weight to experimental data points that were associated with larger changes in pm between successive measurements. A modified Newton’s method was used for minimization of the objective function. The optimization procedure constituted a third iteration loop in which the previously described dual iteration loops were performed in each cycle. For the quasi-steady-state situation, the parameters Rt and Dt always appeared together. Consequently, the best fit estimates of Vmax and the product RtDt were obtained. Results Measurements were made with three batches of rat islets. Two experiments were conducted with batch 1, a single experiment with batch 2, and five experiments with batch 3. Two of the experiments involved the addition of uncoupling and inhibitive agents of oxidative phosphorylation,22 and, therefore, the bestfit Vmax values from those two experiments were not comparable to the other six and are not reported here. Size Distribution. The frequency distribution of islet diameter for batch 3 is plotted in Figure 1, in terms of number and volume fractions. For this batch, equivalent diameters of individual islets ranged from 105 µm to 295 µm. Similar distributions were obtained for batches 1 and 2. For the combined data from all three batches, the number-mean and volume-mean diameters were 173 ( 47 and 214 ( 55 µm, respectively. For modeling purposes, each of the seven islet size groups was represented by a characteristic diameter estimated independently for each islet batch using the volume average of the diameters of individual islets in each size group. Other methods of averaging yielded very similar characteristic diameters for each islet group. Oxygen Consumption Experiments. The results of a typical oxygen consumption experiment are shown in Figure 2. The
Figure 2. Bulk oxygen partial pressure in the medium as a function of time for run 3A; the solid line is predicted by the model using parameter values that gave the best agreement with experimental data.
bulk oxygen partial pressure is plotted as a function of time for the experiment with islets from batch 3, run A. The solid line is predicted by the model using the parameter values that gave the best agreement with the experimental data: Vmax ) 2.30 ( 0.01 × 10-8 mol/(cm3 s) and RtDt ) 1.49 ( 0.04 × 10-14 mol/ (cm mm Hg s). In this experiment, experimental data were recorded every 1.05 s, on average, and every tenth experimental data point was plotted in Figure 2. Theoretical prediction was generated for 1-s time intervals and matched to experimental data by interpolation of three consecutive prediction data points, using a second-order polynomial. The overall agreement between model prediction and experimental data was excellent, although there were regions displaying small systematic deviations in one direction or another. The results from individual experiments are summarized in Table 1. Vmax and RtDt values ranged from 2.30 × 10-8 to 4.29 × 10-8 mol/(cm3 s) and 0.86 × 10-14 to 1.83 × 10-14 mol/(cm mm Hg s), respectively. There were variations between batches. In the five measurements with batch 3, there was no consistent
Ind. Eng. Chem. Res., Vol. 46, No. 19, 2007 6161 Table 1. Best-Fit Values of Vmax and rtDt at 37 °C run
maximum oxygen consumption rate, Vmaxa (10-8 mol/(cm3 s))
oxygen permeability in tissue, RtDta (10-14 mol/(cm mm Hg s))
batch 1 1A 1Bb average
3.03 (n ) 1)
1.66 ( 0.03 1.83 ( 0.04 1.75 ( 0.12 (n ) 2)
batch 2
4.29 ( 0.02
1.23 ( 0.01
2.30 ( 0.00 2.92 ( 0.01 3.46 ( 0.02 2.83 ( 0.01 2.88 ( 0.47 (n ) 4)
1.49 ( 0.04 1.07 ( 0.02 0.90 ( 0.02 0.86 ( 0.01 0.88 ( 0.01 1.04 ( 0.26 (n ) 5)
3.14 ( 0.67 (n ) 6) 3.40 ( 0.77 (n ) 3)
1.24 ( 0.38 (n ) 8) 1.34 ( 0.36 (n ) 3)
batch 3 3A 3B 3C 3D 3Eb average average (by run) average (by batch)
3.03 ( 0.01
a Values tabulated are the mean ( standard error of the estimate for individual runs and the mean ( standard deviation for individual batch and overall averages (by run or by batch). b Uncoupled and inhibited runs. The Vmax value is not comparable to that of other runs.
change in Vmax, but there was a decrease in RtDt for successive measurements. The overall averages for the three batches were Vmax ) 3.40 ( 0.77 × 10-8 mol/(cm3 s) and RtDt ) 1.34 ( 0.36 × 10-14 mol/(cm mm Hg s). Using Rt ) 1.02 × 10-9 mol/(cm3 mm Hg), the batch-average effective diffusion coefficient for oxygen in the islet tissue was Dt ) 1.31 ( 0.36 × 10-5 cm2/s at 37 °C. When averaged for all runs, the best-fit values were similar and are reported in Table 1. We solved the complete transient problem for the experiments in Table 1. In seven out of nine experiments, Vmax increased by less than ∼2% and RtDt decreased usually by less than 1% and 4%-5% in two cases. Thus, for experiments in this study, the quasi-steady-state solution in the spheres was adequate. We also examined the effect of increasing the ratio Vt/Vm by generating artificial experimental data for spheres of uniform size using the full transient solution and fitting parameters with the quasisteady-state model. The error in the fitted values increased as R and Vt/Vm increased. For R ) 150 µm and Vt/Vm ) 0.2, Vmax was ∼12% too low and RtDt was ∼10% too high. Thus, in experiments with high Vt/Vm, the complete transient model is needed. Use of the transient solution required significantly more computing resources than the use of the quasi-steady-state solution. The fitted parameters were not very sensitive to the estimate used for Km. Increasing Km by 50% and decreasing it by 90% caused an increase of 10% and a decrease of 20%, respectively, in the estimate of RtDt. Vmax changed by
