Brain Metastasis Response to Stereotactic Radio Surgery: A Mathematical Approach

: Brain metastases (BMs) are cancer cells that spread to the brain from primary tumors in other organs. Up to 35% of adult cancer patients develop BMs. The treatment of BM patients who have well-controlled extracranial disease and a small number of lesions consists of localized doses of radiation (stereotactic radio surgery (SRS)). Estimating prognosis among BM patients may allow treatments to be chosen that balance durability of intracranial tumor control with quality of life and the side effects of treatment. No mathematical model-based quantitative biomarkers have been determined for estimating prognosis. As a ﬁrst step toward that goal, we describe a mathematical model of growth and response of brain metastasis to stereotactic radio surgery. The mathematical model incorporates some biological mechanisms involved in BM growth and response to SRS and allows the observed dynamics to be accurately described.


Introduction
Cancer is one of the world's major health problems and is the second leading cause of death in industrialized countries. Vast resources have been devoted to cancer research in recent decades, which has resulted in only a small reduction in cancer death rates, of about 1-2% per year [1]. The question has been raised of whether different approaches, complementing classical ones, may be required to achieve more substantial breakthroughs in the war on cancer. One of these may be to increase the use of tools offered by applied mathematics. This is reflected in an increasing number of publications and reviews on the interface between mathematics and cancer [2][3][4][5].
Brain metastases (BMs) are cancer cells that have spread to the brain from tumors in other organs in the body. A substantial number, some 10-35% of adult cancer patients, develop BMs [6]. BMs are a major cancer-related complication and are 10 times more common than malignant primary brain tumors. The incidence of BMs is rising because improved systemic therapies control systemic disease and prolong survival but cross the blood-brain barrier (the highly selective membrane barrier that separates the circulating blood from the brain) too poorly to be able to control BMs.
The treatment of BM patients who have well-controlled extracranial disease and three or fewer measurable lesions typically consists of high localized doses of radiation (stereotactic radio surgery (SRS)) on the visible lesions and sometimes whole-brain radiotherapy (WBRT) to target potentially occult BMs. This therapeutic approach allows metastatic lesions to be controlled in many, but not all BM patients [7].
Very little is known of how to describe mathematically the dynamics of metastatic tumors in the brain and their response to therapies. BMs as seen in MRIs are composed of several compartments: the contrast-enhancing tumor, the necrotic core and an infiltration component.
Estimating prognosis among patients with BMs is clinically relevant, as it may allow clinicians to recommend treatments that balance durability of intracranial tumor control with quality of life and the side effects of treatment. The prognostic index that has become most prominent is the diagnosis-specific Graded Prognostic Assessment (GPA) [39]. This index includes prognostic factors that depend on the type of cancer, age, Karnofsky Performance Scale score, extracranial disease status, number of BMs, etc. For the case of primary brain tumors, mathematical modeling has helped in defining biomarkers of prognosis [40][41][42][43] and response to therapies [44][45][46].
As a first step in constructing prognostic and response biomarkers in BMs, we developed a minimal mathematical model able to describe the longitudinal dynamics of BMs. This was done using BM patient imaging data to feed ordinary differential equations (ODEs) in biologically grounded tumor growth models. First, the extent to which these models could describe the dynamics of untreated BMs was analyzed. Next, the dynamics of response to radiosurgery over time was studied when death mechanisms and damaged cell compartments were included.
The remainder of the paper is organized as follows: First, Section 2 sets out our mathematical model and describes the patient cohort used to validate it and the methodologies used to process the imaging data. Next, Section 3 describes the results. Finally, Section 4 discusses the implications of our study and summarizes our conclusions.

Mathematical Model of Response to Radiosurgery
Our mathematical model to describe growth and response to therapy of brain metastases is based on a set of ordinary differential equations for the different cellular compartments involved in a simplified description of the tumor growth dynamics. The first will be the proliferating tumor cells P(t). Since BMs are often small tumors that are detected at an early phase of growth, far from the brain's carrying capacity, the number of proliferating cancer cells can be described by an exponential model of the form In Equation (1), P(t) represents the number of proliferating cancer cells at time t, and ρ represents the rate of proliferation. A recent paper [47] has shown that the growth dynamics of brain metastasis could be better described by a superlinear growth law. For the sake of simplicity, in this paper we will keep the simpler and more classical growth law given by Equation (1).
Stereotactic radiosurgery is a highly precise form of radiation therapy that lethally damages a fraction 1 − S f of the irradiated proliferating cancer cells. Lethal damage, caused by high doses of radiation, cannot be repaired and occurs via different pathways. Here we will account for the two most important pathways in the response to SRS. The first process is fast and leads to a fraction of lethally damaged cells being incorporated directly into the necrotic cell (N(t)) compartment. The second process leads to a population D(t) of lethally damaged cells that may undergo one or more divisions before dying by mitotic catastrophe [48]. The evolution of this cellular compartment is given by the equation where k is the average number of mitosis cycles completed before death. In practice, k/ρ provides a, typically long, time scale for the death of this population of damaged cells. This type of model have been found to provide a good description for the response of low-grade glioma [35] and prostate cancer [49] to radiation therapy. There is evidence in the literature that SRS destroys tumor vascular beds, thereby deteriorating the tumor microenvironment and leading to indirect tumor cell death [50]. These phenomena can produce regions of hypoxia, tumor necrosis and massive release of tumor antigens, elevating antitumor immune response in a short period of time. Furthermore, tumor hypoxia may persist after vascular injury caused by SRS, and both oxygenated and hypoxic cells are ablated by high-dose radiation [51]. This fact, together with the impossibility of accounting for the number of hypoxic cells in brain metastases through MRI, motivates the non-inclusion of this population within the mathematical model. Thus, we will also account for a population of immune cells I(t) present in the tumor, whose dynamics would be governed by the equation In Equation (3), α is a stimulation parameter that accounts for the immune system activation by the presence of necrosis, and λ I is the decay rate of the immune activation.
Finally, the dynamics of the necrotic cell (or necrosis) compartment will be governed by the equation Thus, we will describe necrosis dynamics by the contribution of damaged cells through mitotic death and the interaction with the immune system, which in turn is activated by the release of a variety of pro-oxidant and proinflammatory cytokines such as tumor necrosis factors [52]. The first term in Equation (4) corresponds to the contribution of damaged cells, and the second describes the elimination of necrotic cells by the action of immune cells. This process is modulated by the constant λ N .
As to the initial data, the tumor is assumed to be composed mostly of proliferating cells before treatment. SRS is performed at a given time t 0 on a tumor cell population T(t 0 ). As stated above, a fraction of tumor cells S f will suffer either no damage or only sublethal damage and will remain viable. A further fraction (1 − S f ) will receive lethal damage of which a fraction will die on a short time scale (i.e., days), and the remainder will move into the compartment of lethally damaged cells. This means that where S f and ∈ (0, 1). Figure 1 summarizes the different compartments included in the model and the effect of radiosurgery. In this paper, we assumed the pretreatment number of immune cells to be very small. Intratumoral areas in lung cancer brain metastases have been reported to contain low numbers of inflammatory cells [53], high levels of immunosuppressive molecules [54] and a suppressed immune microenvironment [55]. In addition to damaging cancer cells, SRS damages healthy tissue cells located in the field of the ionizing radiation. In spite of the increasing spatial precision of current radiosurgery techniques, radiation necrosis of normal cells is a frequent complication. Since brain tissue renewal occurs on very long time scales of many months, this damage appears typically more than a year after SRS [56]. This effect has nothing to do with tumor recurrence but can be confused with it on imaging because of the increase in the inflammatory compartment. To account for that, we will include an explicit additional source term in the necrotic compartment h(t) in Equation (4): Note that T(t) = P(t) + D(t) + N(t) + I(t) is a measure of the number of cellular elements in the tumor and its environment. We assumed it to be related to the observed tumor volume, which appears in MRIs as a combination of active areas (P + D + I) plus a necrotic component (N).
The model was solved and fitted to the available longitudinal volumetric data (see below) using ode45 and fmincon functions included in the scientific software package MATLAB (R2019b, The MathWorks, Inc., Natick, MA, USA). We performed a heuristic procedure to validate that the fitted model attained a global minimum on the set of parameters estimated (see Appendix A for a particular example). All parameters were fitted for each metastasis individually, even in the case of lesions belonging to the same patient.

Patients
Patients included were participants in the study METMATH (Metastasis and Mathematics), a retrospective, multicenter, nonrandomized study approved by the Ethics Committees of Instituto Valenciano de Oncología (code 2011-46, 2011), Hospitales HM (code 18.09.1303-GHM, 2018) and Hospital Universitario Ramón y Cajal (code METMATH, 2018). We reviewed the METMATH records to look for patients satisfying the following inclusion criteria: patients diagnosed with brain metastasis of a primary lung cancer that had undergone at least a T1-weighted MRI examination with contrast before SRS and at least two T1-weighted MRIs after SRS. The time interval between the diagnostic MRI and SRS had to be at most 2 weeks. Patients who received whole-brain radiation therapy within 4 months before SRS, or during the followup period, were excluded. Patients who received surgery were also excluded. In total, 45 patients satisfied the inclusion criteria. Patients with more than one treated brain metastasis were analyzed; 32 patients had only one BM, 11 patients had two BMs and 2 patients had three BMs. Finally, 60 brain metastases were included. Median patient age was 60 years (range 43-80) and sex was 59% male, 41% female.
The dose and fractionation schedule was chosen at the discretion of the treating radiation oncologists and performed using γ-rays. A total of 46 lesions were treated with a single dose (range 17-24 Gy), six lesions with three dose fractions (range 5.5-8.8 Gy) and eight lesions received between four and six dose fractions (range 4-9.5 Gy). The mean follow-up period of the patients studied was 11 months (range 3-30 months).
Contrast-enhanced T1-weighted sequence was gradient echo transformed using 3D spoiled-gradient recalled echo or 3D fast-field echo after intravenous administration of a single dose of gadobenate dimeglumine (0.10 mmol/kg) with a 6-8-min delay.
MRI studies were performed in the axial or coronal plane with either a 1.5 T Siemens scanner (259), a 3 T Philips scanner (45) or a 1 T Philips scanner (13). Imaging parameters were no gap, slice thickness of 0.5-2.0 mm (mean 1.1 mm), 0.4-1.1 mm (mean 0.6 mm) resolutions in the x and y planes and 0.6-2.0 mm spacing between slices (mean 1.5 mm).

Tumor Segmentation
Brain metastasis T1-weighted images were collected in DICOM format and analyzed by the same image expert (OLT, 2 years of experience on tumor segmentation) and reviewed by either a senior radiologist (EA) or an image expert with 5 years of expertise in tumor segmentation (JPB). The tumor volume for each brain metastasis was defined on gadolinium-enhanced magnetic resonance imaging (T1Gd-MRI) as the contrast-enhancing (CE) compartment of the tumor combined with the central non-enhancing (non-CE) compartment enclosed by the contrast (the latter usually represented necrosis).
Segmentation was performed by importing image files into the scientific software package MATLAB (R2019b, The MathWorks, Inc., Natick, MA, USA). Images were manually placed in a 3D box and then semi-automatically delineated using a gray-level threshold chosen to identify the contrast-enhancing volume. Segmentation was corrected manually slice by slice as described in [45]. Figure 2 shows examples of the time evolution of the volume of brain metastases treated with SRS in three patients. The cases chosen provide instances of three different typical behaviors. In the first case (Figure 2a), the patient had a sustained response lasting for at least 17 months. In a second patient (Figure 2b), the tumor decreased in volume for approximately 184 days and relapsed after the initial response. Finally, in the third case (Figure 2c), the tumor continued to grow after radiosurgery treatment. Henceforth, we will denote the lesion behaviors described in the cases shown as monotonically decreasing lesions (MDLs), lesions that are first decreasing and later increasing (DILs) and monotonically increasing lesions (MILs). Our mathematical model was able to describe the different scenarios shown in Figure 2. Figure 3 shows additional examples for 9 of the 60 metastases. The model was able to describe different dynamics of response to treatment in all cases. It also provided good fits for the longitudinal volumetric data of 36 metastases, which showed evolution framed within one of the three behaviors mentioned above: 6 metastases showed an MDL behavior, 19 showed a DIL behavior and 11 showed an MIL behavior. This study included only 13 patients with more than one metastasis treated, and thus a statistical analysis to assess similarities in the behavior of the lesions could not be performed. Four of the 11 patients were observed to have the same post-treatment dynamics in all of their metastases, but the other seven patients had lesions with different longitudinal dynamics after treatment therapy.

The Mathematical Model Describes the Early Inflammatory Dynamics Observed in the Post-SRS Response
Six patients in our dataset had imaging studies performed within the period of 3 months after SRS. These BMs increased in size initially, suggesting treatment failure. However, they decreased in size after a few weeks, which was shown by the second examination after SRS. This is due to an early inflammatory response to the treatment that is assumed to occur in most cases. In fact, this is the reason why the first control post-treatment MRI is typically performed 3 months after SRS once the effect of early inflammation has vanished. Figure 4a,b shows the longitudinal volumetric dynamics of BMs and the best fits using the mathematical model for two of those patients.
It is interesting that when breaking the total tumor volume into the different components (color lines in Figure 4), it was obvious that the initial volumetric growth was actually associated with an increase in immune activation. In both cases, radiation generated a large initial number of necrotic cells that stimulated immune cells. The sum of all the cellular compartments resulted in an increase in the total volume, but in fact, the number of proliferating cells was smaller after SRS treatment.
We calculated T(I max ), representing the time when the immune cell population reaches its peak for the model fits. Figure 4c shows the histogram of the calculated times for the 42 metastases classified within the above response groups. The computed times can be approximated by a gamma distribution with shape parameter a = 1.4 and scale parameter b = 36, and they have a mean value equal to 51 days and a variance equal to 43 days. This means that the peak of inflammation occurs more frequently between months 1 and 2 after SRS. Thus, our model confirms that when MRI scans are performed within 3 months after treatment, the results may be affected by early inflammation events and do not provide a reliable measure of the proliferating tumor component.
It is important to emphasize that all four classes of response dynamics studied so far (monotonically decreasing lesions, lesions that are first decreasing and later increasing, monotonically increasing lesions and early inflammation) were accurately described by the mathematical model given by Equations (1)-(4). Figure 5 shows the mean square errors (MSEs) of the best fits obtained using the mathematical model for this subset of brain metastases.
where n is the total number of follow-up MRIs performed for each metastases, V i is the segmented volume of scan i and V i is the volume estimated with the model at the time of scan i.

Damage to Healthy Tissue Could Lead to Late Inflammatory Response and Radiation Necrosis
Late radiation necrosis is a frequent event after SRS treatment [57,58]. This is typically observed in MRIs as volumetric growth of the lesion, typically between 1 and 2 years after SRS, followed by spontaneous (partial or complete) remission. Radiation necrosis (also called "pseudoprogression") poses a challenge to radiologists since it is very difficult to differentiate it from true tumor progression. In the former case, there is no need for anti-tumoral treatment, while the latter requires a different therapy.
In our dataset, 7 BMs were diagnosed with radiation necrosis and an additional 11 BMs presented late episodes of volume increase followed by volume decrease compatible with that condition. Figure 6 shows two examples of the dynamics and the best fits obtained, using different mathematical models. First, we tried to fit the dynamics using Equations (1)- (4). However, the model did not accurately describe the dynamics (see dashed blue lines in Figure 6a,b) with large MSEs (Figure 6e). To explain the inflammatory response in this group of patients, we included the late damage to the surrounding healthy tissue as described by Equation (9). To do so, the parameters in Equations (1)- (3) and (9) were fitted for each metastasis using a Gaussian form for the damage function h(t) = kNormal(µ, σ). A substantial reduction in the MSE was obtained (Figure 6f) in line with the model dynamics, closely resembling the data (see solid blue lines in Figure 6a,b).

Time to Tumor Progression Can Be Obtained from the Mathematical Model
Finally, our mathematical model allowed for a theoretical estimation of the time of tumor progression from the initial response data.
After therapy, proliferating cells have exponential growth given by the explicit solution P(t) = P 0 e ρt . Similarly, damaged cells have exponential decay given by D(t) = D 0 e −ρt/k . When dT/dt > 0, tumor volume will grow back. In ideal conditions, the immune system is able to efficiently counteract necrotic cells, and both populations would eventually disappear. Under this assumption, tumor progression is determined by the populations of proliferating and damaged cells, i.e., P 0 e ρt − 1 k D 0 e −ρt/k > 0. Taking the most reasonable value, k = 2, i.e., damaged cells dying after two cell cycles on average, and substituting the values of P 0 and D 0 using Equations (5) and (6), the previous inequality is satisfied if  In the first example (see Figure 7a,d), in a tumor without inflammation, the progression time clearly indicates the moment in which the proliferating cells outgrow half of the damaged cell compartment, thus dominating tumor growth. For this patient, there is an increase in tumor volume right after the progression time t p computed from our mathematical model (Equation (10)).
The second and third examples correspond to patients diagnosed with radiation necrosis. The progression time could be used in these cases to distinguish between an increase in the tumor volume caused only by the radiation necrosis (false progression) and an increase caused by real growth. Specifically, the second example (Figure 7b,e) corresponds to a patient with radiation necrosis without progression. In this case, the tumor volume increase can be associated with the radiation necrosis, as it occurs 10 months before progression was expected. In the third example, radiation necrosis occurs after the progression time. In the graph, it is possible to see how even if the tumor volume increases due to the radiation necrosis, it continues growing after the inflammation decreases, suggesting the coexistence of radiation necrosis and tumor progression.

Discussion and Conclusions
In this study, we put forward a mathematical model that describes the effect of SRS on brain metastases. The mathematical model includes only four cellular compartments (proliferating, necrotic, damaged and immune) to account for the different dynamics observed after SRS. Interestingly, the simple model was able to describe the volumetric evolution of the lesions observed in the clinic and the different scenarios related to inflammation.
We wanted to keep the complexity of the model at a minimum. It is easy to construct complex models accounting for many different, and not too well-known, biological processes. However that often results in too many unknown parameters to be fitted using a limited amount of biological data. In our case, the only available information was the MRI data and, more specifically, tumor volume. Retaining only a small number of essential variables/parameters is the best way to avoid overfitting problems. In our mathematical model, we used a single population of proliferating metastatic cells accounting effectively for the tumor dynamics. Brain metastases, and specifically lung cancer ones, are heterogeneous and composed of genetically and phenotypically different subpopulations [59][60][61]. This heterogeneity leads to a differential response to radiation therapy depending on many molecular factors such as EGFR overexpression [62], TopBP1 and Claspin [63], MET [64], CAVEOLIN-1 [65] and many others [66,67]. Indeed, radiation resistance is not only associated with specific individual features of the cancer cells. It is being increasingly recognized that the tumor microenvironment, cell-cell communications and other factors play a role in this complex emergent property [68][69][70]. Here, we intended to develop a minimal model able to describe observed volumetric dynamics, which is the only follow-up information available through standard MRI-imaging.
It is interesting that the model was able to describe qualitatively the observed dynamics without accounting for other relevant biological details of response to radiation therapy such as the differential response of well-oxygenated and hypoxic cells to radiation therapy [48] that has been included in other, more detailed models of radiation therapy of brain tumors by different authors (see e.g., [28,[71][72][73][74]). Another effect not accounted for in this study was the fact that irradiated cells may undergo cell cycle arrest or become quiescent for long time periods after sublethal radiation doses [33]. As imaging techniques progress and provide biological data on metastasis status able to feed more complex mathematical models, it may be relevant to develop models including all those biological processes.
In the context of our simplifying assumptions, we assumed the proliferation rate of the BMs to be the same after treatment as before it. The phenomenon of growth acceleration of residual tumors of certain histologies after fractionated radiotherapy courses, the so-called accelerated repopulation, has been known for a long time [48]. However, repopulation rate has been found to be lower in tumors with increased cell loss, as happens in radiation surgery [75]. In fact, there are no reports of observations of accelerated repopulation after radiation surgery treatments. This is also consistent with the observation in experimental tumors that the triggering of accelerated repopulation requires minimum total treatment durations longer than than 3-4 weeks in human tumors [48].
Other biological processes could also lead to differences in proliferation rate with time. The possibility of actions on the primary tumor influencing the dissemination and growth dynamics of metastasis has been hypothesized on the basis of different sources of biomedical data and studied mathematically [19,76,77]. However, in the case of brain metastases, the blood-brain barrier may provide a chemical firewall for the communication between the primary tumor and metastatic colonies, and the fact is that no solid evidence of such chemical links has been provided. Moreover in the case studied here of stage IV lung cancer, surgical treatment of the primary tumor is not typically a therapeutic option. This is why we did not account for any "communication" between the primary tumor and the BMs in our modeling approach.
Early post-treatment inflammation was very easily accounted for in the model by incorporating early cell death (necrosis) and associated inflammatory response. Radiation necrosis events, assumed to be the result of damage to healthy tissue, were incorporated in a simple and rather ad hoc way by the function h(t). The effect of radiation on healthy tissue as an additional source of necrosis could be analyzed under different forms of this function, but we have chosen the Gaussian form as a first approximation. Induced damage to normal tissue by radiation therapy has been studied using mechanistic mathematical models [78], where the cell population is directly exposed to radiation. Incorporating more details of the SRS-induced damage to healthy tissue could allow us to write mechanistic equations for the normal tissue behavior in response to treatment. As a first approach, we wanted to keep the model simple and with a minimal number of parameters, but more complex mathematical models could account for this effect in a more elegant way.
In this work, we chose S f to be an adjustable parameter instead of a fixed predetermined value. In the context of elementary "static" mathematical approaches to radiation therapy such as the linear-quadratic equation, it is customary to assume a given S f value to solve problems in therapy replanning and other such applications. That given S f value would be characteristic of the tumor histology. However, that approach is of limited use in the context of radiation surgery of brain metastasis because there is a broad variety of potential outcomes to SRS treatment ranging from no response to complete response even for the same primary types. This is probably due to the intra-and interpatient heterogeneity as well as the role of many other elements in the response such as the tumor or immune environments, as discussed previously.
According to the results of the model fits, it seems that this simple model has some limitations in describing the tumor dynamics in some patients undergoing radiation necrosis. The model fitting in these cases showed mean square errors larger than those obtained in the more common scenarios of early response to SRS. This suggests that radiation necrosis is a more complex phenomenon where different biological factors may play a role.
One of the current problems in clinical practice is to differentiate between tumor progression and radiation necrosis, as both display a similar course as observed in the MRI, requiring advanced diagnostic techniques for identification [79,80]. In this study, we obtained an analytical estimate for the progression time due to the growth of remnant proliferating cells after SRS. This estimated time may be the basis for biomarkers helping clinicians to distinguish between progression and radiation necrosis. If the tumor increases in size at a time close to the predicted tumor recurrence, this could be an indication of tumor progression. On the other hand, tumor outgrowths at times much earlier than the estimated by the model could be an indication of the presence of an inflammatory process.
The key point in moving from Equation (10) to a clinically relevant biomarker would be to define the requirements for a reliable estimate of the outgrowth time, as this must be done from MRI data. In principle, only a small number of data points are necessary to obtain the value of this parameter (ρ, , S f ). If there are two MR imaging exams before SRS, they would allow the tumor proliferation rate ρ to be estimated. Interestingly, performing a second MRI immediately before SRS also has advantages from the clinical point of view, since it would allow for a more precise definition of the target volume. It has been reported recently [81,82] that measurable changes occur in brain metastasis over a short period of time, on the order of a week, so a final planning right before SRS would help in achieving greater therapeutic efficacy and provide a second MRI to obtain a baseline growth rate estimate.
The parameters and S f can be obtained using the first few (2-3) MRIs of the standard follow-up after radiosurgery. Thus, patients with two pretreatment MRIs and a 1-year follow-up after SRS could have precise progression time estimates that could be used for comparison with the observed dynamics after the first year. Unfortunately, our database did not have patients fulfilling these requirements (i.e., all patients who presented radiation necrosis or progression had only one pretreatment MRI), so we could not explore the idea in more detail. As future work, we will address this by planning either a retrospective search for such patients or a prospective study.
This work is one of the first mathematical studies describing the different outcomes observed in the response of brain metastasis to radiosurgery using longitudinal tumor imaging data. Brain metastasis is a condition 10 times more frequent than primary brain tumors. It is thus striking that so many mathematical papers have studied primary brain tumors and their time dynamics while BMs have received little attention. Our model is a first step to a mathematical description of that disease. We hope that this work will stimulate further studies addressing this important condition and that richer datasets could be used to feed more elaborate models in the future.
In conclusion, we have developed a mathematical model based on a set of ordinary differential equations that describe the observed longitudinal dynamics of brain metastases after SRS. The model allowed the varying early longitudinal dynamics observed in patients to be accurately described with very few parameters. Radiation necrosis events were described in a simplified way, and in most cases, they were also fitted accurately using the model. We obtained an equation for the expected tumor progression time based on a few parameters that could be the basis for biomarkers to help in discriminating between radiation necrosis and true progression, which currently represents a major challenge in the clinical setting. Informed Consent Statement: Patient consent was waived as the data were obtained from retrospective observational clinical studies that were approved by the corresponding institutional review boards.

Data Availability Statement:
The data presented in this study are available upon request from the corresponding author.

Conflicts of Interest:
The authors declare no conflicts of interest.

BMs Brain Metastases SRS
Stereotactic Radio Surgery RT Radiation Therapy MRI Magnetic Resonance Imaging

Appendix A
To confirm that the optimization procedure implemented in the MATLAB function fmincon was not trapped in a local minimum of the parameter space, we heuristically explored parameter values around the solution. To do so, we plotted the profiles of the objective function by fixing all parameters except one at their optimal values and considered the unset parameter as a variable. The function used is given by the expression MSE = 1 n ∑ n i=1 ( V i (ρ, S f , , λ N , λ I , α, h) − V i ) 2 , where n is the total number of follow-up MRIs performed for each metastases, V i is the segmented volume of scan i and V i (ρ, S f , , λ N , λ I , α, h) is related to the total number of cellular elements in the tumor and its environment T(t i ) given by the model (see Section 2.1). Figure A1 shows two examples of minimum MSE profiles as a function of a single parameter, where the other parameters have been set to the optimal value obtained by fitting the data using fmincon. The simulations suggest that the set of parameter values provided by the MATLAB routine were indeed global minima, although the procedure used cannot guarantee it with absolute certainty.