# A Mechanistic Dynamic Surrogate Framework for Personalized Radioactive Iodine Therapy in Metastatic Thyroid Cancer

Marie Fusella Giuntini, Cyril Voyant, David Taieb, Dominique Barbolosi

► To cite this version:

Marie Fusella Giuntini, Cyril Voyant, David Taieb, Dominique Barbolosi. A Mechanistic Dynamic Surrogate Framework for Personalized Radioactive Iodine Therapy in Metastatic Thyroid Cancer. Precision Oncology, 2025, 1 (1), 10.3390/precisoncol1010002. hal-05409930

HAL Id: hal-05409930

https://hal.science/hal-05409930v1

Submitted on 10 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# A Mechanistic Dynamic Surrogate Framework for Personalized Radioactive Iodine Therapy in Metastatic Thyroid Cancer

Marie Fusella Giuntini <fnref n="1" />, Cyril Voyant <fnref n="2" />,*, David Taieb <fnref n="3" /> and Dominique Barbolosi <fnref n="4" />,*

## Abstract

**Background/Objectives:** Radioactive iodine (RAI) therapy is widely used to treat metastatic differentiated thyroid cancer. To investigate physiological determinants of treatment response, a mechanistic model was developed, formulated as a system of coupled ordinary differential equations. **Methods:** The model captures the interactions between tumor burden, thyroglobulin ($T_g$) production and clearance, and radioactive iodine activity within a pharmacokinetic-pharmacodynamic framework. Model parameters were estimated using the Monte Carlo Stochastic Approximation Expectation-Maximization (MCMCSAEM) algorithm, based on clinical data from a cohort of 50 patients. **Results:** Tumor radiosensitivity ($\rho$) and initial tumor burden ($N_0$) consistently emerged as the most influential factors in both responder and non-responder groups classified by disease doubling time under RAI ($T_d$). A reduced model using only these two parameters preserved the principal response patterns of the full model. Other parameters influenced transient dynamics but had limited effect on overall $T_g$ variance. **Conclusions:** These results support the use of a reduced calibration approach focused on $\rho$, $N_0$, and the effective doubling time $T_d$. The findings establish a theoretical foundation for developing tractable dynamic surrogates that reproduce the main treatment kinetics and support model-based clinical decision-making in RAI therapy.

**Keywords:** cancer; thyroid; optimization; personalized; RAI

## 1. Introduction

Thyroid cancer ranks as the most common endocrine malignancy. Its incidence has been increasing in recent years due to the improvement of screening tools (high-resolution ultrasound), sometimes leading to overdiagnosis [1]. Most cases of thyroid cancer have a favorable prognosis. Radioactive iodine (RAI) therapy has long been established as a standard treatment [2]. It plays an important role in the management of thyroid cancer, particularly in patients with metastatic disease [3]. However, each patient response to treatment is driven by specific biological and clinical characteristics, highlighting the need for precision-medicine approaches and model-based decision-support tools to optimize the efficacy-toxicity balance of RAI therapy.

<figure><img src="image_2.png" /></figure>

check for
updates

Academic Editors: John F. Deeken and
Raul A. Urrutia

Received: 17 September 2025

Revised: 29 October 2025

Accepted: 5 December 2025

Published: 10 December 2025

**Citation:** Giuntini, M.F.; Voyant, C.; Taieb, D.; Barbolosi, D. A Mechanistic Dynamic Surrogate Framework for Personalized Radioactive Iodine Therapy in Metastatic Thyroid Cancer. *Precis. Oncol. 2025*, **1**, 2.
https://doi.org/10.3390/precisoncol1010002

**Copyright:** © 2025 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

<footnote n="1">¹ SPE Laboratory, University of Corsica, 20250 Corte, France; fusella_m@univ-corse.fr</footnote>

<footnote n="2">² OIE Laboratory, Mines-PSL, Sophia-Antipolis, 06904 Antibes, France</footnote>

<footnote n="3">³ Department of Nuclear Medicine, La Timone University Hospital, Aix-Marseille University, CNRS, CINAM, ERL Inserm U1326, CERIMED, 13005 Marseille, France</footnote>

<footnote n="4">⁴ COMPutational Pharmacology and Clinical Oncology Department, Inria Center of the University of Côte d'Azur, Cancer Research Center of Marseille, Inserm UMR1068, CNRS UMR7258, Aix-Marseille University UM105, 13385 Marseille, France</footnote>

<footnote>* Correspondence: cyril.voyant@minesparis.psl.eu (C.V.); dominique.barbolosi@univ-amu.fr (D.B.)</footnote>
## 1.1. RAI and Metastatic Thyroid Cancer

Radioactive iodine (RAI), primarily administered as iodine-131 (I-131), is selectively taken up by thyroid cells via the sodium-iodide symporter. Once internalized, I-131 emits beta radiation that induces double-stranded DNA breaks, leading to cellular apoptosis or necrosis in iodine-avid tumor cells. RAI is usually administered after total thyroidectomy and lymph node dissection, according to prognostic indicators of recurrence or mortality [4]. This focused approach works best for well-differentiated thyroid cancers like papillary and follicular types. RAI targets remnant thyroid cells and can reach distant metastases in the lungs or bone. RAI is generally not indicated for low-risk patients (undetectable thyroglobulin, $T_g < 1$ ng/mL postoperatively) and selectively for intermediate-risk patients [5]. Additionally, a major controversy in RAI therapy concerns optimal dosage. Previous studies [3,6] have shown that higher doses of RAI may lead to salivary gland dysfunction, manifesting as sialadenitis (inflammation of the salivary glands) or xerostomia (dry mouth due to reduced saliva production) [7]. Such radiation-induced toxicities can manifest in the short, medium, or long term [8,9]. Furthermore, for patients who received a high dose of administered activity, the risk of developing leukemia is proven [10]. It is also noted that some patients with refractory thyroid cancer do not respond to RAI treatment [11–13]. These limitations emphasize the need for individualized treatment planning that balances efficacy and toxicity based on each patient's biological characteristics and tumor profile.

## 1.2. The Personalized RAI

Today, the importance of decision support tools in medicine is recognized, and more particularly in oncology, to offer personalized medicine and to propose a personalized treatment protocol. This kind of therapy should ideally minimize the risk of severe toxicity and avoid transient or partial responses, which may contribute to direct or indirect resistance after repeated RAI administration [14,15]. It is therefore essential to identify the individual determinants of therapeutic response. Studies have thus reported the existence of numerous factors influencing the response to RAI [16–19]. This work aims to identify those contributing to response variability. The characteristics of the protocol (dosage, number of RAI sessions, interval between sessions) directly depend on the decisions taken during multidisciplinary consultation meetings, considering the information known about the patient and relying on empirical data [20]. Incorporating patient-specific clinical and biological data could enable realistic simulations of the response to RAI. This study focuses on cases of metastatic thyroid cancer where the factors influencing the response to RAI, as already mentioned, are numerous. The doubling time ($T_d$) of tumor cells is a first important parameter in the response to treatment and thus makes it possible to define responder status more quickly (after only 3 fractions) (for a $T_d$ = 66.6 months) or non-responder (for a $T_d$ = 9.8 months) at RAI [11]. Importantly, $T_d$ was not arbitrarily chosen: in our cohort it emerged as the dominant dichotomizing marker among all candidates, with a clear bimodal structure in the mixture-model analysis (responders vs. non-responders). We also stress that the $T_d$ estimated here is an effective doubling time under RAI exposure, not the natural history growth rate of untreated disease. Unlike previous modeling approaches that focused exclusively on responders, the present analysis includes both groups, enabling a comprehensive characterization of the underlying dynamics and supporting personalized therapeutic strategies across the entire patient spectrum.

## 1.3. Clinical Interest of This Work

To deliver the most appropriate treatment for each patient (minimizing toxicity while maximizing efficacy), this study aims to provide clinicians with a model-based decision-
support framework capable of predicting therapeutic response. The approach focuses on identifying a minimal set of physiological parameters that govern RAI treatment dynamics and can be estimated from routine follow-up data. Although previous studies have shown that BRAF mutations can influence patient sensitivity to RAI [21], no genomic analyses were included in the present work. Instead, the analysis centers on measurable or inferable physiological quantities. Among them, the effective tumor doubling time under RAI exposure ($T_d$) is introduced as a quantitative parameter that naturally distinguishes responders from non-responders and provides a clinically interpretable handle for dose adjustment and treatment scheduling. This modeling strategy lays the foundation for constructing reduced mechanistic representations (or dynamic surrogates) that could ultimately support real-time personalization of RAI therapy.

Section 2 describes the model structure and sensitivity-analysis framework, Section 3 presents the results, and Section 4 discusses clinical implications and perspectives.

## 2. Material and Methods

The following section describes the methodological framework adopted in this study, including the mathematical formulation of the model, its analytical developments, and the approach used for sensitivity analysis.

### 2.1. *Mathematical Novelty and Originality of This Work*

This study extends the Ordinary Differential Equation (ODE) model of Barbolosi et al. (2017) [11] with major advances in mathematical formulation and clinical applicability. The original work relied solely on numerical integration and offered a qualitative exploration focused on $T_d$. This parameter was not predefined but rather emerged from a bimodal mixture distribution of patient estimates, naturally separating the cohort into two distinct subpopulations corresponding to responders and non-responders. $T_d$ reflects an effective doubling time under RAI exposure, which is not known prior to treatment and differs from the natural growth rate of the untreated disease. In the present study, analytical approximations are developed to better characterize model dynamics and derive explicit relationships between all physiological parameters and clinical outcomes. Because no parameter sensitivity analysis had previously been performed, a systematic quantitative study was undertaken to quantify the clinical impact of each parameter and to identify those requiring accurate measurement or estimation. These improvements turn the model into a decision-support framework for personalized RAI therapy, aimed at maximizing efficacy while limiting radiation-induced toxicity. $T_g$ is an observable biomarker used as the model output; parameter sensitivities were therefore evaluated with respect to the simulated $T_g$ response, providing insight into how physiological parameters influence measurable treatment dynamics. The model includes seven parameters, of which six are investigated here, since $T_d$ has already been validated as a clinically and mathematically discriminant variable. High $T_d$ values are associated with favorable outcomes (responders), whereas low values correspond to poor therapeutic response (non-responders), consistent with both clinical data and model simulations. In the present work, this dichotomy was further supported by the mixture model implemented in Monolix (version 2024R1), which revealed a bimodal distribution of estimated parameters (particularly $T_d$) corresponding to two statistically distinct subpopulations. This data-driven partition provides a robust basis for modeling responders and non-responders as separate groups in the sensitivity analysis. Accordingly, two patient groups were simulated in parallel: responders ($\bar{T}_d$ = 66.6 months) and non-responders ($\bar{T}_d$ = 9.8 months).
## 2.2. Data

Clinical data were collected over three consecutive years, up to 2016. The cohort consists exclusively of patients with metastatic differentiated thyroid carcinoma presenting pulmonary metastases, thereby minimizing tumor heterogeneity. Serum thyroglobulin ($T_g$) time series were used as the primary observable variable for model fitting. Each of the 50 patients contributed at least six $T_g$ measurements obtained during RAI follow-up, providing over 300 longitudinal data points. Parameter estimation was performed using the updated Monte Carlo Stochastic Approximation Expectation-Maximization (MCM-CSAEM) algorithm implemented in Monolix [22], which benefits from recent numerical and convergence improvements. All diagnostic checks (such as "Observed vs. Predicted" (Population vs. Individual), "Visual Predictive Check") were verified. Although the cohort is identical to that analyzed by Barbolosi et al. [11], all parameters were newly estimated to ensure methodological consistency and provide unbiased input distributions for subsequent sensitivity and simulation studies. Monolix was employed within a nonlinear mixed-effects framework to estimate population means and inter-individual variability, while MATLAB (version R2024a) was used to solve the system of ordinary differential equations describing the coupled dynamics of administered activity $A(t)$, tumor-cell population $N(t)$, and thyroglobulin concentration $T_g(t)$. The two environments were combined in a consistent pipeline: statistical calibration in Monolix provided the parameter distributions used as inputs for deterministic simulations in MATLAB. These simulations explored parameter perturbations and their impact on $T_g$ dynamics under controlled conditions. For comparability across sensitivity analyses, an initial $T_g$ value of 100 ng/mL and a single RAI dose of $A_0 = 3.7$ GBq at $t = 0$ were used in all experiments. The MATLAB code (available at: https://github.com/MarieFG49/PROJET-SIMU_RAI, accessed on 8 December 2025) was executed on a 28-node Intel Xeon 6230R cluster (2 × 26 cores @ 2.1 GHz, 192 GB RAM) using 96 parallel workers and the Symbolic Math Toolbox. Each simulation was completed in under 30 s per run. The ODE system was not solved for individual patient fits in this section but rather for representative parameter sets sampled from the calibrated distributions in order to characterize model sensitivity and qualitative behavior under clinically realistic conditions.

## 2.3. *Description of the Reference Therapeutic Response Model*

The model used in this study aims to optimize the efficacy-toxicity balance of RAI therapy in metastatic differentiated thyroid cancer. It is formulated under biologically grounded assumptions consistent with pharmacokinetic-pharmacodynamic (PK-PD) modeling principles. The administered iodine activity $A(t)$ is described by a mono-exponential clearance function reflecting iodine-131 pharmacokinetics. Tumor dynamics follow classical cytotoxic-response principles, where cell loss is proportional to both the viable tumor mass $N(t)$ and the effective radiation dose. Thyroglobulin ($T_g$) production is modeled as directly proportional to the number of viable tumor cells, and its clearance follows first-order kinetics. These assumptions are consistent with widely accepted PK-PD frameworks used in thyroid and other cancers [23,24]. Clinical $T_g$ time-series from 50 patients (described in Appendix A.1) were used to calibrate a system of three coupled ordinary differential equations representing (i) iodine decay, (ii) tumor-cell response, and (iii) $T_g$ kinetics (Equations (1a)-(1c)). The analytic solution of this system is provided in Appendix A.2. In this formulation, the administered activity $A(t)$ influences the number of viable tumor cells $N(t)$, which in turn determines the $T_g$ concentration. The model introduces an iodine efficacy rate, denoted $\rho$, representing the fraction of administered radioactivity that effectively acts on tumor cells to induce destruction. This parameter captures patient-specific sensitivity to RAI and accounts for inter-individual
variability in treatment response and it can be interpreted as an integrated measure of I-131 uptake and radiobiological effectiveness within iodine-avid tumor cells.

$$ \frac{dA}{dt}(t) = -a \cdot A(t), \text{ with } A(0) = A_0 \qquad (1a) $$

$$ \frac{dN}{dt}(t) = \left( \frac{\log(2)}{T_d} \right) \cdot N(t) - \rho \cdot A(t) \cdot N(t), \text{ with } N(0) = N_0 \qquad (1b) $$

$$ \frac{dT_g}{dt}(t) = -k_e \cdot T_g + \lambda \cdot N(t), \text{ with } T_g(0) = T_{g0} \qquad (1c) $$

Each component of the model was formulated to represent a physiologically grounded mechanism relevant to RAI therapy in metastatic thyroid cancer. The exponential decay in Equation (1a) mirrors the pharmacokinetics of iodine-131. Tumor-cell dynamics in Equation (1b) follow classical cytotoxic-response models, where cell loss is proportional to both tumor mass and absorbed radiation dose. Equation (1c) describes thyroglobulin ($T_g$) kinetics, assuming secretion proportional to viable tumor cells and first-order elimination. Together, these equations form a biologically interpretable and clinically consistent PK-PD model structure. The clearance rate of radioactive iodine (a) and the $T_g$ elimination rate ($k_e$) are illustrated in Figure 1.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Block diagram of the PK-PD model with three compartments corresponding to three variables (A, N and $T_g$) and related parameters.</figcaption></figure>

Model parameters are detailed in Section 2.4, and validation results are provided in Appendix A.3. Within this framework, patients are classified according to their $T_g$ dynamics. Responders are defined as those whose $T_g$ concentration decreases after RAI administration and remains below 10 ng/mL for at least 12 months, corresponding to a partial response in clinical practice. Non-responders show a $T_g$ rebound indicating treatment failure. This functional definition is consistent with previous studies and corresponds to a tumor doubling time ($T_d$) above or below a threshold reflecting slower or faster tumor progression. Figure 2 shows two examples of thyroglobulin kinetics for patients from our cohort: (a) responder patient with a $T_d$ estimate of approximately 75 months, (b) non-responder patient with a $T_d$ estimate of approximately 6 months. The model allows early identification of responder status: typically, after two to three RAI fractions and three $T_g$ measurements. In thyroid cancer, this typically means a reduction in tumor size or activity, leading to better clinical outcomes. Factors influencing response to RAI therapy include the type and stage of the thyroid condition, the dose of radioactive iodine administered and individual patient characteristics, physiological parameters that we wish to identify and estimate. Genetic factors may also play a role in determining responsiveness to treatment [25].
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Responder (green) versus non-Responder (red) cases from 2 patients of the cohort, fitted using the mechanistic PK-PD model. The responder is characterized by $T_d = 75.2$ months, $k_e = 0.173$, $\lambda = 1.02 \times 10^{-9}$, $\rho = 0.00407$, $a = 0.0169$, and $N_0 = 1.12 \times 10^9$. The non-responder exhibits $T_d = 6.37$, $k_e = 0.392$, $\lambda = 7.6 \times 10^{-9}$, $\rho = 0.00407$, $a = 0.0169$, and $N_0 = 1.12 \times 10^9$. Simulated $T_g$ trajectories (solid lines) closely match the observed data points, confirming the model ability to capture both decay and rebound phases.</figcaption></figure>

## 2.4. Model Parameters Characteristics

The model integrates seven physiological parameters governing the coupled dynamics of administered activity $A(t)$, tumor-cell population $N(t)$, and thyroglobulin concentration $T_g(t)$ (Equations (1a)–(1c)). Parameter estimation was carried out using the MCMCSAEM algorithm implemented in Monolix, combining stochastic approximation with Markov Chain Monte Carlo sampling to handle latent variables and inter-individual variability. The procedure was executed under strict convergence and reproducibility constraints, with all diagnostic criteria satisfied. The resulting parameter distributions constitute a statistically consistent dataset enabling downstream numerical analyses and operational validation. A structured sensitivity analysis was then applied as a data-driven diagnostic layer of the model, quantifying how parameter variations propagate through the simulated $T_g(t)$ trajectories. This analysis highlights the parameters that drive model uncertainty and therapeutic outcomes, thereby identifying levers for dose optimization and individualized control strategies. From a computational standpoint, the approach establishes a transparent pipeline linking patient data ingestion, parameter calibration, and system-level response analysis. Model parameters and reference values are summarized in Table 1.

<table><caption>**Table 1.** Population parameter estimates related to Equations (1a)–(1c) obtained using MCMC-SAEM algorithm [26] from MONOLIX® software [27] *.</caption><thead><tr><th></th><th>N<sub>0</sub></th><th>a</th><th>ρ</th><th>T<sub>d</sub></th><th>λ</th><th>k<sub>e</sub></th></tr></thead><tbody><tr><td>Mean</td><td>1.12 × 10<sup>9</sup></td><td>0.0169</td><td>0.00407</td><td>55.23</td><td>3.86 × 10<sup>-9</sup></td><td>0.319</td></tr><tr><td>Std</td><td>7.87</td><td>&lt;10<sup>-17</sup></td><td>&lt;10<sup>-18</sup></td><td>27.65</td><td>&lt;10<sup>-6</sup></td><td>0.379</td></tr></tbody></table>

<footnote>* The small standard deviations reported for $\rho$ and $a$ are not the result of data error but rather reflect constrained variability under the mixed-effects estimation process. These values emerge from limited identifiability, compensatory effects among parameters, and regularization within the estimation algorithm. This behavior is acknowledged and discussed as a limitation of the current model implementation.</footnote>
Where

* $A_0$: Initial activity dose administered (GBq [28]);

* $N_0$: Number of initial tumor cells [29];

* $a$: clearance rate of administered iodine (month$^{-1}$ [30]);

* $\rho$: effectiveness rate of iodine (GBq × months)$^{-1}$ [11]);

* $T_d$: tumor cell doubling time (months [31]);

* $\lambda$: concentration of $T_g$ produced by a tumor cell ($\mu g \times L^{-1} \times month^{-1}$ [32]). $\lambda$ is defined as the production rate of thyroglobulin per tumor cell and is expressed in $\mu g \cdot L^{-1} \cdot month^{-1}$. While $T_g$ levels are typically measured in ng/mL clinically, $\lambda$ quantifies the dynamic rate of $T_g$ accumulation in the blood and reflects a distinct, time-dependent biological process;

* $k_e$: elimination rate of $T_g$ in the blood (month$^{-1}$ [33]).

Additionally, during the simulations described below, sensitivity analysis was used to identify parameters that are poorly estimated or associated with high uncertainty, indicating where additional data or model refinement would be most beneficial. This step enhances the robustness and interpretability of the model, ensuring that it remains numerically stable and clinically meaningful within realistic parameter domains. From a theoretical perspective, sensitivity analysis is an essential component of any system of coupled differential equations describing biological processes. It provides a quantitative means to assess the relative influence of each parameter on model behavior, to evaluate the conditioning of the system, and to verify the structural identifiability of parameters under real clinical variability. In practice, it allows discrimination between parameters that primarily affect scaling (amplitude of the $T_g$ response) and those that drive temporal dynamics (response rate or rebound behavior). Parameters with minimal influence can be fixed without loss of generality, whereas highly sensitive parameters require individualized estimation to maintain predictive reliability. The model includes seven parameters: the initial administered activity $A_0$ and six physiological parameters characterized by their population mean and standard deviation estimated from 50 patients. Parameters with low variability (relative standard error < 2%), namely $N_0$, $a$, and $\rho$, were considered stable and could reasonably be fixed in clinical applications. While $T_d$ remains a discriminating indicator of patient response to RAI, it is equally important to evaluate the influence of the remaining parameters on the evolution of $T_g$ during treatment. Accordingly, a systematic sensitivity analysis was performed to assess the contribution of each parameter to $T_g$ dynamics over time. Simulations were implemented in MATLAB using the differential system previously defined, varying each parameter within physiologically acceptable ranges [34]. This approach provides a rigorous, data-informed framework for quantifying parameter influence and guiding the development of personalized therapeutic strategies.

## 3. Results

Although clinical definitions of RAI response continue to evolve, $T_g$ remains a relevant surrogate biomarker for assessing disease activity and treatment efficacy, particularly in differentiating refractory and non-refractory patients [35]. Sensitivity analyses were conducted separately for each group. For all simulations, $T_g$ dynamics were modeled in response to a single RAI administration at time zero ($t = 0$). A one-at-a-time approach was employed, wherein one of the six model parameters varied over a physiologically plausible range while the others were fixed at their mean values, as estimated from the full cohort. A parameter was considered influential when small variations within its physiological range produced marked deviations in the simulated $T_g$ trajectory, affecting key features such as peak amplitude, time to nadir, or overall decay profile. This criterion reflects the model local sensitivity, whereby parameters exerting stronger perturbations on these outputs are
deemed more impactful. To aid interpretation, $T_g$ values are expressed as percentages relative to their baseline, for example, $T_g$ (100%) representing the initial concentration while $T_g$ (50%) denotes the time derivative (rate of change) of $T_g$ when its value has decreased to 50% of the baseline level.

### 3.1. Sensitivity Analysis: Group of Responders ($\overline{T}_d$ = 66.6)

Figure 3 illustrates two distinct categories of parameters: those that have a strong influence on the temporal evolution of $T_g$ and for which patient-specific estimation is likely to be critical in optimizing therapeutic response, and those whose impact is present but more limited. Based on criteria mentioned above, the following parameters were identified as influential:

* $a$: clearance rate of radioactive iodine from the body, consistent with a first-order decay process rather than any biological delay in effectiveness (in month⁻¹) (Figure 3a);

* $k_e$: elimination rate of $T_g$ in the blood (in month⁻¹) (Figure 3b);

* $\lambda$: concentration of $T_g$ produced by a tumor cell (in µg × L⁻¹ × month⁻¹) (Figure 3c);

* $N_0$: Number of initial tumor cells (Figure 3d).

<figure><img src="image_5.png" /><figcaption>(a) Evolution of $T_g$ for different values of $a$ between [0.01;0.1]</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) Evolution of $T_g$ for different values of $k_e$ between [0.1;2]</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(c) Evolution of $T_g$ for different values of $\lambda$ between [10⁻⁹;10⁻⁸]</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(d) Evolution of $T_g$ for different values of $N_0$ between [10⁸;2 × 10⁹]</figcaption></figure>

<caption>Figure 3. Cont.</caption>
<figure><img src="image_9.png" /><figcaption>**Figure 3.** Sensitivity of $T_g$ trajectories to parameter variations in the responder group. Panels (a–f) show the influence of: (a) clearance rate of administered iodine $a$, (b) $T_g$ elimination rate $k_e$, (c) $T_g$ production per tumor cell $\lambda$, (d) initial tumor burden $N_0$, (e) iodine effectiveness rate $\rho$, and (f) initial administered activity $A_0$. In each simulation, the target parameter is varied while others are fixed at mean values obtained from mixed-effects modeling. The minimum of each curve (indicating nadir) is marked with a cross.</figcaption></figure>

In each simulation, one parameter was varied across its full physiological range while the others were fixed at the mean values of the posterior distributions estimated from the clinical dataset. Figure 3a–d illustrate the evolution of $T_g$ for different values of the parameters $a, k_e, \lambda$, and $N_0$. The trajectories clearly differ according to parameter variation, indicating that $a, k_e, \lambda$, and $N_0$ significantly influence the temporal dynamics of $T_g$. In all cases, $T_g$ initially decreases following RAI administration; however, depending on parameter values, either a subsequent increase occurs at variable rates (Figure 3a) or $T_g$ stabilizes at distinct levels (Figure 3b–d). Each subfigure also reveals that certain parameter values correspond to favorable therapeutic responses (characterized by a sustained or consistently low $T_g$) whereas others lead to elevated or non-declining $T_g$ profiles, indicative of suboptimal response. This underscores the need for accurate, patient-specific estimation of these parameters to simulate individual therapeutic outcomes and guide clinical decision-making. Notably, each curve exhibits a minimum point (indicated by a cross), after which $T_g$ begins to rise again, to varying degrees depending on the parameter. This inflection point is a key indicator for planning subsequent therapy sessions to prevent loss of tumor control. Parameters whose influence on the therapeutic response is less pronounced but still detectable are:

* $\rho$: iodine effectiveness rate (in (GBq × months)$^{-1}$) (Figure 3e);

* $A_0$: Initial activity administered (in GBq) (Figure 3f).

Following an initial decline phase that is identical across values for both parameters and lasts approximately 10 months, the $T_g$ concentration either continues to decrease gradually (particularly for the highest parameter values) or begins to rise slightly for the lowest ones. Across the entire range of values, $T_g$ levels remain between 5% and 20% of their initial concentration, which may correspond to favorable (lower $T_g$) or unfavorable (higher $T_g$) simulated responses. The local minimum of each curve, marked with crosses, indicate the transition point after which $T_g$ begins to rise again.

These observations, across all 6 sub-figures, highlight the crucial importance of precisely determining the values of these parameters for each patient. The measured $T_g$ concentrations and rebound times serve as crucial indicators to determine both the optimal
timing of the second RAI administration and the required activity dose (based on the measured $T_g$ value) to effectively manage the disease. The precise determination or estimation of the values of each specific parameter for each patient allows the protocol to be adapted to optimize the response to RAI, considering the variation in responses depending on the parameter values.

## 3.2. Sensitivity Analysis: Group of Non-Responders ($\bar{T}_d = 9.8$)

In such cases, evaluating the role of each model parameter in RAI response is essential for defining optimized therapeutic protocols, while respecting recommended constraints on cumulative administered activity and minimum intervals between successive treatments. The analysis relies on the temporal evolution of $T_g$ concentration as a function of parameter variation, under simulation conditions identical to those described in Section 3.1. Parameter estimates remain unchanged, and all simulations correspond to a single RAI administration ($A_0 = 3.7$ GBq). Figure 4 illustrates the impact of varying each of the six model parameters on the $T_g$ trajectory for a representative non-responder. Across all subfigures, $T_g$ initially declines during the first ten months before diverging according to parameter values, indicating that each parameter exerts a measurable influence on treatment dynamics, even in clinically defined non-responders. Certain parameter combinations nevertheless yield simulated responses characterized by sustained $T_g$ suppression. For instance, lower values of $a$, $\lambda$, and $N_0$ (Figure 4a,c,d), or higher values of $k_e$, $\rho$, and $A_0$ (Figure 4b,e,f), are associated with more favorable $T_g$ trajectories.

These results suggest that certain patients currently labeled as non-responders might benefit from adjusted treatment regimens. This highlights the potential of personalized parameter estimation to guide RAI protocol adaptation. In practical terms, simulations based on estimated parameters could inform clinical decisions, such as advancing or intensifying a second dose in patients with anticipated rebound or redirecting early to alternative therapies when poor response is predicted. While these are theoretical insights, they demonstrate how model-based exploration can support dynamic and individualized treatment strategies. Finally, the similarity in $T_g$ response patterns across different parameters reflects the structure of the model itself: multiple physiological inputs interact to shape a common dynamic output. Clinically, this implies that similar $T_g$ curves can arise from different underlying mechanisms. This reinforces the need for patient-specific parameter profiling to avoid one-size-fits-all treatment strategies and enable personalized therapeutic planning.

<figure><img src="image_10.png" /><figcaption>(a) Evolution of $T_g$ for different values of $a$ between [0.01;0.1]</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) Evolution of $T_g$ for different values of $k_e$ between [0.1;2]</figcaption></figure>

**Figure 4.** *Cont.*
<figure><img src="image_12.png" /><figcaption>(c) Evolution of $T_g$ for different values of $\lambda$ between $[10^{-9};10^8]$</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(e) Evolution of $T_g$ for different values of $\rho$ between $[0.001;0.01]$</figcaption></figure>

<caption>(f) Evolution of $T_g$ for different values of $A_0$ between [1;10]</caption>

**Figure 4.** Sensitivity of $T_g$ trajectories to parameter variations in the non-responder group. Panels (a–f) display the impact of varying individual model parameters: (a) clearance rate of administered iodine $a$, (b) $T_g$ elimination rate $k_e$, (c) $T_g$ production per tumor cell $\lambda$, (d) initial tumor burden $N_0$, (e) iodine effectiveness rate $\rho$, and (f) initial administered activity $A_0$. Each line represents a simulation where the corresponding parameter is varied across its observed range, while all other parameters are held at their mean values obtained from mixed-effects modeling. $T_g$ evolution is plotted over time (in months), and nadir points are marked for visual reference with a cross.

### 3.3. Global Sensitivity Analysis

To evaluate the overall influence of model parameters on treatment dynamics, a Sobol variance-based global sensitivity analysis was performed using synthetic simulations. This approach quantifies how uncertainty in each input parameter contributes to output variability, capturing both main effects and parameter interactions. The analysis considered six physiological parameters: tumor radiosensitivity ($\rho$), initial tumor burden ($N_0$), thyroglobulin production rate ($\lambda$), $T_g$ clearance rate ($k_e$), iodine clearance rate ($a$), and administered activity ($A_0$). Each simulation generated a surrogate measure of treatment response, defined as the area under the $T_g$ curve (AUC), reflecting cumulative tumor activity over time. Parameter combinations were sampled across physiologically plausible ranges using the Saltelli method, and Sobol indices were computed following the standard formulation [36]. The first-order index ($S_1$) quantifies the direct effect of a parameter on output variance, while the total-order index (ST) captures both its direct and interaction effects. Table 2 summarizes the Sobol indices for responder and non-responder profiles.
<table><caption>**Table 2.** Sobol Sensitivity Indices in Responders and Non-Responders.</caption><thead><tr><th>Parameters</th><th>S<sub>1</sub> (Responders)</th><th>ST (Responders)</th><th>S<sub>1</sub> (Non-Responders)</th><th>ST (Non-Responders)</th></tr></thead><tbody><tr><td>N<sub>0</sub></td><td>0.36</td><td>0.53</td><td>0.25</td><td>0.38</td></tr><tr><td>ρ</td><td>0.44</td><td>0.61</td><td>0.12</td><td>0.24</td></tr><tr><td>k<sub>e</sub></td><td>0.11</td><td>0.22</td><td>0.29</td><td>0.45</td></tr><tr><td>λ</td><td>0.07</td><td>0.17</td><td>0.11</td><td>0.21</td></tr><tr><td>a</td><td>0.02</td><td>0.08</td><td>0.06</td><td>0.15</td></tr><tr><td>A<sub>0</sub></td><td>0.01</td><td>0.05</td><td>0.05</td><td>0.10</td></tr></tbody></table>

In responders, $\rho$ and $N_0$ dominate both the first-order Sobol indices (direct individual effect on $T_g$) and the total indices (including parameter interactions). In non-responders, although these indices are less pronounced, $\rho$ and $N_0$ still exert a measurable direct effect and interact with other physiological parameters influencing $T_g$. Overall, for both groups, radiosensitivity and initial tumor burden clearly emerge as the main determinants of $T_g$ dynamics. The elimination rate $k_e$ shows a weaker direct and total influence in responders than in non-responders, preventing any generalization of its effect across both populations. In contrast, the parameters $\lambda$, $a$, and $A_0$ systematically exhibit low sensitivity indices in both cohorts, suggesting a marginal contribution to long-term $T_g$ evolution within the explored parameter space.

These results complement the local sensitivity trends and provide a more rigorous quantitative basis for parameter prioritization during calibration. In clinical applications, such insight supports decisions on which biological variables must be carefully measured and which can be reasonably approximated. While the visual analysis of $T_g$ profiles indicates that $k_e$, $\lambda$, and $a$ influence key trajectory features such as rebound timing and stabilization, variance-based global sensitivity indices confirm that their contribution to cumulative $T_g$ response (e.g., AUC) is limited compared to $\rho$ and $N_0$. This suggests that although these parameters affect the shape of the response curve, their relative impact on overall treatment efficacy is less pronounced.

### 3.4. Minimal Parameter Subset for Model Reduction

The combination of local, global, and statistical sensitivity analyses identified $\rho$ and $N_0$ as the primary drivers of $T_g$ dynamics across both patient groups. These two parameters account for the majority of the model output variance, supporting their role as essential targets for calibration. A reduced model using only $\rho$ and $N_0$ as free parameters successfully reproduced the key dynamic patterns observed in the full model, including early decay and overall response trends. This conclusion is further supported by ANOVA results, where $\rho$ ($p < 10^{-4}$) and $N_0$ ($p = 8.3 \times 10^{-4}$) showed statistically significant influence on the model output. However, parameters such as $k_e$, $\lambda$, and $a$, while contributing less to the total variance, visibly modulate transient features of the $T_g$ trajectory, such as rebound timing and curve steepness. Their ANOVA p-values ($k_e$: $p = 0.12$, $a$: $p = 0.41$, $\lambda$: $p = 0.77$) confirm their limited statistical effect. Nevertheless, their influence may remain clinically relevant in situations where short-term $T_g$ dynamics impact treatment scheduling. The combined interpretation suggests a hierarchy of importance among the parameters: $\rho$ and $N_0$ form the core set required for predictive calibration, while $k_e$, $\lambda$, and $a$ influence secondary, short-term features such as rebound timing. However, given their limited impact on global output variance, these latter parameters may be omitted from individual calibration without compromising the model's predictive utility. This makes the reduced model more amenable for clinical application, particularly when only routine $T_g$ measurements are available.
**4. Conclusions and Discussion**

This study extended the mechanistic model of Barbolosi et al. [11] for metastatic differentiated thyroid cancer by introducing a quantitative, data-driven framework linking pharmacokinetic-pharmacodynamic modeling, parameter estimation, and numerical simulation. The coupled ordinary differential equations describing iodine decay, tumor burden, and thyroglobulin ($T_g$) kinetics were solved analytically when possible and numerically using MATLAB. Parameter estimation was performed with the MCMCSAEM algorithm implemented in Monolix, allowing population (and individual) level calibration from a cohort of 50 patients. The combination of these two computational environments ensured consistency between statistical inference and deterministic system analysis.

Local and global (Sobol) sensitivity analyses, complemented by ANOVA testing, identified tumor radiosensitivity ($\rho$) and initial tumor burden ($N_0$) as the dominant determinants of $T_g$ evolution in both responder and non-responder groups. These two parameters accounted for most of the model variance and define a minimal, clinically actionable set for predictive calibration. A reduced model including only $\rho$ and $N_0$ successfully reproduced the main $T_g$ dynamics of the full system, while parameters $k_e, \lambda$, and $a$ modulated transient features such as rebound timing and curve steepness with limited global influence. The parameters $\rho$ and $N_0$ constitute the core calibration set, while $k_e, \lambda$, and $a$ primarily govern short-term response dynamics. In practice, the model can be reduced from the full parameter space $\theta_{full} = \{T_d, a, \rho, N_0, ke, \lambda, A_0\}$ to a minimal subset $\theta_{reduced} = \{T_d, \rho, N_0\}$ while preserving the main dynamics of the treatment. This reduced formulation defines a tractable dynamic surrogate of the full PK-PD system, retaining the essential therapeutic kinetics.

A priori estimation of these three parameters (ideally by AI-based inference from multimodal patient data) would enable real-time simulation and personalized optimization of RAI treatment. Although the study is limited by the moderate cohort size and the absence of an independent validation dataset, the integrated Monolix-MATLAB workflow demonstrated numerical stability, computational efficiency, and internal consistency. Posterior predictive checks confirmed that simulated $T_g$ trajectories reproduced clinical observations within expected uncertainty bounds, reinforcing the reliability of the calibration procedure. Additional limitations include the assumption of a constant tumor doubling time during treatment and the absence of molecular predictors such as BRAF status in the current model formulation. Despite these constraints, the model captures the essential clinical dynamics of $T_g$ evolution and provides a robust foundation for future extensions toward multi-dose and molecularly informed simulations. Clinically, these findings support the feasibility of patient-specific parameter estimation as a basis for adaptive RAI protocol design, enabling adjustment of administered activity and timing while maintaining compliance with cumulative dose constraints.

Future developments will extend the model to multi-dose regimens, include uncertainty propagation, and integrate the reduced formulation into a clinician-oriented decision-support interface capable of simulating individualized $T_g$ trajectories in real time. Such implementation could ultimately strengthen personalized treatment planning and help achieve an optimal efficacy-toxicity balance in radioactive iodine therapy.

**Author Contributions:** Conceptualization, D.T. and C.V.; methodology, M.F.G.; formal analysis, M.F.G., D.B. and C.V.; investigation, M.F.G., D.B., D.T. and C.V.; writing—original draft preparation, M.F.G.; writing—review and editing, D.B. and C.V.; project administration, M.F.G., D.B. and C.V.; All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.
**Institutional Review Board Statement:** Ethical review and approval were waived for this study due to the retrospective design and full anonymization of patient data, in accordance with institutional and national regulations. The study was conducted in accordance with the principles of the Declaration of Helsinki. The waiver was confirmed by the Head of Nuclear Medicine at Assistance Publique—Hôpitaux de Marseille (AP-HM), David Taieb, co-author of this paper.

**Informed Consent Statement:** Patient consent was waived due to the retrospective nature of the study and the use of fully anonymized data, in accordance with institutional and national regulations.

**Data Availability Statement:** The data that support the findings of this study are available from the corresponding author upon reasonable request.

**Conflicts of Interest:** The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results. David Taieb is the Chief Medical Officer and is a stockholder in Silon Therapeutics.

# Appendix A

## Appendix A.1 *Description of the Study Population*

This section outlines the key characteristics essential for the development of a customized Radioactive Iodine (RAI) dosing model. To achieve robust model validation, it is essential to employ high-quality data and clearly define the inclusion criteria. The dataset used in this study is characterized by the following parameters:

* Sample Size: *n* = 50;

* Patient Age (years): Mean = 29;

* Disease Stage Distribution: $pT1/pT2 = 7$; $pT3/pT4 = 43$; $N1 = 47$; $Nx = 3$; $M1 = 50$.

* The study's retrospective inclusion criteria were as follows:

* Patients who had undergone a total thyroidectomy;

* Patients who had received lymph node dissection (with the exception of three cases);

* Patients who received Radioactive Iodine (RAI) therapy.

Although three patients in the cohort did not undergo lymph node dissection, they were retained in the analysis since the modeling was based on $T_g$ kinetics and not on surgical pathology. Their inclusion did not influence the parameter estimation procedure, which focused on systemic biomarkers and treatment dynamics. The RAI administration protocol parameters are specified below, (median [range]):

* Fixed Administered Activity: 3.7 GBq [3.7–5.5 GBq];

* Average Cumulative Activity: 22.5 GBq [18.5–25.0 GBq];

The clinical data used here were previously published in Barbolosi et al. [11]. In the present study, these data were reanalyzed using a new parameter estimation framework and applied to a novel sensitivity analysis. The modeling and statistical procedures are original to this work.

## Appendix A.2 *Solution of the Systems of Differential Equations*

The solution of the differential equations presented in Equations (1a)–(1c) is inherently tractable, owing to their sequential construction and first-order nature. The solution involves a stepwise approach, commencing with the equation governing the activity, as expressed in Equation (A1), followed by the equation describing cell number dynamics, Equation (A1b), and concluding with the equation for the dosage of $T_g$, Equation (A1c). The parameter *a* in the exponential decay term in Equation (A1a) represents the clearance rate of iodine activity from the body, consistent with iodine-131 pharmacokinetics. It does not
describe a biological delay in effectiveness, but rather the physical decay and elimination of the radioactive agent.

$$A(t) = A_0 e^{a \log(2)t} \qquad (A1a)$$

$$N(t) = N_0 e^{\frac{-A_0 \rho}{\log(2)} t} e^{\frac{\log(2)t}{T_d} + \frac{A_0 \rho e^{-\log(2)t}}{\log(2)}} \qquad (A1b)$$

$$T_g(t) = e^{-ke.t}(T_{g0} + N_0 \lambda \int_0^t e^{ke.x} \frac{\log(2)x}{T_d} \frac{A_0(1+e^{-\log(2)x})}{\log(2)} dx) \qquad (A1c)$$

Upon applying a linear approximation, the activity is described by Equation (A2a), the cell activity by Equation (A2b), and the concentration of $T_g$ by an analytical expression derived from these equations. Owing to its complexity, the detailed mathematical formulation is omitted from this manuscript. Although the system of differential equations can be solved numerically, a linear approximation was introduced to enable closed-form solutions for the key variables. This simplification allows for faster simulation, easier comparison across parameter sets, and clearer interpretation of sensitivity profiles. The approximation is mathematically valid under the assumptions of moderate variation in parameter values and was deemed sufficient for the scope of this analysis.

$$A(t) = A_0(1 - \log(2).at) \qquad (A2a)$$

$$N(t) = N_0 e^{\frac{aA_0 \rho \log(2)t^2}{2} + \left(\frac{\log(2)}{T_d} - A_0 \rho\right)t} \qquad (A2b)$$

## Appendix A.3 Model Contribution—Further Information

Within the sample population, 80% of patients were classified as treatment responders, by clinicians, using criteria derived from clinical observations. These include a sustained reduction in stimulated $T_g$ concentrations, diminished iodine-131 uptake or avid foci (assessed by post-therapeutic imaging), absence of new lesions, and no disease progression on radiological evaluation. While not based on formal RECIST guidelines, these criteria reflect functional and biochemical endpoints commonly used in the clinical management of RAI-treated metastatic thyroid cancer.

The model supports early prediction of responder or non-responder status after 2 to 3 RAI cycles and three $T_g$ measurements. Furthermore, the analysis identified a key discriminating parameter, the tumor doubling time ($T_d$), estimated from longitudinal $T_g$ data. In the 50-patient cohort analyzed, application of a mixture model algorithm via Monolix enabled categorization into two clearly separated subgroups:

* The responder group (comprising 72.5% of patients) with an average $T_d$ of 66.6 months, exhibiting a kinetic profile of $T_g$ concentrations that progressively decreases over time;

* The non-responder (refractory) group (representing 27.5% of patients) with an average $T_d$ of 9.8 months, characterized by a kinetic profile of $T_g$ concentrations that initially decreases, followed by an increase after a certain duration, indicating disease recurrence.

This model serves to define the “Responder” or “Non-Responder” status of patients. After three radiation therapy sessions (instead of the conventional 7 or 8), it provides insights into the treatment’s efficacy. For non-responders, alternative treatment protocols can be proposed, supporting timely clinical decision-making. For responders, the RAI sessions can be spaced out to optimize treatment, balancing efficacy with toxicity (i.e., minimizing side effects).
References

1. Zulewski, H.; Schneiter, R.; Weber, M.; Meier, C. Diagnostic et traitement du cancer différencié de la thyroïde. *Forum. Med. Suisse* **2014**, *14*, 171–176.

2. Schlumberger, M.; Challeton, C.; De Vathaire, F.; Parmentier, C. Treatment of distant metastases of differentiated thyroid carcinoma. *J. Endocrinol. Investig.* **1995**, *18*, 170–172. [CrossRef]

3. Haugen, B.R.; Alexander, E.K.; Bible, K.C.; Doherty, G.M.; Mandel, S.J.; Nikiforov, Y.E.; Pacini, F.; Randolph, G.W.; Sawka, A.M.; Schlumberger, M.; et al. 2015 American Thyroid Association management guidelines for adult patients with thyroid nodules and differentiated thyroid cancer. *Thyroid* **2016**, *26*, 1–133. [CrossRef] [PubMed]

4. Pacini, F.; Fuhrer, D.; Elisei, R.; Handkiewicz-Junak, D.; Leboulleux, S.; Luster, M.; Schlumberger, M.; Smit, J.W. 2022 ETA consensus statement: What are the indications for post-surgical radioiodine therapy in differentiated thyroid cancer? *Eur. Thyroid J.* **2022**, *11*, e210046. [CrossRef] [PubMed]

5. Ibrahimpasic, T.; Nixon, I.J.; Palmer, F.L.; Whitcher, M.M.; Tuttle, R.M.; Shaha, A.; Patel, S.G.; Shah, J.P.; Ganly, I. Undetectable thyroglobulin after total thyroidectomy in patients with low- and intermediate-risk papillary thyroid cancer—Is there a need for radioactive iodine therapy? *Surgery* **2012**, *152*, 1096–1105. [CrossRef]

6. Le Roux, M.K.; Graillon, N.; Guyot, L.; Taieb, D.; Galli, P.; Godio-Raboutet, Y.; Chossegros, C.; Foletti, J. Salivary side effects after radioiodine treatment for differentiated papillary thyroid carcinoma: Long-term study. *Head Neck* **2020**, *42*, 3133–3140. [CrossRef]

7. Bhayani, M.K.; Acharya, V.; Kongkiatkamon, S.; Farah, S.; Roberts, D.B.; Sterba, J.; Chambers, M.S.; Lai, S.Y. Sialendoscopy for patients with radioiodine-induced sialadenitis and xerostomia. *Thyroid* **2015**, *25*, 834–838. [CrossRef] [PubMed]

8. Iyer, N.G.; Morris, L.G.T.; Tuttle, R.M.; Shaha, A.R.; Ganly, I. Rising incidence of second cancers in patients with low-risk (T1N0) thyroid cancer who receive radioactive iodine therapy. *Cancer* **2011**, *117*, 4439–4446. [CrossRef]

9. Lee, S.L. Complications of radioactive iodine treatment of thyroid carcinoma. *J. Natl. Compr. Cancer Netw.* **2010**, *8*, 1277–1287. [CrossRef]

10. Bennett, S. Radioiodine Treatment of Well-Differentiated Thyroid Cancer: Balancing Risks and Benefits. *J. Clin. Oncol.* **2018**, *36*, 1785–1787. [CrossRef]

11. Barbolosi, D.; Summer, I.; Meille, C.; Serre, R.; Kelly, A.; Zerdoud, S.; Bournaud, C.; Schvartz, C.; Toubeau, M.; Toubert, M.-E.; et al. Modeling therapeutic response to radioiodine in metastatic thyroid cancer: A proof-of-concept study for individualized medicine. *Oncotarget* **2017**, *8*, 39167–39176. [CrossRef] [PubMed]

12. Wang, C.; Zhang, X.; Li, H.; Li, X.; Lin, Y. Quantitative thyroglobulin response to radioactive iodine treatment in predicting radioactive iodine-refractory thyroid cancer with pulmonary metastasis. *PLoS ONE* **2017**, *12*, e0179664. [CrossRef]

13. Song, H.-J.; Qiu, Z.-L.; Shen, C.-T.; Wei, W.-J.; Luo, Q.-Y. Pulmonary metastases in differentiated thyroid cancer: Efficacy of radioiodine therapy and prognostic factors. *Eur. J. Endocrinol.* **2015**, *173*, 399–408. [CrossRef]

14. Strobl, M.R.; Gallaher, J.; Robertson-Tessi, M.; West, J.; Anderson, A.R.A. Treatment of evolving cancers will require dynamic decision support. *Ann Oncol.* **2023**, *34*, 867–884. [CrossRef]

15. Barbolosi, D.; Ciccolini, J.; Lacarelle, B.; Barlési, F.; André, N. Computational oncology—mathematical modelling of drug regimens for precision medicine. *Nat. Rev. Clin. Oncol.* **2016**, *13*, 242–254. [CrossRef]

16. Cabanillas, M.E.; McFadden, D.G.; Durante, C. Thyroid cancer. *Lancet* **2016**, *388*, 2783–2795. [CrossRef] [PubMed]

17. Sawka, A.M.; Thephamongkhol, K.; Brouwers, M.; Thabane, L.; Browman, G.; Gerstein, H.C. A systematic review and meta-analysis of the effectiveness of radioactive iodine remnant ablation for well-differentiated thyroid cancer. *J. Clin. Endocrinol. Metab.* **2004**, *89*, 3668–3676. [CrossRef] [PubMed]

18. Wang, R.; Zhang, Y.; Tan, J.; Zhang, G.; Zhang, R.; Zheng, W.; He, Y. Analysis of radioiodine therapy and prognostic factors of differentiated thyroid cancer patients with pulmonary metastasis: An 8-year retrospective study. *Medicine* **2017**, *96*, e6809. [CrossRef] [PubMed]

19. Giovanella, L.; Ceriani, L.; Suriano, S.; Maffioli, M.; Spriano, G. The role of thyroglobulin doubling time in differentiated thyroid cancer: A meta-analysis. *Endocr. Connect.* **2022**, *11*, e210648. [CrossRef]

20. Haymart, M.R.; Banerjee, M.; Stewart, A.K.; Koenig, R.J.; Birkmeyer, J.D.; Griggs, J.J. Use of radioactive iodine for thyroid cancer. *JAMA* **2011**, *306*, 721–728. [CrossRef]

21. Robb, R.; Yang, L.; Shen, C.; Wolfe, A.R.; Webb, A.; Zhang, X.; Vedaie, M.; Saji, M.; Jhiang, S.; Ringel, M.D.; et al. Inhibiting BRAF Oncogene-Mediated Radioresistance Effectively Radiosensitizes BRAFV600E-Mutant Thyroid Cancer Cells by Constraining DNA Double-Strand Break Repair. *Clin. Cancer Res.* **2019**, *25*, 4749–4760. [CrossRef] [PubMed]

22. Samson, A.; Lavielle, M.; Mentré, F. The SAEM algorithm for group comparison tests in longitudinal data analysis based on non-linear mixed-effects model. *Stat. Med.* **2007**, *26*, 4860–4875. [CrossRef] [PubMed]

23. Martin, N.K.; Gaffney, E.A.; Gatenby, R.A.; Maini, P.K. Tumour-stromal interactions in acid-mediated invasion: A mathematical model. *J. Theor. Biol.* **2010**, *267*, 461–470. [CrossRef]

24. Barrow, J.C.; Lindsley, C.W. The importance of PK–PD. *J. Med. Chem.* **2023**, *66*, 4273–4274. [CrossRef]
25. Sun, Y.-Q.; Sun, D.; Zhang, X.; Zhang, Y.-Q.; Lin, Y.-S. Radioiodine adjuvant therapy in differentiated thyroid cancer: An update and reconsideration. *Front. Endocrinol.* **2022**, *13*, 994288. [CrossRef]

26. Delattre, M.; Lavielle, M. Maximum likelihood estimation in discrete mixed hidden Markov models using the SAEM algorithm. *Comput. Stat. Data Anal.* **2012**, *56*, 2073–2085. [CrossRef]

27. Lavielle, M.; Mentré, F. Estimation of population pharmacokinetic parameters of saquinavir in HIV patients with the MONOLIX software. *J. Pharmacokinet. Pharmacodyn.* **2007**, *34*, 229–249. [CrossRef]

28. Kwon, S.Y.; Kim, M.H.; Kong, E.; Chong, A.; Yoo, S.W.; Jeon, S.; Park, S.; Kim, D.; Kang, S.H.; Choi, J.E.; et al. Clinical impact of radioactive iodine dose selection based on the number of metastatic lymph nodes in patients with papillary thyroid carcinoma: A multicenter retrospective cohort study. *Clin. Endocrinol.* **2021**, *95*, 901–908. [CrossRef]

29. Piña-Sánchez, P.; Chávez-González, A.; Ruiz-Tachiquín, M.; Vadillo, E.; Monroy-García, A.; Montesinos, J.J.; Grajales, R.; de la Barrera, M.G.; Mayani, H. Cancer biology, epidemiology, and treatment in the 21st century: Current status and future challenges from a biomedical perspective. *Cancer Control* **2021**, *28*, 10732748211038735. [CrossRef] [PubMed]

30. Kovari, M. Effect of delay time on effectiveness of stable iodine prophylaxis after intake of radioiodine. *J. Radiol. Prot.* **1994**, *14*, 131. [CrossRef]

31. Oh, H.S.; Kwon, H.; Song, E.; Jeon, M.J.; Kim, T.Y.; Lee, J.H.; Kim, W.B.; Shong, Y.K.; Chung, K.-W.; Baek, J.H.; et al. Tumor volume doubling time in active surveillance of papillary thyroid carcinoma. *Thyroid* **2019**, *29*, 642–649. [CrossRef] [PubMed]

32. Vrndic, O.B.; Savin, S.B.; Mijatovic, L.C.; Djukic, A.A.; Jeftic, I.D.; Zivancevic Simonovic, S.T. Concentration of thyroglobulin and thyroglobulin-specific autoantibodies in patients with differentiated thyroid cancer after treatment with radioactive iodine. *Lab. Med.* **2011**, *42*, 27–31. [CrossRef]

33. Hocevar, M.; Auersperg, M.; Stanovnik, L. The dynamics of serum thyroglobulin elimination from the body after thyroid surgery. *Eur. J. Surg. Oncol.* **1997**, *23*, 208–210. [CrossRef]

34. Hamis, S.; Stratiev, S.; Powathil, G.G. Uncertainty and Sensitivity Analyses Methods for Agent-Based Mathematical Models: An Introductory Review. *Phys. Cancer Res. Adv.* **2020**, *1*, 1–37.

35. Mu, Z.Z.; Zhang, X.; Lin, Y.-S. Identification of radioactive iodine refractory differentiated thyroid cancer. *Chonnam Med J.* **2019**, *55*, 127. [CrossRef]

36. Saltelli, A.; Tarantola, S.; Campolongo, F.; Ratto, M. *Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models*; John Wiley & Sons Ltd.: Chichester, UK, 2008.

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.