Appendices
Additional information on the compartmental model
The model is based on three key variables: administered activity (A), tumor cell count (N), and thyroglobulin concentration (\(T_g\)), all of which interact dynamically over time (see Eqs. 1a, 1b, and 1c). Several parameters are essential to the model, including the delayed effectiveness rate of iodine (a), expressed in \(\hbox {months}^{-1}\)29; the immediate effectiveness rate of iodine (\(\rho\)), expressed in (GBq \(\times\) months)\(^{-1}\)17; the tumor cell doubling time (\(T_d\)), expressed in months30; the rate of thyroglobulin production per tumor cell (\(\lambda\)), expressed in \(\mu\)g\(\times\) \(\hbox {l}^{-1} \times\) \(\hbox {month}^{-1}\)31; the \(T_g\) elimination rate from the bloodstream (\(k_e\)), expressed in \(\hbox {month}^{-1}\)32; and the initial variables: \(A_0,N_0\) and \(T_{g0}\)33. The average values of these parameters are presented in Table 1 and will be used for subsequent analysis (Fig. 6).
$$\begin{aligned} \frac{dA}{dt}(t)&= -a \times A(t), \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \,\,\text {with} \quad A(0) = A_0 \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \end{aligned}$$
(1a)
$$\begin{aligned} \frac{dN}{dt}(t)&= \frac{\ln (2)}{T_d} \times N(t) – \rho \times A(t) \times N(t), \quad \text {with} \quad N(0) = N_0 \quad \end{aligned}$$
(1b)
$$\begin{aligned} \frac{dT_g}{dt}(t)&= -k_e \times T_g(t) + \lambda \times N(t), \quad \quad \quad \quad \quad \text {with} \quad T_g(0) = T_{g0} \quad \end{aligned}$$
(1c)
Fig. 6The alternative text for this image may have been generated using AI.
Three-compartment pharmacokinetic model for radioiodine therapy. Schematic representation showing the dynamics of radioactive iodine (A), thyroid tissue (N), and serum thyroglobulin (\(T_g\)) with their governing differential equations and interaction pathways. Dashed arrows represent biological interactions while solid arrows indicate decay/clearance processes.
Table 1 Population parameter estimates related to Eqs. (1) Obtained Using the MCMC-SAEM Algorithm18 Implemented in Monolix® Software19.
Description of the study population
A comprehensive sample of 50 patients was used, along with an in-depth comparative analysis between empirical clinical data and model predictions. In this section, we delineate the essential characteristics for the development of a customized RAI dosing model. To ensure robust model verification, it is important to employ high-quality data and clearly define the inclusion criteria. The dataset utilized in this study exhibited the following features: the sample size was \(N = 50\), with a mean patient age of 29 years. Regarding disease stage, there were 7 patients classified as pT1/pT2, 43 patients as pT3/pT4, 47 patients with N1 status, 3 patients with Nx status, and all 50 patients exhibited M1 status. The retrospective inclusion criteria for the study were total thyroidectomy, lymph node dissection (except for three patients), and RAI therapy. The parameters of the RAI administration protocol were as follows: the administered activity ranged from 3.7 to 5.5 GBq, with an average cumulative activity of 22.2 GBq. The frequency of administration was fixed at intervals of every 6 to 8 months. Within the sample, \(\sim 70\%\) of patients were classified as having responded to treatment. This classification was based on specific criteria, including a reduction in stimulated \(T_g\) concentration values, theoretically equivalent to a decrease in tumor uptake or iodine-131 avid foci, the absence of new lesions and the absence of disease progression during radiological assessments. The results of the model can be used to predict a patient responder or non-responder status after two or three RAI therapies and three \(T_g\) measurements. A key discriminating parameter identified in this study is the tumor cell doubling time (\(T_d\)), estimated from patient treatment data, which effectively distinguishes between responders and non-responders to RAI. In the cohort of 50 patients analyzed, the mixture mathematical algorithm implemented in the Monolix software categorized patients into two significantly distinct groups. The responder group comprised 72.5% of patients, with an average \(T_d\) of 66.6 months. These patients exhibited a kinetic profile of \(T_g\) concentrations that progressively decreased over time. In contrast, the non-responder (refractory) group represented 27.5% of patients, with an average \(T_d\) of 9.8 months. This group displayed a kinetic profile where \(T_g\) concentrations initially decreased but subsequently increased, indicating disease recurrence. The model facilitates the determination of a patient responder or non-responder status after just three radiation therapy sessions, compared to the conventional seven or eight sessions. For non-responders, this enables the proposal of alternative treatment protocols, allowing timely adjustments to the therapeutic approach. For responders, the intervals between RAI sessions can be extended to optimize treatment, balancing efficacy and minimizing toxicity, thereby reducing potential side effects.
