Time-Activity data fitting in molecular Radiotherapy: Methodology and pitfalls

Absorbed radiation doses are essential in assessing the effects, e.g. safety and efficacy, of radiopharmaceutical therapy (RPT). Patient-specific absorbed dose calculations in the target or the organ at risk require multiple inputs. These include the number of disintegrations in the organ, i.e. the time-integrated activities (TIAs) of the organs, as well as other parameters describing the process of radiation energy deposition in the target tissue (i.e. mean energy per disintegration, radiation dose constants, etc). TIAs are then estimated by incorporating the area under the radiopharmaceutical ’ s time-activity curve (TAC), which can be obtained by quantitative measurements of the biokinetics in the patient (typically based on imaging data such as planar scintigraphy, SPECT/CT, PET/CT, or blood and urine samples). The process of TAC determination/calculation for RPT generally depends on the user, e.g., the chosen number and schedule of measured time points, the selection of the fit function, the error model for the data and the fit algorithm. These decisions can strongly affect the final TIA values and thus the accuracy of calculated absorbed doses. Despite the high clinical importance of the TIA values, there is currently no consensus on processing time-activity data or even a clear understanding of the influence of uncertainties and variations in personalised RPT dosimetry related to user-dependent TAC calculation. As a first step towards minimising site-dependent variability in RPT dosimetry, this work provides an overview of quality assurance and uncertainty management considerations of the TIA estimation.


Introduction
Radiopharmaceutical therapy (RPT) is a rapidly growing field of nuclear medicine (Fig. 1) with well-documented clinical benefits in specific patient cohorts [1,2].Further improvements based on personalised treatment require individual absorbed dose estimation.Calculating the absorbed dose from RPT is highly complex and has many pitfalls.At a high level, accurate RPT dosimetry requires reliable software to delineate target regions or organs at risk, quantitative measurements of the radiopharmaceutical biokinetics, precise calculation of time-activity curves, information about the physical process of radiation deposition (incl.'S-values' -mean absorbed doses per radioactive decay for a given source-target combination) [3], as well as robust timeintegrated activity estimation.Each step influences the resulting absorbed dose and can be performed as independent actions or a series of Monte Carlo simulations [4].In recent years, coordinated projects have been launched to improve accuracy and harmonise treatment quality with RPT within large-scale international associations.The community has been involving major international associations such as the Society of Nuclear Medicine and Molecular Imaging (SNMMI), the European Association of Nuclear Medicine (EANM), the International Atomic Energy Agency (IAEA), the European Federation of Medical Physics (EFOMP) and the American Association of Medical Physics (AAPM) [5,6], as well as multiple spontaneous research projects (i.e., DosiTest and MRT-Dosimetry) [7,8].However, the first step of evaluating potential variability in dose calculation due to site-specific calculation methodology remains largely untackled.
The SNMMI 177 Lu Dosimetry Challenge 2021 [6] is a recent example of a large-scale project aiming to draw international attention to the fact that RPT dosimetry is overwhelmingly user-and software-dependent.Unfortunately, due to the high complexity of the RPT dosimetry workflow and the limited number of included patient datasets (a total of two) used in this particular challenge, the results mainly emphasise the fact that there is a wide variation in RPT dosimetry methods, without understanding contributions of individual steps within the workflow.Follow-up work focused on understanding and minimising the variations is therefore needed.
The selected function used for time-activity data inter-and extrapolation, together with the number of data samples and uncertainties, significantly affects the calculated area under the TAC (TIA) [9], affecting the calculated absorbed dose.The high user dependency on TIA values hampers the development of dosimetry-guided treatment planning, ultimately leaving the major advantages of RPT over chemotherapy or biological therapy unexploited [3].Examples include the possibility of monitoring local and systemic applications of therapeutic radiopharmaceuticals using imaging and potential treatment optimisation for well-defined patient populations or even individual patients.Despite the importance of the topic, to our knowledge, there is no international general guideline regarding TAC selection and fitting for RPT dosimetry purposes.
To address these limitations, we aim to provide a generic quality framework for TAC modelling in organ-based RPT dosimetry targeted at individuals involved in or interested in setting up patient-specific dosimetry workflows.A condensed summary of steps and procedures involved in TAC modelling, recommendations regarding the determination of the uncertainty of results, selection of the TAC fit models, processing of the data, illustration of typical or expected pitfalls, and suggestions on how to validate the relevance of the result obtained is provided.Although quantification accuracy and precision of the input imaging and biological data affect the accuracy of the TAC, this topic is outside the scope of this work, as it has been extensively covered elsewhere [24].For this work, we assume that images have been appropriately corrected (i.e. for attenuation, scatter, deadtime, and normalisation) and reconstructed according to appropriate  O.V. Ivashchenko et al. manufacturer specifications, local preferences or related quality guidelines [10,11].Therefore, the proposed TAC analysis framework presented in the following sections should be applied to appropriately reconstructed and corrected quantitative image data.

Workflow of TAC fitting
The process of TAC modelling for RPT dosimetry involves many steps (Fig. 2), roughly divided into four phases [12].These are: 1. the collection of the data required for the fit (Data Input), 2. the definition of the model or function used to model the change of organ activity with time (Functions Input), 3. the actual processing of the data and the definition of underlying uncertainties (Data Processing), 4. the selection of the data required for presenting the result of the fit (Data Output).
The following sections will detail these topics independently to simplify understanding of the inputs, processes and outputs.

Data input
Reliable and robust input data are essential for accurately determining TACs.As mentioned above, data must be quantitative (regardless of the method of how this is achieved), and time points must be appropriately defined.In addition, all relevant injection and time parameters must be automatically or manually recorded (i.e.activity administered, residual activity, lag time, radionuclide, etc.).The sections below describe the steps involved once the time series of activity data have been obtained and the relevant organs/tumours have been adequately segmented for further analysis.

Time-Activity Data: How and when to collect it
To perform radionuclide dosimetry, the Medical Internal Radiation Dosimetry (MIRD) scheme requires the input of the time-integrated activity in accumulating tissues [13].Therefore, the activity distribution of an injected radiopharmaceutical at different time points is required, e.g.collected with a gamma camera or PET scanner as a series of data points by acquiring images of organs/tumours at various time points after the radiopharmaceutical is administered to the patient.
For RPT agents, gamma cameras (through planar, Bremsstrahlung or SPECT imaging) are more frequently used because most RPT agents have single gamma photon emissions.Compared to SPECT, PET imaging presents many advantages (e.g., easier quantification, higher spatial & temporal resolution and sensitivity).However, there are only a few therapy nuclides with relevant positron emissions, for example, 90 Y [14,15].
The optimal time points at which to measure the data for calculating the TACs is a highly topical area of research.The MIRD guidelines suggest using at least three time points per exponential term of the fit function.In the case of a bi-exponential function, the measurement times should be the following multiples of the actual effective half-life (T_eff): 0.1, 0.3, 1, 3, and 5 times the T_eff [13].A low number of time points is desirable to improve clinical acceptance of a dosimetry imaging protocol.However, this raises the question of how many time points are required to reliably construct a TAC.
The TAC depends on many physiological kinetic parameters, including but not limited to radiopharmaceutical trapping, reversibility of uptake and half-lives (biological and physical).Biokinetic models can Fig. 2. Schematic representation of the workflows, inputs, and dependencies required for TAC customisation, including the data input, function input, data processing (fitting to data output blocks), and the data output steps.
O.V. Ivashchenko et al. be used to predict the time-dependent distribution of radiopharmaceuticals and their elimination rates along specific excretion pathways.The ICRP has produced reports that give parameters for tissue distribution and retention of a range of elements for age-dependent phantoms, together with data on urinary and faecal excretion [16,17].In actual patients, generic models based on standardised phantoms or volunteers may not deliver good results [18].
Many publications regarding 177 Lu peptide receptor radionuclide therapy (PRRT) illustrated the importance of an appropriate sampling schedule for TAC estimation in personalised dosimetry [19].More recent simulation work investigated sampling schedules for [ 177 Lu]Lu-PSMA therapy using 1 to 4 time points (also with a mix of planar and SPECT imaging), ultimately concluding that 3 time points sampled at 3-4, 96-100 and 192 h performs well in terms of root-mean-square error [20].They also noted that the presence of a late time point (in their case at 192 h) had a positive impact on the accuracy of the TIACs.Some examples of research studies that demonstrate the variation in acquired time points for RPT dosimetry calculations are summarised in Table 1.

A priori knowledge on the Time-Activity data
Even before the time-activity data of a specific patient are measured, useful prior knowledge can be utilised.Firstly, the half-life of the radionuclide used is always known and must be included in the fitting function employed.Secondly, the absolute errors of the time-activity data may be known [24].In this case, these errors should be considered to give less weight to uncertain time-activity data than to more accurate and precise data (see section 2.2.3 for more details).
If the absolute errors of the time-activity data are not known, information about the type of error model (e.g., constant, proportional) may nevertheless be available.If so, a suitable error model can be implemented (section 2.2.2).If the parameters of the error model are known from population data, they can be used for assigning adequate weights to the data during the fitting process.

A priori knowledge about the fit function parameters
Since there can be several local minima during the TAC fitting procedure and thus multiple solutions for the parameter fit, the best solution can be found easier and faster using a priori knowledge about the parameter values.Therefore, adding the following information can increase the probability of finding a good fit: 1. Implementing constraints.For example, the parameter positivity requirement increases the probability of finding a good fit.2. Considering knowledge about the distribution of the fit function parameters in the patient population.For example, assume a parameter is normally distributed, and the population's mean and standard deviations are known.This knowledge can be interpreted as already having a prior measurement of this parameter, effectively representing a Bayesian approach [13,51,52].Consequently, this "prior knowledge", i.e., the population's mean and standard deviations, can be integrated into the fitting procedure (section 2.2.3). 3. Identifying fit function parameters that may not be relevant for the TIA calculation because the values of the parameters are similar for all patients.This is more likely to occur with more complex models.A sensitivity analysis [25] can be used to determine the non-relevant parameters of the used fit function, possibly allowing fixing these parameters.Such an approach can be interpreted as a special case of item (2) with the population's standard deviation being zero.

Input data for the fit algorithm
Due to the small amount of measurement data combined with significant measurement errors, the 'objective function' (definition of the measure of the fit to the data -see section 2.2.3) often does not have a single local minimum.Depending on the initial values for the fit function parameters, the algorithm for minimising the objective function can deliver different results.Therefore, methods for quickly finding suitable starting values for the parameters are advantageous.If guessing appropriate parameters is not an option, random or systematic sampling of the entire parameter space can be used to find good initial values [9].If the distributions of the fit function parameters in the patient population are known, the mean or median values may be used as starting values.
Selecting suitable stopping criteria is essential to ensure the convergence of the fit algorithm.These criteria are usually tolerance or threshold values for changes in the step size and for the change in the objective function.If the step size or the change in the objective function is smaller than these specified tolerance values, the minimisation of the objective function is considered to have been achieved.These tolerance values must therefore be small enough to find a minimum of the objective function with suitable accuracy, yet should not be chosen too small as the search for the minimum of the objective function can require an undesirably large amount of computing time.
Often, the implementations of the fitting algorithms also use a maximum number of iterations constraint, which is not desirable as the algorithm is not converged if the algorithm ceases due to the constraint.

Functions input
The TIA is needed to reliably calculate the absorbed dose.For this integration to be possible, a "curve" from the injection time to a maximum time, which may be infinity, must be defined.Usually, the activity is measured only for a few time points.Therefore, the complete pharmacokinetics must be estimated for all time points by interpolation between the measurement points and extrapolation for times outside the measurement points.
As schematically visualised in Fig. 2, the curve fitting part of the TAC modelling consists of four main components, namely the selection of 1. the fit function to be fit to the time-activity data; 2. the error model function for the time-activity data; 3. the objective function to be minimised during the fit; 4. the actual algorithm used to perform the objective function minimisation.

Fit functions for time-activity data
Many functions (or models) can fit a discrete dataset.Choosing a specific function implies making a specific assumption about the underlying kinetics or model.Consequently, different functions will produce different results for the TAC and absorbed doses.A pharmacokinetic model is a mathematical description of the biodistribution of a known concentration of a substance over time [26].In the following, the models are divided into the following groups, namely *single time points estimates result in different accuracy of the TAC, often lower or (at best) comparable to the initial schedule.
the "empirical", "analytical", and "compartmental" models, as well as the particularly complex whole-body compartment models that contain detailed physiological structures and parameters.All these models, even the "empirical" ones, use assumptions about the underlying physiology.Therefore, there is not always a sharp dividing line between the different groups.
Notably, the integration of prior knowledge into the mathematical models can allow a reduction in the number of measurement data.For example, if it is known that f(t = 0) = 0, and only functions fulfilling this constraint are used, no measurement at time zero (or very near to zero) needs to be performed.

Empirical models.
An example of an empirical model is the trapezoidal rule for the integration between measured time points plus the assumption of physical decay until infinity from the last time point.This model is unphysiological, e.g., due to the non-continuous first derivative.The accuracy of the obtained area under the curve depends on the number of measurement points and their distribution over time.An advantage of this model is that no assumptions about the analytical form (or function) of the TAC are included, which could have led to a systematic error.There are also disadvantages: 1.As typically there are only a few measured time points (possibly not optimally scheduled), numerical integration may lead to errors.2. Assumptions must be made for the times between the injection of the radiopharmaceutical and the first measuring point.Often it is assumed that at the time of injection, the activity in the organ under consideration is zero.Then it is favourable not to choose the first measurement time too late, to minimise the error caused by this assumption [27].3. A further assumption is needed for the time after the last measurement.A conservative estimate assumes that the activity remains in the tissue and thus only decays according to the physical half-life.However, this estimation is often inaccurate, making it unsuitable for use in RPT [28] (liver radioembolisation is one exception).For tumours, this assumption can lead to an overestimation of possible efficacy, and for normal organs may lead to an overestimation of absorbed dose and thus to possible undertreatment.As a result, a final measurement, e.g. after at least five effective half-lives, is desirable to reduce the magnitude of the error.

Analytical models.
Assuming that the relevant biological processes, e.g. the transfer of substances between different compartments, are of first-order kinetics, it can be shown that the kinetics can be described by a sum of exponential functions [26,27]: where A(t) is the activity at time t, λ phys represents the nuclide's physical decay constant, A i and λ i represent the coefficient and rate constant of exponential term i.
In Eq. ( 1), A(t) is a function of time t with the parameters λ i ≥ 0. Integrating λ phys explicitly into the model achieves better results and makes projection for different physical half-lives easier [29,30].The number n of exponential functions needed for a fit depends on the data and can be determined using a model selection procedure (Section 2.3.3.).
The TIA from the injection time t = 0 to infinity using Eq. ( 1) results in where λ phys represents the nuclide's physical decay constant, λ i -biological time constant of term i of the fit function, t -time, A i -the coefficient of respective exponential term i of the fit function, i -number of terms that are integrated in a piecewise manner, N -number of decays of the radiopharmaceutical drugs in the organs.Realistic biological processes always have non-linear components.Therefore, the sums of exponential functions are only approximate solutions.

Compartmental models.
Compartment models are based on grouping "similar" spaces and substances to simplify the mathematical description of physiological processes [26].An organ may be described by a few interconnected compartments.Between these compartments, an exchange of the considered radiopharmaceutical or metabolisation occurs.Physiologically-based pharmacokinetic (PBPK) models reflect (approximately) the actual anatomy and physiology [31].
The mathematical representation of the movement of the radiopharmaceutical within such a compartment model is described by linear or nonlinear differential equation systems.Mostly, investigated models have very few compartments, a consequence of the limited possibilities to measure data from all compartments, the often very small number of measurement data, and the occurring measurement errors.

Whole-Body physiologically based pharmacokinetic models.
A special class of compartment models are PBPK models that encompass the whole body [32,49].In contrast to the compartment models described above, which describe only part of the body and therefore require the measurement of an input function, whole-body models have considerable knowledge about the anatomy and physiological processes represented by the individual compartments.Thus, prior anatomical and physiological knowledge is integrated into the structure of the model.Additional prior knowledge can be integrated by using known values of parameters, e.g.age, body weight or transfer rates between compartments [33][34][35].
It is often incorrectly assumed that the fit function must have fewer parameters than there are measurement data.The fit function may have more parameters, but fewer parameters can be fitted than there are measurement data.Thus, several parameters need to be fixed based on a priori knowledge.For example, physical decay is a known parameter of the fit function that can be factored out [9,27].It is also beneficial to not perform decay correction before fitting a function to the data [29,30].

Error model function for time-activity data
The sources of the uncertainties in the time-activity data are (but are not limited to) the systematic uncertainty corresponding to the volume determination and quantifications, image noise, registration, patient motion, image reconstruction, scatter and attenuation corrections [24].If the uncertainties of the time-activity data are not known, error model functions can be used to estimate the uncertainties from the data.Therefore, appropriate variance models can be designed, for example, as follows [9,12]: Where A is a constant component of the error, B is a multiplicative component of the error, C is the weighting factor of the data point, y i is the data at time point i, f(t i ) is the function value at time point i, σ 2 i is the general form of the variance, ν is the scaling factor of the variance, ϑ 2 i is the variance of the model fitting for data (Eq.( 4)) or model (Eq.( 5)) based weighting is chosen by the user).
Different weightings are created by choosing parameters A, B and C, which can be set fixed (these data are input data, see section 2.1.2) or can be treated as fit parameters in curve fitting.A, B, and C can be chosen based on the suitable error model in their dataset: • For the constant error model and absolute weighting, the values are A = SD, B = 0, and C = 0, with SD corresponds to the standard deviation of the data.• For the proportional error model and absolute weighting, the values are A = 0, B = FSDsquare, and C = 0, with FSD being the assumed or fitted fractional standard deviation of the data points.
Note that for a relative weighting fitting, the variance is estimated from the data and calculated based on equation (3).A fractional standard deviation of 10 % is achieved when setting (A, B, C)=(0, 0.01, 2) and a constant weighting by (A, B, C)=(1, 0, 0).If the parameters of the error model function are also unknown, e.g., ν, they can be fitted as free parameters when minimising the objective function (section 2.2.3).

Objective function
There are infinite possible curves or functions that "fit" a pharmacokinetic dataset to a certain degree.Therefore, a measure of the fit to the data must be defined, often referred to as "objective function".
The objective function usually is the sum of the weighted squared differences of the data points and the fit function [9,36].Also, a priori knowledge (if available) about the parameter values of the fit function can be integrated as additional sums of squares as follows (e.g., based on previous experiments or publications): where P is the likelihood, i denotes the data points, N is the total number of data points, j denotes the parameters, K is the total number of parameters, y i and f(t i ) are the measured data and the function values at time t i , σ i 2 is the variance of data point i, p j is the current value of parameter j, p j is the average population value, and ω j 2 is the associated variance.The Bayes terms (right sum in Eq. ( 6)) correspond to additional measurements of the parameters with the given uncertainty.As it is possible to construct an infinite number of fit functions to be used in Eq. ( 6), it is mandatory to find the "fit function most supported by the data".Usually, several fit functions are tried out and often the bestsuited fit based on experience is selected.The functions should be physiologically sound and the noise in the data should not be fitted.Although such criteria help to reduce possible fit functions, there is still much user interaction, which is why the reproducibility is not optimal.Therefore, model selection criteria should be applied to select the best fit function supported by the data.A prerequisite for the reproducibility of results is the agreement on such a measure.
Very importantly, by choosing a fit function, the user also chooses a priori information to be included in the result.For this reason, it is crucial to use a set of well-suited functions for fitting.

Algorithms for integration
There are a variety of fitting algorithms that solve the fitting problem using different approaches and with different success probabilities [37,38].Therefore, it is essential, on the one hand, to choose a suitable algorithm and, on the other hand, to ensure that the final result presented by the algorithm is a good fit (section 2.3.2).

Data processing
After collecting the input data used for the TAC modelling and making decisions about the type of function to accommodate, one enters a phase of data processing and actionable TAC optimisation (see Fig. 2).During this phase, one must consider how to actually fit the chosen function (e.g.software), what parameters and uncertainties can be calculated/estimated to define the goodness of the fit, and how to make an informed choice about the most appropriate fit results for the data one has.These aspects are covered in chapters 2.3.1 to 2.3.3.

Fitting
There are many software options commonly used for curve fitting, such as simulation analysis and modelling software (SAAMII) [38], Olinda/EXM [39] and general platforms such as MATLAB, Python etc..The fitting process starts with preparing the biokinetic data and choosing an appropriate error model for the data (section 2.2.2.).Inputs of the error model that need to be set are: 1.The type of the error model.The following error models for the fit can be used: constant, proportional, combined, exponential.2. The degree of the error model: whether to fix the degree of error with the absolute error model, e.g. a fixed fractional standard deviation (FSD) to 10 % (ref), or estimate the error from the data with a relative error model.3. The reference of the error model: whether the error model is based on the data or the model.
The next step is to define the computational settings, such as the fitting algorithm and the convergence criterion (section 2.2.4).After passing the function/model structure to the software (section 2.2.1.),the user has to define the starting value of the adjustable parameters and their limits in the software.Good starting values of the parameters can accelerate the convergence of the fitting algorithm and avoid local minima of the objective function during the fitting process.
The fitting process is started by minimising the objective function in Eq. ( 6).If the user has prior information about the adjustable parameters, the objective function with Bayesian terms can be used in the fitting process.
In general, during each iteration of the optimisation process, the algorithm searches for values of the parameters which reduce the value of the objective function.The fitting process stops when the convergence test is satisfied.Convergence is reached when: 1. the relative error between the objective function value at the end of ith iteration ( − 2ln(P)) and the minimised objective function ( − 2ln(P′)) is less than or equal to the convergence criterion δ 2. the absolute difference between the value of the adjustable parameter and the parameter at i-the iteration is smaller than the product of the convergence criterion δ and the absolute difference between upper and lower limit of the parameters.
The convergence criterion δ can be chosen between 10 − 4 to 10 − 7 [40].Function parameters are estimated in the fitting process.As described above (section 2.2.3), the least squares approach with a nonlinear regression technique can be used to minimise the objective function (Eq.( 6)), which leads to a minimum value of the sum of squared error SSE (Eq.( 7)) Here, y i is the measured activity at time ti, and f(x i ) is the output of the function at time ti.As suggested in the literature, the weighting of the activity uncertainty is not included [24].The covariance of the fitted parameters is calculated using the following equation: where COV is the covariance matrix of the fitted function parameters, dof is the degree of freedom calculated as N-K with N as the total number of the activity data, and K is the total number of the fitted function parameters, Jp is a Jacobian matrix of the function used for the fitting (see Supporting information equations s.1-s.2).The standard deviation of the time-integrated activity coefficient (TIAC) is calculated using Gaussian error propagation, taking the covariances of the fitted parameters into account.
A result without a degree of uncertainty devalues, as the uncertainty determines the quality and reliability of the result.Sometimes, the error function for weighting the data accuracy may be unknown.In this case, the fitting may be performed repeatedly for different weightings.If the results do not strongly depend on the chosen weighting, the user can trust the obtained result.At last, performing a fit must always include a goodness-of-fit check, as described in the next section.

Goodness of fits
A check of the quality of the fit function in describing the biokinetic data set is recommended.Therefore, the analysis of the goodness of fit is needed.Verification of the fit accuracy requires at least all low-level checks in the following list, and summarised in more brief form in Table 2 [9]:

Visual inspection of the fitted curves (low level)
The goodness of fit is performed by inspecting the quality of the fitted curves.A good fit corresponds with a fitted curve passing through or near most data points and has random trends [53].This visual test is a powerful diagnostic tool as it is straightforward to understand and often detects patterns or anomalies in the data.

The estimated parameters of the model should be plausible (low level).
For example, if the parameters stand for physiological quantities, they should be within the physiologically possible range and if the parameters stand for the transfer rates or biological decay, they should all be positive.

The value of the coefficient of determination R 2 (low level).
Coefficient of determination or R 2 is the proportion of the variance for the dependent variable that is explained by independent variables.The value of R 2 is between 0 and 1.This value is calculated according to the following equation: where SSR is the sum of squares explained by the model and SST is the total variation of the data.A value of R 2 close to 1 corresponds to a good fit.However, an evaluation of R 2 alone may be an inadequate measure of the goodness of fit [40].

Coefficients of variation (CV) of the fitted parameters (mid-level).
The coefficient of variation of the fitted parameters is a statistical measure of the dispersion of the parameter around the mean.The coefficient of variation is calculated as the ratio of the standard error to the mean value of the fitted parameter.A low number of the CV corresponds to an accurate determination of the function parameters.In contrast, a high number of the CV corresponds to an inaccurate determination of the parameter.Based on the literature, fitting with a CV of the fitted parameters less than 25 % and 50 % is considered precise and acceptable [40], respectively.To pass the goodness of fit criteria, the CV of the fitted parameters of a function should be no higher than 50 %.

Correlation matrix (mid-level)
The correlation matrix is a normalised form of the covariance matrix.The correlation matrix shows the degree of correlation between pairs of parameters.The matrix is squared and symmetric, with the diagonal matrix elements being 1.The matrix elements vary from + 1 (perfect positive correlation) to − 1 (perfect negative correlation).Elements of + 1 and − 1 indicate that the two parameters have opposite and similar effects on the model, respectively.The element value of zero corresponds to two parameters with no correlations.The element of the correlation matrix is calculated as follows: σ X σ Y (10) where Corr(X, Y) is the correlation matrix of parameter X and Y, Cov(X, Y) is the covariance matrix of parameter X and Y, σ X is the standard deviation of parameter X, and σ Y is the standard deviation of parameter Y (see SI). High correlation elements tend to widen the confidence intervals.The ideal value of the off-diagonal element of the correlation matrix is close to zero, which shows that the parameters of a function are uncorrelated or independent of each other.Based on the literature, to pass the goodness of fit criteria, most of the elements of the off-diagonal matrix of the fitted parameters of a function should vary from − 0.8 to 0.8 [40].

The weighted residuals should be randomly distributed (midlevel).
If they show a systematic trend around the zero line, the model has not captured all information in the data.7. To ensure that the final result is reliable, the uncertainty of the calculated TIA value should be verified.Gear et al. [24] introduced uncertainty a guideline for the whole RPT dosimetry pipeline.In this particular case, we are focusing on uncertainty aspects associated with the TAC fitting process only.For example, this includes consideration of an underlying noise model type of the data when the TAC is fitted, or the numbers of degrees of freedom that the fitted model can have.Although TAC fitting-associated uncertainty is often overlooked and not reported consistently, it is a major quantitative factor clearly illustrating how reliable the calculated TIA value is (see Fig. 3,4 and Table 3).
Always several functions are more or less suited to fit given data points.As described in the next section, model selection algorithms can be used to reproducibly select the fit function most supported by the data [42,43].

Model selection
In RPT, it is often unclear which function best approximates the biokinetic data [44].For this, it has been shown that model selection is an essential and critical aspect of scientific data analysis [42].Model selection could increase the reproducibility of the fitting results in contrast to using a rule of thumb by the user.Different criteria exist to perform the model selection, such as the F-test, the Bayesian Information Criterion (BIC) [45], and the Akaike Information Criterion (AIC) [46].It has been shown that the AIC performs better than other criteria for model selection in molecular radiotherapy (Fig. 3).The AIC criterion estimates the Kulback-Leibler distance, quantifying the information loss if a model only approximates the true model [41,42].The corrected AIC (AICc) has been introduced by Hurvich and Tsai [46].It adapts the AIC for the case of N/K < 40 and has been shown as an appropriate method for model selection in RPT.The model selection for all functions passing the goodness-of-fit test with the AICc method is started by calculating the Akaike weight [41,42].The Akaike weight of each function is calculated as follows: Δ i = AICc i − AICc min (12) where P is the estimated objective function minimised for the fitting, AICc min is the lowest value of all fitted functions, Δ is the difference of the AICci of function i and the AICc min , F is the total number of functions that passed the goodness-of-fit test and ω AICci 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 [42].To analyse the stability of the model selection for different data sets, the Jackknife method can be used [42,43].

Data output
As listed in sections 2.1 to 2.3, as the output of the time-activity data fitting, the following information should be provided:

Practical examples and pitfalls
Using methodological information from Sections 2.1 to 2.4, we provide several practical examples, illustrating the importance and actual effect of different TAC fitting steps.
To illustrate these effects, we selected synthetically generated data representing tumour uptake of [ 177 Lu]Lu-PSMA.Noise-free data are used for the examples so that the underlying truth is known.The problems described below are, therefore, even more pronounced when, as in reality, measurement errors are present.The data is available in the Support Information (SI) file and can be used by the reader to follow our calculations independently.The ground truth TAC for the illustrative examples in sections 3.1 and 3.2 was determined based on Eq. ( 17) with parameters given in Table 3.In the following paragraphs, we illustrate how kinetic data, corresponding to a realistic measurement scheme, can be used to select suitable TAC models and how the use of different distributions of measurement points for TAC improves the calculated TIACs and the stated goodness-of-fit criteria from Table 2.

Model selection criteria
First, we show how the most suitable TAC fit model can be selected using the selection criteria listed in section 2.3.3.In this example, we assume that 6 post-injection time points (2,20,44,68,116, and 164 h p. i.) are available to fit the tumour's biokinetics (TAC).These data were collected using SPECT/CT after intravenous injection of ~ 6 GBq of [ 177 Lu]Lu-PSMA-617.
Existing prior knowledge should be considered when defining a set of functions for model selection, as suggested by Burham et al. [41].Therefore, sums of exponential functions are considered for explaining the biodistribution of the radiopharmaceuticals and the constraint f(t = 0) = 0 has to be assumed for the organs.A rapid increase of activity in the organs with a half-life of 1 min due to the blood circulation time in humans could be added as well as suggested in the literature [48].
Given the number of points available for the fit, sums-ofexponentials models were evaluated as plausible candidates, using varying degrees of prior knowledge included in the fit (four models in total).These were: two-parameter (f 2P (t)), three-parameter (f 3P (t)), four-parameter (f 4P (t)), and five-parameter sums of exponentials (f 5P (t)).These functions are listed in Eqs. ( 14) to (17), respectively.Each model includes an exponential factor modelling the known physical decay.
where f nP is the sum of exponentials function with n adjustable parameters, A i are the coefficients of the respective exponential terms with values ≥ 0, λ i are biological rate constants with values ≥ 0, λ bc is the rate of blood circulation approximated as λ bc = ln (2)  1min [47], and λ phys represents the radionuclide physical decay constant.The rate of the blood circulation λ bc reflects the time needed for the administered activity to reach the organ via the blood flow; thus it is not tracer-specific (but may depend on the organ).As the time of the first measurement in RPT is large compared to this blood circulation time constant, the actual value is not relevant.However, the corresponding exponential term is needed to ensure that the functions fulfil the constraint f(t = 0) = 0.
Full information on TAC fit types and prior information incorporated into the models, and their goodness-of-fit performance can be found in Table 3.In Fig. 3 one can visually evaluate the goodness of fit for the models from f 2P (t) to f 5P (t) as well as their residuals' distribution.
In addition to the fit parameters listed in Table 3, a pairwise F-test was performed to verify the best TAC model, as one of the methods listed in chapter 2.3.3.The results of this comparison can be found in Table 4.As expected, the stepwise execution of the F-test confirms that the fiveparameter bi-exponential model f 5P (t), which was used to define the ground truth, best fits the noise-free data.

Distribution of data points
Another important consideration from a logistics perspective of data collection is how selection and distribution of the data points used to fit the TAC affects the fitted curve, ultimately defining the TIA values and uncertainty and the subsequently calculated absorbed dose.To illustrate this effect, we selected the same set of six data points, representing synthetically simulated lesion uptake of [ 177 Lu]Lu-PSMA-617.
First, since obtaining 6 data points is a major challenge in clinical practice, we show how measurements with 4, 3 or even two data points affect the accuracy of calculated TIA values.To this end, we have fitted the same two-parameter mono-exponential model (f 2P ) to measurement points listed in Table 5.Although the model with two parameters (f 2P ) did not yield the best fit results in section 3.1, it is the only type of fit of those tested in this manuscript that can be applied to a dataset containing only two measurements.The results of this comparison can be seen in Fig. 4, A-C.
Second, we want to emphasise that all examples before this point assumed that the data provided for the fit are free of noise.Although noise-free data are useful for mathematical modelling purposes, this is an oversimplification that is not generally valid in clinical practice.Therefore, in Fig. 4D-4F, we show how different realisations of the noise that may be present during the measurement can influence the results of TAC adjustments.For the noise simulations, a 10 % standard deviation proportional error was utilised using Graphpad Prism software (more information is valuable in the SI, noise realisation).The results of this calculation and corresponding goodness of fit data are available in Table 5.

Discussion and Conclusion
In this work, we have shown the importance of selecting accurate fit functions for time-activity data.In terms of steps to be performed in the curve fitting process itself, we demonstrate a detailed workflow with inputs and dependencies, and the mathematical framework to be considered in curve fitting.In addition, we expanded on each of the data inputs that require consideration by analysts, such as the number and timing of data point acquisitions, data required for a fitting algorithm and details on the types of fit functions themselves.We then provided background on the fitting process itself, such as convergence criteria and eventual assessment using goodness-of-fit metrics, and demonstrated the process chain through a series of curve fits of different models and time points to synthetic noise-free and noisy biokinetic data with a shape and %aA/cm 3 timeframe typically encountered from a prostate lesion in a [ 177 Lu]Lu-PSMA-617 dosimetry study.We derive the following conclusions from the application of our workflow to this example dataset: • A five-parameter bi-exponential model (Eq.( 17) fits the 6-time point data best, according to the F-test comparison for 4 physiologically relevant models.• A reduction of sampling time points to 4, 3 and 2 data points (keeping only the later time points) demonstrates that by selecting an appropriate model, a minor reduction in TIA of 1.6 % compared to the ground truth can be achieved even with only retaining the final 2 data points.The uncertainty of the result is, however, high.• Addition of different noise realisations to the same 3 timepoint data shows a markedly different TIA as compared to the ground truth, and an increase in CV and RMSE.
An important consideration within the process, and one that is often overlooked in the collection of imaging timepoint data, is the inclusion of a robust error analysis within dosimetry calculations.In their detailed discussion about the sources of errors in the calculation of absorbed dose from quantified images, EANM dosimetry guidance recommends propagating estimates of the various quantities (count rate, calibration factor and recovery coefficient, etc.) and their associated uncertainties through the calculation process in order to obtain an estimate of mean absorbed dose coupled with an associated uncertainty [24].Within the TAC derivation process specifically, there are error sources such as image reconstruction, patient motion, registration, and attenuation correction to name a few, and error bars for each data point can be calculated from the standard uncertainty of the measured activity.They also recommend that TAC data should be plotted using 95 % confidence intervals due to systematic uncertainty in activity combined with the parameter uncertainties in the fitting algorithm itself.A clinical example demonstrates that up to 40 % fractional error in absorbed dose may be obtained.
As earlier described, it remains important to choose the sampling schedule carefully dependent on a variety of biokinetic factors, such as the biological half-life of the radiotracer.Logistical roadblocks can also be a concern, such as clinician support for a dosimetry schedule, willingness/ability of the patient to return for extra scanning, scanner availability and reimbursement of scanning costs to the healthcare facility.Our results in Fig. 4 demonstrate that although some time points can be removed from the fitting process, down to only 2 datapoints with a slight change in calculated TIA, this comes at the cost of not being able to properly verify the quality of the fit to two points, and relying on visual observation only (Table 5) which is suboptimal for dosimetric calculations.Care should be taken to balance the desired acquisition schedule with the mathematical reliability of the curve-fitting workflow.
In relation to the above point of verification of fit accuracy, our recommendations in Table 2 range from a low qualitative assessment, such as visual observation, up to a medium quantitative assessment by correlation matrix calculation and calculation of the random distribution of weighted residuals, finishing with a high-level calculation of the associated uncertainty.As a general recommendation, we propose that all 7 quality control measures be carried out for each curve fit where possible.However, we understand that this might be challenging to achieve in the routine clinical practice of non-specialized hospitals.We therefore propose that at least all low-level controls, namely the visual inspection, the plausibility of the estimated parameters, and the coefficient of determination, should be verified for each tissue and organ at risk for all RPTs.These controls should be sufficient if a commonly used RPT agent and a generally accepted imaging scheme are used (e.g., the EANM 177 Lu-PSMA guideline).Furthermore, this will not burden the clinical workflow as the main CE-marked software packages used for RPT dosimetry (incl.Hermes, MIM, DosiSoft) already have them integrated into their workflows.At the same time, we would like to emphasise that in studies that are expected to provide new insight into

Table 5
Basic quantitative information for TAC fits listed in Fig. 4. The type of fits, TIA , uncertainty and R 2 parameters are provided.the pharmacokinetics of new RPT agents, we recommend all quality controls listed in Table 2 to be performed.As shown in Table 5 and Fig. 4, TAC fits with apparently identical low-level control values (R2, visual appearance) can have significantly different uncertainty of the results and, as a result, lead to different conclusions in the study.Up to now, not all medium-or high-level controls can be performed in CEmarked software packages.That is why we would like to bring readers' attention to the fact that more elaborative TAC-modeling software should be used for agents that are not yet well studied and that more data analysis options should be incorporated in the CE-marked software.Given that clinical studies using new RPT agents will mainly occur in clinical research settings, it is acceptable to allow longer processing times required to generate TAC curves.Another approach to further improving TAC fit that should be considered, but is not covered in the manuscript, is the ability to include population-based data either in the fit or in the model selection process.In recent work [49,50], population-based model selection approaches have been shown to perform better and result in higher robustness or outcomes compared to individual-based selection approaches.However, coverage of this topic is beyond the scope of the current manuscript.Additionally, we would like to emphasise that, although this work explicitly focused on the good practice of TAC fitting for organ-based dosimetry applications in RPT, more and more centres are using voxelised dosimetry packages.With voxelised dosimetry, however, it will no longer be possible to perform all checks mentioned in Table 2 for every voxel-organ combination.Therefore, new quality measures and good practice guides need to be developed for voxelised RPT dosimetry applications.
In summary, we have presented a detailed workflow to be considered within the mathematical framework of fitting curves (models/functions) to data for dosimetry applications.Each step of the workflow is described, and through example data, demonstrates the pros and cons of a reduced timepoint acquisition scheme and the effect of noise in the measurements.We also provide specific recommendations for the assessment of fitting through quantitative and qualitative control measures.We surmise that this workflow and associated examples could be an invaluable aid to those analysing nuclear medicine dosimetry data in terms of how to improve reliability within a sparsely sampled dataset.

Ethics approval and consent to participate
Not applicable; no patient data were used in the manuscript.

Fig. 1 .
Fig.1.The number of scientific peer-reviewed publications on the subject of RPT for a selection of the most commonly utilised beta particle emitting radionuclides (generated via combined ScienceDirect database and PubMed searches, using radionuclide symbol combinations, and the keyword 'therapy').

Fig. 3 .
Fig. 3.A -visual illustration of four TAC models tested in the model selection example (f2p, f3p, f4p & f5p, corresponding to Eqns 14-17, respectively).B-E -residual plots of the models.Note that f4p and f5p functions are highly similar, and are overlaid in panel A. A systematic deviation in the residual plot of f2p (panel B) demonstrates that the fit function used is inadequate.The same applies to function f3p to a lesser degree.Note the different scales of the residual plots, showing the strong reduction of residuals for functions with more parameters.

1 .
The list of fit functions investigated in the model selection process a. Akaike weight (ideally incl.uncertainty) of the function best supported by the data 2. The fit function best supported by the data b.Parameters of the function with their uncertainties (Section 2.3.1) 3. The goodness of fit of the best fit function c. Visual inspection (Graphs of the data and fitted curves of the best function) d.Visual inspection of the residuals e. Coefficient of variations of the fitted parameters f.Covariance and correlation matrix of the parameters g.Estimated intra-individual variability of the data (Section 2.2.2) if a relative weighting model was chosen.4. Value and uncertainty of the TIA based on the best fit function from the ground truth value.

Fig. 4 .
Fig. 4. Result of TAC fits for lesion uptake of [ 177 Lu]Lu-PSMA-617.Fit curves computed using noise-free data with 4 (A), 3 (B), or 2 (C) data points with the twoparameter fit function (f 2P -solid curves).D-Ftwo-parameter fit function using 4, 3 or 2 data points and three types of noise realisation (see SI).In panels A and D the fit to the three-parameter function (f 3P -dashed curves) is also shown for comparison.

Table 1
Examples of selected clinical studies and the time points utilised for TIAC computation.

Table 2
Recommended quality control measures necessary to verify the accuracy of the TAC fit.At least low-level controls should be performed every time.

Table 3
Parameters of TAC models used in the "model selection" examples.

Table 4 F
-test comparison for four fit models tested in the "model selection criteria" section.
*model listed in the hypothesis suit the data the best.O.V.Ivashchenko et al.