
==== Front
Z Med Phys
Z Med Phys
Zeitschrift für Medizinische Physik
0939-3889
1876-4436
Elsevier

S0939-3889(23)00007-7
10.1016/j.zemedi.2023.01.007
Original Paper
Population-based model selection for an accurate estimation of time-integrated activity using non-linear mixed-effects modelling
Hardiansyah Deni denihardiansyah@ui.ac.id
ab⁎
Riana Ade a
Eiber Matthias c
Beer Ambros J. d
Glatting Gerhard de
a Medical Physics and Biophysics, Physics Department, Faculty of Mathematics and Natural Sciences, Universitas Indonesia, Depok, Indonesia
b Research Collaboration Center for Theranostic Radiopharmaceuticals, BRIN, Bandung, Indonesia
c Department of Nuclear Medicine, Klinikum rechts der Isar, Technische Universität München, München, Germany
d Department of Nuclear Medicine, Ulm University, Ulm, Germany
e Medical Radiation Physics, Department of Nuclear Medicine, Ulm University, Ulm, Germany
⁎ Corresponding author: Dr. sc. hum. Deni Hardiansyah, Physics Department, Faculty of Mathematics and Natural Sciences (FMIPA), Universitas Indonesia, Depok 16424, Indonesia. denihardiansyah@ui.ac.id
20 2 2023
8 2024
20 2 2023
34 3 419427
2 8 2022
13 1 2023
© 2023 The Author(s)
2023
https://creativecommons.org/licenses/by/4.0/ This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).
Purpose

Personalized treatment planning in Molecular Radiotherapy (MRT) with accurately determining the absorbed dose is highly desirable. The absorbed dose is calculated based on the Time-Integrated Activity (TIA) and the dose conversion factor. A crucial unresolved issue in MRT dosimetry is which fit function to use for the TIA calculation. A data-driven population-based fitting function selection could help solve this problem. Therefore, this project aims to develop and evaluate a method for accurately determining TIAs in MRT, which performs a Population-Based Model Selection within the framework of the Non-Linear Mixed-Effects (NLME-PBMS) model.

Methods

Biokinetic data of a radioligand for the Prostate-Specific Membrane Antigen (PSMA) for cancer treatment were used. Eleven fit functions were derived from various parameterisations of mono-, bi-, and tri-exponential functions. The functions' fixed and random effects parameters were fitted (in the NLME framework) to the biokinetic data of all patients. The goodness of fit was assumed acceptable based on the visual inspection of the fitted curves and the coefficients of variation of the fitted fixed effects. The Akaike weight, the probability that the model is the best among the whole set of considered models, was used to select the fit function most supported by the data from the set of functions with acceptable goodness of fit. NLME-PBMS Model Averaging (MA) was performed with all functions having acceptable goodness of fit. The Root-Mean-Square Error (RMSE) of the calculated TIAs from individual-based model selection (IBMS), a shared-parameter population-based model selection (SP-PBMS) reported in the literature, and the functions from NLME-PBMS method to the TIAs from MA were calculated and analysed. The NLME-PBMS (MA) model was used as the reference as this model considers all relevant functions with corresponding Akaike weights.

Results

The function f3a=A1e-λ1+λphyst+A2e-λphyst was selected as the function most supported by the data with an Akaike weight of (54 ± 11) %. Visual inspection of the fitted graphs and the RMSE values show that the NLME model selection method has a relatively better or equivalent performance than the IBMS or SP-PBMS methods. The RMSEs of the IBMS, SP-PBMS, and NLME-PBMS (f3a) methods are 7.4%, 8.8%, and 2.4%, respectively.

Conclusion

A procedure including fitting function selection in a population-based method was developed to determine the best fit function for calculating TIAs in MRT for a given radiopharmaceutical, organ and set of biokinetic data. The technique combines standard practice approaches in pharmacokinetics, i.e. an Akaike-weight-based model selection and the NLME model framework.

Keywords

Akaike weight
Model selection
NLME
PSMA
==== Body
pmc1 Introduction

The estimation of individual absorbed doses in molecular radiotherapy (MRT) is generally helpful for therapy optimisation [1], [2], [3], [4]. For example, the EC Directive 2013/59/Euratom states that absorbed doses (requiring time-integrated activity TIAs as input) in nuclear medicine therapies shall be planned individually [5]. A critical unresolved issue in MRT dosimetry is the question of which fit function to use for the calculation of the TIA, as the obtained TIA value highly depends on the chosen fit function [6], [7]. Therefore, determining the “optimal” fit function is important for accurately calculating TIAs and, subsequently, the absorbed dose. The “optimal” fit function can be selected from a set of suitable models using a model selection procedure [8]. Relevant criteria for the best fit function are that the goodness of fit is satisfactory and that this function is most supported by the data [7]. The latter requires model (or function) selection that could be done using the Akaike Information Criterion (AIC) method [9], [10].

The measurement of biokinetic data is often done at only a few time points in the clinic. The selection of a fit function for an individual patient, i.e. individual-based model selection (IBMS), based on a low number of data, allows only the investigation of fit functions having few parameters, which, therefore, may not correctly reflect the kinetics of the radiopharmaceutical [6]. In addition, model selection has a high uncertainty for a low ratio of the number of data to the number of adjustable parameters. Therefore, population-based model selection is of interest as it can increase the ratio between the number of input data and the number of parameters, which could lead to lower uncertainty in model selection. Recently, we developed a population-based model selection method with a shared function parameter (SP-PBMS) that could lead to a better performance than the IBMS method [6]. The fitting was performed population-wise, including functions with a parameter shared between patients.

Another population-based approach used for TIAs calculation is the non-linear mixed-effects (NLME) modelling. Devasia et al. recently used an NLME model to determine kidney TIAs in PRRT with 177Lu-DOTATATE using a bi-exponential function and biokinetic data measured with SPECT/CT. The NLME model provides valuable information across patients (inter-individual variability) and within patients (intra-individual variability) that could be used to optimise the individual estimation of TIAs [11]. It has been shown that implementing NLME modelling could lead to an acceptable accuracy of calculated TIAs [11]. Although the NLME model has been shown to perform well [11], [12], the bi-exponential function used in the NLME method might not be applicable to biokinetic data of different radiopharmaceuticals, organs, or patient biokinetic data. Therefore, using the best fit function within the NLME model approach may improve the NLME method.

In this work, we present a general method based on population information to improve the calculation of TIAs. It is demonstrated in an example of 177Lu-PSMA therapy using the NLME model and population-based model selection (NLME-PBMS) method. In this NLME-PBMS method, a set of mathematical models (or functions) is defined, an NLME model fit is performed, and the best fit function is selected using the Akaike weight method [13], [14]. In addition, the performance of the NLME-PBMS method is compared to the IBMS and the SP-PBMS methods reported in the literature [6]. In the clinic, the NLME-PBMS method can be used to define the best fit function for a given radiopharmaceutical, organ and population kinetic data set that can be used to fit a “new” patient with the known inter- and intra-individual variability.

2 Materials and methods

2.1 Biokinetic data

Biokinetic data (the time-activity data) of [177Lu]Lu-PSMA-I&T RLT in kidneys were collected from thirteen patients with metastatic castration-resistant prostate cancer [6], [15]. All patients underwent [177Lu]Lu-PSMA-I&T radioligand therapy (RLT) and had post-therapeutic planar whole-body scintigraphies at (1.1 ± 0.7) h, (20.7 ± 2.3) h, and (163.8 ± 2.0) h p.i. One patient had one additional time point at 66.1 h p.i. Three patients had two extra time points at (45.9 ± 1.6) h and (68.7 ± 1.7) h p.i. The activity was quantified from the images by drawing the regions of interest using the geometric mean of anterior and posterior counts with background corrections [6], [15]. There was no scatter correction performed. Calibration factors were calculated based on the whole-body region of interest in the first scan. The institutional review board approved this retrospective analysis (115/18S), and the requirement for informed consent was waived.

2.2 Set of sums of exponential functions

As suggested by Burnham et al., only well-suited functions should be included in the investigated function set [8]. Therefore, sums of exponential functions (one, two and three terms) [6], [7] with increasing complexity were used and analysed to describe the physical and biological processes of the radiopharmaceutical’s pharmacokinetics (Eqs. (1), (2), (3), (4), (5), (6), (7), (8), (9), (10), (11)):(1) f2at=A1e-λ1+λphyst

(2) f3at=A1e-λ1+λphyst+A2e-λphyst

(3) f3bt=A1αe-λ1+λphyst+A11-αe-λphyst

(4) f3ct=A11-αe-λ1+λphyst+A1αe-λphyst

(5) f4at=A1e-λ1+λphyst+A2e-λ2+λphyst

(6) f4bt=A1αe-λ1+λphyst+A11-αe-λ2+λphyst

(7) f5at=A1e-λ1+λphyst+A2e-λ2+λphyst-A3e-λphyst

(8) f5bt=A1e-λ1+λphyst+A2e-λ2+λphyst-(A2+A1)e-λ3+λphyst

(9) f5ct=A1e-λ1+λphyst+A2e-λ2+λphyst+A3e-λphyst

(10) f6at=A1e-λ1+λphyst+A2e-λ2+λphyst-A3e-λ3+λphyst

(11) f6bt=A1e-λ1+λphyst+A2e-λ2+λphyst+A3e-λ3+λphyst

where fi is a fit function with i estimated parameters, the Aj are pre-factors with values ≥0, λphys is the physical decay constant of 177Lu (λphys=ln2/T1/2, with the half-life T1/2=6.6443 d [16]), λj are the biological decay parameters of the radiopharmaceutical. All biological decay constants λj were estimated with a value ≥0 as described in the literature [17], [18].

2.3 Non-linear mixed-effects model

The NLME model was used to estimate the parameters of the SOE functions in Eqs. (1), (2), (3), (4), (5), (6), (7), (8), (9), (10), (11), i.e. Aj and λj. The estimated parameters in the NLME model were defined as a combination of the fixed and random effects. Fixed effects are the mean values of the estimated parameters in the population, while random effects show the variability of the estimated parameters [19]:(12) Pi=TVPi×exp(ETAi)

(13) ETAi=N(0,σi2)

Here Pi is the estimated parameter i, TVPi is the fixed effect of the estimated parameter i in the population, and ETAi is a random number following a Gaussian distribution with mean zero and variance σi2 (random effect) [11], [12], [20].

2.4 Study workflow

The parameters of Eqs. (1), (2), (3), (4), (5), (6), (7), (8), (9), (10), (11) were fitted to the biokinetic data of [177Lu]Lu-PSMA-I&T RLT in kidneys using the NLME method. All fittings were performed in MATLAB software vR2020a. The exponential error model was used for the NLME model fittings [12]. The following computational settings were used for the IBMS and SP-PBMS methods as suggested in the literature [6]: Rosenbrock Algorithm, objective function convergence criterion 10−4, and absolute-based variance model with a fractional standard deviation of 0.15. Starting values of the parameters were chosen based on a trial-and-error approach. The goodness of fit for each SOE function was checked by visual inspection of the fitted graphs and the coefficient of variation CV of the fitted fixed-effects parameters (Table 1 in ref. [7], goodness of fit acceptable for CV < 0.5). The best fit function selection was performed by calculating the corrected Akaike Information Criterion AICc, which is the AIC corrected for the case of N/K < 40 (number of data N and number of parameters K) [8], [9]. The corresponding Akaike weights [7], [8], [9] were calculated as follows:(14) AICc=-2lnP+2K+2K(K+1)N-K-1

(15) Δi=AICci-AICcmin

(16) wAICci=e-Δi2/∑i=iFe-Δi2

where P is the objective function minimized for the fitting, AICcmin is the lowest AICc value of all fitted functions, Δi is the difference of the AICci of function i and the AICcmin, F is the total number of functions that passed the goodness-of-fit test and wAICci is the Akaike weight of function i. The Akaike weights indicate the probability that the model is the best among the whole set of considered models [8], [9]. To analyse the stability of the NLME-PBMS method through model selection for different set of data, the Jackknife method was used as suggested in the literature [9], [21]. Here, the leave-one-out method was applied thirteen times with only twelve patients included for the calculation of the Akaike weights.

The best function in the IBMS method for the same biokinetic data was function f2a (Eq. (1)) [6]. In the SP-PBMS method, the best function describing the biokinetic data of [177Lu]Lu-PSMA-I&T was function f3b¯ (i.e. f3b, Eq. (3) with α= 0.9632) [6]. Therefore, we used functions f2a and f3b¯ to describe the biokinetics of [177Lu]Lu-PSMA-I&T in kidneys to represent IBMS and SP-PBMS methods, respectively. The TIAs were calculated using the functions from IBMS, SP-PBMS and NLME-PBMS methods with integration from t = 0 min to t = 100 000 min. The standard deviations of the TIAs from the IBMS and SP-PBMS method were calculated based on the bootstrap method using the frequency-based sampling for the parameters as suggested by Saltelli et al. [22]:(17) Xis=Gisinωis=12+1πarcsin(sinωis+φi)

where s is a modified scalar variable varying over the range -π/2 < s<π/2, and ωi is the frequency value used for parameter Xi and φiis a random phase-shift chosen uniformly between 0 and 2π. Parameters of the SOE function Xi were sampled from the variability obtained from the IBMS and SP-PBMS methods. The frequencies show a specific number related to the sinusoidal functions assigned to each parameter [23]. The vector of the input parameters of interest Xi was sampled from the log-normal distribution (Eqs. (18), (19)). The sampling of the IBMS and SP-PBMS fitted parameters were generated using Eq. (17). The mean and SD of the parameters in the log-normal distribution were calculated based on the following equations [23]:(18) meanlog-normal=2×ln(meannormal)-0.5×ln(meannormal2+SDnormal2)

(19) SDlog-normal=ln(meannormal2+SDnormal2)-2×ln(meannormal)

where meanlog-normal is the mean of the parameters in log-normal distribution, SDlog-normal is the SD of the parameters in log-normal distribution, meannormal is the mean of the parameters in normal distribution, and SDnormal is the SD of the parameters in normal distribution. The frequency ωi and the number of model evaluations in Eq. (17) were fixed to 1028 and 8193 [23], respectively. The standard deviations of individual TIAs from the NLME-PBMS method was not calculated as the NLME model fitting method only reports the intra-individual variability from the populations [12], [19].

The TIA was also calculated for the model averaging (MA) [8], [10] from the NLME-PBMS method according to the following equation:(20) TIAMA=∑iwAICci×TIAi

where TIAMA is the TIA of MA, wAICci is the Akaike weight of the SOE function i calculated using Eq. (16), and TIAi is the TIA of the SOE function i in the NLME-PBMS method. Relative deviations (RDs) and root-mean-square errors (RMSEs) were used to analyse the accuracy of the calculated TIAs using the IBMS, SP-PBMS, and NLME-PBMS methods (per single fit function) with the TIAs obtained from the NLME-PBMS (MA) method as the reference. The total number of fitted parameters leads to the superiority of the NLME-PBMS method over the IBMS and SP-PBMS methods. For the same data, the total number of parameters in the NLME-PBMS, IBMS and the SP-PBMS method are K = 7 (Table 1), K = 26 (Function f2bin [6]), and K = 27 (Table 1 in [6]) respectively. Based on Eq. (14), with a conservative assumption of equal −2 ln (P) values for all fits and N = 46 data points, the values of the term (2K+2KK+1N-K-1) for the NLME-PBMS, IBMS and SP-PBMS methods are 17, 126, and 138, respectively. These large differences in the AICc leads to Akaike weights (Eqs. (15), (16)) of the best function from the NLME-PBMS, IBMS and SP-PBMS models of about 100%, 2.2 × 10−24% and 5.2 × 10−25%, respectively. This result supports the use of the NLME-PBMS method as the only reference in our study, as including the IBMS and SP-PBMS method in the model averaging will not lead to a relevant change.

The NLME-PBMS (MA) considers all functions that passed the goodness of fit test with their corresponding Akaike weights. The relative deviation RDs and the RMSEs were calculated as follows:(21) RDj=TIAj-TIAMATIACMA

(22) RMSEj=SDRDj2+MeanRDj2

where RD is the relative deviation, TIAj is the TIA value calculated using method j, and TIAMA is the TIA calculated using the MA method. RMSEj is the root-mean-square error over all patients of RDj, SDRDj is the standard deviation of RDj, MeanRDj is the mean of RDj and j corresponds to the IBMS or SP-PBMS or NLME-PBMS methods.

3 Results

The fittings using the NLME method for functions f4b, f5a, f5b, f5c, f6a and f6b did not pass the goodness-of-fit test, i.e. the fitting failed based on the visual inspection of the fitted curves and/or the CV > 0.50 of the fitted fixed effects (Table 1). Function f3a was selected as the function most supported by the data in the NLME-PBMS approach based on the Akaike weight of 51.72 % (Table 1). The Jackknife method was applied to the subset functions f3a, f3c and f4a as these functions had the highest Akaike weights (Table 1). Based on the Jackknife method, the selection of f3a was relatively stable with a median Akaike weight of 50% and range of 38%–74% (Table 1). Function f3a had the highest Akaike weight in 12/13 Jackknife simulations. Therefore, in one case f3c had a higher Akaike weight than f3a (Table 1).Table 1 Goodness of fits and Akaike weights for the NLME model with different SOE functions. The total number of biokinetic data N used in this analysis is 46; the number of parameters of the NLME model for the corresponding SOE function is given in column K.

Equation number	Function name	K	Coefficient of Variation CV (max)Ϯ	Akaike weight (%)‡	Jackknife Akaike weights (% mean (SD); % median [min, max])§	
1	f2a	5	0.07	3.11	–	
2	f3a	7	0.20	51.72	54 (11); 50 [38, 74]	
3	f3b	7	0.16	3.84	–	
4	f3c	7	0.21	22.37	33 (10); 33 [19, 59]	
5	f4a	9	0.27	18.96	13 (7); 14 [2,22]	
6	f4b	9	0.78	–	–	
7	f5a	11	Infinity	–	–	
8	f5b	11	Infinity	–	–	
9	f5c	11	1.26	–	–	
10	f6a	13	1.62E+12	–	–	
11	f6b	13	0.91	–	–	
Ϯ The maximum value of the CVs of the fit parameters (CV is calculated as expω2-1[19], with ω2 being the variance of the fixed effect). Only CV values <0.50 pass the goodness-of-fit check (Table 1 in [7]).

‡ Models passing the goodness-of-fit test are included in the NLME-PBMS (MA) model.

§ The Jackknife analysis was based on functions f3a, f3c and f4a.

The five SOE functions, i.e. f3a, f3c, f4a, f3b, and f2a, were used for the calculation of the TIAs of the NLME-PBMS (MA) method using Eq. (20). Fig. 1 compares the curves obtained from the IBMS (f2a), SP-PBMS (f3b¯), NLME-PBMS (f3a), and NLME-PBMS (MA) methods. The reference model NLME-PBMS (MA) shows a good performance in determining the biokinetic data in all patients (Fig. 1). Visual inspection of the fitted graphs in Fig. 1 shows that the NLME-PBMS (f3a) method has a better or equivalent performance as the SP-PBMS and IBMS methods and is similar to the NLME-PBMS (MA) method. Fig. 2 shows the comparison of the TIAs calculated using the IBMS, SP-PBMS, NLME-PBMS (f3a) and NLME-PBMS (MA) methods. The SDs of the TIAs per-injected activity or time-integrated activity coefficient (TIACs) in Fig. 2 were calculated based on the bootstrap method with sampling according to Eqs. (17), (18), (19). In general, the TIACs from all methods were similar except for P6. Fig. 3 and Table 2 show the RD of the TIAs of the IBMS, SP-PBMS (f3b¯), NLME-PBMS (f3a), NLME-PBMS (f4a), NLME-PBMS (f2a), NLME-PBMS (f3b), and NLME-PBMS (f3c) methods to the TIAs calculated using NLME-PBMS (MA) method. The NLME-PBMS (f3a) has a better performance than the IBMS and SP-PBMS (f3b¯) by a factor of four based on the RMSEs values (Table 2).Figure 1 Time-activity data and fit curves were obtained from the IBMS, SP-PBMS, NLME-PBMS (f3a) and NLME-PBMS (MA) methods.

Figure 2 Time-integrated activity coefficients (TIACs) calculated using the NLME-PBMS (f3a), NLME-PBMS (MA), IBMS (f2a) and SP-PBMS (f3b¯) methods. The standard deviations of the TIACs in the IBMS and PBMS methods were calculated based on the sampling according to Eq. (24).

Figure 3 RDs of TIAs obtained from the IBMS, SP-PBMS (f3b¯), NLME-PBMS (f3a), NLME-PBMS (f4a), NLME-PBMS (f2a), NLME-PBMS (f3b), and NLME-PBMS (f3c) methods.

Table 2 The RDs and RMSEs of different methods in comparison to the NLME-PBMS (MA) reference.

Methods	RD (%)	RMSE of the RD (%)	
	Mean (SD)	Median [min, max]		
IBMS	1 (7)	2 [−11, 13]	7.4	
SP-PBMS	3 (8)	1 [−4, 29]	8.8	
NLME-PBMS (f2a)	−2 (16)	6 [−36, 17]	16.5	
NLME-PBMS (f3a)	−2 (2)	−2 [−4, 3]	2.4	
NLME-PBMS (f3b)	14(19)	21 [−30, 34]	23.5	
NLME-PBMS (f3c)	−1 (1)	−2 [−3, 2]	2.0	
NLME-PBMS (f4a)	4 (3)	3 [0, 8]	4.5	

4 Discussion

Accurate determination of individual TIAs is highly desirable in MRT. Recently, it has been shown that the NLME method could be used for dosimetry in MRT with a relatively accurate determination of TIAs [11], [12]. The advantages of the NLME model are the ability of the method to describe the kinetics of radiopharmaceutical drugs in terms of intra- and inter-individual variability [12], [19]. The information on intra- and inter-individual variability may benefit specific applications such as a single-time-point dosimetry approach [11], [12] in MRT. However, model selection is needed to determine the best fit function with the NLME model approach. Model selection has been shown as an essential and critical aspect of scientific data analysis [8]. Using a model selection procedure for NLME model fitting could increase the reproducibility of the fitting results in contrast to applying the rule of thumbs [24] or simply guessing by the user. The selection of a suitable mathematical model (i.e. function) for calculating TIAs is essential, as using an improper function might invalidate or deteriorate the result [6]. The strong influence of the chosen function on the calculation of RDs can be seen in our analysis results (Fig. 3 and Table 2). Therefore, in this study, we applied the NLME model with model selection to calculate the TIAs. The output of the NLME-PBMS method was then compared with the IBMS and SP-PBMS methods to check if the NLME-PBMS method could be used to determine accurate TIAs in MRT.

Although the functions f3a and f3c are equivalent functions differing only by their parameterisation, function f3a was selected as the best model based on the Akaike weight. This result demonstrates the relevance of model selection as even equivalent functions may lead to a different performance in determining the biokinetic data of radiopharmaceuticals (Fig. 3 and Table 2) due to the assumed difference in the distribution of parameters.

P6 in Fig. 1 demonstrates how the chosen fit function may lead to a considerable difference in TIAs: The fitted curves from the NLME-PBMS (f3a) method seem to be more adequate than the IBMS and SP-PBMS methods in describing the biokinetic data in P6. Accordingly, the NLME-PBMS (f3a) RMSE value is around four times lower than that of IBMS and SP-PBMS (Table 2). However, one disadvantage of the NLME model approach is that the standard deviation of individual TIAs from the NLME-PBMS method is not analysed as the output of the NLME model only reports the intra-individual variability in the population.

The best model selected from the NLME-PBMS method, i.e. f3a, has a similar structure to the best model selected from the SP-PBMS method, i.e. f3b¯ (i.e. f3b, Eq. (3) with α = 0.9632) [6]. Calculating the TIAs using f3b¯ (SP-PBMS, which has 2 adjustable parameters) [6] needs at least three data to obey the constraint of Kmax (maximum number of fitted parameters) = N (number of data) + 1 of the individual fitting of a new patient. In contrast, the fitting of the NLME model uses a population-based fit, where the ratio between the number of data and the number of fitted parameters is higher than in the individual fits [12]. This allows its use in single-time-point dosimetry.

A general problem in MRT is that it is unknown which function fits the available data best. This applies to the case with many data and with few data per organ. Therefore, a different modeller might use other functions for the same dataset, leading to reproducibility issues [6]. Here the proposed method will provide a more reproducible approach as the model selection includes many suitable SOE functions. Available population data in nuclear medicine are usually heterogeneous and sparse. The presented method can be used for this common situation, as pharmacokinetic information of heterogeneous data can be derived from a population and introduced for the individual fit. Selecting the function most supported by the data [9], [10] and using it in the NLME model fitting is very useful also for other applications in molecular radiotherapies, such as in single-time-point imaging dosimetry [11], [12]. Presumably, when a better fit function is used, the accuracy of single-time-point dosimetry also becomes higher.

In the clinical setting, the presented model selection method would be essential in the sequence of steps as follows: (1) collect biokinetic data of a patient population, (2) derive the fit function most supported by the data as presented here using the NLME-PBMS method, and (3) perform dosimetry for new patients using the derived best function and the NLME model fitting with the known inter- and intra-individual variabilities. The fitted parameters from the previously measured population can be used as the starting values for the fitting of the new patients. Another possibility to decrease the total number of the adjustable parameters is by fixing the parameters of the population to the values from the previous population.

In this study, a maximum of five data per patient were available. More data per patient could yield a different best function (with more parameters) for the NLME-PBMS method, even in the same patient population for the same organ and radiopharmaceutical. Therefore, the NLME-PBMS approach needs to have enough data per patient for the patient population used to determine the best fit function. If the aim is the accurate determination of the TIAs, “enough” data can be determined based on the change in the determined TIAs when increasing the number of data per patient (and possibly also the number of parameters of the best fit function).

Compartmental models, including additional physiological information in the investigated function set, could further improve the NLME-PBMS method results. Compartmental models might perform better than the investigated fit functions based on sum of exponential functions, which are most frequently used for curve fitting in dosimetry in molecular radiotherapy.

In this study, function f3c seems to be better than f3a based on the RMSE values (Table 2). This is in contrast to the Akaike weight results which show f3a as the best function obtained from the NLME-PBMS method (Table 1). This reflects that Akaike weight and RMSE values are different measures. The AICc is used to find the model most supported by the data; however, this tells nothing about the difference in the TIAs for the different functions. Therefore, the differences in TIAs were quantified using the RMSE. In addition, both have an associated uncertainty. The differences between the two measures noted here may disappear if more data and patients are used. Another model selection procedure, the Bayesian Information Criterion (BIC) [25], was also used for our data. As a result, the BIC favours the same best function as the Akaike weight.

The method presented in this study is not constrained to the function set used in this manuscript. Thus, “arbitrary” functions, which should be suited to the kinetics [8], can be added to the here investigated function set, and the best fit function can be determined according to the presented algorithm. Therefore, another function could be the best for a different biokinetic data set (e.g. more patient data at different time points, or for different organs and radiopharmaceuticals).

The advantages of the proposed method are achieved by improving both inputs, i.e. (1) the data and (2) the set of models from which the best one is selected. This, in turn, also improves the results for the TIAs:1. Data of a population used in the NLME-PBMS model instead of just a single patient in the IBMS method for estimating TIAs is advantageous. Fitting three biokinetic data using IBMS only allows a mono-exponential function with a maximum of two parameters as the degree of freedom cannot be zero (degree of freedom = N–K [6]). We show that the mono-exponential function could lead to inaccurate determination of TIAs in P6 (e.g. Figure 1, Figure 3). By using the NLME-PBMS model, we could fit an SOE function with a higher number of parameters to patients with a low number biokinetic data.

2. The set of models from which the best one is selected using the Akaike weights is also constrained by Kmax = (N − 2) (Eq. (21)) [7], [8], [9], [26] in both the IBMS and the PBMS methods. Therefore, in our example, the NLME model, in principle, would allow including a much higher number of parameters in the model set functions. However, as seen in Table 1, SOE functions with more than four parameters are unnecessary (based on this biokinetic data set).

The limitations of this study are:1. Quantification of the activities was based on planar images. It has been shown that activity quantification using imaging with three-dimensional approaches such as SPECT and PET is superior to planar imaging [27]. Planar imaging suffers from biases, and the background correction is known to be sensitive to region-of-interest definition [27], [28]. It is recommended that biokinetic data from three-dimensional approaches be used for future studies. As the main aim of this study was to introduce and present an NLME-PBMS method, which can be implemented for any quantitative imaging data, the results obtained with planar imaging from this study should be translatable to the case of SPECT and PET imaging. Improved quantification should, however, improve the population approach. Reducing systematic and random quantification errors will probably result in changed fixed effects, reduced random effects (inter-individual variability), and reduced intra-individual variability.

2. Quantification without scatter correction is a limitation in our study. The presented population method is constructed to use similarities in the biokinetics of the population under investigation. Thus, patients must be measured similarly for the method to be applicable. Therefore, removing scatter from the data by implementing quantitative imaging is expected to improve the results obtained using this population-based method, as discussed in item 1 above.

3. The starting values of the adjustable parameters were chosen based on a trial-and-error approach. This could lead the minimisation algorithm to local minima of the objective function. Therefore, different starting values were investigated. Of these, only fit results showing a good fit based on the goodness of fit criteria having the highest likelihood were used for further evaluation.

5 Conclusions

A critical unresolved issue in MRT dosimetry is the question of which fit function to use for calculating the time-integrated activity. We developed a method to include fitting function selection in a population-based method to determine the best fit function for calculating TIAs in MRT for a given radiopharmaceutical, organ and patient data set. The best fit function can then be used to fit a “new” patient with the known inter- and intra-individual variabilities from the corresponding patient population. The technique combines standard practice approaches in pharmacokinetics, i.e. an Akaike-weight-based model selection and the NLME model framework. The fit function selection accuracy and precision in the NLME-PBMS method are less dependent on noise in the data than in the IBMS method, which relies on only very few data of a single patient. Thus, using the proposed fit function selection and the NLME model will increase the reproducibility of the found best fit function for determining TIAs.

Consent for Publication

All authors read the manuscript and consented to its publication.

CRediT authorship contribution statement

Deni Hardiansyah: Conceptualization, Methodology, Software, Validation. Ade Riana: Software, Validation. Matthias Eiber: Resources. Ambros J. Beer: Resources. Gerhard Glatting: Conceptualization, Methodology, Validation, Resources.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgment

This study was supported by research grant PUTI Q1 from Universitas Indonesia with grant number NKB-475/UN2.RST/HKP.05.00/2022 and mobility grant WCP from Kemenristek Dikti with grant number 3252/E4/DT.04.03/2022.
==== Refs
References

1 Lassmann M. EANM Dosimetry Committee guidance document: good practice of clinical dosimetry reporting Eur J Nucl Med Mol Imaging 38 1 2011 192 200 20799035
2 Glatting G. Lassmann M. Nuclear medicine dosimetry: quantitative imaging and dose calculations Z Med Phys 21 4 2011 246 247 21862300
3 Glatting G. Bardiès M. Lassmann M. Treatment planning in molecular radiotherapy Z Med Phys 23 4 2013 262 269 23597414
4 Hardiansyah D. The role of patient-based treatment planning in peptide receptor radionuclide therapy Eur J Nucl Med Mol Imaging 43 5 2016 871 880 26577941
5 Union, C.o.t.E., European Council Directive 2013/59/Euratom on basic safety standards for protection against the dangers arising from exposure to ionising radiation and repealing Directives 89/618/Euratom, 90/641/Euratom, 96/29/Euratom, 97/43/Euratom and 2003/122/Euratom. Off J EU L13; 2014: 1–73.
6 Hardiansyah D. A population-based method to determine the time-integrated activity in molecular radiotherapy EJNMMI Phys 8 1 2021 82 34905131
7 Kletting P. Molecular radiotherapy: the NUKFIT software for calculating the time-integrated activity coefficient Med Phys 40 10 2013 102504
8 Burnham K.P. Anderson D.R. Model selection and multimodel inference 2002 Springer-Verlag New York
9 Glatting G. Choosing the optimal fit function: comparison of the Akaike information criterion and the F-test Med Phys 34 11 2007 4285 4292 18072493
10 Kletting P. Glatting G. Model selection for time-activity curves: the corrected Akaike information criterion and the F-test Z Med Phys 19 3 2009 200 206 19761098
11 Devasia T. A novel time-activity information sharing approach using nonlinear mixed models for patient-specific dosimetry with reduced imaging time points: application in SPECT/CT imaging post 177Lu-DOTATATE J Nucl Med 62 8 2020 1118 1125 33443063
12 Hardiansyah D et al. Single-time-point estimation of absorbed doses in PRRT using a non-linear mixed-effects model. Z Med Phys; 2022 [in print 20 July 2022].
13 Akaike H. In: Petrov BN, Csaki F, editor. Second international symposium on inference theory. Budapest: Akademiai Kiado; 1973. p. 267–81.
14 Hurvich C.M. Tsai C.-L. Regression and time series model selection in small samples Biometrika 76 2 1989 297 307
15 Kletting P. Modeling and predicting tumor response in radioligand therapy J Nucl Med 60 1 2019 65 70 29748236
16 Kondev F.G. Nuclear data sheets for A=177 Nuclear Data Sheets 159 2019 1 412
17 Williams L.E. On the correction for radioactive decay in pharmacokinetic modeling Med Phys 22 10 1995 1619 1626 8551986
18 Glatting G. Reske S.N. Treatment of radioactive decay in pharmacokinetic modeling: influence on parameter estimation in cardiac 13N-PET Med Phys 26 4 1999 616 621 10227364
19 Mould D.R. Upton R.N. Basic concepts in population modeling, simulation, and model-based drug development-part 2: introduction to pharmacokinetic modeling methods CPT Pharmacometrics Syst Pharmacol 2 2013 e38 23887688
20 Bauer R.J. Guzy S. Ng C. A survey of population analysis methods and software for complex pharmacokinetic and pharmacodynamic models with examples AAPS J 9 1 2007 E60 E83 17408237
21 Shao J. Tu D. The jackknife bootstrap. Springer Series in Statistics 1995 Springer-Verlag New York
22 Saltelli A. Tarantola S. Chan K.P.S. A quantitative model-independent method for global sensitivity analysis of model output Technometrics 41 1 1999 39 56
23 Hardiansyah D. Important pharmacokinetic parameters for individualization of 177Lu-PSMA therapy: a global sensitivity analysis for a physiologically-based pharmacokinetic model Med Phys 48 2 2021 556 568 33244792
24 Siegel J.A. MIRD pamphlet no. 16: techniques for quantitative radiopharmaceutical biodistribution data acquisition and analysis for use in human radiation dose estimates J Nucl Med 40 2 1999 37S 61S 10025848
25 Schwarz G. Estimating the dimension of a model Ann Stat 6 2 1978 461 464
26 Kletting P. The NUKDOS software for treatment planning in molecular radiotherapy Z Med Phys 25 2015 264 274 25791740
27 He B. Frey E.C. Comparison of conventional, model-based quantitative planar, and quantitative SPECT image processing methods for organ activity estimation using In-111 agents Phys Med Biol 51 16 2006 3967 3981 16885618
28 Garkavij M. 177Lu-[DOTA0, Tyr3] octreotate therapy in patients with disseminated neuroendocrine tumors: analysis of dosimetry with impact on future therapeutic strategy Cancer 116 4 Suppl 2010 1084 1092 20127957
