A Mathematical Model of Thyroid Disease Response to Radiotherapy

We present a mechanistic biomathematical model of molecular radiotherapy of thyroid disease. The general model consists of a set of differential equations describing the dynamics of different populations of thyroid cells with varying degrees of damage caused by radiotherapy (undamaged cells, sub-lethally damaged cells, doomed cells, and dead cells), as well as the dynamics of thyroglobulin and antithyroglobulin autoantibodies, which are important surrogates of treatment response. The model is presented in two flavours: on the one hand, as a deterministic continuous model, which is useful to fit populational data, and on the other hand, as a stochastic Markov model, which is particularly useful to investigate tumor control probabilities and treatment individualization. The model was used to fit the response dynamics (tumor/thyroid volumes, thyroglobulin and antithyroglobulin autoantibodies) observed in experimental studies of thyroid cancer and Graves' disease treated with I-131-radiotherapy. A qualitative adequate fitting of the model to the experimental data was achieved. We also used the model to investigate treatment individualization strategies for differentiated thyroid cancer, aiming to improve the tumor control probability. We found that simple individualization strategies based on the absorbed dose in the tumor and tumor radiosensitivity (which are both magnitudes that can potentially be individually determined for every patient) can lead to an important raise of tumor control probabilities.


Introduction
Radioactive iodine, 131 I, therapy (RAI) is a type of molecular radiotherapy (also called targeted radiotherapy) that has been commonly used for the treatment of differentiated thyroid cancer since the 1940s.This therapy, based on the ability of well-differentiated (papillary or follicular) thyroid cancer cells to absorb iodine, uses radiation either to ablate remnant thyroid tissue after surgery or to treat residual, recurrent, or metastatic cancer, such as pulmonary and bone metastases.Radioiodine therapy is suitable for patients with I-avid but surgically unresectable metastases of differentiated thyroid cancer, DTC.In patients with non-avid metastases, due to the expansion of less differentiated cells, RAI shows no benefits, and treatment shifts to other strategies.
One of the most relevant physical quantities affecting the response to therapy is the absorbed radiation dose to lesions [1,2,3].Large differences in the efficacy and cure rates of RAI on metastatic DTC patients are reported in many studies [4].This might be partially due to the large dosimetric uncertainties of RAI.Other factors, such as age, thyroglobulin levels at diagnosis, and stage of the disease have been identified as important prognostic factors [5,6].The systemic nature of this therapy, with a radionuclide that takes several days to disintegrate, makes RAI treatment planning a complex procedure.Due to these complexities, RAI treatment administration is usually guided by non-personalized empirical criteria that establish the activity of 131 I to be administered, as well as cycle schedules.This scenario poses a concern of potentially over/under-dosing patients, and there is a need to develop alternative strategies and tools for the development of individualized treatment planning.Some markers are usually measured in the management of thyroid cancer treatment, as surrogates of response.The most extended tumor marker is thyroglobulin, Tg, which is produced by normal follicular cells, and accumulated as a colloid in an extracellular compartment of the thyroid follicles.Under certain circumstances, such as rapid and disordered thyroid tissue growth, inflammation leading to follicular destruction, or hemorrhage, Tg is also secreted into the systemic circulation.The concentration of Tg in blood is thus proportional to the volume of thyroid tissue, and elevated serum Tg levels can be seen in benign conditions like Graves' disease and thyroiditis, and also in thyroid cancer.Regarding the latter, high serum Tg concentrations are observed in approximately two-thirds of well-differentiated thyroid cancers, including follicular and papillary thyroid carcinomas, up to 50% of anaplastic thyroid carcinomas, and some medullary carcinomas [7].
After partial or total thyroidectomy, serum Tg levels post-surgery reflect the amount of residual thyroid mass, with high levels one month after thyroidectomy indicating incomplete excision of the gland, or recurrent or metastatic disease [8,9].Furthermore, the risk of having recurrent or persistent disease increases as the postoperative Tg rises [10].This is used in the initial follow-up of patients with papillary and follicular thyroid carcinoma, with measurements every 6-12 months, as a guide to initiate RAI for remnant ablation, and as an important surrogate marker of the response of metastatic patients to RAI [11].However, the subset of DTC patients (12%) with low blood Tg levels before thyroidectomy will not show rising serum Tg levels when there is a recurrence of the disease, therefore the importance of knowing the patient serum Tg level before surgery [12].
The measurement of serum Tg is influenced by the coexistence with antithyroglobulin autoantibodies, TgAb, in blood.These antibodies are present in cases of Hashimoto's disease, Graves' disease, 7.5-25% of DTC patients and 10% of the healthy population [13].Their presence impairs the reliability of the most sensitive Tg detection methods currently available (immunometric assays, IMA), leading to levels lower than expected or undetectable [14,15].Therefore, current thyroid cancer management guidelines indicate that TgAb should always be measured in parallel with Tg, to identify potentially unreliable Tg tests [11].The analytical interference is usually attributed to the ability of serum TgAb to block the access to the Tg epitopes, to which the IMA test antibodies bind [7,16].In addition, it has been recently reported that the interference is more significant in vivo than in vitro, suggesting the induction of the enhanced clearance of Tg by TgAb [17,18,19].Since the only source of Tg is thyroid tissue, the persistence of elevated TgAb levels is an indicator of the continued presence of Tg, and, for example, possible metastasis.As TgAb levels change in accordance with Tg concentration, follow-up measurements of these antibodies can be used as a surrogate tumor marker, helping to predict cancer recurrence in patients with disseminated thyroid carcinoma [20,21].
There is an ongoing paradigm shift towards individualized treatment planning in molecular radiotherapy, including RAI.To achieve this goal, there is a need to develop individualized dosimetry methods for molecular radiotherapy [22], but also reliable biomathematical dose-response models that can predict the effect of the treatment [23,24].Mathematical models can be very useful to assist in the analysis and interpretation of experimental results.Once validated, they can be used to investigate "What if?" questions, simulating different treatment strategies to optimize RAI through computational simulation.The use of computer modeling is usually referred to as in silico modeling, as an allusion to the in vivo or in vitro models typically used in biology/biomedicine.Such models might also be useful to guide the management of patients during the course of the treatment, especially if they include the markers usually measured in the clinical practice as surrogates of response.In the case of the thyroid, models of response to cancer (or other diseases) therapy, including RAI, should ideally include the evolution of Tg and TgAb.However, this has not yet been addressed in the literature.
Modeling of thyroid cancer response to therapy was previously studied by Barbolosi et al. [23], who presented a phenomenological model of tumor response to RAI, including dynamics of Tg.Traino et al. [24] provided a differentiated thyroid cancer mass reduction model, based on personalized RAI dosimetry calculations, for post-surgical thyroid remnant or recurrent or metastatic cancer.That study did not consider the evolution of Tg or TgAb.Pandiyan et al. [25] presented a mathematical model of response to methimazole therapy for Graves' disease, which includes many of the compartments discussed above, and models the evolution of TSH, T4 and T3 hormones, as well as other antibodies relevant for Graves' hyperthyroidism, but not Tg nor TgAb.Other mathematical models of Graves' disease relapse after antithyroid drugs [26] or RAI [27] have been also presented.Of special interest is also the article by Degon et al. [28], who presented a mathematical model of Tg production in the thyroid, although the response to RAI was not investigated.Table 1 summarizes these models, presenting their main characteristics and limitations.
In this work, we present a radiobiological model of thyroid response to RAI.We build on the models presented above, while aiming to address some of their limitations.The approach that we Author Summary Limitation Barbolosi et al. [23] Model of metastatic thyroid cancer response to RAI, including dynamics of Tg.Phenomenological model, which may limit the potential for individualization.TgAb is not considered.Traino et al. [24] DTC mass reduction model in response to RAI.Based on personalized dosimetry calculations.
Extension of the phenomenological Linear-Quadratic formalism to model mass reduction.Neither Tg nor TgAb are considered.Degon et al. [28] Normal functioning thyroid model considering TSH, Tg, T3, T4 and TPO.
Neither the response to RAI nor Tg are considered.follow in this work is more mechanistic than those presented in previous studies [23,24,27], and the model specifically includes the dynamics of Tg and TgAb.It consists of a set of differential equations describing the dynamics of different populations of thyroid/tumor cells, and Tg and TgAb concentrations.
The model was fit to published data of patients with metastatic differentiated thyroid cancer treated with several cycles of RAI, and patients with Graves' disease.Direct intercomparison with previous thyroid response models was, however, not possible, as none of them considered together all the variables included in our work: they either missed the dynamics of Tg or TgAb, or did not model the response to RAI.
In addition, we present a stochastic Markov-like version of the model that can be used to compute tumor control probabilities.We used it to investigate RAI treatment individualization strategies aiming to improve the tumor control rates achieved with this therapy.This latter study highlights the potential importance of such models for the individualization of RAI.

The Model
The response of the thyroid cells to RAI is modeled by extending a previous work of the group that describes the response of a population of cells to a continuous radiation dose rate [29].A system of coupled differential equations is used to model the dynamics of different tumor cell populations: undamaged cells, N (t), sub-lethally damaged cells, N s (t), doomed cells, N d (t), and dead cells, N x (t).
The radiation treatment is characterized by a dose rate, r(t), having a time variation that depends on the physical decay of the radioisotope and its biological clearance.Following the rationale behind the classical linear-quadratic model of radiation response [30], radiation damage can be either lethal (∼ ar(t)) or sub-lethal (∼ br(t)).Lethally damaged (doomed) cells cannot recover, and will eventually die, moving to the dead compartment with a rate γ.Sub-lethally damaged cells may recover (we assume a simple exponential repair kinetic term, with characteristic repair rate µ), or may become doomed due to further lethal or sub-lethal damage.Dead cells disappear due to resorption with a resorption rate η.
Cell proliferation is supposed to follow a logistic function, with an exponential rate, λ, and exponential proliferation moderated according to the total number of cells and a saturation limit, N sat .Experimental evidence shows exponential proliferation slowing down with increasing tumor volumes, and therefore growth models must include some sort of saturation mechanism [31].In general, we assume that N d and N s cells may still carry proliferative capacity, which can be different from that of undamaged tumor cells.
The system of equations has the following form: where is the total number of cells at time t.Three compartments are included to model the concentration of serum thyroglobulin and antithyroglobulin antibodies in the patient.The free thyroglobulin compartment, T f , stands for the concentration of thyroglobulin not bound to antithyroglobulin antibodies, A b , while the T b compartment accounts for the concentration of Tg bound to these antibodies: Patients treated with 131 I for thyroid carcinoma distant metastases or Graves' disease may achieve complete remission, including undetectable serum Tg levels, months or even years after treatment [4,32,33].This suggests that, after being damaged by radiation, thyroid cells can still produce detectable amounts of Tg for long periods of time, despite having a limited proliferative capacity [34,35].Based on this, we assume that free thyroglobulin is produced by the cells that keep proliferative capacity (healthy, sub-lethally damaged, and doomed cells).The Tg production rate is modulated by the total thyroglobulin level (the parameter c is a modulation constant).This intends to mimic the Hypothalamus-Pituitary-Thyroid feedback mechanisms that regulate the thyroid hormones pools in the body, as previously considered by other thyroid mathematical models [25,28].
Moreover, Tg can be released from the colloid reservoir into the bloodstream when thyroid cells are damaged by radiation [18].To reflect this, a T f production term is included, proportional to N d , with a characteristic release rate z t .Free Tg is eliminated from blood with a clearance rate k e , and can bind to Tg antibodies (if these are present) with a binding rate k a , moving to a compartment that describes the kinetics of the antibody-bound thyroglobulin serum, T b (t).This binding also removes TgAb, for which we use the same term with a binding rate k ′ a .Bound Tg is assumed to show an enhanced metabolic clearance, with rate k e +k eb , induced by TgAb [17,18,19].
TgAb are intra-follicular antibodies that can bind to immune cells and antigens, not directly linked with thyroid cell destruction per se [20].Thyroid disease is associated with Tg structural changes and leakage, leading to TgAb production [19,20].TgAb are mainly produced by lymphocytes infiltrating the gland and, to a lesser extent, by cervical lymph nodes and bone marrow immune cells [36].Considering this, the TgAb compartment in our model includes two terms of antibody production.The first is proportional to the total serum Tg concentration (through a production rate constant z), therefore being related to the cancer thyroid tissue that produces Tg.The second term is related to lymphocyte production, and, therefore, linked to the population of doomed cells that would release antigens, with an associated production constant z b .A natural elimination term is included to model the disappearance of thyroid antibodies observed in patients after thyroidectomy or thyroid tissue destruction by RAI.The previously mentioned binding of Tg to TgAb also results in an effective elimination for these antibodies, as bound TgAb cannot be detected by the immunometric TgAb assays, in which Tg is the agent used to bind the antibodies in the sample [37].
Figure 1: Sketch of our stochastic model with transfer probabilities between compartments and flux arrows.P i,j is the probability of a cell in compartment i to go to compartment j (i.e., N i → N i − 1, N j → N j + 1).P i,i is the probability of duplication of a cell in compartment i (i.e., N i → N i + 1).(n, d, s, x ) refer to compartments of non-damaged, doomed, sub-lethally damaged, and dead cells, and "out" is used for cells leaving the system.

Modeling Tumor Control Probability-Stochasticity
In radiobiological modeling, one of the main aims is to model the probability that a given treatment may control the tumor in a population of patients, the so-called Tumor Control Probability (TCP).This is usually mechanistically modeled by assuming the clonogenic cell hypothesis: a tumor will be controlled if all cells with clonogenic capacity in it are killed [38].In our model, we associate clonogenic potential to undamaged and sub-lethally damaged cells, N and N s (doomed N d cells may still have some proliferative capacity, but they will eventually die and cannot sustain a tumor long-term).Therefore, a given treatment will lead to tumor control if N (t) + N s (t) = 0 at some time t.
Notice that our model, as presented in the previous subsection, is continuous and deterministic.In order to simulate tumor control probability, we need a discrete (numbers of cells) stochastic model.To do so, and to be able to model TCP, we have identified terms in the equations as transfer probabilities between compartments, and used them to obtain a stochastic Markov model.The model is sketched in Figure 1.Transfer probabilities for each compartment/process are obtained by re-writing the terms in the differential equations as , where i stands for each cell compartment, and each term f j is a transition probability.For each time step, given by ∆t, the discrete probabilities are given by: P i,j is the probability of a cell in compartment i to go to compartment j (i.e., N i → N i − 1, N j → N j + 1).P i,i is the probability of duplication of a cell in compartment i (i.e., N i → N i + 1).n, d, s, and x refer to compartments of non-damaged, doomed, sub-lethally damaged, and dead cells, respectively, and "out" is used for cells leaving the system.
In order to speed up computation time (stochastic modeling is more time-consuming), continuous and stochastic modeling have been combined.The stochastic modeling is activated only when numbers of cells fall below a given threshold (10,000 cells).The combination of continuous and discrete modeling seems adequate, as it is much faster than a purely discrete model, and differences between both approaches only appear for low numbers of cells.

Fit to Experimental Data: Graves' Disease and Tumor Response
There is limited literature providing adequate data to evaluate the capability of our model to fit response to therapy, including the dynamics of cells (mass), Tg, and TgAb.Suitable datasets should include the evolution of thyroid mass/volume, the dynamics of Tg and TgAb after RAI, and the dose to the thyroid/metastases.In a dedicated search, we only found two sets of experimental clinical data, including these variables [18,39].
The study of Latrofa et al. [18] monitored the response to a single administration of 555 MBq RAI in a cohort of 30 patients with Graves' hyperthyroidism.Thyroid tissue volume and serum Tg and TgAb levels were registered at the time of treatment and every 15 days thereafter, up to 3 months.We used reference cell density values to estimate the corresponding thyroid cell numbers [40], and average values of absorbed dose per MBq and effective half-life reported by Traino et al. [27].
We also used the model to fit the experimental data of differentiated thyroid cancer response published by de Keizer et al. [39].In this study, the authors reported the response of a cohort of patients who suffered 131 I-avid metastatic or recurrent disease after near-total thyroidectomy followed by ablative RAI.Patients received recombinant human thyrotropin (rhTSH) before radioiodine administration of 7400 MBq.Tumor response was assessed by comparing serum Tg levels pre-treatment and 3 months post-treatment.In order to have information about tumor mass evolution and at least three follow-up marker measurements, we focused on the patients who received at least two RAI administrations.This study provided personalized dosimetry information, based on 131 I scintigrams, that allows us to define the dose rate profiles in the lesions, r(t), considering variations in 131 I uptake and biodistribution (effective half-life) during the course of the treatment.These are relevant parameters that importantly affect the outcome, but are rarely available from the literature.Uncertainties of tumor mass, Tg and TgAb levels are not reported in that paper, hence we have associated ad hoc uncertainties of 10% (30%) of the first measurement point to tumor masses and Tg levels (TgAb levels).

Evaluation of Tumor Control Probabilities-Modeling Study
We have used the stochastic flavour of the model to investigate TCP in silico and to formulate hypotheses of optimal treatment strategies.In order to generate a population of patients we have taken a set of model parameter values (obtained from fits to experimental data), including uptake and half-life of activities in the tumor, and randomly perturbed them (with normal perturbations with a relative standard deviation of 25%).With this procedure, we have generated 1000 different sets of parameter values, naively representing a cohort of 1000 patients with individualized responses.
The TCP for the population (or subgroup of the population) is computed as: where controls refers to the number of simulations (patients) for which control was achieved (as detailed in Section 2.2), and size refers to the number of simulations (patients), either in the whole cohort (1000), or in subgroups of that cohort.
We have investigated whether treatment individualization can lead to important gains of TCP, i.e., rather than injecting the same activity A 0 to each patient, the injected activity is adjusted according to patient characteristics/model parameters (using parameter values that can be easily measured in the clinic).We have used a single fraction treatment, and we have investigated three hypotheses and individualization strategies, graphically sketched in Figure 2: • Individualization strategy 1 -Patients with low absorbed dose in the tumor (either due to low activity uptake and/or fast activity clearance) will have poor tumor control, while patients with high doses in the tumor will have good control.Patients can be stratified according to the absorbed tumor doses associated with the non-individualized standard treatment and the injected activities adjusted for each group.For high-dose patients, if the TCP is already high, a reduction in activity (and therefore dose) might not have a significant effect on tumor control, but could potentially decrease the likelihood of organ/tissue toxicity.On the other hand, for patients receiving low tumor doses, injected activities could be increased to enhance the tumor doses and improve the likelihood of control.According to this reasoning, patients have been stratified into five groups according to absorbed tumor doses, and the injected activities adjusted for each group.In a clinical scenario, one should consider how such changes in injected activities affect the probability of toxicity.Modeling toxicity is beyond the scope of this work, but we have imposed a constraint in the individualization, requiring that the average injected activity in the population does not change (a sort of populational isotoxicity constraint).
• Individualization strategy 2-The radiosensitivity of tumor cells will also condition tumor response and tumor control.While radiosensitivity cannot be as easily measured as the radiation dose, there are genetic profiling techniques that can infer the radiosensitivity of cells: this is behind the Genomically Adjusted Radiation Dose (GARD) methodology, aiming at individualizing radiation dose according to radiosensitivity profiling [41,42].We have investigated a similar approach.Patients have been stratified into five groups, as in Strategy 1, according to their radiosensitivity (the parameter b of sublethal damage, which will serve as a simple surrogate of radiosensitivity, even if radiosensitivity will depend on more parameters).Those in the groups with lower b values (radioresistant) will have the injected activity scaled up, and those in groups with higher b values (radiosensitive) will receive scaled-down activities, as shown in Figure 2.
• Individualization strategy 3-We will also investigate a strategy that is a combination of 1 and 2 above.The dose groups created in Individualization strategy 1 have been further split into two subgroups, A and B, each containing 50% of the patients in the group, according to the value of parameter b.Those in the subgroup with lower b values (radioresistant) will have the injected activity scaled up, and those in the subgroup with higher b values (radiosensitive) will receive scaled-down activities.

Model Implementation, Resolution, Identifiability, and Data Fitting
The model was implemented in Matlab (The Mathworks Inc.) and it was numerically solved by employing an explicit Euler method [43] with a time step ∆t = 5 min.The same time step was used for the stochastic Markov model.The integration time depends on the problem under study, but can reach several months (for multifraction treatments of DTC).In this situation, only a fraction of the simulated data were saved (typically data are saved every day) in order to prevent unmanageable datasets.Code and data are available from the Harvard Dataverse [44], including a more detailed presentation.
A simulated annealing method [45] was implemented to fit the model to experimental data.The function to be minimized, F , during the optimization was the weighted sum of quadratic differences between model and experimental data: where A exp and A mod refer to experimental and model results, and w to the weights (uncertainties of the experimental data presented in Section 2.3).Indices i and j run over the number of variables (volumes/masses, Tg, TgAb) and time points included in the fit, respectively.Identifiability is an important property when fitting a model to experimental data and refers to whether it is possible to determine the values of the unknown parameters.We have investigated the identifiability of the model here presented by relying on the DAISY package [46].

Identifiability
The model is classified as globally identifiable by DAISY, when considering the following observables (and the rationale for considering them): • r(t): the dose rate could be monitored by using, for example, imaging techniques to characterize the dynamics of activity/dose.Formally, r(t) is an input, rather than an observable.
• N (t) + N s (t) + N d (t) + N x (t): this is related to the dynamics of tumor volumes, which can be again monitored with imaging techniques.• N (t) + N s (t): in our model, N and N s carry clonogenic capacity, and the fraction/ number of clonogenic cells can be obtained by harvesting tumor samples and performing plating experiments.
• T b (t) + T f (t): the dynamics of serum thyroglobulin can be monitored in blood samples.Separating T b from T f may be more difficult, therefore we consider the sum of both variables.
• A b (t): the dynamics of serum thyroglobulin antibodies can be also monitored in blood samples.

Fit to Graves' Disease Data
In Figure 3 we present best-fits to population-averaged response to a single administration of 555 MBq RAI measured in a cohort of 30 patients with Graves' hyperthyroidism.Left panels show fits of the model to volumes (total volumes and fractions corresponding to undamaged, doomed, and dead cells), as well as Tg (total, free, and bound) and TgAb.Best-fitting parameter values are reported in the supplementary materials, Table 2.In general, results obtained with the model match fairly well the trends observed in the measurements.Slow-doomed and dead-cell elimination rates, with half-lives ∼ 1 month, were required to fit the smooth decrease in thyroid volume with time after treatment.Tg and TgAb levels were adequately reproduced, considering the experimental uncertainties.However, the rise observed in TgAb after month 2 (five-fold increase with respect to the minimum value registered) was much smoother in the model (four-fold change).
We have also simulated the effect of interpatient variabilities.In order to simulate a population, best-fitting parameters to the population-averaged data were perturbed with a normal distribution with 25% relative standard deviation (1000 different sets of parameter values to simulate a population of 1000 patients).Confidence intervals reported in Figure 3 (left panels) are similar to the variability of the experimental data.

Fit to Differentiated Thyroid Cancer Response Data
Figures 4 and 5 show best fits of the model to the response of two differentiated thyroid cancer patients treated with RAI after rhTSH stimulation.Best-fitting parameter values are reported in the supplementary materials, Table 2. Patient 1 suffered a local recurrence (one single lesion) and showed detectable TgAb levels (Figure 4).This patient received three 131 I administrations, at months 0, 6, and 12.The tumor mass steadily increased with time after treatment, illustrating the previously reported typical slow response of DTC to RAI.Probably linked to the presence of TgAb, this patient showed rather low Tg levels (panel B), which increased from month 0 to 6, and then decreased to low levels.The model was able to reproduce these trends, although with TgAb changes that are smoother than those experimentally observed (panel C).
Patient 2 presented three metastatic lesions in the lungs, and undetectable TgAb levels.Panels A, B and C in Figure 5 show the evolution of the mass with time for every lesion, which remained stable in the follow-up for lesions 1 and 3, and decreased approximately by 20% at month 6 in lesion 2. Slow-cell elimination rates were again required to fit the evolution of the mass and reproduce the drop in Tg concentration (panel D) associated with the loss of functional, undamaged cells.

In Silico Treatment Individualization
In Tables 3-5 we show the tumor control probabilities calculated with the Markov model, in a cohort of 1000 simulated patients, for individualization strategies 1, 2 and 3, respectively.Statistical uncertainties are not presented for the sake of simplicity, but they can be calculated from the reported TCP values and the number of patients in each group, which, using the binomial distribution, leads to values up to ±5% for 100 patients and TCP∼50%.
Table 3 shows the tumor control probabilities calculated with the Markov model, in a cohort of 1000 simulated patients, for the individualization strategy 1 (a stratification based on the absorbed dose in the tumor under the standard treatment).As expected, the TCP is higher in the groups receiving larger doses: when there is no individualization, the TCP in the lowest dosage group (group 1) is 16% versus 99% in the highest dosage group (group 5, which corresponds to the 20% of the cohort who receive the larger doses).Treatment individualization, with maximum activity perturbation of ±10% and ±20%, leads to an important TCP increase.In groups 1 and 2, who receive the lower doses with the standard activity, the TCP raises from 16% and 68.5% up to 39% and 79.5% (in the scenario with maximum activity perturbation of ±20%).On the other hand, in groups 4 and 5, where the individualization leads to lower injected activities, TCP values decrease slightly, up to 4.5% lower, yet still very close to 100%.When evaluated on the whole cohort, individualization strategy 1 leads to an overall TCP improvement of 3.1% and 5.7% for the 10% and 20% maximum activity boosts, respectively.
Table 4 shows the tumor control probabilities for the individualization strategy 2 (a stratification based on the radiosensitivity of tumors, characterized by the parameter of sub-lethal damage, b).As expected, the TCP is higher in the groups with more radiosensitive tumors: when there is no individualization, the TCP in the least radiosensitive group (group 1') is 48.5% versus 93.5% in most radiosensitive group (group 5').Treatment individualization, with maximum activity perturbation of ±10% and ±20%, leads to a slight TCP increase.In groups 1' and 2', the TCP raises from 48.5% and 71% up to 67.5% and 79.5% (in the scenario with maximum activity perturbation of ±20%).On the other hand, in groups 4' and 5', where the individualization leads to lower injected activities, TCP values decrease by a ≃ 8%.When evaluated on the whole cohort, individualization strategy 2 leads to an overall TCP improvement of 1.7% and 2.4% for the 10% and 20% maximum activity boosts, respectively.These improvements are lower than those obtained with strategy 1, probably because, as mentioned, radiosensitivity has a complex dependence on several parameters, while only parameter b has been considered for this stratification.
Table 5 presents the results of individualization strategy 3, which further splits the 5 patient groups of strategy 1 into two subgroups each, attending to their radiosensitivity.Row 1 corresponds to non-individualized treatment and reflects the effect of dose andradiosensitivity on TCP: subgroups A comprise the more radioresistant patients, who present lower TCP values than patients in subgroups B. This individualization strategy is able to further improve the control probability of the whole population, improving the TCP of the whole cohort by up to 8% (versus 5.7% for strategy 1 ).

Discussion
In this piece of work, we have presented a mechanistic biomathematical model of thyroid response to RAI.The model includes different populations of thyroid cells with varying degrees of radiationinduced damage, and the dynamics of Tg and TgAb, which are important surrogates of treatment response.The model can adequately fit the response to RAI of DTC and Graves' disease observed in experimental studies, as shown in Figures 3-5.
The ultimate goal of in silico biomathematical response models is to assist in the design of optimal therapeutic strategies, including treatment individualization.In this regard, we used a stochastic version of the model to investigate the treatment gain (measured in terms of tumor control probability) that could be achieved by DTC RAI individualization.Administered activities were tailored considering patient-specific dosimetry, as well as radiosensitivity features.In a simulated patient cohort, we have found that individualization with moderate variations of the injected activity may lead to an improvement of the TCP (up to ≃+8% for the whole cohort, and up to ≃+20% for specific groups of patients, see Tables 3-5), exemplifying how the model can be used to study novel therapeutic strategies.The model can be easily used to analyze other treatment parameters such as the time between RAI administrations in multiple fraction treatments.
Given the large number of parameters involved in our model, it is convenient to use experimental data to set fixed values for some of them, or at least to establish ranges within which the parameter values would be expected to lie, based on population-averaged values or clinical/observational studies.This is of paramount importance to ensure that the model can fit experimentally observed dynamics with parameters that are biologically sound and to avoid overfitting, and merits further discussion: • Radiobiological experiments with tumor cell cultures irradiated at different dose rates are useful to determine radiosensitivity parameters and repair rate ranges [47].These data served to fix the repair half-life to 2 hours, within the range of values reported in Steel et al. [48], and to set a threshold to the radiosensitivity parameters a and b equal to 0.45.The parameter values associated with tumor cell elimination are more uncertain: the mechanisms by which radiation induces the different modes of cell death, their relative relevance, and variability between cancer cell types are still a matter of debate [49].Thyroid tumors are mostly of epithelial origin [50], which has been associated with low levels of early apoptosis.Mitotic catastrophe, observed to peak five days after irradiation, may contribute more importantly to cell-killing in these tumors [51,52].Other pathways such as senescence or radiation-induced stimuli of the innate and adaptive immune cells would extend the kinetics of cell death (the time when lethally damaged cells die).Although the relative relevance of these pathways remains uncertain [53], the delayed effect observed in the response of thyroid to RAI [33,35] suggests that they may play a significant role.As a result, doomed and dead cell elimination rates are expected to be low, with effective half-lives of the order of several weeks/months.The experimental data used in this work for model fitting confirmed this, as masses/volumes decreased slowly, requiring cell elimination half-lives that ranged from 15 days to values as high as 2 years, for the dead cells of the DTC local recurrence.
• Regarding thyroid tumor markers, the information available to estimate value ranges for many of the parameters associated with the dynamics of these compounds is very scarce, or non-existent in some cases.The natural elimination of serum Tg from the body is an exception to this.The follow-up of Tg levels after thyroidectomy in non-metastatic DTC patients allows us to characterize this process, with measured Tg half-lives ranging from 3.7 h to 4.3 days in humans [8,54,55].In our data fitting process, Tg elimination rates tended to take values lower than those clinically measured, and we set a maximum Tg half-life of 10 days.
• A relation between serum Tg level and tumor burden is commonly suggested, but there are very few quantitative studies on the subject.Bachelot et al. measured serum Tg levels after thyroid hormone withdrawal and tumor mass in patients with follicular and papillary thyroid cancer with lymph node metastases and negative TgAb tests [56].Patients had previously undergone resection of the thyroid gland and showed low (0.5%) or no significant uptake in the tumor bed, indicating that these metastases were the only source of serum Tg.Lymph node metastases were localized by a preoperative 131 I whole-body scan and completely removed for accurate mass determination.We used the mass and Tg data reported in this study to estimate an average value for bulk Tg production of our model.The large reported variability of tumor masses and Tg levels led us to set a range for the Tg production rate ([4.05 × 10 −4 , 1.1 × 10 −3 ] ng 2 /(ml 2 day −1 )).
Best-fitting parameter values for each set of experimental data are presented in the supplementary materials, Table SM1.It is observed that the cell proliferation rates required to fit the response of Graves' disease to RAI are lower than those obtained when fitting thyroid cancer data, which is an expected result due to the typically high proliferative capacity of tumor cells.Best-fitting parameters indicate that the radiosensitivity of cells in Graves' disease is lower than the radiosensitivity of tumor cells in DTC, which is expected due to the typically higher radiosensitivity of rapidly-proliferating tumor cells.

Conclusions
In this work, we have presented a mathematical model of thyroid response to radiotherapy.The model includes several populations of cells, undamaged (viable), sub-lethally damaged, doomed (non-repairable), and dead cells, as well as two tumor markers, the serum concentration of thyroglobulin and antithyroglobulin autoantibodies.These two substances are relevant in the management of thyroid cancer patients, being used as surrogates of treatment response in the follow-up after RAI, to evaluate whether tumor control has been achieved or make the decision to deliver further treatment.The model considers several processes such as cell proliferation, sub-lethal and lethal radiation damage, repair, elimination of dead cells, and the release, interaction and elimination of Tg and TgAb.To our knowledge, no previous radiotherapy thyroid response model has considered these two thyroid tumor markers together.Previous response models are restricted to the tumor mass [24,27], consider mass and Tg [23], or other antibodies relevant for Graves' hyperthyroidism [25].However, the interaction of Tg and TgAb is an important feature that is worth to be modeled, as a good understanding of it may contribute to a better understanding of treatment response and improve treatment design.
Models like this can assist in the interpretation of clinical data, being useful to interpret patient follow-up data, for example, they may help to determine the presence of subclinical burden from Tg/TgAb dynamics, or when to deliver additional treatment in multi-fraction therapies.A more ambitious role for this type of model is to assist in the design of optimal therapeutic strategies.Aiming precisely at this objective, we have developed a stochastic version of our model that can be used to compute tumor control probabilities, and have used it to investigate in silico treatment individualization according to patient-specific dose/radiosensitivity features.
Among the most significant findings of our work we highlight that (a) the model adequately fits the response to RAI of DTC and Graves' disease observed in experimental studies, and we hypothesize that (b) treatment individualization involving moderate variations of the injected activity can lead to important improvements of TCP.
The model has limitations though, which must be considered for clinical application.In particular, the model contains multiple fitting parameters, and parameter values should be taken with caution, as mass and dosimetry uncertainties strongly affect the values associated with radiosensitivity, interpatient variability is large and there is scarce literature about TgAb-associated parameters.We have not performed a direct intercomparison of fits obtained with our model and previous models (discussed in Section 1), as none of them considered all the variables included in our work and, therefore, a direct comparison was not possible.
Molecular radiotherapy, including RAI, is undergoing an important expansion, due to the good clinical results achieved.Further validation of the presented model will be possible once more RAI response data (including dose, thyroid mass/volume, and Tg and TgAb dynamics) become available.In addition, the same modeling methodology could be of interest to investigate the dynamics of different tumors treated with different molecular radiotherapies.Equations ( 1)-( 4), modeling tumor response, could be of direct application to different tumors, while new equations should be developed to account for the dynamics of relevant biomarkers.
Data Availability Statement: Data and code supporting the reported results are available from https://doi.org/10.7910/DVN/PBLNH1.

Figure 2 :
Figure 2: Sketch of stratification and therapy individualization.Three strategies are presented.In Strategy 1, simulated patients are split into 5 groups according to absorbed dose in the tumor.Dose thresholds (D 1 , D 2 , D 3 , D 4) are defined to have 1/5 of patients in each group.Strategy 2 is similar, but patients are stratified into five groups according to their radiosensitivity (characterized through the parameter of accumulation of sub-lethal damage, b).Sensitivity thresholds (S 1 , S 2 , S 3 , S 4 ) are defined to have 1/5 of patients in each group.Strategy 3 is a combination of both previous strategies: the same five groups of Strategy 1 are created (according to absorbed dose in the tumor), but each of them is further split into two subgroups (A and B) according to the radiosensitivity of tumor cells, as in Strategy 2. The threshold to split each group is defined as the median of the group, in order to have 1/2 of patients in subgroup A and subgroup B. For individualization, administration of 131 I activity A is prescribed differently in each group (Strategy 1 and Strategy 2) or subgroup (Strategy 3), as a function of activity perturbations ∆A and δA.Individualized activities are defined in such a way to verify that the mean injected activity in the whole cohort does not change, no matter what individualization strategy is used.

Table 2 :
List of parameters associated to the model, symbols used in equations, units and values calculated in the fitting process to: average population Graves' hyperthyroidism data (column 3), and individual DTC response data, one local recurrence and one metastatic patient, columns 4 and 5 respectively.(*) The same proliferation rate was considered for undamaged and sublethally damaged cells for the sake of simplicity.(**) The same value was used for ka and ka'.

Figure 3 :
Figure 3: Evolution of population-averaged thyroid markers with time after 131 I treatment for Graves' disease.Experimental measurements, asterisks with standard deviation error bars, and best fits obtained with the model presented in this work: (A) thyroid volume: total, thick solid line; and volumes corresponding to undamaged cells, dot line; doomed cells, dashed line; and dead cells, dash-dot line.(B) Tg concentration: total, solid line; free Tg compartment, dashed line; and bound Tg compartment, dash-dot line.(C) TgAb concentration.Panels (D-F) show 68% confidence intervals for thyroid volume, Tg, and TgAb, respectively (dashed lines), as well as population-averaged fits.In order to simulate a population, best-fitting parameters to the population-averaged data were perturbed with normal distribution with 25% relative standard deviation (1000 different sets of parameter values to simulate a population of 1000 patients).

Figure 4 :
Figure 4: Evolution of DTC local recurrence with time for a 131 treatment consisting of three administrations of 7400 MBq each.Experimental measurements, represented by asterisks with error bars, and best fits obtained with the model.(A) Tumor mass: total indicated by a thick solid line; masses corresponding to undamaged cells by a thin solid line; doomed cells by a dash-dot line; and dead cells by a dashed line.(B) Total Tg concentration.(C) TgAb concentration.

Figure 5 :
Figure 5: Evolution of three DTC lung metastases with time for a 131 I treatment with two administrations of 7400 MBq each.Experimental measurements, asterisks with standard deviation error bars, and best fits obtained with the model.Panels (A-C) show tumor masses corresponding to every lesion: total by a thick solid line; masses corresponding to doomed cells by a dash-dot line; and dead cells by a dashed line.Panel (D) shows total thyroglobulin concentration.

Table 1 :
List of reviewed mathematical models of thyroid/thyroid cancer response.

Table 3 :
Overall and group-specific tumor control probability when the injected activity is not adapted to each group (T CP ni ), and when it is adapted according to Strategy 1 with ∆A = 0.1 * A 0

Table 4 :
Overall and group-specific tumor control probability when the injected activity is not adapted to each group (T CP ni ), and when it is adapted according to Strategy 2 with ∆A = 0.1 * A 0

Table 5 :
Overall (whole cohort, WC) and group-specific tumor control probability when the injected activity is not adapted to each group (T CP ni ), and when it is adapted according to Strategy 3 with ∆A = 0.1 * A 0 (T CP ind1 ) or ∆A = 0.2 * A 0 (T CP ind2 ).The stratification in groups is performed according to absorbed dose in the tumor and radioresistance of tumor cells.