LymphoDose: a lymphocyte dose estimation framework—application to brain radiotherapy

Lymphocytes play a central role in the response to radio-chemotherapy and immuno-radiotherapy (Deutsch et al 2019). However, one of the iatrogenic effects commonly observed in radiotherapy (RT) is the appearance of severe radiation-induced lymphopenia (sRIL), which is defined as a significant decrease in the absolute lymphocyte count (ALC) under 500 lymphocytes per mm³ of blood (‘Common Terminology Criteria for Adverse Events (CTCAE),’ 2017). Approximately 30%–40% of sRIL are observed during or after the end of RT treatment in patients with gliomas (Mallick et al 2022, Zhang et al 2022). sRIL can persist for several months after the end of RT and represents an independent risk factor for reduced overall survival (OS) for both radio-chemotherapy and immuno-radiotherapy treatments (de Kermenguy et al 2023), with an estimated hazard ratio of 1.86 for patients with glioblastoma (El Houat et al 2023). Several independent variables have been identified as bad prognostic of sRIL such as a large volume of irradiation, a normo-fractionated treatment (compared with hypo-fractionated), photon irradiation (compared with proton) or the use of concomitant chemotherapy (de Kermenguy et al 2023). Recently, dosimetric constraints for immunological organs have been compiled (Venkatesulu et al 2022), with the goal to move to lymphocyte-sparing RT strategies to improve treatment outcomes by enhancing the anti-tumor immune response (Galluzzi et al 2023).

Several research groups have been interested in developing models to specifically assess doses received by blood lymphocytes (Hammi et al 2020, Sung and Cho 2022, Beekman et al 2023). This could be particularly relevant for brain tumors that contain few lymphoid structures and progenitor cells (van Hooren et al 2021), and where many of the irradiated lymphocytes in the RT field are found in circulating blood. However, none have yet considered the out-of-field (OOF) dose to the head and neck (H&N) region rich in lymphoid structures, or the slow recirculation and homing of lymphocytes between lymphoid organs such as the spleen or lymph nodes and blood (which only contains about 2% of the total lymphocyte pool at any given time (Ganusov and De Boer 2007, Sender et al 2023)). For these reasons, it is not trivial to extrapolate the doses received by circulating blood to the doses received by recirculating lymphocytes.

In this study, we propose a new in silico framework with highly efficient numerical implementation based on stochastic continuous time semi-Markov chain process to evaluate doses received by the lymphocyte pool after volumetric modulated arc therapy (VMAT) irradiation on 33 patients with primary brain tumors. Three cases of the model, corresponding to different physical and biological scenarios, were tested: (H1) no recirculation of lymphocytes outside of blood; (H2) recirculation between several lymphoid organs; (H3) recirculation between several lymphoid organs with OOF dose to lymphoid structures of the H&N region.

Modelling lymphocyte pool dynamics

To assess the doses received by the lymphocyte pool during irradiation, we modelled their dynamics within the organs. We distinguished two distinct but interdependent processes: (i) the slow recirculation of lymphocytes between lymphoid organs (Ganusov and Tomura 2021), and (ii) the rapid circulation of lymphocytes through all organs via the bloodstream. These 2 phenomena were modelled in this work by implementing 2 interconnected compartmental models which we will refer to as ${{\mathscr{M}}}_{1}$ and ${{\mathscr{M}}}_{2}$ respectively in the remainder of the manuscript. As illustrated in figure 1, both models implement a stochastic continuous-time semi-Markov chain process (Çinlar 1969) to simulate lymphocyte paths between compartments during the whole treatment. A continuous-time semi-Markov chain is a stochastic process with finite state space that changes states according to an arbitrary random variable and a transition matrix. In practice, the implementation of a continuous-time semi-Markov chain process consists of associating 2 parameters with each compartment of the ${{\mathscr{M}}}_{1}$ and ${{\mathscr{M}}}_{2}$ models: (1) a continuous distribution of lymphocyte residence times, which enables a residence time to be stochastically determined for each lymphocyte passage through that compartment via inverse transform sampling (ITS), and (2) a vector of transition probabilities to the connected compartments, in order to stochastically determine which compartment the lymphocyte will migrate to.

Zoom In Zoom Out Reset image size

Model ${{\mathscr{M}}}_{1}$ (figure 1 and Supplementary Data 1) is based on the work of Ganusov et al (Ganusov and Auerbach 2014), who proposed a model of lymphocyte recirculation via circulating blood between lungs, liver, spleen, subcutaneous lymph nodes (SCLNs), mesenteric lymph nodes (MLNs) and Peyer’s Patches (PPs), fitted to experimental data on lymphocyte migration in rats (Smith and Ford 1983). In this model, lung and liver compartments have extremely short mean residence times (0.46 min and 0.88 min respectively) compared with the other lymphoid compartments but captured most blood lymphocytes (78.2% and 17.4% respectively). Based on the hypothesis that lymphocytes circulating in these 2 organs do not a priori leave the bloodstream, the model was modified in our study by grouping 3 compartments {circulating blood, lung, liver} into a single compartment called Blood circulation. In accordance with (Ganusov and Auerbach 2014), residence times were distributed according to an exponential law with a parameter ${\lambda }_{M1}=\tfrac{1}{\mathrm{Mean}{\rm{\unicode{x02008}}}\mathrm{residence\; time}}$ for the Spleen, PPs and Blood Circulation compartments, and according to a gamma distribution with a shape ${k}_{M1}^{{\prime} }=2$ and a scale ${\theta }_{M1}=\tfrac{\mathrm{Mean}{\rm{\unicode{x02008}}}\mathrm{residence}\,\mathrm{time}}{{k}_{M1}^{{\prime} }}$ for the SCLNs and MLNs compartments.

Model ${{\mathscr{M}}}_{2}$ (figure 1 and Supplementary Data 1) was based on the ICRP Publication 89 (2002), which describes a transition matrix for blood flow between 28 organs in the adult male and female human body. Mean residence time values were adapted from the HEDOS model (Beekman et al 2023), considering a total blood volume of $5.3$ l and $3.9$ l and a cardiac output of $6.5$ l.min⁻¹ and $5.9$ l.min⁻¹ for male and female patients respectively. In accordance with the data given in (Shin et al 2021), residence times were distributed according to a Weibull distribution with a shape ${k}_{M2}=2$ and a scale ${\lambda }_{M2}=\tfrac{\mathrm{Mean}\,\mathrm{residence}\,\mathrm{time}}{{\rm{\Gamma }}(1+\tfrac{1}{{k}_{M2}})},$ with ${\rm{\Gamma }}$ the gamma function.

The code was developed in Julia language version 1.9.3 (Bezanson et al 2017) and simulations were performed considering ${10}^{5}$ particles.

Cohort of patients with glioblastoma

A cohort of patients with glioblastoma treated after 2018 with VMAT irradiation was created. The cohort consisted of 33 patients with glioblastoma, whose characteristics and treatment parameters are described in Table 1 and detailed in supplementary data 2. The cohort was composed of 20 male and 13 female patients with a prescribed dose of 60 Gy, 50 Gy and 40.5 Gy for 25, 5 and 3 patients respectively. The planning treatment volumes (PTVs) had a mean volume of 278 ± 121 cm³ and were irradiated at a mean dose rate of 1.14 ± 0.47 Gy min⁻¹.

Table 1. Table describing the cohort of 33 patients with glioblastoma. std: standard deviation.

	Patients with glioblastoma (total = 33)
Age (years)	Median	60.0
	Range (min–max)	38–80
Sex (n (%))	Male	20 (60.6%)
	Female	13 (39.4%)
Tumor localisation (n (%))	Temporal	13 (39.4%)
	Frontal	5 (15.1%)
	Others	15 (45.5%)
PTV volume (cm3)	Mean ± std	278 ± 121
Dose rate to PTV (Gy/min)	Mean ± std	1.14 ± 0.47
Prescribed dose (Gy)	60	25 (75.8%)
	50	5 (15.1%)
	40.5	3 (9.1%)
Dose per fraction (Gy)	2.0	25 (75.8%)
	2.5	5 (15.1%)
	2.7	3 (9.1%)
Number of fractions NF (n (%))	30	25 (75.8%)
	20	6 (18.2%)
	15	2 (6.0%)
Number of arc per fraction, NA (n (%))	1	4 (12.1%)
	2	18 (54.5%)
	3	11 (33.4%)
Arc duration ${{\boldsymbol{t}}}_{{{\boldsymbol{A}}}_{{\boldsymbol{i}}{\boldsymbol{,}}{\boldsymbol{j}}}}$(s)	Mean ± std	55 ± 17
Number of segments per arc NS (n (%))	= 91	21 (63.6%)
	< 91	9 (27.3%)
	> 91	3 (9.1%)

Simulation of patients with glioblastoma treatment

Out-of-field dose evaluation

Model comparison and sensitivity analysis

Lymphocyte dynamics

Out-of-field dose

Estimated doses to blood and lymphocytes

Case H1 (figure 2, H1) resulted in the irradiation of almost all the circulating blood (V_{>0 Gy} = 99.8 + 0.7%) at the end of VMAT treatment, at a mean dose of 309.9 ± 74.7 mGy (figure 1). The irradiated lymphocytes pass through the field at an average of 10.0 ± 4.9 times during the ${N}_{F}$ treatment fractions (figure 2, H1). Considering lymphocyte recirculation in H2 case (table 2 and figure 2, H2), the fraction of lymphocytes receiving a non-zero dose dropped to 40.4 ± 10.2%, with a lower mean dose estimated at 52.6 ± 21.1 mGy, especially because irradiated lymphocytes passed through the irradiation field only 1.58 ± 0.91 times on average. The most realistic scenario, including the OOF dose component (case H3) (table 2 and figure 2, H3), showed that almost all lymphocytes (97.6 + 2.5%) were irradiated (in the field or OOF in the H&N lymph nodes) at the end of the VMAT treatment, mostly due to the OOF dose (mean dose equal to 265.6 ± 48.5 mGy).

Zoom In Zoom Out Reset image size

Table 2. Different metrics summarizing the volume of lymphocytes irradiated and the dose received by the blood (H1) and of the recirculating lymphocyte (H2 and H3) during VMAT irradiation, averaged over the whole cohort. ${V}_{\gt 0{Gy}}$ (volume of lymphocytes having received a non-null dose), ${V}_{\gt 0.125{Gy}}$ (volume of lymphocytes having received a dose greater or equal to 0.125 Gy), the average, median and 98th percentile doses are reported. The last 3 metrics are calculated on all blood/lymphocytes (with or without irradiation) and on irradiated blood/lymphocytes only (${V}_{\gt 0{Gy}}$). For H3, ${V}_{\gt 0{Gy}}$ represents lymphocytes that have been irradiated at least once either in the field or out of the field in the lymph nodes of the H&N region.

Model	${{\boldsymbol{V}}}_{\gt 0{Gy}}$ (%)	${{\boldsymbol{V}}}_{\gt 0.125{Gy}}$ (%)	Average dose (mGy)		D50% (mGy)		D2% (mGy)
			To all lymphocytes	To ${{\boldsymbol{V}}}_{\gt 0{Gy}}$	To all lymphocytes	To ${{\boldsymbol{V}}}_{\gt 0{Gy}}$	To all lymphocytes	To ${{\boldsymbol{V}}}_{\gt 0{Gy}}$
H1	99.8 ± 0.7	87.6 ± 10.4	309.9 ± 74.7	309.9 ± 74.7	285.4 ± 71.3	285.9 ± 71.1	699.1 ± 215.0	699.4 ± 215.3
H2	40.4 ± 10.2	3.5 ± 2.1	19.8 ± 4.8	52.6 ± 21.1	0.4 ± 0.9	34.3 ± 15.6	161.8 ± 55.6	223.7 ± 116.5
H3	97.6 ± 2.5	82.2 ± 9.5	259.5 ± 48.7	265.6 ± 48.5	255.8 ± 49.9	262.7 ± 50.5	559.8 ± 103.2	561.6 ± 103.7

Sensitivity analysis

To our knowledge, this framework is the first to simulate radiation doses received by lymphocytes during standard of care brain irradiation considering (i) that lymphocytes are mostly found in lymphoid organs and not only in the peripheral circulation and (ii) that out of the field doses may significantly contribute to radio-induced lymphopenia. We used two interconnected compartmental models that mimic the presumed two-speed average journey of lymphocytes between blood, where they migrate fast, and lymphoid organs, where they may reside for a long time. Considering in- and out-field irradiation, the simulation revealed that 97.6% of total body lymphocytes may have passed through a radiation-receiving area during the full course of radiation therapy in patients treated with chemoradiation for glioblastoma, with 82.2% ± 9.5% of total body lymphocytes receiving doses >0.125 Gy. On average on the patient cohort, irradiated lymphocytes would receive a mean cumulative radiation dose of 265.6 mGy ± 48.5 mGy when treated with brain-directed conventional VMAT for glioblastoma.

Other teams previously intended to simulate radiation doses received by circulating lymphocytes upon brain irradiation, without considering recirculation or OOF doses (Hammi et al 2020, Shin et al 2021). When considering that all lymphocytes are kept in the peripheral circulation as a closed system during radiation treatment of glioblastoma, our H1 model estimated that 99.8% ± 0.7% of all circulating lymphocytes received radiation at any dose during treatment, with a mean dose of 309.9 mGy ± 74.7 mGy. This is consistent with the findings obtained by Shin et al who modelled a dose to circulating blood cells after 30 fractions of brain-directed intensity modulated radiation therapy (IMRT) treatment (8). Hammi et al obtained a blood mean dose of 133 mGy after 30 fractions of brain-directed IMRT when precisely modelling the cerebral vasculature using a sophisticated 4D blood flow model (Hammi et al 2020). However, accurate estimation of the radiation dose delivered to lymphocytes is a complex task considering lymphocytes’ recirculation and homing process, mediated by specific biological mechanism such as L-selectin (or CD62L) receptor and Sphingosine-1-Phosphate (S1P) pathways. Lymphocytes can traffic along lymph vessels and through lymphoid organs, thus transiently exit the blood circulation and reintegrate it later, implying that the whole blood-circulating lymphocyte pool can be renewed up to 11 times a day (GOWANS 1957, Gough and Crittenden 2022). Thus, peripheral lymphocytes present in the blood at the time of irradiation can be different from a radiation fraction to another, and models restricted to circulating lymphocytes may significantly misestimate the lymphocyte dose.

Thus, considering these recirculation mechanisms within the H2 and H3 cases, the lymphocyte dose values estimated by our model (figure 2, H2 and H3) can be compared with the radiosensitivity of lymphocytes. A recently published review by Paganetti (Paganetti 2023) synthesized several in vivo studies and estimated that the average T lymphocytes dose–response $\alpha$ values was ∼0.6 Gy⁻¹, which corresponds to a 50% lethal dose (LD50) of 1.15 Gy with a linear model of cell survival. Moreover, lymphocytes radiosensitivity has been shown to depend on multiple factors, including subtype lineage (B cells being the most sensitive, followed by T cells and NK cells), activation status (naïve cells being more sensitive than educated ones), function and most importantly, local microenvironment at the time of radiation, with tumor-infiltrating lymphocytes being more radioresistant than circulating ones (Trowell 1952, Geara et al 1992, Arina et al 2019, Heylmann et al 2021).

Despite the notable radiosensitivity of circulating lymphocytes, our results (regardless of the H1, H2 or H3 model used) hardly explain the 30 to 40% severe lymphopenia rate often observed in series of patients with brain tumors treated with chemoradiation (de Kermenguy et al 2023), especially as the total estimated dose is usually delivered over a 6-week period, which presumably gives sufficient time for DNA damage repair between fractions in cells exposed to sub-lethal low-dose radiation. This calls into question the relative contribution of the direct cytotoxic effect of lymphocyte irradiation on lymphopenia incidence after standard of care radiotherapy for glioblastoma (figure 3). Importantly, independent series showed that glioblastoma and more broadly, intracranial tumors are intrinsically characterized by a S1P1-mediated sequestration of T cells in the bone marrow, which could contribute to the worsening of lymphopenia upon lymphotoxic treatment (Mahaley et al 1977, Chongsathidkiet et al 2018).

Zoom In Zoom Out Reset image size

In addition, lymphocytopenia can be influenced by the concomitant exposure to systemic medications, including temozolomide and corticosteroids. Temozolomide is an alkylating agent that rapidly passes through the blood–brain barrier and synergies with irradiation to kill brain tumor cells. As a myelosuppressive agent, it can cause killing with concurrent irradiation, can also induce a lymphopenia (Wick and Weller 2005, Brandes et al 2006) that preferentially affects the CD4+ T compartment (Su et al 2004). In a study that enrolled 562 elderly patients with glioblastoma receiving hypofractionated radiotherapy with or without temozolomide, the rates of severe lymphopenia were 27.3% and 10.3% respectively in the groups with or without temozolomide (Perry et al 2017). In addition, patients suffering from advanced glioblastoma often require long-term steroid therapy, which can impair immune functions by suppressing pro-inflammatory cytokine release pathways (Goodman et al 2023) or by redistributing blood circulating lymphocytes to the bone marrow (Fauci 1975). Some brain tumors, such as meningiomas, generally do not require chemotherapy or steroids as part of standard treatment (Maggio et al 2021). Studying patients with such tumors could make it possible to separate the effects of chemotherapy and/or steroids from those of radiotherapy on the incidence of lymphopenia.

Then, some experimental results tend to support the hypothesis of other underlying biological effects to explain RIL, such as radiobiological effects at a distance (Kapoor et al 2015)(spatial and/or temporal) indirectly impacting lymphocyte homeostasis (Piotrowski et al 2018) or recirculation (Vianello et al 2013). Another hypothesis relies on the role of immune suppressive cells in response to radiation. Ghosh et al showed that in patients with glioblastoma, chemoradiation-induced lymphopenia was a direct consequence of the increased representation of myeloid-derived suppressive cells (MDSC) in the peripheral blood after treatment, with radiation therapy triggering an aberrant myelopoiesis in the bone marrow (Ghosh et al 2023). MDSC directly suppress T cells activation and proliferation, at least partly by inhibiting pro-inflammatory cytokine release and by driving the differentiation of CD4+ T cells towards immunosuppressive regulatory T cells (Kumar et al 2016, Cassetta et al 2020). In mice models grafted with glioblastoma cells, blocking MDSC using pharmaceutical compounds was sufficient to prevent radiation-induced lymphopenia, suggesting that this mechanism may have a greater impact than direct lymphocyte cytotoxicity induced by radiation (Ghosh et al 2023).

Several perspectives are open for future research. Firstly, clinical and pre-clinical data are essential to validate our results and to complete the modelling, for example to take into account lymphocyte death and homeostasis (Pham et al 2023). Secondly, the recirculation model used is based on rats data. To our knowledge, there are no global studies that exist today on human data. There may be important differences between humans and rats in terms of lymphocyte recirculation and homing, especially because the two species have partially different immune systems (den Braber et al 2012, Shay et al 2013, Wildner 2019). Moreover, even if the estimated mean transit time of ${{\mathscr{M}}}_{1}$ Blood circulation was consistent with results observed in other articles for naive lymphocytes (Smith and Ford 1983, Ager 1994, Ganusov and Tomura 2021), several different lymphocyte subpopulations with different recirculation patterns coexist in the blood, and therefore models different from the ${{\mathscr{M}}}_{1}$ model are required. Thus, the ${{\mathscr{M}}}_{1}$ model was established considering naive lymphocytes, but recirculation patterns and residence times can vary greatly when considering activated lymphocytes or tissue resident lymphocytes (Ganusov and Tomura 2021). As an example, $\gamma \delta$ T cells in sheep had significantly higher levels of L-selectin, a surface receptor important for homing in lymph nodes, and therefore recirculated more actively than other subpopulations such as $\alpha \beta$ T cells or B cells (Young et al 2000). Similarly, a large B cell subpopulation excluded from lymphatic recirculation and showing increased migration to the spleen has been identified (Young et al 1997). This sub-population, which a priori spends more time in the circulating blood, may be more exposed to radiation during treatment and may therefore be negatively selected. This selection of cell subpopulations according to their recirculation features could participate to explain the ALC ‘plateau’ observed in some studies (Ellsworth et al 2020, de Kermenguy et al 2023).

In addition to the specific dynamics of these subpopulations, it would be interesting to study the effect on recirculation of the irradiation or the presence of a tumor, the latter can strongly influence immune system by creating tertiary lymphoid structures or tumor-associated high endothelial venules (Blanchard and Girard 2021). All these considerations call into question the relevance of the characterisation of RIL currently used. Indeed, simple blood sampling does not seem to be able to describe all the complexities of the processes involved, and therefore imposes major obstacles on the interpretation of these data. There is no guarantee that the blood sample taken at a given time is representative of the whole immune system, either quantitatively or qualitatively (lymphocyte sub-populations and their functionality). A same effect on ALC could arise from different causes, depending on the irradiation site for example. Thus, in-depth studies of lymphocyte recirculation and homing in humans using recent imaging methods (Meyblum et al 2023) could shed light on the question of RIL, but also in the much broader context of radio-immunotherapeutic combinations.

Thirdly, average residence times used in the ${{\mathscr{M}}}_{2}$ model, for the Brain compartment in particular, are those of circulating blood. However, it could be significantly higher for the lymphocytes it contains, in particular because of their ability to adhere to endothelial cells, or simply because their diameter (∼7 μm)(Downey et al 1990) is of the same order of magnitude as diameters of some cerebral capillaries (Marín-Padilla 2012). Moreover, the dose-volume histogram used to assess the dose in the irradiation field (the ${{\mathscr{M}}}_{2}$ Brain compartment in this study) only takes into account the tumor volume and contains no information on tumor location. The underlying assumption is that blood is distributed relatively homogeneously in the brain. A further development would be to simulate the vascularization of the irradiated area (as in (Hammi et al 2020)), to better take into account tumor location and potential inhomogeneity of blood distribution in the encephalon. In addition, the modelling of the OOF dose is currently approximate. The fraction of lymphocytes in H&N SCLNs was not patient-specific and assumed a homogeneous distribution of lymphocytes between all SCLNs, and the code used to calculate the OOF dose was trained on paediatric OOF dose data (161 of whom were over 15 years old) from old-generation accelerators and probably overestimates OOF doses. Even if the neural network showed good performance in inferring OOF dose from unknown machines (irradiators not included in the training set), especially for close to the field regions (Benzazon et al 2024), experimental measurements are required to validate these results. Finally, it would be interesting to consider the whole-body OOF dose to lymphoid structure, and not just in the H&N lymph nodes.

We have developed an in-silico framework that estimates the radiation dose received by recirculating lymphocytes, applied to brain VMAT irradiation. The model considers both the lymphocytes recirculation flow and the OOF dose to H&N lymphoid structures. It is an important preliminary step towards a better understanding of the radiobiological mechanisms underlying RIL, and the model could be applied to other irradiation locations in the future. From a clinical point of view, reducing the incidence of RIL could not only improve the response to radiotherapy and chemoradiation, but also enhance the potential of radioimmunotherapy combinations or abscopal effect.

This study was supported by the grant N°ANR-21-RHU5-0005 within the FRANCE2030 investment plan and 2018-1-PL BIO-06-IGR-1².

Main work and first draft of the manuscript: FDK, Project initiation and supervision: CR and ED, Generation of dose map for each irradiation segment: RV and CM, Creation of the patient cohort: EL and VJ, Development of the OOF dose calculation code and generation of OOF dose maps: NB, Advice on the analysis and manuscript: PM, MM’h, ID, DM, CC and MM.

All data that support the findings of this study are included within the article (and any supplementary information files).

The utilization of this retrospective training cohort was performed in accordance with the General Data Protection Regulation (GDPR) and approved by the local ethical committee (IRB number 2023-199).

Not needed for this work.

The authors declare no conflict of interest.