Abstract
We proposed an alternative to a monoexponential model of radioiodine kinetics to obtain a more accurate estimate of absorbed doses to postsurgical thyroid remnants. We suggested that part of the difference between the predicted and the actually absorbed therapeutic doses of 131I, usually explained by radiation damage of thyroid cells, can be attributed to errors resulting from inadequate sampling of data and oversimplified modeling. Methods: A standard monoexponential model and alternative biphasic model (incorporating both radioiodine uptake and clearance) were used on 2 sets of patient data to fit time–activity measurements after administration of diagnostic and therapeutic activities of radioiodine. One set of data consisted of 633 records of routine measurements, and the second set consisted of 71 prospectively collected records with measurements performed more frequently and for a longer time. The time–activity curves derived from the 2 models were used to calculate residence times for diagnostic and therapeutic activities of 131I, and the respective residence times were compared using the paired t test. Errors of fitting and prediction of therapeutic time–activity data were also calculated. Results: With both models, a statistically significant difference (P < 0.01) was found between residence times after diagnostic administration of 131I and residence times after therapeutic administration of 131I. However, the effects of biphasic modeling and of improved sampling substantially reduced the difference (P < 0.01). Errors of fitting and prediction were smaller with the biphasic model than with the monoexponential model (P < 0.01). Conclusion: The biphasic model more accurately predicts 131I kinetics when applied to measurements in the short interval after diagnostic administration of radioiodine. The minimum requirement for the biphasic model is measurement twice a day at intervals > 6 h for at least 3 d after administration.
Radioiodine therapy of differentiated thyroid cancer is intended to ablate thyroid tissue or thyroid cancer metastases remaining after surgical treatment by administration of sufficiently large doses of radioiodine. There are 2 commonly used algorithms for prescribing the therapeutic administered activity, fixed or individualized (1). Considering individual variations in radiopharmaceutical kinetics and accumulating mass of the target tissue, patient-specific radiation dosimetry can improve the accuracy of dose estimates and increase the safety and cost-effectiveness of radioiodine therapy (1–3). The patient-specific treatment- planning paradigm, initially proposed by Benua et al. (4), consists of serial time–activity measurements over selected regions after the administration of a tracer activity of 131I. The kinetic data are then used to determine the cumulated activities in target tissues and to predict the therapeutic dose. After the administration of therapeutic activity, serial time–activity measurements and cumulated activity calculations are repeated and the absorbed dose is compared with the predicted dose. Experience has shown that in individual patients, the therapeutic dose actually absorbed in the target tissue is less than would have been predicted using only the kinetics of the pretreatment tracer radioiodine study (5–8). The explanation is that the therapeutic dose results in acute radiation damage that causes rapid leakage of iodide from the damaged thyroid cells (1). More recently, the concept of thyroid stunning was introduced, and it is being investigated (9).
The radiation dose to a target tissue is proportional to the area under the time–activity curve for that tissue, assuming (as with radioiodine and thyroid tissue) that the dose contribution of other tissues is negligible. The curve is fitted to analytic models, and the area under the curve is calculated by numeric or analytic integration. Accurate estimation of the area requires adequate sampling over a sufficiently long period (10,11). In clinical practice, however, only a limited number of measurements is usually available over a short interval of several days. In part because the required number and frequency of activity measurements depend on the number of adjustable parameters of the model, the function often used in practice is a simple monoexponential.
In this study, we evaluated an alternative biphasic model incorporating both radioiodine uptake and clearance. Although the complexity (the number of estimated parameters) is the same as in the monoexponential model, the biphasic model has the potential to provide more accurate estimates using data limited to initial periods after iodine administration. We suggested that a significant part of the difference between the predicted and actually absorbed therapeutic doses of 131I can be attributed to errors from inadequate sampling of data and oversimplified modeling.
MATERIALS AND METHODS
Patient-specific radiation dosimetry is based on measurement of 131I kinetics and on estimation of the accumulating mass of the target tissue. This study focused on examination of 131I kinetics. The parameter selected for comparison of monoexponential and biphasic kinetic models was residence time τ (in days):
Eq. 1
that is, the ratio of cumulated activity Ac (in MBq · d) and administered activity A0 (in megabecquerels). The cumulated activity is the area under the time–activity curve and represents the sum of all nuclear transitions in the region and interval of interest. In the MIRD schema, the residence time is used to calculate the mean absorbed dose per unit administered activity (10,11):
Eq. 2
where D (in grays) is the mean absorbed dose and S (Gy/MBq · d) is the mean absorbed dose per unit cumulated activity (for simplicity, arguments indicating source and target regions are omitted).
Patients
Radioiodine treatment was performed 6 wk after total or nearly total thyroidectomy for differentiated thyroid carcinoma. Patients were on a low-iodine diet and were hypothyroid, with the level of serum thyroid-stimulating hormone exceeding 30 mIU/L. Thyroid hormone substitution started 24 h after therapeutic administration of 131I. For the study, only the patients with neck lesions receiving the first radioiodine treatment were considered. Two groups of patients were analyzed. In the first group (n = 633), time–activity measurements after 131I administration were performed according to a routinely used standard protocol. In the second group (n = 71), the measurements were performed more frequently and for a longer time with respect to the standard protocol.
The first group of patients consisted of 506 females (age range, 9–81 y; mean age ± SD, 47.0 ± 14.3 y) and 127 males (age range, 11–78 y; mean age ± SD, 47.8 ± 14.9 y). Three hundred forty-one patients (53.9%) had papillary carcinoma, 162 (25.6%) had follicular carcinoma, and 130 (20.5%) had the mixed forms. The second group of patients consisted of 54 females (age range, 16–77 y; mean age ± SD, 45.3 ± 16.4 y) and 17 males (age range, 19–72 y; mean age ± SD, 46.6 ± 16.4 y). Fifty-seven patients (80.3%) were treated for papillary and 14 (19.7%) for follicular thyroid carcinoma.
Diagnostic activity (70–75 MBq) and therapeutic activity (3–8 GBq; mean, 4.2 ± 1.4 GBq) of 131I were administered in the form of Na131I solution. If the uptake in thyroid remnants exceeded 5%–9% of diagnostic administered activity, a surgical revision was considered. Otherwise, therapeutic activity was applied immediately after evaluation of diagnostic measurements (i.e., 3–6 d after the diagnostic application).
Measurement Protocols
According to the standard protocol, the counting rate was measured over the remnants of thyroid gland and other known lesions in the neck using a collimated scintillation probe with a NaI(Tl) detector connected to a multichannel analyzer. Measurements were performed in the constant distance (20 cm) between the detector and the inspected region at 3, 24, 48, and 72 h after the diagnostic administration. After the administration of therapeutic activity, the measurements were performed daily for 10–20 d.
To compare both curve-fitting models and to validate the biphasic model in detail, a prospective study with a modified protocol was performed on 71 consecutive patients. According to the modified protocol, measurements were performed twice a day at an interval of approximately 6 h for an extended period (i.e., up to 6 d after administration of diagnostic activity and up to 20 d after administration of therapeutic activity).
A vial with an iodine solution of known activity was used for calibration. The vial was placed in the thyroid phantom, consisting of a polyethylene cylinder of diameter 150 mm and height 160 mm with an eccentric hole of diameter 35 mm and depth 100 mm for the vial. The shortest distance between the hole and phantom surface was 5 mm. Calibration was performed before measurement and repeated every 2 h. For measurements after administration of therapeutic activities of 131I, a lead diaphragm (12 mm) was inserted in front of the probe to reduce the incident counting rate. All measured counting rates were corrected for background activity.
Monoexponential Model
In routine practice, the activity A(t) (in megabecquerels) of 131I in thyroid remnants at time t (in days) after administration is modeled (for t ≥ t1 ≅ 1 d) by the monoexponential function:
Eq. 3
where A(0) (in megabecquerels) is a constant and Tef (in days) is the effective half-life of radioiodine. The corresponding model describing measured values Am(t) of the activity A(t) in logarithmic form is:
Eq. 4
where n is normal noise. There are 3 unknown patient-specific parameters in the model to be estimated: the time t1 of maximum activity in the lesion, the constant A(0), and the effective half-life Tef. Obtaining the correct value of the time t1 (often substituted by the peak time of measured counting rate or fixed to 1 d after iodine administration) is critical for accurate curve fitting using Equation 4, especially when the number of measurements is low. Therefore, time t1 was estimated by shifting t1 between the highest counting rate values until the minimum error of the fit was achieved. Noise n reflects measurement and approximation errors as well as the random character of the observed biophysical processes. Assumed normality of the noise implies that the remaining 2 unknown parameters, A(0) and Tef, are optimally estimated by the least squares method. The residence time was then calculated as:
Eq. 5
Alternative Biphasic Model
In the biphasic model, the activity A(t) at time t measured since the time of administration (t = 0) is expressed as:
Eq. 6
where C1, C2, and C3 are patient-specific constants and Tp (in days) is the physical half-life of 131I (8.04 d). The corresponding model describing measured values Am(t) of the activity A(t) in logarithmic form is:
Eq. 7
where the normal noise n has the same justification as in the traditional monoexponential model. The assumed normality of additive noise again implies that the constants C1, C2, and C3 are optimally estimated by least squares. The residence time is calculated using Equation 1, in which the cumulated activity Ac is substituted by the area under the curve modeled by Equation 7. The area is obtained by a simple numeric integration.
Statistical Evaluation
Using monoexponential and biphasic models, residence times were estimated for 633 records with routine (relatively infrequently sampled) time–activity measurements and for 71 records with more frequent time–activity measurements. All estimated residence times were converted to natural logarithms to reflect their approximately lognormal distribution. The values of ln τd and ln τt were obtained from data measured after administration of diagnostic and therapeutic activities, respectively, of 131I. Comparisons were made pairwise in individual patients, and the differences were evaluated by a paired t test (12). Mean values and variances between the groups of patients were compared by a standard t test and F distribution (12).
Errors of fitting and prediction for monoexponential and biphasic models were calculated for the group of patients with more frequently performed measurements after therapeutic administration of 131I. Errors were calculated as the sums of squared deviations between measured values of activity and the values modeled by Equations 4 and 7. Errors of fitting were summed over all data points available for the curve fitting. Errors of prediction were summed over the points not used for the curve fitting after the parameters of the respective models were estimated from measurements in the short initial interval and the curves extrapolated to a longer time (Fig. 1). Occasionally, the number of measurements was close to the number of parameters and the fit was almost errorless. Seven of 71 data records with errors close to 0 were excluded from evaluation to avoid bias.
Example of curve fitting in 1 patient using monoexponential (A and C) and alternative biphasic (B and D) models applied to all (A and B) and initial (C and D) data. In each plot, only black points are used for curve fitting. Dashed lines show fitted curves. Figure shows ability of respective models to fit data and to predict course of time–activity curve when fitted to limited number of initial data points. Values of estimated parameters are t1 = 0.89 d, ln A(0) = 6.29, Tef = 2.29 d, τ = 0.22 d (A); C1 = 0.41, C2 = 5.80, C3 = −0.24, τ = 0.37 d (B); t1 = 0.89 d, ln A(0) = 5.74, Tef = 7.07 d, τ = 0.69 d (C); and C1 = 0.34, C2 = 5.72, C3 = −0.21, τ = 0.36 d (D). Average ln errors of fitting are 0.06 (A) and 0.03 (B). Average ln errors of prediction are 2.96 (C) and 0.08 (D). When only 1 measurement a day is evaluated, the errors increase to 0.07 (A), 0.06 (B), 3.88 (C), and 0.43 (D), respectively. exp. = exponential.
RESULTS
An example of time–activity curves fitted to measured data points in an individual patient is given in Figure 1.
Logarithms of residence times τd and τt estimated after the administration of diagnostic and therapeutic activities of 131I are compared in Tables 1 and 2. Mean values and ratios of τd and τt are summarized in Figure 2. With diagnostic records, the effect of the model was not significant in either of the 2 groups of patients. The effect of frequent sampling, however, resulted in significantly longer estimates of residence times with both models. After the administration of therapeutic activity, the effects of both the model and the sampling frequency were found to be significant. The longest estimates of residence times τd and τt were found with the biphasic model and frequently sampled data.
Mean values and ratios of τd and τt obtained by monoexponential (M) and biphasic (B) models for 633 routine records with infrequent sampling (M633, B633) and for 71 experimental records with frequent sampling (M71, B71).
Logarithms of Estimates of Residence Times τd After Diagnostic Administration of 131I
Logarithms of Estimates of Residence Times τt After Therapeutic Administration of 131I
Paired differences ln τd − ln τt are summarized in Table 3. Ideally, if there is no thyroid stunning or other possible effects causing the differences between the residence time estimated from diagnostic and therapeutic records, the differences should be 0 (i.e., the residence times estimated from diagnostic and therapeutic records should be identical). For the same reason, ratios τd/τt shown in Figure 2 should ideally equal 1. The condition is best approximated by the biphasic model when applied to frequently sampled data. However, even in the best approximation, the differences between residence times of 131I after diagnostic and therapeutic administrations remain significant.
Paired Differences Between Logarithms of Estimates of Residence Times ln τd − ln τt
Errors of fitting and prediction in logarithmic form are given in Tables 4 and 5. In comparison with the monoexponential model, the alternative model provides estimates with significantly smaller errors of fitting and prediction when 2 measurements a day are used for curve fitting. When only 1 measurement a day is used, the differences become insignificant.
Errors of Fitting of Measured Time–Activity Curves
Errors of Prediction of Measured Time–Activity Curves
DISCUSSION
The activity of radioiodine used for the ablation of thyroid remnants is not standardized. Tumor size, lymph node involvement, and distant metastases seem to have implications for the level of administered activity of 131I (13). Therapeutic activities used in practice vary widely from <1 to >12 GBq (13–15). The frequency of complications has been reported to increase with the dose (13), and in many patients very low doses have been shown to be sufficient to prevent recurrence of the disease (14,15). Prediction of therapeutic doses based on the determination of patient-specific radioiodine pharmacokinetics and lesion uptake thus could help to standardize protocols and improve both the safety and the cost-effectiveness of the treatment. Besides other factors (unknown accumulating mass of the thyroid remnants), the problem of dose prediction is that the residence time estimated after administration of therapeutic activity of 131I is substantially shorter than that estimated after administration of diagnostic activity (5–8) and that the resulting actual dose absorbed in the target tissue is thus lower than predicted. One of the possible explanations currently discussed in the literature is the effect of thyroid stunning: relatively low activities of 131I used either for diagnostic scintigraphy or time–activity measurements may result in a time-dependent reduction of thyroid 131I uptake (9). Recent studies on the possible stunning effect of low doses of radioiodine are, however, controversial, finding both its presence (16–18) and its absence (19,20).
To analyze thyroid stunning and other biologic effects of radiation in detail, all other possible effects that may contribute to the differences between the true and the predicted therapeutic doses should carefully be excluded. The aim of this study was to analyze the effects of modeling and curve fitting in time–activity measurements and to propose a practical and useful method to improve the planning of radioiodine therapy. All other factors affecting the calculation of therapeutic doses, especially the mass of the target tissue, were ignored.
A simple monoexponential model fitted to the initial part of the elimination curve does not properly predict the complete time–activity curve. Such a model neglects both the uptake and the long-term retention phases of 131I kinetics. Reliable fitting of more realistic models (i.e., models with more adjustable parameters) requires a higher number and frequency of measurements over a longer interval than is feasible in clinical practice.
The motivation of the alternative biphasic model is described in our previous reports (21,22). Briefly, the aim was to approximate the results of more complex multicompartmental models that, in practice, cannot successfully be estimated because of a lack of data. An additional important requirement was linearity in parameters that would guarantee both simplicity and robustness of estimation. The resulting biphasic model emerged from a series of experiments testing various simple analytic descriptions that produced time–activity curves similar to those provided by more complex compartmental models. The model takes into account the uptake phase and, using a limited number of well-distributed samples over the initial period, also predicts the long-term retention well. The biphasic model thus represents a compromise between simple but inadequate and complex but demanding models.
The biphasic model yields estimates of therapeutic residence time averaging approximately 50% higher than those yielded by the standard monoexponential model. As can be seen in Figures 1 and 2 and in Tables 3 and 5, the main advantage of the new model is its significantly better predictive ability than that of the monoexponential model. The difference between both models decreases with increasing number of measurements and the length of the interval used for curve fitting.
Besides the type of the model, the number, frequency, and duration of time–activity measurements was found to be important for accurate estimation of dose. According to experience based on this study, reliable curve fitting requires that time–activity measurements be performed after the administration of a diagnostic activity of 131I twice a day for at least 72 h. Daily measurements should be separated by at least 6 h. Measurement after the therapy should last longer, but the sampling frequency of data points can be lower.
The results obtained with the biphasic model show that radiation damage to thyroid cells may not be the only explanation for the short residence time of 131I found after administration of therapeutic activities of radioiodine. The differences between diagnostic and therapeutic kinetics can also be explained by the different lengths of time during which the measurements are performed, the different frequencies and regularities of data sampling, and the different models used to fit data. However, even with well-sampled data and the biphasic model, the difference between residence times estimated after the administration of diagnostic and therapeutic activities of 131I remains significant and is likely related to the radiobiologic effect of the 131I radiation.
CONCLUSION
The important feature of analytic models used for fitting of time–activity curves is their ability to predict, that is, extrapolate, the course of activity behind the interval of measurement. In comparison with the traditional monoexponential model, the alternative biphasic model better predicts the kinetics of 131I when applied to the short initial interval after diagnostic administration of radioiodine. After therapeutic administration of 131I, the biphasic model provides estimates of residence time that are closer to those predicted from diagnostic time–activity measurements. At the same time, the biphasic model requires estimation of only the same number of unknown patient-specific parameters as does the monoexponential model. Even with the biphasic model, both diagnostic and therapeutic time–activity measurements have to be performed as long and as frequently as is practical. The minimum for the biphasic model seems to be 2 measurements a day separated at an interval of at least 6 h for at least 72 h. Later measurements can be performed once a day.
APPENDIX
The parameters of both monoexponential and biphasic models are estimated by the regression C = (X′X)−1X′Y, where matrix C contains the estimated parameters of the respective model, matrix X contains an independent variable (time), and matrix Y contains a dependent variable (activity). X′ is the transpose of X.
In the monoexponential model, the content of the matrices Y, X, and C is Y = [ln A(t1), ln A(t2),…, ln A(tk)]′, X = [x1, x2], x1 = [1, 1,…, 1]′, x2 = [−t1, −t2,…, −tk]′, and C = [ln A(0), ln 2/Tef]′, where A(ti) (i = 1, 2,…, k) is the activity measured at time ti, k is the number of measurements, and A(0) and Tef are the estimated parameters of the monoexponential model (Eqs. 3 and 4). Time t1 is the peak time of the count rate estimated from measurements.
In the biphasic model, the content of the matrices Y, X, and C is Y = [ln A(t1) + t1 ln(2)/Tp, ln A(t2) + t2 ln(2)/Tp,…, ln A(tk) + tk ln(2)/Tp]′, X = [x1, x2, x3], x1 = [1, 1,…, 1]′, x2 = [ln t1, ln t2,…, ln tk]′, x3 = [t12/3 ln t1, t22/3 ln t2,…, tk2/3 ln tk]′, and C = [C1, C2, C3]′, where Tp is the physical half-life of 131I and C1, C2, and C3 are the estimated parameters of the biphasic model. Practical implementation of the estimation procedures requires introduction of physical limits (e.g., 0 initial value of cumulated activity) and use of the method of restricted least squares.
Acknowledgments
This study was supported in part by IGA grants 4581-3 and NN 5382-3/99 from the Ministry of Health of the Czech Republic and by grants 102/99/1564 and 102/00/D072 from the Grant Agency of the Czech Republic.
Footnotes
Received Nov. 1, 2000; revision accepted Mar. 5, 2001.
For correspondence or reprints contact: Jindřiška Heřmanská, PhD, Department of Nuclear Medicine, FN Motol, V úvalu 84, CZ-150 00 Praha 5, Motol, Czech Republic.
