A Novel Approach to Parametric Studies in Nuclear Medicine

Abstract

Finding the rate constants (rates of transition) of radiopharmaceuticals between organs or different regions of interest (ROIs) in nuclear medicine quickly and accurately provides rich physiological data that is essential for determining drug dosages and performing precise radiation dosimetry. This article describes a novel exact analytical approach to achieving this goal. This can be done in a very short time period at the beginning of treatment with application of tracer levels of radiopharma, and using activities in the ROIs construct the rate constants. Here, using randomly chosen rate constants, time activity curves (TACs) for a four-compartment model are created, and using different sampling regimens of these TACs, the rate constants are reconstructed in a large number of simulations. The least square error (LSE) and relative LSE in the reconstructed parameters are shown to be better than $1\times {10}^{-4}$ . While studies based on optimization methods have been described previously, this paper achieves the rate constants accurately and robustly using analytical methods, thus improving on the clinical applicability of these methods.

Introduction and Motivation

With the continued refinement of radiological science, the reliability of quantitative data has improved markedly. However, the results of nuclear imaging studies are principally considered qualitatively, with quantitative assessment typically limited to measurements of standard uptake value within regions of interest (ROIs) which correspond to organs and neoplasms. This has not been due to lack of effort by scientists, engineers, and clinicians; underlying mathematical and statistical complexities have limited progress, particularly for modalities that require time-dependent data acquisition in the clinical setting, relating the time activity curves (TACs) for different organs or ROIs in real time.

The physiologically based pharmacokinetic (PBPK) models with which we are concerned are a mainstay of the literature in physiology and nuclear medicine.^1,2 In the language of differential equations, such models provide the equations of motion à la Newton, for the transport and distribution of drugs or other compounds throughout the body. Modeling approaches identify organs or tissues as interconnected compartments that exchange substances, namely the drug or compound of interest, between them. The parameters of such models, called rate constants, encode physiologically meaningful information that is specific to the individual. These rate constants are not accessible by clinical examination or qualitative radiological evaluation. This leads to an inverse problem that asks if one can determine a unique set of rate constants from a dynamic nuclear medicine scan. In Steiner et al., ³ a model to achieve exactly these rate constants is developed by performing a “bottom-up” optimization method over the parameter space determined by the rate constants. That article also references the literature extensively and compares the results to those that are reported in the literature using bi-exponential and multiexponential models. While an excellent fit to TAC is achieved, the optimization problem is nonconvex and can suffer from getting trapped in local minima, not to mention that gradient descent optimization methods are computationally intensive.

Here, we present a novel approach to solving this inverse problem by applying a new mathematical result that produces highly accurate results in a very short time period and completely analytically. The mathematical theory and algorithm presented in Ellingson and Jafari ⁴ side-steps the need to perform an elaborate optimization and almost instantaneously solves the problem for the rate constants, based on sampling of the TAC from multiple ROIs. Furthermore, the proposed exact analytical method allows the determination of the rate constants from a minimal set of data acquired over a short time span from possibly small activities administered to the patient; thus, this method has potential for clinical applications. This method provides many advantages over existing techniques in allowing the estimation of physiological parameters in a broad set of scenarios, making the required dosage and the subsequent dosimetry more accurate than is presently available.^2,5–8

Theory

In nuclear medicine and radiopharmaceutical therapy (RPT), the injected radiopharmaceutical (RP) distributes over time throughout tissues, commonly modeled as flow from one compartment to another throughout the body. A multicompartment model is often used to mathematically describe these transitions (see Steiner et al., ³ Hobbs, ⁵ Siegel et al., ⁶ Nestorov, ⁹ Gospavic et al., ¹⁰ and the references therein) in terms of a system of differential equations.

The solution of this system represents the activity in each compartment as a linear combination of exponential functions of the form:

$$f(t)=\sum _{n=1}^N{c}_n{e}^{{\lambda }_{n}t},$$

(1)

for an $N$ compartment system, where ${\lambda }_n$ are the eigenvalues of the transition matrix modeling the dynamics of the flux between the different compartments (also see Steiner et al. ³ ). In particular, the problem presented by Figure 1 can be described by the following initial value problem:

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

(2)

Figure 1.

An Example Compartmental Model.

where $\mathbb{A}$ represents the time-dependent vector of activities in the different compartments and the matrix $\mathbb{M}$ is given by

$$\mathbb{M}=\left[\begin{array}{cccc}-(k_{12}+k_{13}+k_{14})& k_{21}& k_{31}& k_{41}\\ k_{12}& -k_{21}& 0& 0\\ k_{13}& 0& -(k_{31}+{\gamma }_3)& 0\\ k_{14}& 0& 0& -(k_{41}+{\gamma }_4)\end{array}\right].$$

(3)

Mathematically, $\mathbb{M}$ is the transpose of the negative out-degree Laplacian of a weighted digraph. In equation (1), $c_n$ are related to the components of the eigenvectors, and ${\lambda }_n$ are the eigenvalues of $\mathbb{M}$ . To determine the TACs and time integrated activity (TIA) in each compartment, one needs to determine the $c_n$ and ${\lambda }_n$ in equation (1) to a relatively high degree of precision. Therefore, one may ask if a pharmacokinetic (PK) “topogram” can be constructed relatively quickly from the administration of small levels of RP to provide the rate constants between ROIs dynamically and in real time. This data may then be used to determine the exact dosage to be used to achieve the desired outcomes, achieve appropriate tumoricidal dosage, while having optimal tissue-sparing effects. Eventually, this may also aid in developing novel RP with improved physiologic and PK characteristics.

This problem, which is also known as the de Prony problem (Gaspard de Prony, 1755–1839), has received significant attention over the past 50 years by well-known numerical and computational analysts (see, e.g. Golub and Pereyra, ¹¹ Pereyra, ¹² Sarkar and Pereira, ¹³ Paluszny et al., ¹⁴ Hussen and He, ¹⁵ and the many references therein). Although not fully incorporated into practice, currently these parameters can be determined using numerical optimization methods that fit the data acquired on gamma cameras, single photon emission tomography (SPECT), or positron emission tomography (PET) machines. The limiting factors seem to be that data acquisition is required to be over long periods of time, frequently over equally spaced time points, and can be quite time consuming. A good fit may be achieved using gradient descent and similar methods ³ ; however, this involves long post-processing steps and the numerical methods are complex. Exponential curve fitting routines are also used and are based on least square error (LSE) minimization (regression)^7,16; however, the latter methods suffer from isolation of each compartment and do not view the system as a whole.

Here, we present a complete analytical solution to this problem and demonstrate how these analytical methods may be implemented on nuclear medicine platforms to provide nearly instantaneous results. In particular, this solution eliminates the uniform sampling requirements of the de Prony/Sarkar ¹³ methods (matrix pencil method (MPm)) and can be achieved over short time periods up to statistical concerns. Having attained the $c$ s and $\lambda$ s, we relate these to the rate constants/transfer coefficients of the compartmental models. In the discussion, we discuss the implicit time-limiting assumptions on the data acquisition. Several examples will be shown to demonstrate the utility of the algorithms and their validation.

Existence and Uniqueness of the Solution

For nondegeneracy and as is the case in fully dissipative systems, as in RPT, we assume that the eigenvalues ${\lambda }_n$ are strictly negative and distinct. The special case of ${\lambda }_n=0$ (for some $n$ ) may be dealt with via a simple translation. Clearly, if the distinctness assumption is removed, equation (1) reduces to a linear combination of $N-1$ or fewer terms. We also note that $f(t)$ in equation (1) is nonnegative for all $t\ge 0$ , as it designates the activity at time $t$ in a given compartment. Thus, the main problem can be reduced to finding $c_n$ and ${\lambda }_n$ in equation (1). The mathematical methods developed are far more general, and do not depend on these particular assumptions, as can be viewed in Ellingson and Jafari, ⁴ and they readily extend to the case of complex exponentials. For RPT applications, we focus on the case where the parameters are real.

With $f(t)$ given by equation (1), define

$$I_{(k)}(t):={\int }_0^t\cdots {\int }_0^{t^{(k+1)}}f(t^{(k+1)})\text{d}{t}^{(k+1)}\cdots \text{d}{t}^{(1)},k=1,2,\dots ,$$

(4)

where $t^{(k)}$ is used to denote the $k$ th variable of integration, and there are $k+1$ integrals. The following theorem (see Ellingson and Jafari ⁴ for proof and detailed discussion), which is the main result from which this paper is derived, may be stated in terms of $I_{(k)}(t)$ , or equivalently in terms of the time-dependent moments of $f(t)$ .

Theorem 1

Let $f(t)=\sum _{n=1}^N{c}_n{e}^{{\lambda }_{n}t}$ , assume that $f(t)\ge 0$ , ${\lambda }_n$ are distinct, and ${\lambda }_n\ne 0$ for all $n$ . Given $\{f(t_i):i=1,\dots ,2N\}$ , and $\{I_{(k)}(t_i):k=1,\dots ,N,$ $i=1,\dots ,2N\}$ , then $c_n$ and ${\lambda }_n$ are uniquely determined for $n\in \{1,2,\dots ,N\}$ . Furthermore, the solution is elementary and relies on the eigenvalues of a particular companion matrix associated with $\{I_{(k)}(t_i):k=1,\dots ,N\}$ .

What is occurring here is a change of basis in the expression of the original equation from a basis consisting of $I_{(k)}(t)$ and powers of $t$ to a basis consisting of exponential functions and powers of $t$ . The change of basis from the first basis to the new basis gives the values of the vector $({\lambda }_1,\dots ,{\lambda }_N{)}^{\text{T}}$ . Once the $\lambda$ s are determined, equation (1) is a linear equation in $c$ ’s and they can be readily determined. The uniqueness of the $c$ ’s and $\lambda$ ’s is expressed by the following:

Proposition 2

If $r\ge N$ and $t_{2N+1},\dots ,t_{2r}$ are additional points where $f(t_i)$ are specified, then $c_n$ and ${\lambda }_n$ will be the same for $n=1,\dots ,N$ and $c_{N+1},\dots ,c_r$ will be zero.

For the detailed proofs of these results and several motivating examples, the reader is referred to Ellingson and Jafari. ⁴ Please note that the proofs presented in Ellingson and Jafari ⁴ extend far beyond applications of this method to RPT, and the numerical considerations of that article delve deeper into the algorithmic nuances of this method. Relating the $c$ ’s and $\lambda$ ’s to the rate constants will be described below.

Determining the Rate Constants

Returning to the multicompartmental models presented in the introduction for PBPK models, and generalizing them appropriately to the case of $N$ compartments, the transfer of pharmaceutical between these compartments is given by a system of $N$ differential equations, similarly to those described in equation (2). We make the following assumptions for our model:

Homogeneity. Within each compartment, there is no variation in pharmaceutical concentration.

Linear transfer. The rate of pharmaceutical transfer between compartments is proportional to the concentration in the originating compartment.

Conservation. When the pharmaceutical moves between compartments, no additional pharmaceutical is created or destroyed. However, each compartment can have an excretion term which has its rate determined in the same linear fashion as for transfers. Furthermore, we assume that the compartment will not produce the pharmaceutical natively (This assumption implies that the transfer and excretion coefficients are nonnegative).

Mixing. The RP under consideration is thoroughly mixed with its endogenous analogs.

Since we are approaching this through the lens of nuclear medicine, we will no longer consider the concentrations of a pharmaceutical and instead consider the activity of an RP. These have a proportional relationship, so such a change is merely a rescaling of the system.

Let $A_i$ , $i=\{1,\dots ,N\}$ , be the activity in compartment $i$ , and let the transfer coefficient for activity transfer from compartment $i$ to compartment $j$ be denoted by $k_{ij}$ , $i\ne j$ . The excretion coefficient from compartment $i$ is ${\gamma }_i$ . Using this notation, equation (2) can be written as

$$\frac{\text{d}}{\text{d}t}{A}_i(t)=-({\gamma }_i+\sum _{j=1,j\ne i}^N{k}_{ij})A_i(t)+\sum _{m=1,m\ne i}^N{k}_{mi}{A}_m(t)+F_i(t),i=1,\dots ,N,$$

(5)

with the initial condition

$$A_i(0)=A_{i0}\ge 0,$$

(6)

where $F_i$ represents the forcing function acting on compartment $i$ . Representing this in vector form, we have

$$\frac{\text{d}}{\text{d}t}\mathbb{A}(t)=-\mathrm{\Gamma }\mathbb{A}(t)-(K\mathbb{1})\circ \mathbb{A}(t)+K^{\text{T}}\mathbb{A}(t)+\mathbb{F}(t),$$

(7)

where $\mathrm{\Gamma }$ is the diagonal matrix with ${\gamma }_i$ as its nonzero entries, $K$ is the hollow matrix where the entry in the $i$ th row and $j$ th column is $k_{ij}$ and has zero entries on the diagonal, $\mathbb{1}$ is the column vector consisting of all $1$ s, and $\circ$ is the element-wise multiplication operator. This can be written more simply as

$$\frac{\text{d}}{\text{d}t}\mathbb{A}(t)=\mathbb{M}\mathbb{A}(t)+\mathbb{F}(t),$$

(8)

where $\mathbb{M}$ is an $N\times N$ matrix with entries $m_{ii}=-({\gamma }_i+\sum _{j\ne i}{k}_{ij})$ and $m_{ij}=k_{ji},$ where $i\ne j$ . It is worth remarking that if all excretion coefficients are zero ( $\gamma =0$ ), then $\mathbb{M}$ is singular. This represents the case of whole body nuclear imaging since it is atypical for a patient to void during imaging. Non-singularity can be restored by not making the typical correction for nuclear decay. A particular case of (8) for a four-compartment model was given in (2) and (3).

For simplicity of presentation, we first consider the case in which $\mathbb{F}(t)=[\delta (t),0,\dots ,0{]}^{\text{T}}$ and return later to the case of more general forcing functions. For this, we have the solution form

$$A_i(t)=\sum _{j=1}^N{c}_{ij}{e}^{{\lambda }_{j}t},$$

(9)

where ${\lambda }_j$ are the eigenvalues, ${\left[\begin{array}{c}{c}_{1j}\dots c_{Nj}\end{array}\right]}^{\text{T}}$ are associated with the eigenvectors of $\mathbb{M}$ , and $A_i(t)$ is the time evolution of the activity in the $i$ th compartment. Notice that we assume that there are no repeated eigenvalues in this formulation. This is true if and only if the discriminant of the characteristic polynomial of $\mathbb{M}$ is nonzero. Since the set of matrices that have repeated eigenvalues has measure zero, in biological contexts, we can disregard this possibility.

Given the activity measurement, $\mathbb{A}(t)$ , and the algorithm discussed, we numerically arrive at the $c_{ij}$ s and ${\lambda }_j$ s. Next, we have before us a new inverse problem, specifically that of determining $k_{ij}$ s and ${\gamma }_i$ s from the $c_{ij}$ s and ${\lambda }_j$ s. Returning to (8) and substituting the solution arrived at from the nuclear imaging study, we have

$$\left[\begin{array}{c}\sum _{j=1}^N{c}_{1j}{\lambda }_j{e}^{{\lambda }_{j}t}\\ ⋮\\ \sum _{j=1}^N{c}_{Nj}{\lambda }_j{e}^{{\lambda }_{j}t}\end{array}\right]=\mathbb{M}\left[\begin{array}{c}\sum _{j=1}^N{c}_{1j}{e}^{{\lambda }_{j}t}\\ ⋮\\ \sum _{j=1}^N{c}_{Nj}{e}^{{\lambda }_{j}t}\end{array}\right].$$

(10)

At first glance, one might mistakenly think that this is a system of $N$ equations for $N^2$ unknowns. However, using the linear independence of the exponentials, we arrive at $N^2$ equations by separating each equation into $N$ equations by isolating terms bearing the same exponential. Thus, we have

$$c_{ij}{\lambda }_j{e}^{{\lambda }_{j}t}=\sum _{p=1}^N{m}_{ip}{c}_{pj}{e}^{{\lambda }_{j}t}\Rightarrow c_{ij}{\lambda }_j=\sum _{p=1}^N{m}_{ip}{c}_{pj}.$$

(11)

Recalling how $m_{ij}$ are defined, we have

$$c_{ij}{\lambda }_j=-c_{ij}({\gamma }_i+\sum _{p=1,p\ne i}^N{k}_{ip})+\sum _{q=1,q\ne i}^N{c}_{qj}{k}_{qi},\mathrm{\forall }i,j=1,\dots ,N.$$

(12)

By uniqueness of ${\lambda }_j$ s, we have a unique solution for our transfer and excretion coefficients $k_{ij}$ and ${\gamma }_i$ .

Nonzero Forcing Functions

We consider here forcing functions that are relevant to this application within nuclear medicine. Forcing functions that are being considered share the property that, at some point in time, uptake stops. That is to say, there exists some $T$ such that $\mathbb{F}(t)=\mathbf{0}$ , $t>T$ . Gravity infusion, and finite time interval infusion methods have been dealt with in Steiner et al. ³

Only using measurements taken at $t>T$ , we treat the prior dynamics of the system during the time of pharmaceutical uptake (when $\mathbb{F}(t)⪰\mathbf{0}$ ) as merely determining the initial conditions for the system going forward, where the algorithm can be implemented to numerically solve for ${\lambda }_j$ and $c_{ij}$ , which yield $k_{ij}$ and ${\gamma }_i$ by the algebra above. There is an additional step we can take in analyzing the activity function now that the ${\lambda }_j$ and $c_{ij}$ are known:

$$\mathbb{A}(t)-\left[\begin{array}{c}\sum _{j=1}^N{c}_{1j}{e}^{{\lambda }_{j}t}\\ ⋮\\ \sum _{j=1}^N{c}_{Nj}{e}^{{\lambda }_{j}t}\end{array}\right]={\mathbb{A}}_p(t),$$

(13)

where ${\mathbb{A}}_p(t)$ is the particular solution for which the forcing function is responsible. This has application in reducing wait times for nuclear imaging studies since uptake times can be described for various drug delivery systems.

If we wish to also consider forcing functions which include an exponential decay factor in lieu of being zero-valued after some amount of time, this is also acceptable so long as the decay occurs on a shorter time scale than the dynamics of the system to allow for measurements to be taken without significant influence by $\mathbb{F}(t)$ . By the nature of nuclear imaging methodology relying upon counting nuclear decay events, there is a lower bound on activity level contribution from the particular solution where it becomes indistinguishable from quantum noise inherent to counting systems.

Results and Validation

To demonstrate the validity of this technique, we propose an example instance of the forward problem and, from the forward problem solution alone, reconstruct the system, that is, solve the inverse problem.

Suppose that the following randomly chosen $K$ and $\gamma$ for a four-compartment system in Figure 1 are given.

$$K=[k_{ij}]=\left[\begin{array}{cccc}0& 0.2& 0.3& 0.3\\ 0.1& 0& 0.4& 0.4\\ 0.4& 0.3& 0& 0.2\\ 0.1& 0.3& 0.2& 0\end{array}\right],\gamma =[{\gamma }_i]=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right]$$

(14)

with initial condition $\mathbb{A}(0)={\left[\begin{array}{cccc}1& 0& 0& 0\end{array}\right]}^{\text{T}}$ in place of a forcing function.

Using these values in (8), we have

$$\frac{\text{d}}{\text{d}t}\mathbb{A}(t)=\left[\begin{array}{cccc}-0.8& 0.1& 0.4& 0.1\\ 0.2& -0.9& 0.3& 0.3\\ 0.3& 0.4& -0.9& 0.2\\ 0.3& 0.4& 0.2& -0.6\end{array}\right]\mathbb{A}(t)+\mathbb{F}(t).$$

Solving the forward problem (see, e.g. Steiner et al. ³ ), this yields the model TACs in Figure 2.

Figure 2.

Time Activity Curves for Four-Compartment Model.

Sampling these curves on the interval $[0,1]$ with step size $0.05$ , we create estimations for $I_{(k)}(t)$ using the trapezoid rule. Following the method above, we reconstruct the rate constants as $\hat{K}$ and $\hat{\gamma }$ .

$$\hat{K}=\left[\begin{array}{cccc}0& 0.2000& 0.3000& 0.3000\\ 0.1004& 0& 0.3996& 0.3995\\ 0.4001& 0.3001& 0& 0.1999\\ 0.0997& 0.2996& 0.2003& 0\end{array}\right],\Vert \hat{\gamma }{\Vert }_2=6.1642\times {10}^{-5}$$

(15)

The reader should compare the values in the matrices (14) and (15) entry-wise. The resulting relative error from $\hat{K}$ to $K$ is

$$\frac{\Vert \hat{K}-K{\Vert }_2}{\Vert K{\Vert }_2}=9.8693\times {10}^{-4}.$$

This agreement is encouraging and repeatable with randomized $K$ s and $\gamma$ s. In 10,000 trials for the four-compartment model with $k_{ij}$ chosen as independent random variables uniformly distributed between $0$ and $1$ and $\gamma ={\left[\begin{array}{cccc}0& 0& 0& U(0,1)\end{array}\right]}^{\text{T}}$ , the error and relative error histograms are tabulated in Figure 3(a) and (b), respectively. This result can be improved by increasing the sampling density of the activity curves and by regularization of the methods described above.

Figure 3.

Comparison of $\hat{K}$ Errors and Relative $\hat{K}$ Errors Over 10,000 Simulations, with the Time Activity Curve (TAC) Sampled Uniformly 100 Times Over $[0,1]$ .

The choice of sampling interval also affects the agreement between the model parameters in an evocative manner. Returning to the model parameters we first considered in (14), we once again sample on the interval $[0,1]$ but with variable step sizes. Considering uniform sampling on this interval beginning with eight samples ( $2N$ ) and increasing to 100. Figure 4 shows how the relative error in determination of $\hat{K}$ depends on the number of the activity data points in $[0,1]$ sampled uniformly, and Table 1 shows the corresponding errors where the sampled activity data points are chosen randomly in $[0,10]$ using uniform random distribution of data points. Figure 5 shows $\hat{K}$ relative error histogram in 100 simulations where 250 random uniformly distributed time data points on the interval $[0,10]$ are used.

Figure 4.

Convergence of $\hat{K}$ Relative Error with Increasing Number of Samples on the $[0,1]$ Interval.

Figure 5.

$\hat{K}$ Relative Rrror Histogram in 100 Simulations Where the Activity is Sampled at 250 Random Uniformly Distributed Time Data Points on the Interval $[0,10]$ .

Table 1.

$\hat{K}$ Relative Rrror on Nonuniformly Sampled Time Points on the Interval $[0,10]$ Using Uniform Random Variable Distributions.

Sampled Time Points	Relative Error
20	$0.035483$
50	$0.00064424$
100	$0.00023338$
250	$1.3029\times {10}^{-6}$
500	$7.2717\times {10}^{-9}$

An important point arises here that needs clarification. In the nuclear medicine setting, acquiring the activity in an ROI at different times often requires the patient to return to the lab and have the activity re-measured. Hence, at best, one may expect to have three to five data points along the TACs. Table 1 gives the relative errors in $\hat{K}$ for 20–500 data points! That is certainly not practical. The reader should note that our samples in the different figures are distributed in time intervals $[0,1]$ , $[0,4]$ , or $[0,10]$ , and the TAC is generated over 1-min, 4-min, and 10-min intervals. Theorem 1 proves that $2N$ data points are sufficient for the complete reconstruction of the parameters in (1) and the rate constants $K$ ; however, as the integrals $I_{(k)}(t)$ in (4) may be somewhat inexact with very few data points, basing these integrations on more data points improves the estimations. Extracting the TAC acquired in a short time period (e.g. 5 min) in list mode and sampling it as desired will permit estimation of the rate constants very accurately. Thus, there is no need for the patient to return and cause an interruption in the workflow of the clinic.

Discussion

Comparison of the parameters in (14) and their estimation in (15) shows that the described algorithm faithfully determines the rate constants from activity measurements. The relative LSEs in the 16 parameters are $9.8\times {10}^{-4}$ . Figure 3(a) and (b) shows this fidelity repeated in 10,000 numerical simulations where the $k_{ij}$ and ${\gamma }_i$ are randomized. Figure 4 shows the convergence of the method as the sampling of the modeled TAC is improved. These numerical tests demonstrate the accuracy and robustness of the method.

Data acquisition remains a limiting factor for diagnostic accuracy and, therefore, clinical utility. This is illustrated in Table 1. These limitations arise primarily from estimations of the integrals $I_{(k)}(t)$ when few data points are used, but they can be significantly improved using cubic splines or other interpolation methods. An example of this is shown in Figure 6. While Figure 6 shows the reconstruction of $\lambda$ s and $c$ s, the error due to sampling is in the reconstructions of these parameters. Deriving the rate constants from these parameters involves solving a well-conditioned system of linear equations, and is only limited by the machine precision (e.g. $1\times {10}^{-16}$ ) and the noise in the data. Figure 6 compares interpolation methods, which are a necessary step for accurate estimation of $I_{(k)}(t)$ s. With a large number of samples of the function, cubic spline interpolation yields superior agreement, while in the low sampling domain ( $⪅20$ ), linear interpolation provides better results.

Figure 6.

Relative Error in Parameter Estimates Using Linear Interpolation (Trapezoidal Rule) (a) and Using Cubic Spline Interpolation (b).

The novel mathematical approach presented here has numerous immediate and potential applications in imaging modalities. The radiologist may feel that the computational methods currently utilized are sufficient for clinical purposes, particularly with relatively simple interpolations or standard reconstructions that provide both a semi-quantitative and overall qualitative characterization of a particular pathology when interpreted within a clinical picture. However, this result is of immediate benefit in the acquisition of nuclear medicine images. Until now, the acquisition was dependent on photon counts reaching a sufficient threshold as a function of time, balancing radiation dose and signal-to-noise ratios with a reasonable (physiological and comfort) duration under the camera. Simple linear or spline interpolation is utilized to assist in data processing. The method presented here allows for a reduction in RP dose by tailoring the dosage to the individual patient in a quantitatively precise manner.

Nuclear medicine provides an additional dimension of patient physiology by utilizing receptor- or function-specific ligands or molecules bound with radioactive agents for imaging, treatment, or both. PET and SPECT technology also provide higher clinical certainty with localization that will benefit from the aforementioned methods. In addition, whole-body imaging in PET and SPECT acquisitions can show the early flow phases as the radiotracer is distributed, reaches the target and critical organs, and decays or is excreted physiologically. Each patient’s physiology may be unique and the method can profile their specific physiologic distribution to optimize the study protocol to answer the clinical question.

The next logical step in this effort toward personalized medicine is to create a model with parameters unique to the patient that determines the RP kinetics of a small dose delivered to the patient to optimize the dose of a theranostic agent to the index tumor. Further administrations and overall protocols can be determined based on radiologic response, with further research required to evaluate treatment protocols and long-term outcomes. Furthermore, the parameters from this patient-specific model can constitute a new set of quantitative imaging biomarkers that might allow for improved characterization of both disease and natural variation. This high-dimensional data can be studied using machine-learning methods and incorporated into the features of neural networks.

PKs within other modalities, including photon-counting CT, MRI, and contrast-enhanced mammography, can further characterize a pathology by relating to and confirming with histological studies. Specific ROIs, such as index tumors, can be targeted for this purpose under each modality to which the above method can be applied to provide subcharacterizations with opportunity for improved prognostic accuracy.

Applying such a model to tracer kinetics may require some additional discussion. The problem of radioactive decay is a nonissue; nuclear medicine systems already make corrections for the physical half-life of the tracer. One might reasonably point out how the transfer coefficients may vary according to physiological conditions. However, within small time intervals, it is reasonable to assume that the transfer coefficients are relatively stable. We make the critical observation of an application where the relative stability gives us the opportunity to describe how these transfer coefficients vary according to physiological conditions. This could yield a phenomenological account of how transfer coefficients vary under differing physiological conditions. Examples of such phenomena that might alter these transfer coefficients include tissue saturation and altered metabolism due to disease. A further generalization that will be considered in future research is relaxing the homogeneity assumption by accounting for spatial distribution within compartments.

A concern some readers might have regarding the clinical viability of the technique described is the number of data points required. For the analytical method to determine a unique solution, $2N$ points are required, where $N$ is the number of compartments considered. In the context of nuclear medicine, this is not a problem. As data are collected and recorded in list mode, thereby time-stamping detection events, within a single scan, there are many data points collected, thus easily satisfying the requirement.

The limitations of this method may be reduced to uncertainties in the parameter estimations due to system noise. Additional mathematical work is being performed to characterize the outcome of this method in noisy systems and the sensitivity of the estimations to noise.

Advantages Over Existing Methods

Optimization methods for parameter fitting must contend with the nonconvexity and high-dimensionality of the problem posed. Convergence to local minima would yield high agreement between the observed and simulated TAC but disagreement between the ground truth parameter values and those found by the optimization procedure. In this event, the physiological significance of these parameters would be lost. Additionally, optimization methods are computationally intensive.

The MPm of Sarkar et al. ¹³ requires uniform sampling in the time domain. This makes the method impractical for certain long-term dynamics in this setting. In addition, MPm involves singular value decomposition of the Hankel matrix that is constructed from the dataset, which can be computationally intensive for large datasets. In addition, in the case of short data length, MPm can suffer from spectral leakage that would affect the accuracy of the parameter estimations. It is well-known that this method is also sensitive to noise and is poorly conditioned in certain settings; however, we will postpone a comparison of the two methods in the noisy setting to an upcoming article.

The analytical method described overcomes many of these challenges by directly solving for the parameters. This allows for the study of multicompartmental models with many more compartments than were able to be previously considered.

Variable Transfer and Excretion Coefficients

We alluded to the notion of variable transfer and excretion coefficients earlier. Under Michaelis–Menten kinetics, the flux rate of activity between containers $i$ and $j$ becomes

$$k_{ij}{A}_i(t)⟶\frac{{\varphi }_{ij}^{max}}{{\mu }_{ij}+A_i(t)}{A}_i(t),$$

(16)

where ${\varphi }_{ij}^{max}$ is the maximum flux rate from containers $i$ to $j$ and ${\mu }_{ij}$ is the Michaelis constant for the flux. Notice that for $A_i(t)\ll {\mu }_{ij}$ flux rate is almost linear with $A_i(t)$ . Therefore, in this lower limit, the methodology for the strictly linear system also applies here, but instead of determining ${\varphi }_{ij}^{max}/{\mu }_{ij}$ as $k_{ij}$ .

To the radiologist, the problem solved is one that initially might seem clinically unimportant since PET data are sufficiently abundant such that linear or spline interpolation provides all the qualitative information they need. However, the accuracy and precision of the PK estimates described allow for a confident approach to patient care and creation of individualized protocols.

Conclusion

We conclude with some immediate observations. The methods discussed here emphasize that the choice of sampling the TAC is completely arbitrary in determining the rate constants associated with compartmental models. The algorithm and the data presented show that the accuracy and precision of recovering the rate constants are excellent, whether the data are sampled randomly or uniformly, and how it relates to the number of data points sampled. The appearance of the iterated integrals or the moments of $f(t)$ in intervals $[0,t_i]$ , $i=1,\dots ,2N$ , introduces sampling and statistical considerations in applied problems, where the data may be discrete or quantized. Expressing the iterated integrals as polynomials in terms of the moments (see Ellingson and Jafari, ⁴ Proposition 2) opens a new connection to be exploited much further. The vast literature on moments (see Akheizer, ¹⁷ Dragotti et al., ¹⁸ and the references therein) will undoubtedly prove useful in extending these methods.

We have already performed simulations by adding statistical noise to the TAC readings and reconstructed the rate constants from these noisy data. Interestingly, where the sampling is relatively dense, the results are consistent with the results presented here and fit the randomly chosen rate constants quite well; not unexpectedly, for coarser sampling regimens, the errors become somewhat larger. The analysis of the statistical uncertainties associated with the noisy activity data is proving to be quite extensive. Once the results on the statistics of the problem are completed, applications of these models to patient data will be center stage to bring these results into clinical practice. Without these models, applications to patient data are unverifiable.

Finally, a natural extension of this technique is dealing with the situation where we do not know the exact number of exponential basis functions in the sum. The over-determined problem has already been considered in Proposition 3 in Ellingson and Jafari, ⁴ while the under-determined problem will require optimization over the set of parameters $c_n$ and ${\lambda }_n$ . Many interesting questions arise in this context, whose treatment will require significant mathematical analysis and will be highly valuable in patient care.

Several topics for future work have already been discussed above and in Section “Dicussion.” It can be useful to collect these topics here. Of course, this is by no means an exhaustive list.

The homogeneity assumption within each ROI will directly impact the calculation of RP dosage and the organ radiation dose. Introducing diffusion and other partial differential equation-based models that account for RP distribution in each ROI would improve the dose estimates. This is an interesting topic and one that is worthy of pursuing.

Introducing the patient-specific PK model parameters as a physiologically meaningful quantitative imaging biomarker into machine-learning methods for diagnostics, therapeutics, and prognostics.

Using Michaelis–Menten kinetics, the rate constants are dependent on the concentration of activity in each compartment. These models lead to saturation of activity/concentrations and to nonlinear systems of differential equations, whose solutions may not be sums of exponentials. Fractional calculus models have been used by various authors to model biological phenomena. The reader is referred to Khan et al.,^19,20 and the references therein for inquiries in this direction.

Finally, the study of sensitivity and stability of the estimates of the rate constants to noise and other nonstochastic parameters (e.g. patient motion) is interesting to achieve full clinical birth of these models.

Abbreviations

LSE

Least square error

PK

Pharmacokinetic

PET

Positron emission tomography

SPECT

Single photon emission computed tomography

TAC

Time activity curve

TIA

Time-integrated activity

PBPK

Physiologically based pharmacokinetic

SUV

Standard uptake value

RP

Radiopharmaceutical

RPK

RP kinetics

RPT

RP therapy

ROI

Region(s) of interest

MPm

Matrix Pencil method