Radiation necrosis after radiation therapy treatment of brain metastases: A computational approach

Metastasis is the process through which cancer cells break away from a primary tumor, travel through the blood or lymph system, and form new tumors in distant tissues. One of the preferred sites for metastatic dissemination is the brain, affecting more than 20% of all cancer patients. This figure is increasing steadily due to improvements in treatments of primary tumors. Stereotactic radiosurgery (SRS) is one of the main treatment options for patients with a small or moderate number of brain metastases (BMs). A frequent adverse event of SRS is radiation necrosis (RN), an inflammatory condition caused by late normal tissue cell death. A major diagnostic problem is that RNs are difficult to distinguish from BM recurrences, due to their similarities on standard magnetic resonance images (MRIs). However, this distinction is key to choosing the best therapeutic approach since RNs resolve often without further interventions, while relapsing BMs may require open brain surgery. Recent research has shown that RNs have a faster growth dynamics than recurrent BMs, providing a way to differentiate the two entities, but no mechanistic explanation has been provided for those observations. In this study, computational frameworks were developed based on mathematical models of increasing complexity, providing mechanistic explanations for the differential growth dynamics of BMs relapse versus RN events and explaining the observed clinical phenomenology. Simulated tumor relapses were found to have growth exponents substantially smaller than the group in which there was inflammation due to damage induced by SRS to normal brain tissue adjacent to the BMs, thus leading to RN. ROC curves with the synthetic data had an optimal threshold that maximized the sensitivity and specificity values for a growth exponent β* = 1.05, very close to that observed in patient datasets.


Introduction
Brain metastases (BM) are the most common intracranial malignancies in adults, with around 25% of patients with cancer developing brain metastases during the course of their diseases [1,2].Most BMs originate from primary cancers from lung, breast or from melanoma [2].Improvements in detection and primary tumor treatment leading to longer survival of cancer patients have resulted in an increase in the incidence of BMs in the last years [3].
Stereotactic radiosurgery (SRS), a form of external beam radiation therapy, is the preferred treatment option for the management of BMs in patients with less than five lesions because of its excellent local control rates [4].
Radiation necrosis (RN) is an inflammatory reaction that appears between 6 and 24 months following SRS in 5% to 25% of treated patients [5].It appears as contrast-enhancing lesions, resembling the appearance of tumor progressions on routine magnetic resonance images (MRIs) what poses a critical problem for physicians.RN may resolve spontaneously, and does not require further work-up, while progressive BMs require a prompt therapeutic action, typically in the form of open brain surgery.Thus the distinction between the two conditions is a problem with direct clinical implications.
There are two main pathophysiological theories for the development of RN in the brain, and it seems likely the combination of both is the best explanation for the process.One explanation suggests that radiation damages different types of healthy brain cells, resulting in apoptosis [6], while another theory suggests that radiation damages the vasculature [7].Damage to endothelial cells, which causes cell death in a delayed manner [8] and the release of inflammatory mediators [9], is a critical factor in the emergence of RN.Astrocytes, oligodendrocytes, and oligodendrocyte progenitor cells are also damaged by radiation [10,11] resulting in the release of hypoxia-inducible factor 1α and VEGF [12].Finally, VEGF induces ICAM-1, increasing the inflammatory response and edema, together with the enhancement seen on CE images.Necrotic cells release a variety of damage-associated molecular patterns (DAMPs) that act as signals to attract pro-inflammatory innate immune cells (macrophages and neutrophils) and lymphocytes.This influx of macrophages and neutrophils is considered an inflammatory response, as seen during infection.Once clearance of the dead cells is complete, the macrophages switch to an anti-inflammatory, pro-resolving phenotype that supports tissue regeneration [13,14].
Although biopsy sampling is considered the gold standard for distinguishing between tumor progression and radiation necrosis, it is subject to certain limitations, particularly in cases where biopsy sampling poses a safety concern, such as in patients with brainstem metastases [4].Conventional MRI is not appropriate for discriminating between both conditions since they share similar characteristics.Several alternative methods have been employed to differentiate tumor progression and RN, including perfusion MR imaging [15], positron emission tomography [16], the use of artificial intelligence [17] and radiomic methods [18,19], or changes in anatomical boundaries [20].However, many of these methods lack independent validation, require complex computational methodologies, or are in need of further study.Consequently, there is a need for diagnostic techniques based on routinely-used imaging methods, such as post-contrast T1-weighted images (T1-WI) [20].
Recently, a simplified version of the Von Bertalanffy [35] growth equation was applied to observational longitudinal volumetric data on lesion growth in patients [36,37].This type of analysis revealed its usefulness in discriminating between conditions, showing that RNs exhibit quicker growth dynamics compared to recurrent BMs [38].
The aim of this study was to understand the longitudinal dynamics of RN development and how different are its growth patterns from those of progressive lesions.To do so, two mathematical models were built.The first, a simple compartmental model, allowed the growth of RN events to be understood.The second, based on a discrete stochastic simulator, allowed the response to SRS to be studied, including more details of the biological processes behind post-SRS inflammation.

Materials and methods
In the planning of radiation therapy and despite improvements in the accuracy of radiation therapy treatments, low doses of radiation are administered to the tissue surrounding the lesion to ensure that not only visible tumor cells but also those that are not visible in medical images are eliminated, to reduce the likelihood of tumor recurrence.However, this approach has the unavoidable consequence of damaging healthy tissues around the metastatic lesion.
While radiation kills tumor cells and reduces tumor volume, it also affects healthy cells, which die by mitotic catastrophe when they are trying to renew.Renewal time for tumor cells is fast, whereas brain cells take months to reach the renewal point.This damage to healthy cells and the subsequent inflammation process governed by the immune system can lead to an apparent increase in tumor volume, which is typically observed several months or even years after SRS [39] as depicted in Fig 1.

Modelling the growth dynamics of radiation necrosis at a cellular level using a compartmental model
To simulate the development of radiation necrosis over time, a simple mathematical model was first developed, describing the evolution of cells taking into account the impact of SRS on both the tumor cells and the healthy cells in the tissue surrounding the lesion, as illustrated in Fig 2 .To calibrate the model, we used MRI data, which resulted in a distribution of model parameters.By using these distributions, our goal was to ensure that the growth exponent obtained was consistent with the expected values, demonstrating consistency regardless of initial volume and parameter variations.
It was assumed that the dynamics of the three leading populations involved in radiation necrosis events, tumor cells T(t), necrotic cells N(t) and immune cells I(t) are governed by the following equations: where T + N + I represents the total volume of the lesion.The initial time t = 0 corresponds to the first available MRI, usually around 3 months after SRS.Eq (1a) describes the growth of the tumor cells, which is assumed to be exponential.The first term in (1b) accounts for the dynamics of healthy cells that are killed by the late effects of radiation when trying to renew as given by which accounts for the fact that cells have the same probability of reaching the renewal point during a given time period t N , assumed to be 2 months, as it has been identified as the duration wherein the effects of irradiation are significant [40].The value of k provides information about the number of healthy cells that die at each time.This means that the total number of cells killed by radiation is given by HðtÞ dt: ð3Þ H will depend on the initial BM volume, since typically, a region of 1 mm width surrounding the tumor is irradiated to kill potentially infiltrating tumor cells.This is the typical irradiation margin used in BM radiation therapy treatments.Hence, we can estimate k � H T /60 by assuming that an equal number of cells will die each day during the 2-month period, t N (60 days).
The second term in (1b) represents necrotic cells cleared away by immune cells I(t).Eq (1c) assumes that the immune response is stimulated by the presence of necrotic cells, while they have a growth rate θ and a death rate λ I .
Parameters estimation for model Eq (1).The growth rate ρ was computed for each lesion from the volumes in the first two points at which growth was observed before any treatment was applied, i.e.
This choice implicitly assumes that the tumor preserves the same growth dynamics after treatment, an assumption taken here for simplicity, although the real dynamics would be more complicated [37].
The death rate λ I was taken to be 0.07 days −1 , leading to a half-life of around two weeks, that is the typical lifetime of activated effector immune cells [41,42].
Parameters λ N , γ and θ were estimated at a lesion level by fitting the model to available longitudinal volumetric data from 20 patients in [43].The data consisted of three volumetric time points for each BM, displaying longitudinal volumetric growth.The MATLAB functions fminsearch and ode45 were used to perform the fitting, returning values in the range [1.85 − 2.36] � 10 −11 days −1 , [1.84 − 2.06] � 10 −7 days −1 and [0.139 − 0.249] days −1 , for λ N , γ and θ, respectively.Thus, the parameter values obtained were very consistent for the different BMs.The immune growth rate θ was found to be about 5 days, which seems a reasonable value.In contrast, the necrotic cell elimination rate λ N and the rate of activation of immune cells γ were found to be very long.This finding may be influenced by the delayed effect of healthy cell death and the activation of the immune system associated with the appearance of RN.Table 1 shows a summary of all parameters and variables used in the compartmental model.
Calculation of the β exponent.To compute β for the simulations of model ( 1), the simulated time interval was split into three time slots, and for each slot, a random time value was chosen.This method was designed to mimic the clinical imaging follow-up, where scans are typically performed in three-months time windows but with a substantial variability due to real-world constraints.Those times t 0 , t 1 , t 2 and the computed volumes V 0 , V 1 , V 2 were used to compute the growth exponent β using the simplified Von Bertalanffy growth equation [35] Solving Eq (5), gives Since there is information on the dynamics at three-time points (t 0 , V 0 ), (t 1 , V 1 ) and (t 2 , V 2 ), the two parameters α and β can be completely determined by evaluating (6) at the times t 0 , t 1 , t 2 , giving Eq (6) is an algebraic equation for β that was solved by using the MATLAB function fzero (which returns the root of a nonlinear function) for each set of known values A dataset of 500 virtual BMs was generated for each choice of initial volume, ranging from 0.5 cm 3 to 3.0 cm 3 .The parameters λ N , γ, and θ were allowed to take uniformly distributed random values within ranges chosen to include the parameters obtained when fitting volumetric data.Specifically, λ N was assigned a value within the range [1.8 − 2.4] � 10 −11 days −1 , γ within [1.8 − 2.1] � 10 −7 days −1 , and θ within [0.14 − 0.25] days −1 .By varying the initial volumes and parameter values, it was possible to obtain a wide range of β values and to analyze the resulting trends.

A discrete stochastic simulator for the response to SRS and radiation necrosis classification
The mathematical model used for the prior analysis included only a limited number of relevant biological elements.A fuller study was carried out using a discrete stochastic brain metastasis simulator (DSBMS) based on Ref. [44].The model included six types of cellular population: normal and damaged healthy cells, normal and damaged tumor cells, immune cells and necrotic cells.For these populations, the main biological processes were implemented at the cellular level: mitosis, migration and cell death.
The discrete model focuses on describing cell populations instead of individual cells.The spatial domain was set as a 3D grid discretized in cubic compartments (voxels) of side length 4x, fixed at 1 mm.Each voxel has a specific dynamic that depends on its occupation and surroundings and can contain cells from each subpopulation with an upper limit indicated as local carrying capacity K.The different cell populations attempt to perform all basic processes at each time step.These processes can be described by a binomial distribution with a probability associated with the process.That is, the number of cells successfully undergoing division, death, migration or transition to another population within the discrete model are calculated voxel-wise and state-wise at each time step.These numbers are thus calculated by randomly sampling the corresponding binomial distribution, whose N will be the number of cells in each population within a voxel, and whose probability will be the rate of the process modulated by the time step length 4x.All processes, their probabilities and associated binomial distributions are described below in detail.
Tumor growth in the discrete model.To simulate the growth of the metastatic lesion before SRS, the biological processes of cell division, death and migration were implemented similarly to the model in [44].Let n t , n n and n h be the total numbers of tumor, necrotic and healthy cells inside a given voxel, respectively.The total numbers of proliferating and dying tumor cells are drawn from the binomial distributions Bðn t ; P t p Þ and Bðn t ; P t d Þ, respectively.
The number of migrating tumor cells is drawn from the respective binomial distribution Bðn t ; P t m Þ.Then these cells are distributed around a neighborhood of 26 voxels (Moore neighborhood) according to a multinomial distribution Y * Mult(X, P Moore ), where Y is a vector of 26 components giving the number of cells that migrate to each voxel, X is the number of migrating cells that has been previously computed and P Moore is a vector of 26 components giving the probabilities of migration to each surrounding voxel [44].The probability P Moore is proportional to 1, 1= ffi ffi ffi 2 p or 1= ffi ffi ffi 3 p depending on whether voxels share a face, edge or vertex with the central voxel respectively.Performing the migration in this way reproduces a diffusive process, wherein migration depends on cell density gradients and distances.For simplicity, it was assumed that all tumor cells belong to the same clonal population without including mutation events or phenotype changes.
For healthy cells, it was assumed that the levels of cell division and death remain balanced due to the ability of these cells to self-regulate.The biological process of migration is the only one that is affected by the evolution of tumor cells.Therefore, the number of migrating healthy cells is drawn from the binomial distribution Bðn h ; P h m Þ being displaced by the pressure exerted by the tumor cell colonization when the total number of cells in the voxel exceeds 45% of its maximum capacity.Fig 3A shows a slice of an actual simulation, where the colors indicate voxel occupation.Note that each voxel can contain a different number of cells.
To construct tumor growth rules in the DSBMS the probabilities of tumor cells reproduction, migration and death were first considered to be given by These probabilities were modulated by the relationship between the time step and the characteristic time of the process.The reproduction probability decreases with voxel occupation while the migration probability increases, simulating competition for space and resources.To mimic the effect of local vascular damage and hypoxia on tumor cell apoptosis signaling, it was activated once the voxel occupation reached a cell number greater than 75% of its carrying capacity.This was incorporated into the probability of death P t d by the hyperbolic tangent term as a function of the voxel capacity.
As to motility, the probability of healty cells migration was given by All probabilities were calculated at each time step and voxel, providing the number of cells undergoing mitosis, apoptosis and migration for every population.Parameter τ represents the characteristic times of each process.K is a carrying capacity.
Response to radiosurgery.A single dose of SRS was simulated in-silico.SRS leads to a fraction of tumor cells S f surviving without damage and remaining viable, a fraction (1 − S f ) receiving lethal damage.Of those, a fraction � dies on a short time scale (i.e.days), and the remaining fraction 1−� moves into the compartment of damaged cells.Radiation therapy induces immune cells and lethal damage to a fraction (1 − S n ) of healthy cells surrounding the Let n d , n i and n hd be the total number of damaged tumor cells, activated immune cells and damaged healthy cells inside a given voxel, respectively.The numbers of damaged tumor and healthy cells that die by radiation are drawn from binomial distributions Bðn d ; P d d Þ and Bðn hd ; P h d Þ, respectively.Then, the number of necrotic cells increases by adding these populations.Additionally, the immune system is activated and immune cells move to the irradiated region to remove necrotic cells.The number of necrotic cells eliminated by interaction with immune cells and the number of immune cells activated are drawn from binomial distributions Bðn n ; P n e Þ and Bðn i ; P i a Þ, respectively.Furthermore, the activated immune cells are removed naturally and this process is simulated from the binomial distribution Bðn i ; P i d Þ.In an analogous way to tumor cells migration, the number of migrating immune cells is drawn from the binomial distribution Bðn i ; P i m Þ using the same algorithm.
Different probabilistic events where included to develop this part of the DSBMS.First, the probabilities of death of damaged cells are given by Damaged tumor cells die by mitotic catastrophe after k cycles of mitosis while trying to repair the damage caused [45].The probability of this event (Eq (9a)) is modulated by the relationship between the time step, the division time rate t h r and the number k of the mitosis cycle.Similarly, damaged healthy cells die by mitotic catastrophe while trying to renew themselves.The time frame of this death event is longer since, as stated above, healthy cells have a low proliferation rate.For this reason, P h d described by Eq (9b) has a time dependence and a similar structure to the standard logistic function.This probability increases when damaged healthy cells are closer to their k-th division.Here t represents the time elapsed from radiosurgery to the evaluation step and 1/2σ is the compression parameter.
It was assumed that the probability of immune cells activation to be given by The immune system activation event takes into account the proportion of necrotic and immune cells within each voxel, and the characteristic activation time when there is at least one immune cell inside.P i a decreases with the voxel occupancy, simulating competition for space and resources.It was assumed that the size of immune cells is different from other cell populations, being q times the size of a healthy or tumor cell.
The probability of necrotic cells elimination is described in the DSBMS This probability was modulated by the ratio between the two populations, inversely to the activation process.In this case, P n e is also affected by the carrying capacity of the voxel.
Finally, the probabilities of immune cells death and migration were assume to have the form These probabilities were modulated by the relationship between the time step and the characteristic time of the process.Additionally, migration events were affected by voxel occupancy, making it harder to move in crowded contexts.The notations introduced to define each population of cells are summarized in Table 2.
In all the previous formulations, the parameters τ included in each probabilistic event represent the characteristic times of each process (see Table 3).
Estimation of parameters for the DSBMS.We choose domain simulation sizes in line with those found in the clinical setting for metastatic lesions pre-SRS treatment, which are around 0.5 − 2 cm 3 with maximum tumor sizes up to 10 cm 3 .Since high-resolution MRI voxel size is around 1 mm 3 , we chose that voxel size 1 mm 3 .Hence, space is discretized as a hexahedral mesh consisting of L × L × L spatial units (voxels), where L = 60.The time step was fixed to 4t = 4 hours.From typical cell sizes [46], the carrying capacity of a single voxel was estimated to be K = 2 × 10 5 cells, assuming similar sizes for tumor and healthy cells (* 10-12 μm).Basal rates of tumor cell division, death and migration were chosen using Bayesian criteria based on the doubling times estimations [47] and imaging data from real BMs in [28].The value k = 2 was taken, i.e. damaged tumor cells dying by mitotic catastrophe after two cell cycles on average [45].
Parameters related to the immune system were obtained from data of the microglial cell population [48,49], which constitutes 0.5%-16.6% of the total number of cells in the human brain and the most abundant type of immune cell in the brain.Then, the initial number of immune cells was set at 10% of the healthy cells surrounding the tumor [50].Immune cells are activated if within voxels are cells damaged or necrotic by SRS.Furthermore, it is known that the immune cells increase in size by 50% (q = 3/2) after activation [51].It was assumed that the mean lifetime of immune cells is around two months and the activation time is in the range of 12 to 20 hours [50].
All the proposed parameters are associated with cellular processes, which combined result in whole-tumor rates.Cellular traits were randomly sampled from the range of allowed basal rates for each simulation.While the division time for a single cell may seem extensive, it's important to note that in our model, the utilized division time represents an average across all tumor cells.Notably, not all cells within a tumor undergo division simultaneously.Taking this into account, the term 'division rate' actually refers to the 'averaged division rate'.This concept extends to parameters such as death time and migration coefficient.The parameters were calibrated using the ABC rejection algorithm.Several simulations were performed using input parameters sampled from prior distributions.Simulations that satisfied realistic brain metastases features, such as size and growth speed, were accepted, and the input parameters that produced them were retrieved to build refined posterior distributions.This process was repeated iteratively until a convergence criterion was met.
Virtual BMs simulations.To simulate tumor growth dynamics after SRS, a set of simulations of BMs was run, starting from 10 3 tumor cells, allowing them to grow until they reached diagnostic volumes in the range 0.5 − 2 cm 3 .Radiosurgery events were simulated and posttreatment tumor evolution continued as previously described.Each simulation had a different set of basal rates, sampled randomly from the ranges specified in Table 3.This study focused its attention on cases with either late inflammation or progression, but always in cases with initial response.Thus, simulations displaying longitudinal volumetric growth in the first four months after SRS were rejected.
To simulate the effect of radiosurgery, the voxel survival fraction was defined where S^f is the maximum survival fraction, n n is the number of necrotic cells in a voxel and K is the carriyng capacity.S^f ranges from 0 to 1, with values approaching zero indicating a heightened susceptibility to damage from radiosurgery, while values close to 1 result negligible damage.In this manuscript we explored the whole range of values [0,1] for this parameter in order to consider both the case of radiosensitive and radioresistant tumors.The well-oxygenated cells are less resistant to radiation (more radio sensitive).Thus, cells that are farthest away and that do not get enough oxygen and nutrients to survive are those that are found in the voxels with the highest number of necrotic cells.Two sets of simulations.were carried out.A first group of 200 BM simulations was performed under the ideal condition of no damage to the healthy tissue surrounding the tumor (S n = 1).This group would correspond to an idealized scenario, that is highly unlikely to happen in clinical practice.The second group of 200 BM simulations accounted for the damage induced by SRS to healthy tissue next to the lesion, between 30% and 90% cell death per voxel due to lethal damage (0.1 � S n � 0.7).The latter corresponded to the scenario that would be expected to happen in the clinical setting.Besides, to establish the criterion for distinguishing between recurrence + inflammation (R & I) and the inflammation group (I) within this second set of simulations, it was necessary to set a limit for the survival fraction of tumor cells.We assumed that below 10% maximum survival of tumor cells, the volumetric growth of the lesion would be primarily associated with inflammation and would thus be classified as inflammation, while above 10% survival, it would be classified as relapse + inflammation.
To calculate the lesion volume, the set of voxels reaching more than 45% of their carrying capacity was considered, looking at all cell types including necrotic cells.This set was denoted by ν and N ν = |ν| was calculated as the number of elements in ν.Since the volume of a voxel was taken to be equal to 1 mm 3 , the tumor volume is equal to the number of occupied voxels N ν in mm 3 .Blood vessels start suffering damage when the cellular density is high enough, thus leading to the gadolinium contrast being released into the brain tissue.This leads to a positive signal in the postcontrast T1-weighted images that is then visualized as being tumor tissue.Thus, the inclusion of ν intends to reproduce this real-world feature and make our simulations comparable with available MRI studies.This is a way to determine which voxels should be included in the calculation of the total volume.

Calculation of β exponent for virtual BMs.
To study the dynamics of post-treatment volumetric growth for simulated tumors, the exponent β was calculated from the volumetric simulation data using Eq (7).
To perform the calculation of β exponent for a virtual tumor, attention was paid to the dynamic behavior of tumors after the second follow-up (six months after radiosurgery).Taking into account that the average time between follow-ups in the clinic is three months, for each case, we calculated the following three volume measurements were calculated from six months post-SRS, that is, t = 6, 9 and 12 months.After that, the volumes were fitted to the growth law and the value of the exponent β was calculated.
Finally, we obtained the estimated b value was obtained for the corresponding simulated tumor as the median value of all the calculated values.
Computational implementation.The DSBMS was coded both in Matlab (R2020b, Math-Works, Inc., Natick, MA, USA) and Julia (version 1.5.3).The main workspace and simulation sections were coded in Julia, while data analysis and plotting were coded in Matlab.Simulations were performed on a 4-core 16 GB 2.7 GHz MacBook Pro.Computational cost per simulation was of the order of 1-3 minutes.

Human data
The human data utilized in this study were sourced from the dataset provided in [43].A dataset of BMs collected as part of a retrospective, multicenter, nonrandomized study, treated with stereotactic radiotherapy and followed up with magnetic resonance imaging (MRI).Twenty cases with at least 3 consecutive imaging studies available, and classified as radiation necrosis, were used, of which 6 correspond to breast primary tumor, 13 to non-small cell lung cancer, and 1 to small cell lung cancer.When comparing the three subgroups, no differences were found in the growth exponents β.The volumetric data used in this paper was calculated from volumetric contrast-enhanced T1-weighted MRI sequences, and the main interest was on total volume measurements.

Compartmental mathematical model reproduces the accelerated growth dynamics of RN events
The compartmental model provides information on the evolution of cell populations and fits patient data.To simulate the RN event, we conducted simulations using the mathematical model described in Eq (1).This model allows us to track the evolution of each type of cell over time, providing insight into how the populations of different cell types change during the development of RN.Fig 4A shows the results of one of those simulations.Initially, only tumor cells are present, and they grow exponentially.Meanwhile, necrotic cells steadily increase as damaged healthy cells die when they attempt to renew.Next, it was observed that the primary drivers of the fast-growing phenomenon associated with RN were the immune cell populations.These immune cells contributed the most to the total number of cells and played a critical role in the RN process.
Fittings were also performed of the solutions of Eq (1) to actual data from patients diagnosed with RN.This led to finding individual parameter values providing excellent fits to the data.Fig 4B shows five examples of these fitted lesions.(see Table 4 for the particular values for each fit).
Analysis of model behaviour reveals high growth exponents in terms of volume growth.To further explore the behavior of the mathematical model, we conducted simulations with different initial volumes ranging from 0.5 to 3.0 cm 3 .Specifically, we simulated 500 RN events per volume using a randomized approach where the model parameters took values within a given range.The exponents obtained in these simulations are shown in Fig 4C The findings revealed that the median β values obtained from the mathematical model are consistently larger than 1.5, regardless of the initial volumes or model parameters used in the simulations.

Stochastic mesoscopic model reflects the observed dynamics in BM recurrences and RN events
The stochastic mesoscale brain metastasis growth simulator was set up using the biological rules described in 'Methods'.A total of 400 simulations were performed of virtual BMs treated with SRS displaying growth after treatment and the exponent β was calculated.In each case, 20 β growth exponents were calculated as explained in 'Methods'.Additionally, the median b was obtained for each simulation.In two of the cases, sublinear growths ( b < 1) were obtained for recurrences.These simulations were generated with little or no damage to healthy tissue, i.e, S n = 1 (Fig 6A and 6B) and S n = 0.7 (Fig 6E and 6F).On the other hand, when there was substantial damage to healthy tissue S n = 0.1, the volumetric evolution displayed a superlinear growth ( b > 1), as shown in Fig 6C and 6D.
Post-SRS BMs are characterized throught growth exponents.Thus, in order to characterize the volumetric growth post-SRS of the BMs using the scaling exponent, the values of β were calculated for the set of 400 virtual BMs.The results obtained for the first group (S n = 1) are shown in Fig 7A .The values of β were grouped according to the value of Sf used in the simulation.Medians and quartiles of the box plots were mostly below 1, although for small values of Sf , there were a set of outliers with high estimates of β.These exponents described the dynamics of relapsing lesions.
For the second group, which included damage to healthy tissue, there were two behaviors observed in silico.Growth exponents allow differentiation between radiation necrosis and recurrence.The computational results suggest that the β value could be used to distinguish the inflammatory response from tumor progression.The ANOVA test for the comparison with virtual BMs leads to significant differences between the inflammatory response group and relapses groups (p = 1.85×10 −12  We obtained AUC = 0.97 and the optimal threshold calculated to maximize the sensitivity and specificity values was β threshold = 1.05.This means that inflammatory events show faster growth dynamics than relapses.Furthermore, the sensitivity obtained at β threshold was equal to 0.96, indicating its ability to accurately identify true positive cases.
We investigated how the AUC changes based on the chosen S f bound that defines the R+I and I groups.However, it is important to observe that as simulated BMs with a greater proportion of tumor cell survival after SRS are incorporated into group I, the resulting β* threshold tends to be less than 1, what is related with the fact that in presence of tumor cells growing, the corresponding β is smaller than 1.

Discussion and conclusion
Radiation necrosis is a common adverse effect associated with radiation treatment.To gain a deeper understanding of the growth dynamics of RN events, computational frameworks were developed, based on mathematical models of increasing complexity.These models provide mechanistic explanations for the observed growth dynamics of RN, and offer insight into the underlying mechanisms of this phenomenon.
Recent research has revealed that RNs tend to exhibit faster growth dynamics than recurrent BMs, which can help distinguish between the two entities [38].However, until now, no mechanistic explanation has been provided for these observations.The mathematical models developed in this work offer a promising way to fill this gap and provide a mechanistic understanding of the differences in growth dynamics between RN and BM.
Although RN is thought to be the result of synergy between vascular damage and immune activation, the precise mechanism underlying this phenomenon remains unclear [52,53].The first mathematical model presented in this study was a compartmental model.The hypothesis was that healthy tissue surrounding the lesion becomes damaged due to radiation, leading to necrosis and the activation of immune cells.This mechanism effectively reproduces the processes that take place during RN events.Additionally, the simulations consistently yield growth exponents larger than 1.5, regardless of the initial values.Even though the model might appear to be overfitted, in fact the model is basically a two-compartment model, since the equation for tumor cells is not coupled with the other two and the parameter is calculated from the total volume of the first two time points.
However, there are limitations to this model.While it provides a basic understanding of the mechanisms involved in RN, it may be too simplistic for a more detailed understanding of the biological processes involved in the response to SRS.It is important to have a more comprehensive understanding of these processes, not only for RN but also for recurrent events.Differentiating RN from tumor progression is one of the main challenges in the management of BMs, since the treatment approach for each of these varies greatly, and accurate differentiation is essential for effective treatment [54].An additional complication arises from the fact that RN is usually heterogeneous, where necrotic cells coexist with tumor cells [55].This complicates the differentiation process and highlights the need for a more detailed understanding of the underlying biological processes.
To the authors' knowledge, only two studies have implemented mathematical models to date to differentiate RN from progression.Both studies used a biomechanically coupled tumor growth model that describes a logistic growth with a diffusive term in the presence of an external force.To estimate the model parameters, the first study [33] used ten lesions, with five exhibiting tumor recurrence and five showing RN.The second study [34], based on the same model, investigated 78 BMs and achieved an area under the ROC curve of 0.75 when differentiating progression from RN.In contrast, our DSBMS model achieved a superior ROC of 0.97, using the growth exponent as the classification parameter.However, it is crucial to note that these results are not directly comparable due to the utilization of synthetic data in our study.
To account for the complexities of the system, a stochastic mesoscale simulator was also developed, providing a more detailed representation of the complex interactions occurring within BMs.In this simulator, BMs are described as a composition of different clonal populations seen at the voxel level, with three populations taken into account before treatment: healthy cells, tumor cells, and necrotic cells, and three more to account for the response after SRS: immune cells, damaged tumor cells, and damaged healthy cells.
Brain metastases growth was described on a mesoscopic scale in the simplest possible way, taking into account only one clonal population.Although it is known that these tumors present cellular heterogeneity in their composition, it is not fully characterized.For this reason, our results are not focused on describing the growth dynamics of BMs optimally.On the other hand, it has been known that tumor treatments induce a reduction of the clonal complexity at the point of maximal response due to the action of selective pressures of the drugs [56], which supports the simplification of a single clonal population to describe the response to SRS.
Continuous models can potentially neglect spatial correlations between the locations of individuals, especially when different species or subpopulations are taken into account.Therefore, the model developed here accounts for these relationships by combining discrete, spatial, and stochastic dynamics.This made it possible to analyze the dynamic behavior of the tumor in terms of volume through a more complete and precise approach.The model was used to simulate the different possible scenarios according to the damage caused to the healthy tissue surrounding the tumor, and the survival fraction of the tumor cells after therapy.Depending on whether the healthy tissue surrounding the lesion is damaged or not, this leads to different values of the growth exponents.
The discrete model explored the coexistence of necrotic cells with tumor cells, resulting in heterogeneous RNs.However, the interaction between the tumor and immune cell population was not taken into account, although it has been hypothesized that radiation induces damage signals within the tumor, making it for the immune system to detect [57].There is only limited data available on the dynamics of immune response after radiation to a tumor lesion, but it could be interesting to include this in future work.Despite this, an increase in the necrotic cell population in the stochastic simulator indirectly affects the tumor cell population by saturating the voxels in question.In these cases, the probabilities of tumor cell proliferation and death are decreased and increased, respectively, and slower growth dynamics of this population are observed.
The stochastic simulator was tested to determine its effectiveness in distinguishing between recurrent lesions and RN events.The model classified both entities at 97% accuracy (AUC = 0.97), which demonstrated its potential to differentiate between these two conditions.The threshold exponent β found, despite being larger than 1, was smaller than those computed by the continuous model and that observed in human beings.This can be explained by the inclusion of virtual simulations with small but not zero residual diseases (S f > 0) after SRS in the inflammation group, which may affect overall volumetric growth dynamics.Another drawback is the low spatial resolution when cell densities decrease in response to therapy, which may affect the lesion volume measurement.In the case of this work, the resolution scale of the MRI data was maintained to facilitate integration with image data and allow computationally feasible simulations up to the macroscopic whole-tumor scale.However, enhancing the resolution may improve patient classification.This is a consideration for the future, especially as technologies continue to advance, providing greater precision in magnetic resonance imaging.
However, these results also suggest that the value of the exponent β could have a direct clinical application.The model supports that when three MRI studies are available satisfying the inclusion criteria, computing β for a particular patient could allow us to diagnose the growth as inflammation or tumor progression.Furthermore, this result is consistent with the study carried out in [38] on a group of real patients with and without radiation necrosis.In the present study, the area under the curve (AUC) associated with the receiver operating characteristic (ROC) curve was found to be 0.97, which is significantly higher than the AUC of 0.74 reported in the previously cited study that used actual patient data.One potential explanation for this difference could be the presence of confounding effects and a higher variability between patients in the actual patient data.In a previous study [37], it was emphasized that while segmentation influenced the computation of the growth exponent β in certain cases, especially for the smallest lesions, the majority of lesions demonstrated a consistent β even in the presence of added variability.
In conclusion, this paper has presented two models to explain the phenomenon of RN and to differentiate it from tumor progression.The first model, a simplistic continuous model, provides a plausible explanation for the high growth exponents observed in human beings during RN events.The second model, a stochastic mesoscale simulator, provided a more detailed representation of the complex interactions occurring within BMs, and explains the growth of BMs with varying degrees of damage.The simulator was tested to distinguish between recurrent lesions and RN events and proved to be a robust theoretical support to the use of the growth exponent β in differentiating between these two conditions.The development of accurate and detailed models could help improve the management of BMs and ultimately improve patient outcomes.

Fig 1 .Fig 2 .
Fig 1. Longitudinal dynamics of the 3D reconstructions of a brain metastasis treated with SRS showing radiation necrosis.Dots correspond to the measured volumes together with their 3D reconstructions and the solid line is the result of interpolating longitudinal volumetric data (shown only to guide the eye).The patient was a 50-year-old male with a non-small-cell lung cancer primary, who underwent a single session SRS with a dose of 20 Gy.In this case, the inflammatory lesion exhibited its peak volumetric expansion around 12 months after treatment.The solid line is a cubic spline interpolation shown to guide the eye.https://doi.org/10.1371/journal.pcbi.1011400.g001

Fig 3 .
Fig 3. A Slice of a tumor simulation before SRS with DSBMS, mapping the distribution of cell density.A representation of the populations of healthy (blue), tumor (red) and necrotic (black) cells within a voxel according to their location is shown.B. Example of a single-shot treatment plan for a virtual simulation of SRS.The target is outlined in yellow, and it is the area most affected by SRS.The green line encloses another area affected with less intensity.C. Spatial distribution of cell populations just before and after SRS.Voxels may be occupied by more than one cell population but the dominant populations per voxel are shown.https://doi.org/10.1371/journal.pcbi.1011400.g003 Fig 5 presents two simulation examples, capturing distinct scenarios of recurrence and radiation necrosis (RN) events, offering insights into their volumetric changes over time and the distribution of the diverse cell populations.DSBMS correctly reproduces the volumetric dynamics of BM growth after SRS.The volumetric growth dynamics after therapy of the virtual BMs generated was studied.Fig 6 shows three examples of these in silico simulations.The first column (A, C, E) shows the dynamics of the different cell populations present: proliferating tumor cells, damaged cells, necrotic cells, immune cells and total tumor cells for the three cases.The second column (B, D, F) shows the longitudinal volumetric dynamics of the simulation displayed in the first column.
Fig 7B shows the scatter plot of the b median calculated for the virtual BMs

Fig 4 .
Fig 4. Results of simulations of the compartmental model Eq (1).a. Example of a simulation of model Eq (1)showing the typical dynamics observed post-SRS in the context of RN events.The simulation illustrates the evolution of each cell type, with the initial number of tumor cells set to N T = 5 � 10 7 , and parameter values ρ = 0.07 days −1 , λ N = 2.3 � 10 −11 days −1 , γ = 1.9 � 10 −7 days −1 and θ = 0.17 days −1 .b. Fittings using Eq (1) (solid lines) for the longitudinal volumetric growth data (circles) for five patients diagnosed with RN.The parameters used for the fitting process are shown in Table4.c. Distribution of growth exponents, β values obtained for simulations with varying initial volumes ranging from 0.5 to 3.0 cm 3 by simulating 500 RN events per volume using a randomized approach where model parameters took values in the predefined range.Note that values of β between 1 and 2 were obtained systematically.https://doi.org/10.1371/journal.pcbi.1011400.g004

Fig 5 .Fig 6 .
Fig 5.The dynamics of brain metastases (BMs) recurrences and radiation necrosis (RN) events are captured by the stochastic mesoscopic model.a. Recurrent BM. b.RN event.White points represent the total tumor volume at different time points, and black lines depict interpolations between the points (provided as visual guidance).Additionally, spatial distribution of the four cell populations illustrate the dynamic shifts in dominance among cell types for both scenarios.The dominant population per voxel is shown.https://doi.org/10.1371/journal.pcbi.1011400.g005 https://doi.org/10.1371/journal.pcbi.1011400.g006according to the different values of (Sf , S n ) simulated.Values of b > 1 were obtained for cases where SRS eliminates most of the tumor cells (values of S^f � 0:1).Here, the volume re-growth was due to the inflammatory component.Otherwise, for larger tumor remnants re-growth simulations (values of 0:1 < Sf < 1) the calculated b exponents were less than 1.This behavior corresponds to tumor relapse.In addition, the β values calculated for the cases simulated with 0.1 � S n � 0.7 and Sf � 0:1 can be seen in Fig 7C.Despite the observed variability, the β values obtained were typically greater than 1.
). Box plots for the different subgroups are shown in Fig 8A The area under the ROC curve (AUC) in Fig 8B illustrates the ability of the exponent b to discriminate between response groups.

Fig 7 .
Fig 7. Comparison of box plots for the growth exponents β calculated for virtual BMs performed with the DSBMS.A. Growth exponents β values computed for the group of 200 simulations with S n = 1 and 0 � Sf � 1. B. Growth exponents β values computed for the group of 200 simulations with 0.1 � S n � 0.7 and 0 � Sf � 0:1.C.Scatter plot that shows the β median calculated for the virtual BMs which were simulated with different values of (Sf , S n ).Basal rate parameters for the simulations are shown in Table3.https://doi.org/10.1371/journal.pcbi.1011400.g007 Fig 9 displays the outcomes from this exploration, ranging the S f bound from 0.04 to 0.4.Additionally, the optimal value of the β* threshold, obtained from the ROC curve, is indicated for each case.As we can see, the classification power of β slightly decreases as the value of S f increases.The result provide insights into the robustness of β as a classifier.

Table 2 . Notation and description of the cell populations in the DSBMS model.
i numbers of immune cells inside a voxel n hd numbers of damaged healthy cells inside a voxel https://doi.org/10.1371/journal.pcbi.1011400.t002