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

Abstract

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.

Citation:Ocaña-Tienda B, León-Triana O, Pérez-Beteta J, Jiménez-Sánchez J, Pérez-García VM (2024) Radiation necrosis after radiation therapy treatment of brain metastases: A computational approach. PLoS Comput Biol 20(1): e1011400. https://doi.org/10.1371/journal.pcbi.1011400

Editor:James Gallo, University at Buffalo - The State University of New York, UNITED STATES

Received:July 31, 2023; Accepted:January 21, 2024; Published: January 30, 2024

Copyright: © 2024 Ocaña-Tienda et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability:The original code for the mesoscopic model has been made available in the repository https://github.com/JuanJS117/MesoscopicModel.

Funding:This research has been supported by the Spanish Ministerio de Ciencia e Innovación (grants PID2019-110895RB-100 and PDC2022-133520-I00), Junta de Comunidades de Castilla-La Mancha (grant SBPLY/21/180501/000145) to VMPG and BOT is supported by the Spanish Ministerio de Ciencia e Innovación (grant PRE2020-092178). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

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].

Mechanistic mathematical models have been used to study the growth of untreated metastasis in different scenarios [21–23], the interaction between primary and secondary tumors [24], cancer metastasis networks [25–27], and, and, in recent years, the response to treatment [28–34]. However only a few of them have studied specifically brain metastases [23, 28, 29, 33, 34].

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.

[Figure omitted. See PDF.]

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

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.

[Figure omitted. See PDF.]

Fig 2. Schematic representation of the populations involved in the response to SRS.

The post-treatment scenario represents the growth of radiation necrosis events as defined in Eq (1). After SRS healthy cells around the tumor become damaged and die by necrotic catastrophe. The appearance of necrosis stimulates immune cells leading to an apparent growth of the lesion.

https://doi.org/10.1371/journal.pcbi.1011400.g002

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:(1a) (1b) (1c)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(2)which accounts for the fact that cells have the same probability of reaching the renewal point during a given time period tN, 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(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 ≈ HT/60 by assuming that an equal number of cells will die each day during the 2-month period, tN (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.(4)

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.

[Figure omitted. See PDF.]

Table 1. Summary of parameters for the compartmental model.

https://doi.org/10.1371/journal.pcbi.1011400.t001

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 t0, t1, t2 and the computed volumes V0, V1, V2 were used to compute the growth exponent β using the simplified Von Bertalanffy growth equation [35](5)

Solving Eq (5), gives(6)

Since there is information on the dynamics at three-time points (t0, V0), (t1, V1) and (t2, V2), the two parameters α and β can be completely determined by evaluating (6) at the times t0, t1, t2, giving(7)

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 V0, V1, V2, t0, t1, t2.

A dataset of 500 virtual BMs was generated for each choice of initial volume, ranging from 0.5 cm3 to 3.0 cm3. 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 △x, 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 △x. 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 nt, nn and nh 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 and , respectively. The number of migrating tumor cells is drawn from the respective binomial distribution . Then these cells are distributed around a neighborhood of 26 voxels (Moore neighborhood) according to a multinomial distribution Y ∼ Mult(X, PMoore), 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 PMoore is a vector of 26 components giving the probabilities of migration to each surrounding voxel [44]. The probability PMoore is proportional to 1, or 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 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.

[Figure omitted. See PDF.]

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

To construct tumor growth rules in the DSBMS the probabilities of tumor cells reproduction, migration and death were first considered to be given by(8a) (8b) (8c)

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 by the hyperbolic tangent term as a function of the voxel capacity.

As to motility, the probability of healty cells migration was given by(8d)

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 Sf surviving without damage and remaining viable, a fraction (1 − Sf) 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 − Sn) of healthy cells surrounding the tumor, where Sn represents the survival fraction of healthy cells remaining undamaged and viable. Fig 3C shows an example of the spatial distribution of cell population before and after SRS, as described above.

Let nd, ni and nhd 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 and , 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 and , respectively. Furthermore, the activated immune cells are removed naturally and this process is simulated from the binomial distribution . In an analogous way to tumor cells migration, the number of migrating immune cells is drawn from the binomial distribution 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(9a) (9b)

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 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, 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(9c)

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. 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(9d)

This probability was modulated by the ratio between the two populations, inversely to the activation process. In this case, 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(9e) (9f)

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.

[Figure omitted. See PDF.]

Table 2. Notation and description of the cell populations in the DSBMS model.

https://doi.org/10.1371/journal.pcbi.1011400.t002

In all the previous formulations, the parameters τ included in each probabilistic event represent the characteristic times of each process (see Table 3).

[Figure omitted. See PDF.]

Table 3. Summary of parameter values used for the stochastic model.

https://doi.org/10.1371/journal.pcbi.1011400.t003

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 cm3 with maximum tumor sizes up to 10 cm3. Since high-resolution MRI voxel size is around 1 mm3, we chose that voxel size 1 mm3. 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 △t = 4 hours. From typical cell sizes [46], the carrying capacity of a single voxel was estimated to be K = 2 × 105 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 103 tumor cells, allowing them to grow until they reached diagnostic volumes in the range 0.5 − 2 cm3. Radiosurgery events were simulated and post-treatment 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(10)where is the maximum survival fraction, nn is the number of necrotic cells in a voxel and K is the carriyng capacity. 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 (Sn = 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 ≤ Sn ≤ 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 mm3, the tumor volume is equal to the number of occupied voxels Nν in mm3. 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.

The three-time instances were taken within the following ranges: t0 ∈ (180, 180+ 15) days, t1 ∈ (t0+ 80, t0+ 100) days and t2 ∈ (t1+ 80, t1+ 100) days. Then, β was estimated for each simulation for the volume and combination of t0, t1 and t2. This procedure was repeated 20 times. Finally, we obtained the estimated 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, MathWorks, 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.

Results

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.

[Figure omitted. See PDF.]

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 NT = 5 ⋅ 107, 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 Table 4. c. Distribution of growth exponents, β values obtained for simulations with varying initial volumes ranging from 0.5 to 3.0 cm3 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

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).

[Figure omitted. See PDF.]

Table 4. Model Eq (1) parameters best fitting the longitudinal dynamics of RN events. Volumetric data were taken from Ref. [43] and lesions 1 to 5 correspond to patients ID numbers: 30031, 40024, 40001, 40176, and 40042, respectively.

https://doi.org/10.1371/journal.pcbi.1011400.t004

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 cm3. 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. 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.

[Figure omitted. See PDF.]

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

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. In each case, 20 β growth exponents were calculated as explained in ‘Methods’. Additionally, the median was obtained for each simulation. In two of the cases, sublinear growths () were obtained for recurrences. These simulations were generated with little or no damage to healthy tissue, i.e, Sn = 1 (Fig 6A and 6B) and Sn = 0.7 (Fig 6E and 6F). On the other hand, when there was substantial damage to healthy tissue Sn = 0.1, the volumetric evolution displayed a superlinear growth (), as shown in Fig 6C and 6D.

[Figure omitted. See PDF.]

Fig 6. DSBMS simulations of longitudinal tumor growth dynamics after SRS.

The first column (A,C,E) shows the dynamics of the total populations of proliferating cells (red lines), damaged cells (orange lines), necrotic cells (violet lines), immune cells (green lines) and total tumor cells (blue lines). The second column (B,D,F) shows the longitudinal tumor volumetric dynamics. Blue lines inside the zoom box represent different fits of β exponent solving Eq (7) and the mean is shown. Subplots (A-B) correspond with a tumor simulation with no damage to healthy tissue (Sn = 1, ), subplots (E-F) correspond with small damage to healthy tissue (Sn = 0.7, ) and subplots (C-D) correspond with high damage to healthy tissue (Sn = 0.1, ). Basal rate parameters for simulations in this figure are h, h, h, h, h, h, h and h.

https://doi.org/10.1371/journal.pcbi.1011400.g006

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 (Sn = 1) are shown in Fig 7A. The values of β were grouped according to the value of used in the simulation. Medians and quartiles of the box plots were mostly below 1, although for small values of , there were a set of outliers with high estimates of β. These exponents described the dynamics of relapsing lesions.

[Figure omitted. See PDF.]

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 Sn = 1 and . B. Growth exponents β values computed for the group of 200 simulations with 0.1 ≤ Sn ≤ 0.7 and . C. Scatter plot that shows the β median calculated for the virtual BMs which were simulated with different values of (, Sn). Basal rate parameters for the simulations are shown in Table 3.

https://doi.org/10.1371/journal.pcbi.1011400.g007

For the second group, which included damage to healthy tissue, there were two behaviors observed in silico. Fig 7B shows the scatter plot of the median calculated for the virtual BMs according to the different values of (, Sn) simulated. Values of were obtained for cases where SRS eliminates most of the tumor cells (values of ). Here, the volume re-growth was due to the inflammatory component. Otherwise, for larger tumor remnants re-growth simulations (values of ) the calculated exponents were less than 1. This behavior corresponds to tumor relapse. In addition, the β values calculated for the cases simulated with 0.1 ≤ Sn ≤ 0.7 and can be seen in Fig 7C. Despite the observed variability, the β values obtained were typically greater than 1.

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). 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 to discriminate between response groups. 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.

[Figure omitted. See PDF.]

Fig 8.

A. Box plots showing the comparison of the growth exponents between the different simulated BMs: relapse group (R), whose response is characterized by tumor progression (, Sn = 1), relapse and inflammation group (R & I), whose response is characterized by tumor progression and inflammation (, 0.1 ≤ Sn ≤ 0.7) and inflammation group (I), whose response is characterized by inflammation (, 0.1 ≤ Sn ≤ 0.7). B. ROC curve for the discrimination between tumor progression (R and R&I groups) and inflammatory response (I group) according to the growth exponent .

https://doi.org/10.1371/journal.pcbi.1011400.g008

We investigated how the AUC changes based on the chosen Sf bound that defines the R+I and I groups. Fig 9 displays the outcomes from this exploration, ranging the Sf 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 Sf increases. The result provide insights into the robustness of β as a classifier. 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.

[Figure omitted. See PDF.]

Fig 9. Exploring the correlation between the AUC of the ROC curve and ranging the Sf bound from 0.04 to 0.4.

The colorbar illustrates the corresponding calculated β* threshold.

https://doi.org/10.1371/journal.pcbi.1011400.g009

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 (Sf > 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.

References

1. 1. Kaal EC, Niel CG, Vecht CJ. Therapeutic management of brain metastasis. Lancet Neurology. 2005; 4(5): 289–298. pmid:15847842

* View Article

* PubMed/NCBI

* Google Scholar

2. 2. Lamba N, Wen PY, Aizer AA. Epidemiology of brain metastases and leptomeningeal disease. Neuro Oncol. 2021; 23(9):1447–1456. pmid:33908612

* View Article

* PubMed/NCBI

* Google Scholar

3. 3. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer statistics, 2022. CA Cancer J Clin. 2022; 72(1): 7–33. pmid:35020204

* View Article

* PubMed/NCBI

* Google Scholar

4. 4. Suh JH, Kotecha R, Chao ST, Ahluwalia MS, Sahgal A, Chang EL. Current approaches to the management of brain metastases. Nat Rev Clin Oncol. 2020; 17(5): 279–299. pmid:32080373

* View Article

* PubMed/NCBI

* Google Scholar

5. 5. Kohutek ZA, Yamada Y, Chan TA, Brennan CW, Tabar V, Gutin PH, et al. Long-term risk of radionecrosis and imaging changes after stereotactic radiosurgery for brain metastases Journal of Neuro-Oncology. 2015; 125(1): 149–156. pmid:26307446

* View Article

* PubMed/NCBI

* Google Scholar

6. 6. Kolesnick R, Fuks Z. Radiation and ceramide-induced apoptosis. Oncogene. 2003; 22(37):5897–5906. pmid:12947396

* View Article

* PubMed/NCBI

* Google Scholar

7. 7. Abdulla S, Saada J, Johnson G, Jefferies S, Ajithkumar T. Tumour progression or pseudoprogression? A review of post-treatment radiological appearances of glioblastoma. Clin Radiol. 2015; 70(11):1299–1312. pmid:26272530

* View Article

* PubMed/NCBI

* Google Scholar

8. 8. Peña LA, Fuks Z, Kolesnick RN. Radiation-induced apoptosis of endothelial cells in the murine central nervous system: protection by fibroblast growth factor and sphingomyelinase deficiency. Cancer Res. 2000; 60(2):321–327. pmid:10667583

* View Article

* PubMed/NCBI

* Google Scholar

9. 9. Baker DG, Krochak RJ. The response of the microvascular system to radiation: a review. Cancer Invest. 1989; 7(3):287–294. pmid:2477120

* View Article

* PubMed/NCBI

* Google Scholar

10. 10. Fike JR, Rosi S, Limoli CL. Neural precursor cells and central nervous system radiation sensitivity. Semin Radiat Oncol. 2009; 19(2):122–132. pmid:19249650

* View Article

* PubMed/NCBI

* Google Scholar

11. 11. Kudo S, Suzuki Y, Noda SE, Mizui T, Shirai K, Okamoto M, et al. Comparison of the radiosensitivities of neurons and glial cells derived from the same rat brain. Exp Ther Med. 2014; 8(3):754–758. pmid:25120594

* View Article

* PubMed/NCBI

* Google Scholar

12. 12. Nordal RA, Nagy A, Pintilie M, Wong CS. Hypoxia and hypoxia-inducible factor-1 target genes in central nervous system radiation injury: a role for vascular endothelial growth factor. Clin Cancer Res. 2004; 10(10):3342–3353. pmid:15161688

* View Article

* PubMed/NCBI

* Google Scholar

13. 13. Sachet M, Liang YY, Oehler R. The immune response to secondary necrotic cells. Apoptosis. 2017; 22(10):1189–1204. pmid:28861714

* View Article

* PubMed/NCBI

* Google Scholar

14. 14. Proskuryakov SY, Konoplyannikov AG, Gabai VL. Necrosis: a specific form of programmed cell death?. Exp Cell Res. 2003; 283(1):1–16. pmid:12565815

* View Article

* PubMed/NCBI

* Google Scholar

15. 15. Shah R, Vattoth S, Jacob R, Manzil FF, O’Malley JP, Borghei P, et al. Radiation necrosis in the brain: Imaging features and differentiation from tumor recurrence. Radiographics. 2012; 32(5):1343–1359. pmid:22977022

* View Article

* PubMed/NCBI

* Google Scholar

16. 16. Chao ST, Suh JH, Raja S, Lee SY, Barnett G. The Sensitivity and Specificity of FDG PET in Distinguishing Recurrent Brain Tumor from Radionecrosis in Patients Treated with Stereotactic Radiosurgery. Int J Cancer. 2001; 96(3):191–197. pmid:11410888

* View Article

* PubMed/NCBI

* Google Scholar

17. 17. Kim HY, Cho SJ, Sunwoo L, Baik SH, Bae YJ, Choi BS, et al. Classification of true progression after radiotherapy of brain metastasis on MRI using artificial intelligence: a systematic review and meta-analysis. Neurooncol Adv.2021; 3(1):vdab080. pmid:34377988

* View Article

* PubMed/NCBI

* Google Scholar

18. 18. Zhang Z, Yang J, Ho A, Jiang W, Logan J, Wang X, et al. A predictive model for distinguishing radiation necrosis from tumour progression after gamma knife radiosurgery based on radiomic features from MR images. Eur Radiol. 2018; 28(6):2255–2263. pmid:29178031

* View Article

* PubMed/NCBI

* Google Scholar

19. 19. Larroza A, Moratal D, Paredes-Sánchez A, Soria-Olivas E, Chust ML, Arribas LA, Arana E. Support vector machine classification of brain metastasis and radiation necrosis based on texture analysis in MRI. Journal of Magnetic Resonance Imaging. 2015; 42(5):1362–1368. pmid:25865833

* View Article

* PubMed/NCBI

* Google Scholar

20. 20. Lasocki A, Sia J, Stuckey SL. Improving the diagnosis of radiation necrosis after stereotactic radiosurgery to intracranial metastases with conventional MRI features: a case series. Cancer Imaging. 2022; 22(1):33. pmid:35794677

* View Article

* PubMed/NCBI

* Google Scholar

21. 21. Barbolosi D, Benabdallah A, Hubert F, Verga F. Mathematical and numerical analysis for a model of growing metastatic tumors. Math Biosci. 2009; 218(1):1–14. pmid:19121638

* View Article

* PubMed/NCBI

* Google Scholar

22. 22. Hartung N, Mollard S, Barbolosi D, Benabdallah A, Chapuisat G, Henry G, et al. Mathematical modeling of tumor growth and metastatic spreading: Validation in tumor-bearing mice. Cancer Res. 2014; 74(22):6397–6407. pmid:25217520

* View Article

* PubMed/NCBI

* Google Scholar

23. 23. Bilous M, Serdjebi C, Boyer A, Tomasini P, Pouypoudat C, Barbolosi D, et al. Quantitative mathematical modeling of clinical brain metastasis dynamics in non-small cell lung cancer. Sci. Rep. 2019; 9(1):13018. pmid:31506498

* View Article

* PubMed/NCBI

* Google Scholar

24. 24. Álvarez-Arenas A, Souleyreau W, Emanuelli A, Cooley LS, Bernhard JC, Bikfalvi A, Benzekry S. Practical identifiability analysis of a mechanistic model for the time to distant metastatic relapse and its application to renal cell carcinoma. PLoS Comput Biol. 2022; 18(8):e1010444. pmid:36007057

* View Article

* PubMed/NCBI

* Google Scholar

25. 25. Chen LL, Blumm N, Christakis NA, Barabasi AL, Deisboeck TS. Cancer metastasis networks and the prediction of progression patterns. Br J Cancer. 2009; 101(1):749–758. pmid:19707203

* View Article

* PubMed/NCBI

* Google Scholar

26. 26. Scott JG, Basanta D, Anderson AR, Gerlee P. A mathematical model of tumour self-seeding reveals secondary metastatic deposits as drivers of primary tumour growth. J R Soc Interface. 2013; 10(1): 20130011. pmid:23427099

* View Article

* PubMed/NCBI

* Google Scholar

27. 27. Newton PK, Mason J, Bethel K, Bazhenova LA, Nieva J, Kuhn P. A stochastic Markov chain model to describe lung cancer growth and metastasis. PLoS ONE. 2012; 7:e34637. pmid:22558094

* View Article

* PubMed/NCBI

* Google Scholar

28. 28. León-Triana O, Pérez-Beteta J, Albillo D, Ortiz de Mendivil A, Perez-Romasanta L, et al. Brain Metastasis response to stereotactic radio surgery: A mathematical approach. Mathematics. 2021; 9(7), 716.

* View Article

* Google Scholar

29. 29. Watanabe Y, Dahlman EL, Leder KZ, Hui SK. A mathematical model of tumor growth and its response to single irradiation. Theor. Biol. Med. Model. 2016; 136. pmid:26921069

* View Article

* PubMed/NCBI

* Google Scholar

30. 30. Camacho A, Jerez S. Bone metastasis treatment modeling via optimal control. Journal of Mathematical Biology. 2018; 78(1-2):497–526. pmid:30132065

* View Article

* PubMed/NCBI

* Google Scholar

31. 31. Benzekry S, Tracz A, Mastri M, Corbelli R, Barbolosi D, Ebos JM. Modeling Spontaneous Metastasis following Surgery: An In Vivo-In Silico Approach. Cancer Res. 2016; 76(3):535–547. pmid:26511632

* View Article

* PubMed/NCBI

* Google Scholar

32. 32. Barbolosi D, Summer I, Meille C, Serre R, Kelly A, Zerdoud S, et al. Modeling therapeutic response to radioiodine in metastatic thyroid cancer: A proof-of-concept study for individualized medicine. Oncotarget. 2017; 8:39167–39176. pmid:28389624

* View Article

* PubMed/NCBI

* Google Scholar

33. 33. Narasimhan S, Johnson HB, Nickles TM, Miga MI, Rana N, Attia A, et al. Biophysical model-based parameters to classify tumor recurrence from radiation-induced necrosis for brain metastases. Med Phys. 2019; 46(5):2487–2496. pmid:30816555

* View Article

* PubMed/NCBI

* Google Scholar

34. 34. Dohm AE, Nickles TM, Miller CE, Bowers HJ, Miga MI, Attia A, et al. Clinical assessment of a biophysical model for distinguishing tumor progression from radiation necrosis. Med Phys. 2021; 48(7):3852–3859. pmid:34042188

* View Article

* PubMed/NCBI

* Google Scholar

35. 35. Von Bertalanffy L. Quantitative laws in metabolism and growth. Q Rev Biol. 1957; 32(3):217–231. pmid:13485376

* View Article

* PubMed/NCBI

* Google Scholar

36. 36. Pérez-García VM, Calvo GF, Bosque JJ, León-Triana O, Jiménez J, Pérez-Beteta J, Belmonte-Beitia J, et al. Universal scaling laws rule explosive growth in human cancers. Nature Phys. 2020; 16(12):1232–1237. pmid:33329756

* View Article

* PubMed/NCBI

* Google Scholar

37. 37. Ocaña-Tienda B, Pérez-Beteta J, Jiménez-Sánchez J, Molina-García D, Ortiz de Mendivil A, Asenjo B, et al. Growth exponents reflect evolutionary processes and treatment response in Brain Metastases NPJ Syst Biol Appl. 2023; 9(1):35. pmid:37479705

* View Article

* PubMed/NCBI

* Google Scholar

38. 38. Ocaña-Tienda B, Pérez-Beteta J, Molina-García D, Asenjo B, Ortiz de Mendivil A, Albillo D, et al. Growth dynamics of brain metastases differentiate radiation necrosis from recurrence. Neurooncol Adv. 2022; 5(1):vdac179. pmid:36726366

* View Article

* PubMed/NCBI

* Google Scholar

39. 39. Sheline GE, Wara WM, Smith V. Therapeutic irradiation and brain injury. Int J Radiat Oncol Biol Phys. 1980; 6(9):1215–1228. pmid:7007303

* View Article

* PubMed/NCBI

* Google Scholar

40. 40. Mineyeva OA, Bezriadnov DV, Kedrov AV, Lazutkin AA, Anokhin KV, Enikolopov GN. Radiation Induces Distinct Changes in Defined Subpopulations of Neural Stem and Progenitor Cells in the Adult Hippocampus. Front Neurosci. 2019; 12:1013. pmid:30686979

* View Article

* PubMed/NCBI

* Google Scholar

41. 41. Mkango SB, Shaban N, Mureithi E, Ngoma T. Dynamics of breast cancer under different rates of chemoradiotherapy. Comput Math Methods Med. 2019; 2019:5216346. pmid:31611927

* View Article

* PubMed/NCBI

* Google Scholar

42. 42. Villasana M, Radunskaya A. A delay differential equation model for tumor growth. J Math Biol. 2003; 47(3):270–294. pmid:12955460

* View Article

* PubMed/NCBI

* Google Scholar

43. 43. Ocaña-Tienda B, Pérez-Beteta J, Villanueva-García JD, Romero-Rosales JA, Molina-García D, Suter Y, et al. A comprehensive dataset of annotated brain metastasis MR images with clinical and radiomic data. Sci Data. 2023; 10(1):208. pmid:37059722

* View Article

* PubMed/NCBI

* Google Scholar

44. 44. Jiménez-Sánchez J, Martínez-Rubio A, Popov A, Pérez-Beteta J, Azimzade Y, Molina-García D, et al. A mesoscopic simulator to uncover heterogeneity and evolutionary dynamics in tumors. PLoS Comput Biol. 2021 17(2):e1008266. pmid:33566821

* View Article

* PubMed/NCBI

* Google Scholar

45. 45. Joiner MC, van der Kogel A. Basic Clinical Radiobiology (4th ed.). CRC Press. 2009.

46. 46. Milo R, Jorgensen P, Moran U, Weber G., BioNumbers-the database of key numbers in molecular and cell biology. Nucleic Acids Research. Nucleic Acids Res. 2009; 38:D750–D753. pmid:19854939

* View Article

* PubMed/NCBI

* Google Scholar

47. 47. Kobets AJ, Backus R, Fluss R, Lee A, Lasala PA, Evaluating the natural growth rate of metastatic cancer to the brain. Surg Neurol Int. 2020; 11: 254. pmid:33024592

* View Article

* PubMed/NCBI

* Google Scholar

48. 48. He B, Wang X, He Y, Li H, Yang Y, Shi Z, et al. Gamma ray-induced glial activation and neuronal loss occur before the delayed onset of brain necrosis. FASEB J. 2020; 34(10):13361–13375. pmid:32851721

* View Article

* PubMed/NCBI

* Google Scholar

49. 49. Wagner S, Lanfermann H, Wohlgemuth WA, Gufler H. Effects of effective stereotactic radiosurgery for brain metastases on the adjacent brain parenchyma. Br J Cancer. 2020; 123(1):54–60. pmid:32362656

* View Article

* PubMed/NCBI

* Google Scholar

50. 50. Askew K, Li K, Olmos-Alonso A, Garcia-Moreno F, Liang Y, Richardson P, et al. Coupled proliferation and apoptosis maintain the rapid turnover of microglia in the adult brain. Cell Rep. 2017; 18(2): 391–405. pmid:28076784

* View Article

* PubMed/NCBI

* Google Scholar

51. 51. Davis BM, Salinas-Navarro M, Cordeiro MF, Moons L, De Groef L. Characterizing microglia activation: a spatial statistics approach to maximize information extraction. Sci Rep. 2017; 7(1):1576. pmid:28484229

* View Article

* PubMed/NCBI

* Google Scholar

52. 52. Hopewell JW, van der Kogel AJ. Pathophysiological mechanisms leading to the development of late radiation-induced damage to the central nervous system. Front Radiat Ther Oncol. 1999; 33(1):265–75. pmid:10549496

* View Article

* PubMed/NCBI

* Google Scholar

53. 53. Miyatake S, Nonoguchi N, Furuse M, Yoritsune E, Miyata T, Kawabata S, Kuroiwa T. Patho-physiology, diagnosis, and treatment of radiation necrosis in the brain. Neurol Med Chir. 2015; 55(1):50–9.

* View Article

* Google Scholar

54. 54. Lee D, Riestenberg RA, Haskell-Mendoza A, Bloch O. Brain Metastasis Recurrence Versus Radiation Necrosis: Evaluation and Treatment. Neurosurg Clin N Am. 2020; 31(4):575–587. pmid:32921353

* View Article

* PubMed/NCBI

* Google Scholar

55. 55. Furuse M, Nonoguchi N, Yamada K, Shiga T, Combes JD, Ikeda N, et al. Radiological diagnosis of brain radiation necrosis after cranial irradiation for brain tumor: a systematic review. Radiat Oncol. 2019; 14(1):28. pmid:30728041

* View Article

* PubMed/NCBI

* Google Scholar

56. 56. Dagogo-Jack I, Shaw A. Tumour heterogeneity and resistance to cancer therapies. Nat Rev Clin Oncol. 2018; 15:81–9. pmid:29115304

* View Article

* PubMed/NCBI

* Google Scholar

57. 57. Alfonso JCL, Papaxenopoulou LA, Mascheroni P, Meyer-Hermann M, Hatzikirou H. On the Immunological Consequences of Conventionally Fractionated Radiotherapy. iScience. 2020; 23(3):100897. pmid:32092699

* View Article

* PubMed/NCBI

* Google Scholar

Citation: Ocaña-Tienda B, León-Triana O, Pérez-Beteta J, Jiménez-Sánchez J, Pérez-García VM (2024) Radiation necrosis after radiation therapy treatment of brain metastases: A computational approach. PLoS Comput Biol 20(1): e1011400. https://doi.org/10.1371/journal.pcbi.1011400

References

1. Kaal EC, Niel CG, Vecht CJ. Therapeutic management of brain metastasis. Lancet Neurology. 2005; 4(5): 289–298. pmid:15847842

2. Lamba N, Wen PY, Aizer AA. Epidemiology of brain metastases and leptomeningeal disease. Neuro Oncol. 2021; 23(9):1447–1456. pmid:33908612

3. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer statistics, 2022. CA Cancer J Clin. 2022; 72(1): 7–33. pmid:35020204

4. Suh JH, Kotecha R, Chao ST, Ahluwalia MS, Sahgal A, Chang EL. Current approaches to the management of brain metastases. Nat Rev Clin Oncol. 2020; 17(5): 279–299. pmid:32080373

5. Kohutek ZA, Yamada Y, Chan TA, Brennan CW, Tabar V, Gutin PH, et al. Long-term risk of radionecrosis and imaging changes after stereotactic radiosurgery for brain metastases Journal of Neuro-Oncology. 2015; 125(1): 149–156. pmid:26307446

6. Kolesnick R, Fuks Z. Radiation and ceramide-induced apoptosis. Oncogene. 2003; 22(37):5897–5906. pmid:12947396

7. Abdulla S, Saada J, Johnson G, Jefferies S, Ajithkumar T. Tumour progression or pseudoprogression? A review of post-treatment radiological appearances of glioblastoma. Clin Radiol. 2015; 70(11):1299–1312. pmid:26272530

8. Peña LA, Fuks Z, Kolesnick RN. Radiation-induced apoptosis of endothelial cells in the murine central nervous system: protection by fibroblast growth factor and sphingomyelinase deficiency. Cancer Res. 2000; 60(2):321–327. pmid:10667583

9. Baker DG, Krochak RJ. The response of the microvascular system to radiation: a review. Cancer Invest. 1989; 7(3):287–294. pmid:2477120

10. Fike JR, Rosi S, Limoli CL. Neural precursor cells and central nervous system radiation sensitivity. Semin Radiat Oncol. 2009; 19(2):122–132. pmid:19249650

11. Kudo S, Suzuki Y, Noda SE, Mizui T, Shirai K, Okamoto M, et al. Comparison of the radiosensitivities of neurons and glial cells derived from the same rat brain. Exp Ther Med. 2014; 8(3):754–758. pmid:25120594

12. Nordal RA, Nagy A, Pintilie M, Wong CS. Hypoxia and hypoxia-inducible factor-1 target genes in central nervous system radiation injury: a role for vascular endothelial growth factor. Clin Cancer Res. 2004; 10(10):3342–3353. pmid:15161688

13. Sachet M, Liang YY, Oehler R. The immune response to secondary necrotic cells. Apoptosis. 2017; 22(10):1189–1204. pmid:28861714

14. Proskuryakov SY, Konoplyannikov AG, Gabai VL. Necrosis: a specific form of programmed cell death?. Exp Cell Res. 2003; 283(1):1–16. pmid:12565815

15. Shah R, Vattoth S, Jacob R, Manzil FF, O’Malley JP, Borghei P, et al. Radiation necrosis in the brain: Imaging features and differentiation from tumor recurrence. Radiographics. 2012; 32(5):1343–1359. pmid:22977022

16. Chao ST, Suh JH, Raja S, Lee SY, Barnett G. The Sensitivity and Specificity of FDG PET in Distinguishing Recurrent Brain Tumor from Radionecrosis in Patients Treated with Stereotactic Radiosurgery. Int J Cancer. 2001; 96(3):191–197. pmid:11410888

17. Kim HY, Cho SJ, Sunwoo L, Baik SH, Bae YJ, Choi BS, et al. Classification of true progression after radiotherapy of brain metastasis on MRI using artificial intelligence: a systematic review and meta-analysis. Neurooncol Adv.2021; 3(1):vdab080. pmid:34377988

18. Zhang Z, Yang J, Ho A, Jiang W, Logan J, Wang X, et al. A predictive model for distinguishing radiation necrosis from tumour progression after gamma knife radiosurgery based on radiomic features from MR images. Eur Radiol. 2018; 28(6):2255–2263. pmid:29178031

19. Larroza A, Moratal D, Paredes-Sánchez A, Soria-Olivas E, Chust ML, Arribas LA, Arana E. Support vector machine classification of brain metastasis and radiation necrosis based on texture analysis in MRI. Journal of Magnetic Resonance Imaging. 2015; 42(5):1362–1368. pmid:25865833

20. Lasocki A, Sia J, Stuckey SL. Improving the diagnosis of radiation necrosis after stereotactic radiosurgery to intracranial metastases with conventional MRI features: a case series. Cancer Imaging. 2022; 22(1):33. pmid:35794677

21. Barbolosi D, Benabdallah A, Hubert F, Verga F. Mathematical and numerical analysis for a model of growing metastatic tumors. Math Biosci. 2009; 218(1):1–14. pmid:19121638

22. Hartung N, Mollard S, Barbolosi D, Benabdallah A, Chapuisat G, Henry G, et al. Mathematical modeling of tumor growth and metastatic spreading: Validation in tumor-bearing mice. Cancer Res. 2014; 74(22):6397–6407. pmid:25217520

23. Bilous M, Serdjebi C, Boyer A, Tomasini P, Pouypoudat C, Barbolosi D, et al. Quantitative mathematical modeling of clinical brain metastasis dynamics in non-small cell lung cancer. Sci. Rep. 2019; 9(1):13018. pmid:31506498

24. Álvarez-Arenas A, Souleyreau W, Emanuelli A, Cooley LS, Bernhard JC, Bikfalvi A, Benzekry S. Practical identifiability analysis of a mechanistic model for the time to distant metastatic relapse and its application to renal cell carcinoma. PLoS Comput Biol. 2022; 18(8):e1010444. pmid:36007057

25. Chen LL, Blumm N, Christakis NA, Barabasi AL, Deisboeck TS. Cancer metastasis networks and the prediction of progression patterns. Br J Cancer. 2009; 101(1):749–758. pmid:19707203

26. Scott JG, Basanta D, Anderson AR, Gerlee P. A mathematical model of tumour self-seeding reveals secondary metastatic deposits as drivers of primary tumour growth. J R Soc Interface. 2013; 10(1): 20130011. pmid:23427099

27. Newton PK, Mason J, Bethel K, Bazhenova LA, Nieva J, Kuhn P. A stochastic Markov chain model to describe lung cancer growth and metastasis. PLoS ONE. 2012; 7:e34637. pmid:22558094

28. León-Triana O, Pérez-Beteta J, Albillo D, Ortiz de Mendivil A, Perez-Romasanta L, et al. Brain Metastasis response to stereotactic radio surgery: A mathematical approach. Mathematics. 2021; 9(7), 716.

29. Watanabe Y, Dahlman EL, Leder KZ, Hui SK. A mathematical model of tumor growth and its response to single irradiation. Theor. Biol. Med. Model. 2016; 136. pmid:26921069

30. Camacho A, Jerez S. Bone metastasis treatment modeling via optimal control. Journal of Mathematical Biology. 2018; 78(1-2):497–526. pmid:30132065

31. Benzekry S, Tracz A, Mastri M, Corbelli R, Barbolosi D, Ebos JM. Modeling Spontaneous Metastasis following Surgery: An In Vivo-In Silico Approach. Cancer Res. 2016; 76(3):535–547. pmid:26511632

32. Barbolosi D, Summer I, Meille C, Serre R, Kelly A, Zerdoud S, et al. Modeling therapeutic response to radioiodine in metastatic thyroid cancer: A proof-of-concept study for individualized medicine. Oncotarget. 2017; 8:39167–39176. pmid:28389624

33. Narasimhan S, Johnson HB, Nickles TM, Miga MI, Rana N, Attia A, et al. Biophysical model-based parameters to classify tumor recurrence from radiation-induced necrosis for brain metastases. Med Phys. 2019; 46(5):2487–2496. pmid:30816555

34. Dohm AE, Nickles TM, Miller CE, Bowers HJ, Miga MI, Attia A, et al. Clinical assessment of a biophysical model for distinguishing tumor progression from radiation necrosis. Med Phys. 2021; 48(7):3852–3859. pmid:34042188

35. Von Bertalanffy L. Quantitative laws in metabolism and growth. Q Rev Biol. 1957; 32(3):217–231. pmid:13485376

36. Pérez-García VM, Calvo GF, Bosque JJ, León-Triana O, Jiménez J, Pérez-Beteta J, Belmonte-Beitia J, et al. Universal scaling laws rule explosive growth in human cancers. Nature Phys. 2020; 16(12):1232–1237. pmid:33329756

37. Ocaña-Tienda B, Pérez-Beteta J, Jiménez-Sánchez J, Molina-García D, Ortiz de Mendivil A, Asenjo B, et al. Growth exponents reflect evolutionary processes and treatment response in Brain Metastases NPJ Syst Biol Appl. 2023; 9(1):35. pmid:37479705

38. Ocaña-Tienda B, Pérez-Beteta J, Molina-García D, Asenjo B, Ortiz de Mendivil A, Albillo D, et al. Growth dynamics of brain metastases differentiate radiation necrosis from recurrence. Neurooncol Adv. 2022; 5(1):vdac179. pmid:36726366

39. Sheline GE, Wara WM, Smith V. Therapeutic irradiation and brain injury. Int J Radiat Oncol Biol Phys. 1980; 6(9):1215–1228. pmid:7007303

40. Mineyeva OA, Bezriadnov DV, Kedrov AV, Lazutkin AA, Anokhin KV, Enikolopov GN. Radiation Induces Distinct Changes in Defined Subpopulations of Neural Stem and Progenitor Cells in the Adult Hippocampus. Front Neurosci. 2019; 12:1013. pmid:30686979

41. Mkango SB, Shaban N, Mureithi E, Ngoma T. Dynamics of breast cancer under different rates of chemoradiotherapy. Comput Math Methods Med. 2019; 2019:5216346. pmid:31611927

42. Villasana M, Radunskaya A. A delay differential equation model for tumor growth. J Math Biol. 2003; 47(3):270–294. pmid:12955460

43. Ocaña-Tienda B, Pérez-Beteta J, Villanueva-García JD, Romero-Rosales JA, Molina-García D, Suter Y, et al. A comprehensive dataset of annotated brain metastasis MR images with clinical and radiomic data. Sci Data. 2023; 10(1):208. pmid:37059722

44. Jiménez-Sánchez J, Martínez-Rubio A, Popov A, Pérez-Beteta J, Azimzade Y, Molina-García D, et al. A mesoscopic simulator to uncover heterogeneity and evolutionary dynamics in tumors. PLoS Comput Biol. 2021 17(2):e1008266. pmid:33566821

45. Joiner MC, van der Kogel A. Basic Clinical Radiobiology (4th ed.). CRC Press. 2009.

46. Milo R, Jorgensen P, Moran U, Weber G., BioNumbers-the database of key numbers in molecular and cell biology. Nucleic Acids Research. Nucleic Acids Res. 2009; 38:D750–D753. pmid:19854939

47. Kobets AJ, Backus R, Fluss R, Lee A, Lasala PA, Evaluating the natural growth rate of metastatic cancer to the brain. Surg Neurol Int. 2020; 11: 254. pmid:33024592

48. He B, Wang X, He Y, Li H, Yang Y, Shi Z, et al. Gamma ray-induced glial activation and neuronal loss occur before the delayed onset of brain necrosis. FASEB J. 2020; 34(10):13361–13375. pmid:32851721

49. Wagner S, Lanfermann H, Wohlgemuth WA, Gufler H. Effects of effective stereotactic radiosurgery for brain metastases on the adjacent brain parenchyma. Br J Cancer. 2020; 123(1):54–60. pmid:32362656

50. Askew K, Li K, Olmos-Alonso A, Garcia-Moreno F, Liang Y, Richardson P, et al. Coupled proliferation and apoptosis maintain the rapid turnover of microglia in the adult brain. Cell Rep. 2017; 18(2): 391–405. pmid:28076784

51. Davis BM, Salinas-Navarro M, Cordeiro MF, Moons L, De Groef L. Characterizing microglia activation: a spatial statistics approach to maximize information extraction. Sci Rep. 2017; 7(1):1576. pmid:28484229

52. Hopewell JW, van der Kogel AJ. Pathophysiological mechanisms leading to the development of late radiation-induced damage to the central nervous system. Front Radiat Ther Oncol. 1999; 33(1):265–75. pmid:10549496

53. Miyatake S, Nonoguchi N, Furuse M, Yoritsune E, Miyata T, Kawabata S, Kuroiwa T. Patho-physiology, diagnosis, and treatment of radiation necrosis in the brain. Neurol Med Chir. 2015; 55(1):50–9.

54. Lee D, Riestenberg RA, Haskell-Mendoza A, Bloch O. Brain Metastasis Recurrence Versus Radiation Necrosis: Evaluation and Treatment. Neurosurg Clin N Am. 2020; 31(4):575–587. pmid:32921353

55. Furuse M, Nonoguchi N, Yamada K, Shiga T, Combes JD, Ikeda N, et al. Radiological diagnosis of brain radiation necrosis after cranial irradiation for brain tumor: a systematic review. Radiat Oncol. 2019; 14(1):28. pmid:30728041

56. Dagogo-Jack I, Shaw A. Tumour heterogeneity and resistance to cancer therapies. Nat Rev Clin Oncol. 2018; 15:81–9. pmid:29115304

57. Alfonso JCL, Papaxenopoulou LA, Mascheroni P, Meyer-Hermann M, Hatzikirou H. On the Immunological Consequences of Conventionally Fractionated Radiotherapy. iScience. 2020; 23(3):100897. pmid:32092699