A Pharmacokinetic Model Determination of Time Activity Curves in Radiopharmaceutical Therapy

Abstract

Introduction and Purpose

Radiopharmaceutical therapy (RPT) dosimetry can be challenging to perform due to sparse data measurements and variations in how the time activity curve (TAC) is determined. In this work, a single system of equations was theoretically derived to estimate the TAC.

Methods

A pharmacokinetic (PK) model was developed to estimate patient specific rate constants for a given set of body compartments. The PK model and an optimizer were numerically implemented to determine the rate constants and, using these physiologic data, to generate TACs and time integrated activities (TIAs) for 3 tissue systems from clinical data gathered in 5 patients. A fourth (aggregate) tissue compartment is added using conservation of activity considerations.

Results

Feasibility of the PK model was demonstrated by successfully generating TACs and TIAs for all patients in a manner comparable to existing methods in the literature. The data are compared to smaller sampling regimes. Differences between the 3- and 4-compartment models show that conservation of activity considerations should be part of TAC estimations.

Conclusion

The results here suggest a new paradigm in RPT in using the rate constants so identified as a diagnostic tool and as a vehicle to achieving individualized tumorcidal dose and/or the maximum tolerable dose to normal tissues.

Introduction

Radiopharmaceutical therapy (RPT) is a rapidly evolving modality that uses internally committed radiopharmaceuticals to treat disease through delivery of radiation dose. ¹ The benefit of RPT is that for some disease sites (eg, oncologic), the radiopharmaceutical has a high target to non-target ratio (also known as the tumor-to-normal ratio or TNR) that spares normal tissue while delivering tumorcidal dose to malignant lesions. ²

Despite high TNRs for many RPTs, normal tissues often still receive substantial dose. One currently accepted paradigm for RPT with substantial normal tissue dose is to deliver a total activity that will result in a dose: (1) known to be tumorcidal to a given tumor; and/or (2) that does not exceed the maximally tolerable dose (MTD) for any normal tissue structure. ³ Normal tissue structures that receive relatively high dose during RPT are known as organs-at-risk (OAR) and the MTD is the dose threshold below which deterministic effects are not expected to be present in a given OAR.

The dose to the tumor(s) and/or OAR(s) can be estimated through sparse measurements of activity using serial imaging and/or blood draws at known time points following a given RPT cycle.^4–9 The activity contained within each tumor and/or OAR is plotted as a function of time and is fit with a known function to create a time activity curve (TAC); the TAC can be integrated to determine the time integrated activity (TIA), which can be used to estimate tissue dose using a variety of formalisms and softwares.^8,10

The TAC fit can be performed using any user defined function, but is most often performed using standard functions such as combinations of exponential functions; ¹¹ these methods have been widely implemented. ¹⁰ These fits are typically performed for each tissue of interest independent of other tissues (eg, the tumor TAC is created independent of the kidney TAC for a given patient). Variability in activity measurement, TAC generation, and TIA calculation has demonstrably large variations in estimated radiation dose to tumors and/or normal tissue. ¹² This approach can be considered as a “top down” approach to modeling the TACs and, while useful, do not provide additional information that may be of significant clinical relevance.

The objective of this article is to develop and demonstrate feasibility of a novel pharmacokinetic (PK) model that uses a single system of equations to simultaneously determine the rate constants for all compartments as a function of one another for a given patient, and produce the TACs of interest. As the patient by definition is a single system with all compartments interconnected, this PK model may provide an initial step in advancing the determination of the TAC in patient specific RPT dosimetry. This is a “bottom up” model, where the rate constants, which are related to many physiologic parameters, are determined and they are used to determine the TACs. The determined rate constants may provide additional diagnostic/disease stage classifications that have not been pursued at this point. PK models have been widely used in pharmacology,^13,14 but their use in RPT has thus far been relatively limited.^15–17

To accomplish this objective, the PK model was theoretically developed in the “Theoretical considerations: PK Model” and “Conservation of Activity and a 4-Compartment” sections to provide the TAC for typical RPT infusion regimens. An optimization algorithm based on a hybrid gradient descent method is presented over the rate constant parameter space in the “Rate Constants: Optimization” section and its numerical implementation is subsequently described in the “Numerical Methods Implementation: Optimization” section. This optimization allowed for rapid convergence of the PK model to the known patient data. The patient data used in this study consisted of activity in 3 organ systems (blood, tumor, and kidney) measured at multiple time points using the clinical data provided by the publicly available American Association of Physicists in Medicine (AAPM): Understanding Time-Activity Curve and Time Integrated Activity Variations in Radiopharmaceutical Therapy (TACTIC) Challenge (https://www.aapm.org/GrandChallenge/TACTIC/default.asp) for patients undergoing Lutetium-177 Prostate Specific Membrane Antigen (Lu-177 PSMA) Therapy. The PK model and optimization resulted in rate constants that were used to create patient specific TACs for each tissue of interest which were compared to TACs determined using typical single and multiple exponential single-organ models. Using a conservation of activity scheme, this model is extended to include an additional separate aggregate compartment that includes all other tissues, excreted activity, and decayed activity. While measurements from this compartment were not available, conservation considerations lead to determination of activity in this compartment and thus additional tuning of the TACs. The data from the 3-compartment and 4-compartment models are compared, and further, numerical fits to smaller subsets of the data (2- and 3-observations) are computed and compared to the 5-data point computations. These comparisons provide some additional evidence of the goodness-of-fit between the multiple-time point (MTP) models in the full PK setting. The addition of the fourth compartment from conservation considerations suggests that this should be a significant step in TAC estimations. Finally, it is important to emphasize that this is not a fit to the TACs; the fit to the TAC provides a validation of how well the determined rate constants model the different patient data.

Methods

Theoretical Considerations: PK Model

Let $A_1^i(t)$ , $A_2^i(t)$ , and $A_3^i(t)$ be the activities in a given patient’s blood pool (blood), tumor(s), and a given OAR (kidneys) for patient i at time t. Figure 1 demonstrates the PK model that will be solved to determine the rate constants and TIAs in this problem. We use

$${\mathbb{A}}^i(t)=\left[\begin{array}{c}{A}_1^i\\ A_2^i\\ A_3^i\end{array}\right](t)$$

to denote the vector-valued function, and note that ${\mathbb{A}}^i(0)=[1,0,0{]}^T$ . The initialization of $\mathbb{A}$ at time 0 assumes that the activity is introduced as a delta function, normalized to 1. For alternative regimens of introduction of activity, we use ${\mathbb{A}}^i(0)=[0,0,0{]}^T$ and use the forcing function $\mathbb{F}(t)$ to model these infusion models. Clearly, this problem is linear under scaling, and multiplying ${\mathbb{A}}^i(0)$ by a constant multiple propagates through the calculations.

Figure 1.

The 3 compartment PK model.

In typical scenarios and as is the case here (the patient data provided by the TACTIC challenge), specification of activities can be somewhat vague. The theoretical results are normalized to incorporate the corrections made in the TACTIC data, which was provided as % activities within each organ system at each time point. Rate constants will have units of $t^{-1}$ , so since time is specified in $\mathrm{hr}$ in the present setting, the rate constants have units of ${\mathrm{hr}}^{-1}$ .

The problem presented in Figure 1 can be described mathematically as

$$\{\begin{array}{c}\frac{\mathrm{d}{\mathbb{A}}^{\mathrm{i}}}{\mathrm{d}\mathrm{t}}=\mathbb{M}{\mathbb{A}}^i+\mathbb{F}(t)\\ {\mathbb{A}}^i(0)=[0,0,0{]}^T\\ \mathbb{F}(t)=[\delta (t),0,0{]}^T\end{array}$$

(1)

It is important to notice that a physical solution of this problem requires

$$\underset{t\to \mathrm{\infty }}{lim}\mathbb{A}(t)=0$$

due to biologic excretion and physical decay of the radioisotope. The matrix $\mathbb{M}$ is given by

$$\mathbb{M}=\left[\begin{array}{ccc}-({\alpha }_1+{\beta }_1+{\gamma }_2)& {\alpha }_2& {\beta }_2\\ {\alpha }_1& -{\alpha }_2& 0\\ {\beta }_1& 0& -({\beta }_2+{\gamma }_1)\end{array}\right]$$

It is assumed that the rate constants ${\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,{\gamma }_1$ , and ${\gamma }_2$ are all strictly equal to or greater than 0, and are independent of time and activity levels in each of the compartments.

Solving for the eigenvalues of the matrix $\mathbb{M}$ , and noting that

$$det(\mathbb{M})=-{\beta }_1{\alpha }_2{\gamma }_1-{\alpha }_2{\beta }_2{\gamma }_2-{\alpha }_2{\gamma }_1{\gamma }_2<0$$

and the trace of $\mathbb{M}=-({\alpha }_1+{\alpha }_2+{\beta }_1+{\beta }_2+{\gamma }_1+{\gamma }_2)$ , $\mathbb{M}$ is nonsingular $(detM\ne 0)$ , and since it has a strictly negative trace, the system is fully dissipative. Thus the matrix of eigenvectors of $\mathbb{M}$ has full rank and is invertible. Setting ${\beta }_2=0$ will not change the invertibility or dissipativity of $\mathbb{M}$ ; this was considered as a special case.

Denoting the eigenvalues of $\mathbb{M}$ by ${\lambda }_1,{\lambda }_2,{\lambda }_3$ , and their corresponding eigenvectors by ${\zeta }_1,{\zeta }_2,{\zeta }_3$ , the solution of the homogeneous differential equation associated with (1) is

$$\begin{array}{rl}\mathbb{A}(t)& =c_1{\zeta }_1{e}^{{\lambda }_{1}t}+c_2{\zeta }_2{e}^{{\lambda }_{2}t}+c_3{\zeta }_3{e}^{{\lambda }_{3}t}\\ & =[{\zeta }_1{e}^{{\lambda }_{1}t}{\zeta }_2{e}^{{\lambda }_{2}t}{\zeta }_3{e}^{{\lambda }_{3}t}][c_1{c}_2{c}_3{]}^T\\ & =\zeta \mathrm{diag}(e^{{\lambda }_{1}t},e^{{\lambda }_{2}t},e^{{\lambda }_{3}t})[c_1{c}_2{c}_3{]}^T\end{array}$$

(2)

Here $\zeta$ is the $3\times 3$ matrix of eigenvectors of $\mathbb{M}$ , and $[c_1{c}_2{c}_3{]}^T$ is the transpose of the row vector of coefficients. Introducing the change-of-variable

$$\mathbb{A}=\zeta \mathbb{B}\mathrm{or}\mathbb{B}={\zeta }^{-1}\mathbb{A}$$

(3)

and substituting into (1), we have

$$\zeta \frac{\mathrm{d}\mathbb{B}}{\mathrm{d}\mathrm{t}}=\mathbb{M}\zeta \mathbb{B}(t)+\mathbb{F}(t)$$

Multiplying both sides on the left by ${\zeta }^{-1}$ and noting that ${\zeta }^{-1}\mathbb{M}\zeta =\mathrm{diag}({\lambda }_1,{\lambda }_2,{\lambda }_3)=:\mathrm{\Lambda }$ results in

$$\frac{\mathrm{d}\mathbb{B}}{\mathrm{d}\mathrm{t}}=\mathrm{\Lambda }\mathbb{B}(t)+{\zeta }^{-1}\mathbb{F}(t)=\mathrm{\Lambda }\mathbb{B}(t)+\mathbb{H}(t)$$

(4)

where $\mathbb{H}(t):={\zeta }^{-1}\mathbb{F}(t)$ . Equation (4) is a diagnonalization of the system in (1) and can be solved easily. If $B_j$ and $h_j$ are used to denote the j-th component of the vectors $\mathbb{B}$ and $\mathbb{H}$ respectively, then

$$B_j(t)={\int }_0^t{e}^{{\lambda }_j(t-s)}{h}_j(s)ds$$

(5)

Combining (5) and (3) gives the desired activities at time t. Of course, the integrated total activity (which is the TIA) denoted here by the vector $\mathbb{T}$ is obtained by

$${\mathbb{T}}^i={\int }_0^{\mathrm{\infty }}{\mathbb{A}}^i(t)dt$$

The superscript i refers to the i-th patient. This superscript is shown here to emphasize that the total committed activity for each patient may be different. The total activity in each organ are the components of the vector $\mathbb{T}$ .

Three special cases are worthy of consideration to flesh out the computations in the implementation section, compare the evolution of the activities and to partially (and implicitly) determine the rate constants on infusion models. Each of these cases has been and/or is clinically used to administer RPT (noting that very short injection times can be considered a delta function).

Case 1: Suppose that the isotope is introduced intravenously as a delta function at time $t=0$ . Then $\mathbb{F}(t)=(\delta (t),0,0{)}^T$ , and by the above $\mathbb{H}(t)={\zeta }^{-1}\mathbb{F}(t)$ . That is, if the first column of ${\zeta }^{-1}$ have components $(w_1,w_2,w_3{)}^T$ , then $\mathbb{H}(t)=\delta (t)(w_1,w_2,w_3{)}^T$ . Putting these into (5), the components of the vector $\mathbb{B}$ will be

$$B_j(t)=e^{{\lambda }_{j}t}{w}_j$$

(6)

and $\mathbb{A}(t)=\zeta \mathbb{B}(t)$ .

Case 2: If the isotope is introduced intravenously more slowly but at a steady (constant) rate, then $\mathbb{F}(t)=(f_1(t),0,0)$ , where

$$f_1(t)=\{\begin{array}{cc}1& 0\le t\le t_d\\ 0& t>t_d\end{array}$$

As in Case 1 above, $\mathbb{H}(t)=f_1(t)(w_1,w_2,w_3{)}^T$ . Putting this into (5) and integrating each component for $0<t<t_d$ , we have

$$B_j(t)=\frac{w_j}{{\lambda }_j}{e}^{{\lambda }_{j}t}(1-e^{-{\lambda }_{j}t})$$

(7)

and for $t\ge t_d$ ,

$$B_j(t)=\frac{w_j}{{\lambda }_j}{e}^{{\lambda }_{j}t}(1-e^{-{\lambda }_j{t}_d})$$

(8)

and, as before, $\mathbb{A}(t)=\zeta \mathbb{B}(t)$ . We set $t_d=10\min$ in our simulations.

Case 3: For gravity infusion, saline is added via gravity to vial of isotope that is intravenously connected to the patient. The saline drip increases the volume of the vial and results in the liquid vial contents entering the patient at a constant rate (mL/s) with an activity concentration that decreases exponentially as saline replaces the isotope. The amount of isotope infused can be described by an exponential

$$f_1(t)=11.88e^{-(\ln (2)/ti)t}$$

where $ti$ is the half-life of activity remaining in the vial and t is the time after starting the infusion. A previously reported $ti=3.5$ min ¹⁸ for a similar Lutetium-177 therapy was assumed for the Lu-177 PSMA therapy gravity infusion injection due to similarities in the infusion methods for both therapies and a lack of literature providing $ti$ for the Lu-177 PSMA therapy. This $ti$ results in the introduction of 95% of the activity in the first 15 min. Clearly the total infused activity at time t is the integral of this function from 0 to t, namely $\frac{11.88ti}{\ln 2}[1-e^{-(\ln (2)/ti)}t]$ . The coefficient $11.88$ normalizes the time integrated activity from $[0,\mathrm{\infty })$ to 1.

Putting this infusion regimen into (4) and (5) and integrating, similarly to Case 2, the following result is arrived at

$$B_j(t)=11.88w_j{e}^{{\lambda }_{j}t}\left[\frac{1-e^{-({\lambda }_j+\frac{\ln 2}{t_i})t}}{{\lambda }_j+\frac{\ln 2}{t_i}}\right]$$

(9)

Since the infusion function $f_1$ is exponentially decaying, the expressions in $[\cdot ]$ above will be a function of t for all $t\ge 0$ .

Before moving forward, there is a subtle point that requires re-visiting and will be important for the determination of the rate constants. There are 2 ways to introduce the initial activity in (1). The initial activity can be inserted as a delta (forcing) function (vector) $\mathbb{F}(t)$ or the problem may be treated as an initial-value problem (IVP) associated with a homogeneous system of equations, that is, $\mathbb{A}(0)=\mathbb{A}(t=t_0)$ . As an IVP, (2) describes the solution of the system, where the coefficient vector $[c_1,c_2,c_3{]}^T$ will be given by ${\zeta }^{-1}\mathbb{A}(0)$ , and thus

$$\mathbb{A}(t)=\zeta \mathrm{diag}(e^{{\lambda }_{1}t},e^{{\lambda }_{2}t},e^{{\lambda }_{3}t}){\zeta }^{-1}\mathbb{A}(0)$$

As a forced system of differential equations, where $\mathbb{A}(0)=0$ , Equations (4) and (5) describe the solution. Since the forcing function $\mathbb{F}=\mathbb{A}(0)\delta (t)$ is introduced as a delta function at time $t=0$ , $\mathbb{H}(t)={\zeta }^{-1}\mathbb{A}(0)\delta (t)$ . Once integrated,

$$\mathbb{B}(t)=\mathrm{diag}(e^{{\lambda }_{1}t},e^{{\lambda }_{2}t},e^{{\lambda }_{3}t}){\zeta }^{-1}\mathbb{A}(0)$$

and $\mathbb{A}(t)=\zeta \mathbb{B}(t)$ will agree with the solution of the IVP. Stated differently, the system in (1) is Markovian, and at any time the activity may be used as the initial activity for future times. At a cursory glance the reader may wonder if the units in this problem are consistent. The left-hand side of (1) is a rate, and $\mathbb{F}$ on the right-hand side may not appear as a rate. However, if $\mathbb{F}$ is chosen as a vector of delta functions, the 2 formulations of the problem are equivalent. If the problem is discretized over the time axis, in each discrete time interval, the clock is “reset to zero,” and the initial activity of that time interval is the activity that is evolved throughout that interval.

Conservation of Activity and a 4-th Compartment

The 3-compartment model shown in Figure 1 was necessitated as the activities in only those compartment were measured clinically. However, using a conservation scheme, a 4-th compartment may be introduced to account for significant communication between blood and all other organs throughout the body, and possible elimination from those areas. While the measured activity in the 4-th compartment is not available, theoretical activity in the 4-th compartment and the resulting eliminated activities can be calculated. As will be seen, this will allow significant improvement in the fitting of the theoretical results to the clinically measured data.

To maintain brevity and as the solution here follows the 3-compartment model closely, only the changes will be highlighted here. Let $A_1^i(t)$ , $A_2^i(t)$ , $A_3^i(t)$ , and $A_4^i(t)$ be the activities in patient i’s blood, tumor(s), OAR (kidneys), and all the other “reachable” organs at time t, respectively. As before, we suppress the exponent i, henceforth. Figure 2 demonstrates the 4 compartments and the rate constants of communication between these reservoirs. We use

$$\mathbb{A}(t)=[A_1,A_2,A_3,A_4{]}^T(t)$$

(10)

to denote the vector-valued function of the activities, and set $\mathbb{A}(0)=[1,0,0,0{]}^T$ . Now (1) holds with the $4\times 4$ matrix $\mathbb{M}$ given by

$$\mathbb{M}=\left[\begin{array}{cccc}-({\alpha }_1+{\beta }_1+{\kappa }_2)& {\alpha }_2& {\beta }_2& {\kappa }_1\\ {\alpha }_1& -{\alpha }_2& 0& 0\\ {\beta }_1& 0& -({\beta }_2+{\gamma }_1)& 0\\ {\kappa }_2& 0& 0& -({\kappa }_1+{\gamma }_2)\end{array}\right]$$

Similarly to the 3-compartment model, it is easily verified that the trace of this matrix $-({\alpha }_1+{\alpha }_2+{\beta }_1+{\beta }_2+{\gamma }_1+{\gamma }_2+{\kappa }_1+{\kappa }_2)$ and $det(\mathbb{M})<0$ . We note that ${\gamma }_1$ is the rate at which the activity leaves the kidney and ${\gamma }_2$ is the corresponding rate at which the activity leaves the “other body compartments” (eg, nonurinary excretion). So if $A_{R_1}$ is the activity eliminated from the kidney, and $A_{R_2}$ is the activity eliminated from the other body parts,

$$\begin{array}{r}{A}_{R_1}={\gamma }_1{A}_3\\ A_{R_2}={\gamma }_2{A}_4\end{array}$$

(11)

So, $A_4$ will be

$$\begin{array}{rl}{A}_4& =1-(A_1+A_2+A_3+A_{R_1}+A_{R_2})\\ & =1-(A_1+A_2+(1+{\gamma }_1)A_3++{\gamma }_2{A}_4)\\ & =\frac{1}{1+{\gamma }_2}[1-(A_1+A_2+(1+{\gamma }_1)A_3)]\end{array}$$

(12)

Strictly speaking, once the activities are corrected for their physical half-life, clinical measurements of the given activities $A_1$ , $A_2$ , and $A_3$ should provide an estimate of the activity elsewhere. Hence a good estimate for measured $A_4$ , ${\hat{A}}_4$ , may also be given by

$$(1+{\gamma }_2){\hat{A}}_4+{\gamma }_1{\hat{A}}_3={\hat{A}}_4+{\hat{A}}_{R_1}+{\hat{A}}_{R_2}=1-({\hat{A}}_1+{\hat{A}}_2+{\hat{A}}_3)$$

(13)

which after simplifying will be

$$\widehat{A_4}=\frac{1}{1+{\gamma }_2}[1-({\hat{A}}_1+{\hat{A}}_2+(1+{\gamma }_1){\hat{A}}_3)]$$

(14)

The identity of (12) and (14) is not an accident and demonstrates the matching of the theoretical estimations and the clinically measured activities. Comparing $A_4$ derived from (12) and ${\hat{A}}_4$ calculated from (14) will allow us to perform the optimization step in the numerical simulations over all the 8 parameters.

Figure 2.

The 4-compartment PK model.

Rate Constants: Optimization

The above compartmental PK model provides an exact solution of the activities as a function of time in each compartment. For specificity, we shall describe our numerical simulations for the 3-compartment model; the 4-compartment model descriptions are exactly analogous, with 2 additional parameters ${\kappa }_1$ and ${\kappa }_2$ introduced. The results for both 3 and 4 compartment models will be presented to highlight the differences.

To provide agreement with experimentally obtained data, the following steps were performed: The error is minimized between the experimentally determined activities $\hat{\mathbb{A}}$ at the specified times and the activities at those times $\mathbb{A}$ predicted by the PK model by varying the 6 parameters (the rate constants) of the PK model, that is,

$$\underset{{\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,{\gamma }_1,{\gamma }_2}{min}\Vert \mathbb{A}(t)-\hat{\mathbb{A}}(t){\Vert }_2$$

(15)

where $\Vert \cdot {\Vert }_2$ is the standard Euclidean norm (root mean square of its components). The gradient descent approach was applied in the 6-dimensional parameter space to achieve this optimization efficiently. As the activity levels in the different organs are significantly different and the problem is over multiple scales, a set of weights are introduced that are proportional to the ratio of the maximum activities in each of the organs. These weights ensure that the significantly smaller errors in magnitude, but not as a percentage of the activity level in those organs, are not ignored by the minimization routine, and the errors in each of the organs are accounted for as a percentage of the activity in those organs.

Numerical Methods Implementation: Optimization

Implementation of this PK model can be described methodically as follows, with code snippets performing the most relevant of these operations provided in Appendix A. This numerical solution was tested using the data available from the TACTIC Challenge data, which is reproduced in Appendix B. The data provided in Appendix B was described by the TACTIC challenge as [verbatim] “25 synthetic datasets of peri-therapeutic time-activity curves in kidneys, blood, and tumor after [Lu-177] PSMA therapy. These data are generated based on randomly selected high-quality clinical datasets, with realistic uncertainties added. While a sample of data for 2 patients is shown throughout the results section, in Appendix C the corresponding data for entire five patient cohort are shown for comparison sake. Since the clinically provided data here relies on the AAPM data, a full discussion of the statistical uncertainties inherent in the data due to photon counting, finite time interval data acquisition, patient variations and other factors is not feasible. Also, no additional information such as body weight, patient age and other parameters were available to quantify variations between the different patients.

The implementation follows the following steps:

As each set of calculations is patient specific, the following calculations will be done for each patient separately and the patient specific indices will be omitted.

Define the 6-dimensional parameter orthant $(0,\mathrm{\infty }{)}^6$ corresponding to the parameters ${\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,$ ${\gamma }_1,{\gamma }_2$ , and start with an initial random choice in this orthant.

Construct the $3\times 3$ matrix $\mathbb{M}$ and using a solver determine the eigenpairs $({\lambda }_i,{\zeta }_i)$ , $i=1,2,3$ for $\mathbb{M}$ .

Construct the $3\times 3$ matrix $\zeta$ with columns ${\zeta }_i$ , $i=1,2,3$ , and find the inverse matrix ${\zeta }^{-1}$ . $\zeta$ will always be invertible, as the matrix $\mathbb{M}$ is nonsingular.

Using Equation 5, compute the components of the vector $\mathbb{B}(t)$ for each of the given times $t=1.4,$ $18.5,42.2,65.8,138.4$ h. Assemble this data into a $3\times 5$ matrix.

The $3\times 5$ matrix $\mathbb{A}$ , where each column of $\mathbb{A}$ correspond to the activities in $(\mathrm{blood},\mathrm{tumor},\mathrm{kidney}{)}^T$ , respectively, and the different columns are the activities at the 5 different times $1.4,18.5,42.2,65.8,138.4$ hr is then obtained by multiplying $\mathbb{A}=\zeta \mathbb{B}$ .

Find the $L^2$ -error E as follows:

$$E=\sqrt{\sum _{j=1}^5(\sum _{i=1}^3|w_i({\mathbb{A}}_{i,j}-{\hat{\mathbb{A}}}_{i,j}){|}^2)}$$

E is the error between the PK model activities $\mathbb{A}$ at time $t_j$ and the experimentally measured activities $\hat{\mathbb{A}}$ at those times in all organs, and $w_i$ are the weights for the 3 organs that corrects for the different scales of measured activities in these organs. In particular

$$w_1=\frac{1}{{\hat{\mathbb{A}}}_{1max}},w_2=\frac{1}{{\hat{\mathbb{A}}}_{2,max}},w_3=\frac{1}{{\hat{\mathbb{A}}}_{3,max}}$$

Clearly, the PK model activities are a function of the parameters.

Partitioning a bounded subset of the parameter orthant, say $(0,1{]}^6$ , into $n^6$ subcubes, for each choice of the parameter, we calculate E at the $3^6-1=728$ adjacent nodes of this lattice. The smallest error at these 728-points gives us the next lattice point.

A random starting point and random step sizes for each parameter are generated, then the above step is repeated until the process converges at a certain number of significant digits. The number of significant digits can be increased by increasing the number of iterations while decreasing the step sizes whenever E for any of the adjacent nodes is not less than the E of the current rate constants.

Since the matrix $\mathbb{M}$ is not symmetric, the dependence of the eigenvalues of this matrix on the matrix entries (ie, the rate constants) is not a convex problem. Furthermore, the eigenvectors will also depend on the matrix entries. Lack of convexity suggests that there may be local minima in the optimization step leading to local minima. Since convergence to a global minimum depends strongly on the choice of the initial guess, the errors are initially calculated on a coarse grid in the parameter space to choose the best starting guess. An elementary calculation of the characteristic polynomials of the matrix $\mathbb{M}$ in (1) gives

$$p(\lambda )=-{\lambda }^3+\mathrm{tr}(\mathbb{M}){\lambda }^2-[({\alpha }_1+{\alpha }_2)({\beta }_2+{\gamma }_1)+({\beta }_1+{\gamma }_2)({\alpha }_2+{\gamma }_1)]\lambda +det(\mathbb{M})=0$$

Note that all the coefficients of this polynomial are strictly negative (the trace and determinant are strictly negative), and the solutions are real and negative. While

$$\begin{array}{rl}{\lambda }_1+{\lambda }_2+{\lambda }_3& =\mathrm{tr}(\mathbb{M})\\ {\lambda }_1{\lambda }_2{\lambda }_3& =\det (\mathbb{M})\end{array}$$

the dependence of the eigenvectors on the matrix entries does not allow one to reduce the minimization problem to one solely based on the eigenvalues of the matrix $\mathbb{M}$ . The calculation and minimization is repeated a hundred times (with random initial guesses) to guarantee the global minimum is achieved.

Finally, the rate constants are used to estimate activity in each organ at time increments of $1/60$ hr from $t=0$ to $t=\mathrm{\infty }$ (practically, $\mathrm{\infty }$ in this case can be assumed to be a time point where the activity of the radiopharmaceutical is expected to go to 0 Bq, such as at $t=10$ physical half-lives). Simpson’s rule is used to integrate the resultant curve to determine the TIA, where the TIA has units hrx % activity, where % activity can be converted to a unit of activity (eg, Bq) through multiplication by the total injected activity.

Results

Figures 3 to 5 show the experimental patient data and the resultant PK model predictions of activity with rate constants optimized to minimize the LSE for the 3-compartment model. These 3 figures each show the kidney, blood, and tumor logarithm of % activities, respectively, for Patient 10 for the 3 infusion methods described above. Patient 10 was chosen as a representative example of the data set; all other patients demonstrated similar results. Note that to retain the same patient descriptors used in the AAPM TACTIC data, the 5 patients are numbered 6 to 10.

Figure 3.

Comparison of PK model prediction to known patient data for the kidney in Patient 10 for all 3 infusion methods using the 3-compartment model. The rate constants, estimated time integrated activity, and LSE are provided in Table 1.

Figure 4.

Comparison of PK model prediction to known patient data for the blood in Patient 10 for all 3 infusion methods using the 3-compartment model. The rate constants, estimated time integrated activity, and LSE are provided in Table 1.

Figure 5.

Comparison of PK model prediction to known patient data for the tumor in Patient 10 for all 3 infusion methods using the 3-compartment model. The rate constants, estimated time integrated activity, and LSE are provided in Table 1.

The optimized rate constants for the 3 infusion regimens for Patient 10 using the 3-compartment model are shown in Table 1 as a representative example. A comparison of the optimized rate constants for steady-state infusion using the 3-compartment model for all 5 patients is shown in Table 2. The TIA and LSE for each organ system are also shown in these tables. It is important to point out that this model is not a fit to the activities. Rather, the rate constants between the different compartments are determined, and those are used to infer the activities in the different compartments.

Table 1.

Comparison of the Rate Constants, Time Integrated Activities, and LSEs for Different Infusion Models in Patient 10 Using the 3-Compartment Model.

Rate constant	Infusion method
or value	Steady state (5 min)	Instantaneous	Gravity
${\alpha }_1$	$6.61\times {10}^{-4}$	$5.59\times {10}^{-4}$	$5.74\times {10}^{-4}$
${\alpha }_2$	$2.38\times {10}^{-2}$	$2.51\times {10}^{-2}$	$2.40\times {10}^{-2}$
${\beta }_1$	$5.07\times {10}^{-2}$	$4.12\times {10}^{-2}$	$4.38\times {10}^{-2}$
${\beta }_2$	$0.00\times {10}^0$	$0.00\times {10}^0$	$0.00\times {10}^0$
${\gamma }_1$	$2.79\times {10}^{-2}$	$2.79\times {10}^{-2}$	$2.79\times {10}^{-2}$
${\gamma }_2$	$2.19\times {10}^0$	$1.77\times {10}^0$	$1.88\times {10}^0$
Blood TIA	$4.47\times {10}^{-1}$	$5.52\times {10}^{-1}$	$5.19\times {10}^{-1}$
Kidney TIA	$8.10\times {10}^{-1}$	$8.15\times {10}^{-1}$	$8.13\times {10}^{-1}$
Tumor TIA	$1.24\times {10}^{-2}$	$1.23\times {10}^{-2}$	$1.24\times {10}^{-2}$
LSE	$3.87\times {10}^{-3}$	$3.92\times {10}^{-3}$	$3.92\times {10}^{-3}$

Abbreviations: LSEs, least square errors; TIA, time integrated activity.

Table 2.

Comparison of the Rate Constants, Time Integrated Activities, and LSEs for the Steady-State Infusion Method in All 5 Patients Using the 3-Compartment Model.

Patient	6	7	8	9	10
${\alpha }_1$	$9.76\times {10}^{-4}$	$5.78\times {10}^{-4}$	$1.57\times {10}^{-3}$	$1.00\times {10}^{-4}$	$6.61\times {10}^{-4}$
${\alpha }_2$	$9.00\times {10}^{-3}$	$3.51\times {10}^{-3}$	$9.49\times {10}^{-3}$	$7.17\times {10}^{-3}$	$2.38\times {10}^{-2}$
${\beta }_1$	$3.30\times {10}^{-2}$	$3.95\times {10}^{-2}$	$3.84\times {10}^{-2}$	$8.94\times {10}^{-2}$	$5.07\times {10}^{-2}$
${\beta }_2$	$0.00\times {10}^0$	$0.00\times {10}^0$	$0.00\times {10}^0$	$0.00\times {10}^0$	$0.00\times {10}^0$
${\gamma }_1$	$2.37\times {10}^{-2}$	$2.23\times {10}^{-2}$	$3.49\times {10}^{-2}$	$2.36\times {10}^{-2}$	$2.79\times {10}^{-2}$
${\gamma }_2$	$1.72\times {10}^0$	$1.94\times {10}^0$	$1.73\times {10}^0$	$2.03\times {10}^0$	$2.19\times {10}^0$
Blood TIA	$5.70\times {10}^{-1}$	$5.06\times {10}^{-1}$	$5.67\times {10}^{-1}$	$4.73\times {10}^{-1}$	$4.47\times {10}^{-1}$
Kidney TIA	$7.91\times {10}^{-1}$	$8.93\times {10}^{-1}$	$6.23\times {10}^{-1}$	$1.78$	$8.10\times {10}^{-1}$
Tumor TIA	$5.46\times {10}^{-2}$	$4.73\times {10}^{-2}$	$8.42\times {10}^{-2}$	$5.41\times {10}^{-3}$	$1.24\times {10}^{-2}$
LSE	$3.19\times {10}^{-2}$	$1.95\times {10}^{-2}$	$1.26\times {10}^{-2}$	$1.48\times {10}^{-2}$	$3.87\times {10}^{-3}$

Abbreviations: LSEs, least square errors; TIA, time integrated activity.

As a comparative measure, Figure 6 was generated with the 5 data points from Patient 6 using the steady-state PK 3-compartment model and also by fitting the known patient data to:

a single exponential model with 2 parameters of the form $p_1\exp (-p_{2}t)$ , constrained to the maximum value at $t=0$ ;

a bi-exponential model with 3 parameters of the form $p_1(\exp (-p_{2}t)-\exp (-p_{3}t))$ , constrained to be 0 at $t=0$ ; and

a bi-exponential model with 4 parameters of the form $p_1\exp (-p_{2}t)-p_3\exp (-p_{4}t)$ .

These functions were fit to the data using the Solver add-on (minimization of objective function) in Microsoft Excel (Microsoft Corporation, Redmond WA). The LSE found using the Excel solver addon by minimizing the objective function was comparable to the LSE presented in Tables 1 and 2 ( ${10}^{-2}$ to ${10}^{-3}$ ). These LSE limits of our solver easily explain the dispersion between the clinically measured values and the model predictions. As it is demonstrated in later figures improvements have been achieved by accounting for the 4-compartment.

Figure 6.

Comparison of pharmacokinetic (PK) model prediction (steady state) using the 3-compartment model and best-fit lines with error minimized using the Microsoft Excel Solver add-on to known data (blue circles).

Due to various clinical limitations, often it is hard to obtain measurements at 5 different time points. The fit here has also been carried out using 2 time-point and 3 time-point to compare to the given 5 time-point measured data. Figures 7 to 10 show these results, where different combination of 2-points and 3-data points are used for the optimization process. The optimized rate constants and blood–kidney–tumor TIAs and the LSE are shown on each of these figures for comparison.

Figure 7.

Comparison of pharmacokinetic (PK) model prediction with 2 time-points and measured data using the 4-compartment model. The PK model was optimized using only with the patient data acquired at 1.4 and 18.5 h; the fit only utilizes these 2 data points and is compared to all clinically measured data.

Figure 8.

Comparison of pharmacokinetic (PK) model prediction with 2 time-point measured data using the 4-compartment model. The PK model was optimized using only with the patient data acquired at 18.5 and 42.2 h; the fit only utilizes these 2 data points and is compared to all clinically measured data.

Figure 9.

Comparison of pharmacokinetic (PK) model prediction with 2 time-point measured data using the 4-compartment model. The PK model was optimized using only with the patient data acquired at 65.8 and 138.4 h; the fit only utilizes these 2 data points and is compared to all clinically measured data.

Figure 10.

Comparison of pharmacokinetic (PK) model prediction with 3 time-point measured data using the 4-compartment model. The PK model was optimized using only with the patient data acquired at 1.4, 18.5, and 42.2 h; the fit only utilizes these 3 data points and is compared to all clinically measured data.

Using the conservation scheme developed in the Conservation of Activity and a 4-Compartment section, Figures 11 and 12 were generated. Figure 11 shows the activities in the blood, kidney, and tumor in Patient 10 using the 4-compartment model, the corresponding 8-parameter rate constants, the different organ TIAs and the LSE. Figure 12 is to be compared to Figure 13, the 3-compartment model results initially presented in Figures 3 to 5, but here put in a similar format as in Figure 12. The rate constants, TIAs, and LSEs for these 2 plots are summarized in Tables 1 and 3. Finally, Figure 14, to be compared to Figure 6, compares the 4 compartment model activities to the exponential and biexponentials fits of those data.

Figure 11.

Comparison of pharmacokinetic (PK) model prediction where conservation is used to provide a complete 4-compartment model.

Figure 12.

Complete 4-compartment model with different infusion regimens.

Figure 13.

Complete 3-compartment model with different infusion regimens (Figures 3 to 5 in semilog format).

Figure 14.

Comparison of pharmacokinetic (PK) model prediction (steady state) using the 4-compartment model and best-fit lines with error minimized using the Microsoft Excel Solver add-on to known data (blue circles).

Table 3.

Comparison of the Rate Constants, Time Integrated Activities, and LSEs for Different Infusion Models in Patient 10 Using the 4-Compartment Model.

Rate constant	Infusion method
or value	Steady state	Instantaneous	Gravity
${\alpha }_1$	$5.32\times {10}^{-5}$	$5.51\times {10}^{-5}$	$5.83\times {10}^{-5}$
${\alpha }_2$	$3.89\times {10}^{-2}$	$3.52\times {10}^{-2}$	$3.69\times {10}^{-2}$
${\beta }_1$	$7.89\times {10}^{-1}$	$7.79\times {10}^{-1}$	$7.59\times {10}^{-1}$
${\beta }_2$	$5.18\times {10}^{-3}$	$5.37\times {10}^{-3}$	$5.22\times {10}^{-3}$
${\gamma }_1$	$2.69\times {10}^{-2}$	$2.93\times {10}^{-2}$	$2.92\times {10}^{-2}$
${\gamma }_2$	$0.00\times {10}^0$	$0.00\times {10}^0$	$0.00\times {10}^0$
${\kappa }_1$	$2.07\times {10}^{-3}$	$2.01\times {10}^{-3}$	$2.37\times {10}^{-3}$
${\kappa }_2$	$3.41\times {10}^{-2}$	$3.37\times {10}^{-2}$	$3.09\times {10}^{-2}$
Blood TIA	$1.46\times {10}^0$	$1.47\times {10}^0$	$1.51\times {10}^0$
Kidney TIA	$3.59\times {10}^1$	$3.29\times {10}^1$	$3.32\times {10}^1$
Tumor TIA	$1.99\times {10}^{-3}$	$2.29\times {10}^{-3}$	$2.38\times {10}^{-3}$
LSE	$7.94\times {10}^{-3}$	$6.60\times {10}^{-3}$	$7.80\times {10}^{-3}$

Abbreviations: LSEs, least square errors; TIA, time integrated activity.

Discussion

In RPT, patient-specific dosimetry is the calculation of radiation dose delivered to malignant lesions, tissue systems, and/or the whole body of individual patients. Often patient-specific dosimetry can be challenging. This is due to several factors. One factor is the sparse collection of data following RPT; the number of data points collected is generally clinically limited to a maximum of 5 post-RPT imaging and/or blood draw data points over a period of 5 to 10 days, with fewer data points (eg, 1 or 2) being most preferable for clinical workflow and patient convenience. ⁸ Another factor is that every patient is heterogeneous; patient anatomy (eg, given tissue size) and physiology (eg, uptake, distribution, metabolism, and excretion of the radioactive drug by a given tissue) can very greatly per patient. ¹

Because of these factors, historically patient-specific dosimetry has been performed by fitting simple functions such as exponentials or bi-exponentials to the collected data. These functions are individually fit to data from each tissue of interest; the parameters of the functions and thus the final TIA for each organ are independent of all other organs.

Physiologically, the patient is a single system; the rate constants between various tissues and thus the resultant TIA are interconnected. In this article, a PK model was developed to simultaneously determine the rate constants (which can be used to determine the TAC and TIA) for all organ systems of interest in a given patient. Extraction of the rate constants using the PK model may provide a new diagnostic tool that would have considerable clinical applications. One potential application may be the estimation of activity in organ systems where the activity is not measured through a priori knowledge and sparse measurements of other patient organs. An additional application may be in the development of new radiopharmaceuticals with improved rate constants that would target specific sites.

However, the main purpose of this article was to demonstrate feasibility of the PK model using clinical patient data. Figures 3 to 5 demonstrate the PK model developed in this article can be utilized to model the TAC in several organ systems using sparse clinical measurements of activity in patients. Importantly, these figures have scales on the y-axis that are different by several orders of magnitude, that is, $1\times {10}^{-1}$ , $1\times {10}^{-3}$ , and $1\times {10}^{-4}$ , and the dispersion in the data is somewhat related to roundoff errors and the mesh grid in gradient descent implementation. The mesh in the gradient descent is adaptively optimized in real-time not to miss the optimal rate constants up to $1\times {10}^{-6}$ scale.

In semi-log plots, while the fit for the blood compartment may not appear as good as the fit for the tumor and kidney compartments, the calculated percent errors for each of the organs (blood, kidney, and tumor) were 0.13%, 0.29%, and 6.06% for the blood, kidney, and tumor, respectively. The different order of magnitude change in scales for the blood compartment compared to the other organs is the cause of this apparent dispersion.

Figure 6 demonstrates the results of the PK model are comparable to existing single exponential and bi-exponential models (3 and 4 parameters) described in the literature.^8,19–21 This comparison is significantly improved in the 4 compartment case (to be discussed shortly) and it further validates the feasibility of the PK model compared to other existing methods. While the superiority of the model described was not assessed other than through analysis of LSE, comparability of the described PK model to existing models shows promise. Intuitively utilizing a PK model provides a co-dependency of the different compartments and this would be superior to evaluating each compartment as a separate system.

Figure 6 shows the PK model provides excellent agreement with the single exponential model and the 2 bi-exponential models for the kidney. In the tumor, there was less agreement, but this plot allows for some interesting interpretations. The first is that the single exponential model, while reasonably fitting the data, is not valid from a biological sense for an organ system such as the tumor, because the single exponential model assumes the maximum activity concentration in each organ occurs at $t=0$ . Second, while not shown in the plot, for this patient at $t=60$  s the value of the 4 parameter bi-exponential fit is $-1.75\times {10}^{-2}$ ; negative activity concentrations are not physically possible and the magnitude of this value is approximately 2 orders of magnitude larger than the maximum activity contained in the lesion at any given time ( $6\times {10}^{-4}$ ); indeed one of the great challenges of using the 4 parameter bi-exponential fit is that the activity can be negative.

When reviewing the activity in blood as shown in Figure 6, the single exponential model and the 2 bi-exponential models cannot be differentiated from one another and fit the known data very well. The PK model shows a rapid increase in blood activity due to the infusion, which is what occurs physically (particularly when considering the delta forcing function). While the PK model appears to model the late washout phase of blood well, the initial washout phase appears relatively poorly modeled with the blood activity concentration decreasing much faster than demonstrated by the patient data. This is most likely due to the use of 3 compartments and with the excretion pathway ${\gamma }_2$ exiting the system with no return. The analysis using the 4-compartment model (with the fourth compartment being the whole body with return) as shown in Figures 11 and 12 has demonstrated significant improvements in using the PK model to estimate blood activity. Our conservation scheme described in the Conservation of Activity and a 4-Compartment normalizes the data provided by the AAPM TACTIC Challenge in the 3 given compartments, accounts for the excretion of activities and thus allows us to account for the “whole body” activities. These results are shown in Figure 12, which may be compared to Figures 5 to 7, or more straight-forwardly to Figure 13. The comparison between the 4 compartment results and exponential and bi-exponential fits are shown in Figure 14. This plot should be compared to Figure 6. The improvement in the 4-compartment fit for the blood compartment in Figure 12 compared to that of the 3-compartment model in Figure 13 and improved matching to the exponential fits adds additional relevance to the conservation scheme developed.

Additionally, the challenge of fitting all compartments with a PK model were not unexpected, as shown in a prior publication ¹⁵ . PK models often have difficulty in fitting all data sets. Thus, while the results of this article demonstrate feasibility of using a PK model, it is clear the model can be further refined and calibrated for each patient, even if relatively few data points are collected (< 5) or if data is not collected for all organ systems. This PK model may be utilized to fill in missing gaps and provide a more accurate patient dose calculation. It is hoped that as this PK model matures, future studies will demonstrate improvements of the PK models over existing methods for producing TACs and TIAs.

The PK model additionally demonstrated some interesting trends as shown in Table 1, which, while somewhat illustrative in nature (since the optimization is done relative to the given data set), shows close agreement between the estimated TIA and rate constants in the blood, kidney and tumor for all infusion methods. The variations in the TIA were within 1% of each other, and the variability in the extracted rate constants were about 5%.

Interestingly, Table 2 demonstrates the variability in the TIA and rate constants among the Patients 6 to 10 can be by as much as a factor of 4. The tumor TIA in Patient 10 was $1.24\times {10}^{-2}$ compared to $5.46\times {10}^{-2}$ for Patient 6. Variability between patients is well-known, ¹ and thus this result is not surprising but perhaps demonstrates that the utilization of a standard model such as an exponential or bi-exponential may not be as patient specific as a single PK model.

As one final point of discussion, the PK model allows for modeling of the actual administration of the radiopharmaceutical used in RPT (instant bolus injection, steady-state infusion over a known time period, and gravity infusion). To the knowledge of the authors, the currently utilized curve fitting methods (exponential and bi-exponential) do not model the time dependent administration. For some organs, such as blood (and thus bone marrow), that may result in underestimation of dose. Of course, the PK model as demonstrated in the above results also appears to underestimate the activity in blood; this will be investigated through the implementation of more compartments and also the incorporation of more physiological data.

The fitting problem has also been carried out using 2 time-points and 3 time-points to compare to the given 5 time-point measured data. Figures 7 to 10 show these results. The rate constants, TIAs, and the LSEs are also shown in these figures. Interestingly and quite surprisingly, optimization with smaller data points leads to quite small LSEs, almost comparable to the 5-data point results. This demonstration may serve as a possible validation of using very small data point sets and having confidence that the optimized results for the rate constants and activities are relatively close to those obtained with more data points. Figure 7 to 9 also vary the 2 time points chosen and show comparable results.

Brosch-Lenz et al^22,23 compare single-time point (STP) to MTP imaging protocols for Lu-177 PSMA therapy and show data for a single photon computed tomography (SPECT) scan at 48 or 72 h scans. While their results are comparable to MTP data, their methods provide a posteriori dosimetry in Lu-177 PSMA studies. A closely related article ²⁴ provides similar data for Lu-177 PSMA therapies for kidney and tumor based on the SPECT images at 2 imaging time points. Peters et al ²⁵ results are consistent with the above studies and their results show that they can obtain lesion dosimetry that is within 20% (compared to 14% for MTP studies) using one time point SPECT imaging at 168 h. Jackson et al ²¹ provided similar results for Lu-177 PSMA RPT and investigate the STP scans that would provide the “best” dosimetry estimates for patients. In a technical report, Theresa Ida Götz ²⁶ provides a deep analysis of Lu-177 therapies TIA and provides the relationship of the TIA with sampling, voxel sizes, and numerical modeling performed for a typical patient. Several additional recent studies on Lu-177-PSMA SPECT imaging may also be cited.^27–29

Returning to the analytical modeling of the input function, if a therapy is designed to provide certain tumorcidal and MTD doses, the administered activity cannot be performed retrospectively and changing future therapies may not be effective due to tissue stunning or other similar effects. Thus, utilization of a PK model may allow for more patient-specific dosing using optimized input functions to achieve the rate constants and dose levels that will result in the greatest likelihood of normal tissue sparing and disease control. This is feasible since the PK model can provide the rate constants for a given tissue system from activity measurements obtained by tracer studies and sampling at only a few distinct time points. The resulting rate constants can then be used to create the entire TAC and thus calculate the activity needed for the desired outcome (tumorcidal dose or reaching a given MTD). This method is analogous to standard practice in both radiation oncology (eg, intensity modulated radiation therapy to provide a certain dose distribution) and radiology (eg, automatic exposure control to provide a certain image quality). This can be inferred from Figures 3 to 5 and is also discussed in the literature. ¹¹ As there is an expected significant variation across patients, this model can lead to developing individualized therapies and better long-term outcomes for patients.

Conclusion

In this article, a theoretical PK model framework was developed to estimate rate constants and TIAs from post-therapy activity measurements in patients undergoing RPT. This PK model used one system of equations to describe the flow of the radiopharmaceutical in the patients body through various compartments. The PK model determined the rate constants and used these physiologic parameters to determine the TACs. The fit to clinical data provided validation, and comparison with existing methodologies demonstrated feasibility.

This PK model based determination of the rate constants provides a diagnostic tool that deserves a closer look to categorize and classify normal tissues, tumors, and patients. Looking at populations of patients and establishing statistics for these parameters would allow physicians to determine if individual patients are significantly different from the population. In addition, and much more importantly, this or similar type PK models may be used to create steady-state activities and TIA in the tumor, by prescribing appropriate infusion regimens to achieve those patterns. This inverse problem is the subject of the immediate next research effort aimed to refine the above findings.

Finally, a thorough analysis of the dependence and variation of the rate constants on the counting statistics, regional nature of measured TACs, and dosimetry using Monte Carlo simulations is the subject of ongoing research.