Dosimetric analysis of a compartmental model for radioligand uptake in tumor lesions

Radioligand therapy is a targeted cancer therapy that delivers radiation to tumor cells based on the expression of specific markers on the cell surface. It has become an important treatment option in metastasized neuroendocrine tumors and advanced prostate cancer. The analysis of absorbed doses in radioligand therapies has gained much attention and remains a challenging task due to individual pharmacokinetics. As an alternative to the often used sum of exponential functions in intra-therapeutic dosimetry, a basic compartmental model for the pharmacokinetics of radioligands is described and analyzed in this paper. In its simplest version, the model behavior is determined by the uptake capacity and the association constant and can be solved analytically. The model is extended with rates for excretion from the source compartment and externalization from the lesion compartment. Numerical calculations offer an insight into the quantitative effects of the model parameters on the absorbed dose in the tumor lesion. This analysis helps understanding the importance of clinically relevant factors, e.g. the effect on absorbed doses of modified radioligands that bind to albumin. Using clinical data, the potential application in intra-therapeutic dosimetry is illustrated and compared to the bi-exponential function which lacks a mechanistical basis. While the compartmental model is found to constitute a feasible alternative in these examples, this has to be confirmed by further clinical studies.


Introduction
In recent years, radioligand therapy has become an important treatment option for metastasized neuroendocrine tumors and advanced prostate cancer [1,2]. This treatment is based on a radiolabeled ligand that selectively binds to markers on tumor cells and is internalized into the cell. It thereby enables the selective transportation of a therapeutic radioisotope to the tumor cells. The target depends on the tumor entity and is represented by the somatostatin receptor for neuroendocrine tumors and by the prostatespecific membrane antigen (PSMA) for prostate cancer. The physical properties of the β-emitter 177 Lu have made it a widely used radioisotope in radioligand therapies [3]. As the therapeutic effect of radioligand therapy is affected by the absorbed dose in the tumor lesions, the analysis of absorbed doses in radioligand therapy is of particular interest [4]. The absorbed dose is determined by several factors, e.g. the radioligand uptake capacity of the tumor cells, the association kinetics of the radioligand with the target and the externalization rate from the tumor cells. In intratherapeutic dosimetry, the absorbed dose in tumor lesions is usually estimated by fitting a time-activity curve to a time series of scintigraphic imaging after administration of the radioligand [5]. Due to requirement of substantial effort and resources, the time series is generally limited to a few measurements [6]. The curve fitting is usually performed with a sum of exponential functions which lack a biological basis [5]. On the other hand, several detailed models with many compartments and interactions have been proposed to study effects of model parameters on absorbed doses [7][8][9][10]. This will be outlined in detail in section 4.3 . These models are in need of comprehensive knowledge about the pharmacokinetics of the respective radioligand which often limits their applicability in Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. intra-therapeutic dosimetry. Furthermore, the analysis is often hindered by the complexity of the model and the effect of a single parameter on the model behavior can be difficult to determine. Based on this disparity between the sum of exponential functions and the comprehensive models, the aim of this paper is to describe and systematically analyze a simplified model that is reduced to the basic pharmacokinetic mechanisms of radioligand therapy and could be of help in intra-therapeutic dosimetry.
This analysis is structured as follows: In section 2, a basic compartmental model for radioligand uptake is described and analyzed with respect to the absorbed dose in tumor lesions. In section 3, the effects of an excretion rate from the source compartment and externalization rate from the lesion compartment on the absorbed dose are analyzed. In section 4, potential clinical applications and limitations are discussed.

Model development
The model consists of two compartments. Compartment X is the radioligand source compartment and can be thought of as the surrounding compartment of the tumor lesion or the blood compartment, depending on the specific situation to be simulated. Compartment Y represents the tumor lesion. The compartment variables x and y represent the amount of radiopharmaceutical in the X and Y compartments, respectively. Homogenous distributions of the radioligand in the compartments are assumed. The variables x and y do not account for radioactive decay, and it is therefore assumed that the pharmacokinetics of the radiopharmaceutical is not affected by the radioactive decay. The compartment variables are normalized to the total amount of added radioligand. For all simulations, initial conditions are set to x(0) = 1 and y(0) = 0, i.e. in the beginning (time t = 0) all of the added radioligand is located in compartment X. For this reason, compartment X is also referred to as source compartment. It is worth noting that using variables describing the amounts of radioligand and not the concentration, as well as scaling the variables and the model parameters to the amount of added radioligand does not lead to a loss of generality. It simplifies numerical calculations and understanding in this analysis, which focuses on relative changes in model parameters. The variable r represents the amount of free receptors for the ligand on the surface of the tumor lesion cells normalized to the total amount of radioligand in the system. The association of the radioligand with the receptor is assumed to be a bimolecular reaction following second order kinetics. The resulting association rate is therefore proportional to the product of the receptor density on the tumor surface and the radioligand concentration in compartment X. Since both tumor surface and volume are held constant, the association rate is proportional to the total amount of receptors and ligands and can be expressed as kxr with an association constant k. Dissociation of the radioligand-receptor complex is not considered in this model. It is assumed that every receptor-ligand complex is instantaneously internalized into the cell, i.e. that the amount of receptor-radioligand complexes equals the internalized amount of radioligand y. The instantaneous internalization is certainly a debatable assumption, but as will be shown in subsection 4.1, it is a reasonably good approximation for some radioligand therapies. Due to internalization, the amount of receptors on the tumor decreases linearly with the amount of internalized radioligand-receptor complex equaling y. Receptor reproduction or recycling are not part of this model, as these mechanisms could be of limited importance on the relatively short timescale of radioligand uptake [11]. By defining m as the maximum amount of receptors on the tumor lesion normalized to the added amout of radioligand, the expression for the free receptors is given by r = m − y.
The parameter m will also be referred to as uptake capacity. In this first model, externalization from the tumor compartment and excretion from the source compartment are not taken into consideration. Summing up, the above considerations result in the following model with two compartments and a rate from the source compartment X to the tumor compartment Y: In this simplest version, the model is similar to a typical receptor-ligand model [12]. Once a solution y (t) is obtained, the time-activity curve for the tumor lesion a(t) can be calculated by multiplying y(t) with a radioactive decay factor and a constant factor α that is needed for calibration and conversion into units of activity: where λ is the exponential decay constant.

Solution
It is clear from equations (1) and (2) that the sum x + y is a constant, which equals 1 according to the initial conditions. By replacing y in equation (1) with the expression 1 − x, the two-dimensional system of differential equations can be transformed into a onedimensional equation: In the case of m = 1, the non-trivial solution meeting the initial condition x(0) = 1 equals 1/(1 + kt). Otherwise, the solution to this Bernoulli differential equation equals: To meet the initial condition x(0) = 1, the constant c is set to −m resulting in the following solution for the two-dimensional system: As expected, the compartment variable y approaches 1 if m 1, meaning that all added radioligand will accumulate in the tumor compartment if the uptake capacity is larger than the total amount of radioligand. Otherwise, the tumor compartment variable y approaches the uptake capacity: This effect of the uptake capacity on the uptake in the tumor compartment can also be seen in figure 1. Figure 2 shows how a higher association constant leads to faster uptake in the tumor compartment.

Calculation of the absorbed dose
The absorbed dose can be calculated according to the MIRD formalism [13]. If cross-absorbed dose from the source compartment is considered negligible and the S-value, defined as absorbed dose rate per unit of activity with units Gy/(Bq×s), assumed to be time independent, the expression for the absorbed dose  results in an expression proportional to the integrated time-activity curve: As this analysis focuses on relative changes of the absorbed dose and not on the absolute value, a new quantity D r ≔ D/(αS) is introduced. This quantity is proportional to the absorbed dose. To simplify reading flow, D r will hereafter be referred to as absorbed dose. The subscript will indicate the relative character of the quantity. It follows that: To evaluate this improper integral, it can be split into two parts: where b is a positive real number. The second term represents an additive error term and can be minimized by choosing a sufficiently large value for b. The upper limit of the uptake in the tumor lesion equals the total amount of added radioligand. Consequently, the error term fulfills the following inequality: Thus, to evaluate the absorbed dose with an error smaller or equal to f, the upper limit of the integral can be set to . An f of 0.1 will be used for all following calculations. Simpson's rule is applied for approximation of definite integrals. Figure 3 shows a contour plot visualizing the dependence of the absorbed dose D r on the uptake capacity and the association constant k. It is clear that the absorbed dose increases with a higher value of k, which signifies a faster uptake of radionuclide into the cells. However, the increase of the absorbed dose due to a faster uptake is limited by the fact that it cannot exceed the absorbed dose obtained in the case of instantaneous uptake. Thus, variations of the association constant have no substantial impact on the absorbed dose if the time scale of radioligand uptake is much faster than the time scale of radioactive decay, which applies to the example in subsection 4.1. A similar relation can be observed for the uptake capacity and the absorbed dose. The relation between the absorbed dose and the uptake capacity is nearly linear in the range from m = 0 to m = 1. In the rather unlikely case that the uptake capacity m is larger then the total amount of added radioligand, the absorbed dose is limited by the finite amount of radioligand.
3. Extending the model 3.1. Excretion from the source compartment Several studies have found that a large fraction of the administered radiopharmaceutical in radioligand therapies is excreted from the body during the first hours after administration [14,15]. Thus, it is worth analyzing the quantitative effect of an excretion rate from the source compartment on the absorbed dose in the tumor lesion. Assuming a first-order elimination from the source compartment, the extended model is described by the following equations: A solution to this system of differential equations cannot be given as a simple analytical expression. Therefore, a numerical approach was used applying the 4th order Runge-Kutta method [16]. The step size was set to 0.1 h which allows for sufficient accuracy given that radioligand uptake dynamics usually occurs within a time frame of multiple hours [17,18]. Some properties that affect the absorbed dose in the tumor compartment can be derived from the contour plots 4 and 5. First of all, it is clear that the excretion rate from the source compartment reduces the absorbed dose. Furthermore, it is intuitively clear that an increased excretion rate can be compensated by increasing the association constant to reach the same absorbed dose. Figure 4 illustrates this relation between the excretion and association constant and the absorbed dose in the case of an uptake capacity of 0.5. Interestingly, the isolines of this contour plot are almost straight lines indicating that a nearly fixed ratio of increments in u and k is necessary to stay on a given isoline. The isolines in figure 5, which visualizes the effect of the excretion constant and the uptake capacity on the absorbed dose in the case of k = 0.15 h −1 , can be explained as follows: If the uptake capacity is smaller than the added radioligand, i.e. m < 1, slow excretion from the source compartment will not prevent the uptake in the tumor compartment from approaching the uptake capacity m. Thus, the excretion will not have a large impact on the absorbed dose if the absorbed dose is mainly determined by the amount of radioligand uptake. But once the excretion is fast enough to significantly reduce the amount of radioligand uptake, the absorbed dose decreases substantially as indicated by the curves of the isolines.

Externalization from the tumor lesion
Many previous analyses indicate that the externalization of the radiopharmaceutical from the tumor lesion is an important determinant of the absorbed dose [17,19,20]. As the exact molecular mechanism of externalization is not fully known, the mathematical formalism has to rely on some assumptions, and the presented model is therefore just one possibility of including externalizaton. Based on the results in [17], it is assumed that the externalization follows a firstorder kinetics. If, in addition, the simplification that the externalized radioligand does not recirculate to the  source compartment is applied, the model is described by the following differential equations: Taking into account the model development in section 2.1, it is clear that equations (15) and (16) include a mechanism of receptor recycling, as receptors are represented on the lesion surface if the amount of internalized radioligand decreases. Figure 6 shows that the quantitative effect of the externalization constant s is heavily dependent on the uptake capacity m. If the uptake capacity is low, the externalization can be compensated to some extent by the uptake of the remaining radioligand from the source compartment. Otherwise, even a small externalization rate cannot be compensated and leads to a decreased absorbed dose.

Comparison with an in-vitro experiment
Before proceeding to possible clinical applications, the presented model in its simplest version as defined by equations (1) and (2) is compared with results from a cell-based experiment published in 2021 by Mansi et al [17]. In this experiment, the radioligand binding and internalization of human embryonic kidney (HEK)-293 cells stably transfected with human somatostatin receptor 2 (SST 2 ) and seeded in 6-well plates were analyzed. For this purpose, the cells were incubated at 37°C in a medium containing an SST 2 peptide agonist, namely [ 177 Lu]Lu-DOTA-TATE (2.5 nM). A blocking ligand was used to distinguish between specific and non-specific binding and uptake. After 0.5, 1, 2 and 4 hours, the amounts of specific cell surface-bound and internalized radioligand were determined and expressed as a percentage of the total applied activity.
As can be seen in figure 7, the surface-bound fraction of radioligand was minimal (< 5 %). The assumption of instantaneous internalization of receptor-radioligand complexes is therefore a reasonable approximation of reality in this case. Moreover, it shows that the proposed model fits the amount of internalized radioligand well. Obviously, this is just one example and comparisons with other radioligands and receptors are necessary, but it is an important one as [ 177 Lu]Lu-DOTA-TATE is the first approved radiopharmaceutical for treatment of neuroendocrine tumors [21]. Time-activity curves obtained from posttherapeutic imaging in animal studies and in humans indicate similar uptake duration of radiolabeled somatostatin analogues or PSMA ligands in tumor lesions compared to the presented cell-based experiment [22,23].

Intra-therapeutic dosimetry
The goal of intra-therapeutic dosimetry is to determine the absorbed dose in tumor lesions and organs at risk during radiopharmaceutical therapy, which is of fundamental importance in the assessment of treatment efficacy and safety. In radioligand therapies, the estimation of absorbed doses in tumor lesions is usually based on scintigraphic imaging at several time points after administration of the radiopharmaceutical [5]. Based on the imaging data, the radionuclide uptake in tumor lesions is quantified. The absorbed dose is then derived from an estimated time-activity curve. Although intra-therapeutic dosimetry is recommended according to several guidelines, routine protocols for the exact procedure in determining the absorbed doses in radioligand therapies are missing [24][25][26]. In particular, several methods for curve fitting of time-activity curves have been suggested.
Most suggested methods are based on fitting a sum of exponential functions, but lack a biological basis [5].
The frequently used bi-exponential function with the biologically essential requirement a(0) = 0 is given by and has three adjustable parameters. Based on the model described in sections 2 and 3, an alternative approach is presented in the following. For this purpose, the model with excretion rates from the source compartment and lesion compartment is used: k m y x sy. 1 9 = + --According to equation (3), the time-activity curve is obtained by multiplying y(t) with the radioactive decay factor and adjusting to units of activity. To overcome overfitting due to few time points in intra-therapeutic dosimetry, it might be useful to reduce the number of adjustable parameters. For example, the excretion rate from the source compartment could be estimated from population pharmacokinetic studies to reduce the number of adjustable parameters from four to three. In what follows, this procedure is by example applied to intra-therapeutic dosimetry data from a study about [ 177 Lu]Lu-PSMA therapy carried out by Peters et al [23]. In this study, ten patients with lowvolume hormone-sensitive metastatic prostate cancer underwent two cycles of [ 177 Lu]Lu-PSMA therapy. Lesion dosimetry was based on SPECT/CT scans at five time points after each therapy. It is known that unbound [ 177 Lu]Lu-PSMA is rapidly cleared from the body via renal excretion. Based on the finding that four hours after administration approximately half of the administered activity is excreted, the rate u is set to a corresponding value of 0.17 h −1 [14]. Using a nonlinear least-squares algorithm, the resulting model can be fitted to the intra-therapeutic dosimetry data [27]. Obviously, the comparison with the uptake dynamics of two index lesions during two therapy cycles does not constitute a proper model evaluation. In particular, more measurements in the first few hours after administration are necessary to evaluate the model in the important initial phase and compare it to other models. In this phase, figure 8 illustrates a notable difference between the presented model and the biexponential curve leading to a larger calculated absorbed dose when using the bi-exponential curve. Furthermore, figure 8 clearly demonstrates the potential advantage of using a mechanistic model. Assuming a correct description of the underlying mechanisms, the fit results can give an insight into physiological parameters. This could be of particular interest in the given situation to study therapy-induced changes in physiological parameters between the first and second therapy cycle. Understanding these mechanisms could be useful to optimize therapy strategies, e.g. by adjusting administered activities.
Another strategy that has gained much attention in recent years is the development of albumin-binding radioligands. The binding to albumin was found to increase the retention time in blood substantially. For example, an analysis of [ 177 Lu]Lu-DOTA-TATE therapy reported a half-life of 4.7 hours corresponding to a decay constant of approximately 0.15 h −1 in the initial fast elimination phase attributed to renal clearance [15].
For the modified radioligand [ 177 Lu]Lu-DOTA-EB-TATE, which binds to albumin, the analysis of the fast initial excretion phase from the blood compartment resulted in a half-life of 9.5 hours  [18]. As illustrated, the presented models can serve as simple tools to study the effects of therapy strategies affecting pharmacokinetic parameters.

Limitations
The simplicity of the presented model comes along with several limitations that could restrict applications. An important limitation is based on the assumption of homogeneous distribution of the radioligand in the compartments. This assumption does not account for important aspects of cellular dosimetry like cell morphology and sub-cellular radioactivity distributions. Several experimental and theoretical studies have illustrated the effect of these aspects on absorbed doses in radioligand therapies [29][30][31]. Another limitation results from the assumption of instantaneous internalization. This assumption does certainly not hold for every radioligand. For example, studies have reported a low internalization rate for somatostatin receptor antagonists [32]. These radiolabeled antagonists have been proposed in recent years as an alternative to radiolabeled agonists [33]. A simulation study incorporating a variable internalization rate in PSMA radioligand therapy showed that the higher the internalization rate, the larger the required radioligand amount for optimal tumor-to-kidney ratios was [7]. Another limiting factor for the proposed model is the lack of a radioligand input function to the source compartment. The simplified, homogeneous source compartment does not account for the complex process of application and circulation of the radioligand through the body in the initial phase. Several studies regarding positron emission tomography (PET) have underpinned the importance of input functions for pharmacokinetic models [34,35]. A further limitation of the proposed model is the restriction to one source compartment and one tumor compartment. Multiple tumor lesions, other important organs in radioligand metabolism or crossirradiation between compartments are not part of the model. In contrast, various models with numerous compartments have been developed to study the pharmacokinetics in radioligand therapies. For example, Hardiansyah et al performed a global sensitivity analysis of a pharmacokinetic model for [ 177 Lu]Lu-PSMA therapy including several compartments and physiological relevant mechanisms [8]. In this model, the receptor density was found to be the most important parameter for absorbed dose in tumor lesions. Similarly, Ittaqa et al performed a global sensitivity analysis of a comprehensive pharmacokinetic model for radioligand therapy in patients with meningioma and found that receptor measurement is important for treatment improvement [9]. These findings about the importance of receptor density are in agreement with the analysis in sections 2.3 and 3. Based on another multicompartmental model, Zaid et al analyzed the pharmacokinetics in therapies with alpha-emitter-labeled somatostatin analogues and compared two radionuclides [10]. In another study, Hardiansyah et al found that a pharmacokinetic model together with simulated dynamic PET could be useful for treatment planning in therapies with radiolabeled somatostatin analogues [36]. Other proposed models aim at predicting the tumor response after radioligand therapy [37]. One example is a study by Kelk et al that found that voxel-based dosimetry could help optimize treatment with [ 177 Lu]Lu-PSMA [38]. Based on these limitations and the increasing knowledge about cellular mechanisms of radioligand metabolism, the proposed model could be extended in several ways. However, it should be noted again that an extension of the proposed model leads to much more complex models that are of limited mathematical tractability and in need of comprehensive knowledge about physiological parameters.

Conclusion
In this paper, a basic compartmental model for radioligand pharmacokinetics including a finite uptake capacity in the tumor lesion was analyzed. The simplified model allowed for a thorough analysis of the model dynamics and the effects of parameters on the absorbed dose in the tumor lesion. Based on examples, the possible application in intra-therapeutic dosimetry as an alternative to the often used sum of exponential functions was illustrated. However, for a proper evaluation of the model and its clinical applications this theoretical analysis has to be followed by clinical studies. The analyzed effects of the model parameters on the absorbed dose can be of particular interest for optimization of radioligand therapy strategies, e.g. in understanding the change in absorbed doses when using modified radioligands that bind to albumin. This knowledge will become even more valuable with the increasing use of radioligand therapies. Moreover, with the introduction of dynamic whole-body PET the value of kinetic models for radioligands has substantially increased in diagnostic procedures according to the theranostic principle [39].

Data availability statement
No new data were created or analysed in this study.