In brief, biokinetic data of [111In]In-DOTATATE in kidneys from eight patients with either meningioma (n = 4) or neuroendocrine tumours (n = 4) were used in this study [15, 21]. An activity of (140 ± 14) MBq of [111In]In-DOTATATE was administered intravenously to the patients as a (51 ± 8) min infusion. Planar whole-body scintigraphies using a double-head gamma camera (ECAM, Siemens, Erlangen, Germany) were performed at T1 = (2.9 ± 0.6) h, T2 = (4.6 ± 0.4) h, T3 = (22.8 ± 1.6) h, T4 = (46.7 ± 1.7) and T5 = (70.9 ± 1.0) h p.i. [21] using a medium energy collimator with energy windows A1 = 171 keV (width 15%), A2 = 245 keV (15%), B1 = 142 keV (18%), and B2 = 205 keV (18%). Background correction and self-attenuation were included in the measurement of organ activity as a function of time according to the MIRD pamphlet number 16 [22]. The percentage of the administered activity in kidneys was used in this study [21]. Biokinetic data of [111In]In-DOTATATE were used as a surrogate for predicting the kinetics of [90Y]Y-DOTATATE used for peptide-radionuclide therapy, as suggested in the literature [21, 23].

In this study, the following sums of exponential (SOEs) functions with 3, 4, and 5 parameters and different parameterisations were used to fit the biokinetics of [111In]In-DOTATATE in kidneys (Eqs. 1–11). The different parameterisations were investigated to demonstrate that the NLME modelling yields different results for different parameterisations.

$$f_ \left( t \right) = A_ e^ + \lambda_}}} } \right)t}} - A_ e^ + \lambda_}}} } \right)t}}$$

(1)

$$f_ \left( t \right) = \frac + \lambda_}}} } \right) \times \left( + \lambda_}}} } \right)} \right)}} \left( + \lambda_}}} } \right) - \left( + \lambda_}}} } \right)} \right)}}\left[ + \lambda_}}} } \right)t}} - e^ + \lambda_}}} } \right)t}} } \right]$$

(2)

$$f_ \left( t \right) = A_ \frac + \lambda_}}} } \right) \times \left( + \lambda_}}} } \right)} \right)}} + \lambda_}}} } \right) - \left( + \lambda_}}} } \right)} \right)}}\left[ + \lambda_}}} } \right)t}} - e^ + \lambda_}}} } \right)t}} } \right]$$

(3)

$$f_ \left( t \right) = \frac \times \lambda_ } \right)}} \left( - \lambda_ } \right)}}\left[ } \right)t}} - e^ } \right)t}} } \right]$$

(4)

$$f_ \left( t \right) = A_ e^ + \lambda_}}} } \right)t}} - A_ e^ + \lambda_}}} } \right)t}} - \left( - A_ } \right) - e^ + \lambda_}}} } \right)t}}$$

(5)

$$f_ \left( t \right) = A_ \alpha e^ + \lambda_}}} } \right)t}} - A_ \left( \right)e^ + \lambda_}}} } \right)t}} - A_ \left( \right)e^ + \lambda_}}} } \right)t}}$$

(6)

$$f_ \left( t \right) = A_ \left( \right)e^ + \lambda_}}} } \right)t}} - A_ \alpha e^ + \lambda_}}} } \right)t}} - A_ \left( \right)e^ + \lambda_}}} } \right)t}}$$

(7)

$$f_ \left( t \right) = \frac }} + \lambda_}}} }}} \right) - \left( + \lambda_}}} }}} \right) - \left( + \lambda_}}} }}} \right)} \right\}}}e^}}} t}} \left\ \right)e^ t}} - \alpha e^ t}} - \left( \right)e^ t}} } \right\}$$

(8)

$$f_ \left( t \right) = \frac }} + \lambda_}}} }}} \right) - \left( + \lambda_}}} }}} \right) - \left( + \lambda_}}} }}} \right)} \right\}}}e^}}} t}} \left\ t}} - \left( \right)e^ t}} - \left( \right)e^ t}} } \right\}$$

(9)

$$f_ \left( t \right) = A_ e^ + \lambda_}}} } \right)t}} - A_ e^ + \lambda_}}} } \right)t}} - \left( - A_ } \right)e^ + \lambda_}}} } \right)t}}$$

(10)

$$f_ \left( t \right) = A_ e^ + \lambda_}}} } \right)t}} + A_ e^}}} } \right)t}} - A_ e^ + \lambda_}}} } \right)t}} - \left( + A_ - A_ } \right)e^ + \lambda_}}} } \right)t}}$$

(11)

where \(f_\) is a fit function, \(i\) is the total number of the estimated parameters, \(A_ \left( \right)\) are the prefactors of the fit function with values \(\ge 0\), \(\lambda_}}}\) is the physical decay constant of 111In (\(\lambda_}}} = \ln \left( 2 \right)/T_ = 1.72 \times 10^ \min^\) [21]), \(\lambda_\) is the rate of blood circulation of 1 min \(\left( = \frac}} \right)\), \(\lambda_\) are the biological decay constants of the radiopharmaceutical with values \(\ge 0\), and the \(\alpha\) values are the fractional contributions of the corresponding exponentials with values between 0 and 1. As described in Burnham et al. [24], existing prior knowledge should be taken into account when selecting the functions to be used for model selection. Therefore, on the one hand only sums of exponential functions were considered [19, 25] and on the other hand the constraint \(f_ \left( \right) = 0\) was implemented. In addition, for functions with 4 and 5 parameters, a rapid increase in activity in the kidneys with a half-life of 1 min was added, which is caused by the blood circulation time in humans. SOE functions with less than three parameters did not pass the goodness-of-fit test and were not included in the analysis.

Parameters in the NLME model consist of the fixed and random effects (Eqs. 12–13) as reported in the literature [8, 15, 26]. Fixed effects describe the mean values of the estimated parameters in the population, while random effects describe the inter-patient variability of the estimated parameters between subjects in the population [27].

$$P_ = }_ \times \exp \left( }_ } \right)$$

(12)

$$}_ = N\left( ^ } \right)$$

(13)

where \(P_\) is the estimated parameter \(j\) in an SOE function, \(}_\) is the fixed effect of the estimated parameter \(j\), and \(}_\) the random effect. \(}_\) is a random number following a Gaussian distribution with mean zero and variance \(\sigma_^\). Parameters of the exponential functions in Eqs. (1–11) were modelled as the combination of a fixed effect and an inter-patient variability (random effect) plus the intra-patient variability.

The parameters of the SOE functions (Eqs. 1–11) were fitted to the biokinetic data of [111In]In-DOTATATE in kidneys (“Biokinetic data” Section) using the NLME method. All NLME model fittings and simulations were performed in MATLAB software vR2020a. As suggested in the literature, an exponential error model with log transformation was used [15]. The MS–NLME method is performed using the Akaike weight. The SOE function with the highest Akaike weight was selected as the fit function most supported by the data. The Akaike weights indicate the probability that the model is the best among the analysed models [24, 28]. The Akaike weights [19, 24, 28] of the SOE functions were calculated as follows:

$$} = - 2\ln \left( P \right) + 2K + \frac \right)}}$$

(14)

$$\Delta_ = }_ - }_$$

(15)

$$w_}_ }} = e^ }}}} /\mathop \sum \limits_^ e^ }}}}$$

(16)

where \(}\) is the corrected Akaike Information Criterion value, \(P\) is the obtained minimum objective function, \(}_\) is the lowest \(}\) value of the SOE functions, \(\Delta_\) is the difference of the \(}_\) of SOE function \(j\) and \(}_\), \(F\) is the number of SOE functions in the model set, and \(w_}_ }}\) is the Akaike weight of function \(j\). The stability of the best SOE function obtained from the MS–NLME method was tested using the Jackknife method [28, 29]. In this method, the leave-one-out method was applied eight times with only seven patients included for the calculation of the Akaike weights.

The performance in determining the TIAs for STP dosimetry of the best SOE function obtained from the MS–NLME method was compared to the performance of the often used bi-exponential function \(f_\) [8, 19, 20]. The parameters of the bi-exponential function (\(f_\) Eq. (4)) were fitted to a patient with only STP biokinetic data by simultaneously fitting within the NLME model framework this new patient’s limited measurement with all data points of all other patients. Biokinetic data at time points T3 = (22.8 ± 1.6) and T4 = (46.7 ± 1.7) p.i. were used for the STP fitting as suggested in the literature [15].

TIAs from the STP fitting using the bi-exponential function (\(f_\)) were calculated by integrating the individual simulated time-activity curves from t = 0–100,000 min (TIASTP_f3d). The STP NLME model fittings were repeated using the best model obtained from MS–NLME method, followed by calculating the corresponding TIAs (TIASTP_MS–NLME).

Relative deviations (RDs) and root-mean-square errors (RMSEs) were used to analyse the accuracy of the calculated TIASTP_f3d and TIASTP_NLME–PBMS with the TIAs obtained from the all-time-point fittings using the best model obtained from MS–NLME (TIAATP_MS–NLME) as the reference. The relative deviation RDs and the RMSEs were calculated according to

$$}_ = \frac}_ - }_} - },m}} }}}_} - },m}} }},$$

(17)

$$}_ = \sqrt }_ } \right)^ + \left( }_ } \right)^ } ,$$

(18)

where \(}_\) is the relative deviation of STP method \(k\) of patient \(m\),\(}_\) is the root-mean square over all patients of \(}_\), \(}_\) is the standard deviation of \(}_\), \(}_\) is the mean of \(}_\), and \(k\) determines the function used for the NLME modelling.