Mechanistic model of radiotherapy-induced lung fibrosis using coupled 3D agent-based and Monte Carlo simulations

Background Mechanistic modelling of normal tissue toxicities is unfolding as an alternative to the phenomenological normal tissue complication probability models. The latter, currently used in the clinics, rely exclusively on limited patient data and neglect spatial dose distribution information. Among the various approaches, agent-based models are appealing as they provide the means to include patient-specific parameters and simulate long-term effects in complex systems. However, Monte Carlo tools remain the state-of-the-art for modelling radiation transport and provide measurements of the delivered dose with unmatched precision. Methods In this work, we develop and characterize a coupled 3D agent-based – Monte Carlo model that mechanistically simulates the onset of the radiation-induced lung fibrosis in an alveolar segment. To the best of our knowledge, this is the first such model. Results Our model replicates extracellular matrix patterns, radiation-induced lung fibrosis severity indexes and functional subunits survivals that show qualitative agreement with experimental studies and are consistent with our past results. Moreover, in accordance with experimental results, higher functional subunits survival and lower radiation-induced lung fibrosis severity indexes are achieved when a 5-fractions treatment is simulated. Finally, the model shows increased sensitivity to more uniform protons dose distributions with respect to more heterogeneous ones from photon irradiation. Conclusions This study lays thus the groundwork for further investigating the effects of different radiotherapeutic treatments on the onset of radiation-induced lung fibrosis via mechanistic modelling.


Introduction
Precision, efficacy and non-invasiveness have made radiotherapy (RT) a first-choice treatment option nowadays for a large portion of cancer patients.Nevertheless, the risk of developing radiation-induced pathologies remains significant and the number of reported injuries increases concurrently with the number of treatments.Notably, more than 50% of the cancer patients are treated with RT 1 and previous studies reported that thorax irradiation led to radiation-induced lung injuries (RILI -namely, pneumonitis and fibrosis) in up to 30% of the cases 2 .These pathological conditions are thought to be triggered mainly by the damaged alveolar epithelium and the resulting inflammation, if not resolved within a few weeks' time, can lead to lung stiffening and eventually death 3 .Consequently, greater efforts are required to improve our understanding of the mechanisms underlying the paths that link radiation damage and toxicity in normal (i.e.non-tumoural) tissues, with the ultimate goal of widening the therapeutic window.
To this aim, normal tissue complication probability (NTCP) models have been implemented that estimate the probability of developing new pathologies as a function of the delivered dose for a given organ [4][5][6] .However, NTCP models widely employed in clinical practice restrict the set of inputs to the delivered dose distribution and a few macroscopic organ-specific parameters that fail to enfold a mechanistic description of the pathways involved.Moreover, the lack of spatial information in the dose-volume histograms (DVHs) that encode the dose distributions' data results in incomplete descriptions of the treatment plans and subsequent erroneous shifts of the NTCP curves 7,8 .By combining Agent-Based (AB) modelling and a Monte-Carlo (MC) simulator we propose an innovative approach in an attempt to address the aforementioned shortcomings (see figure 1).AB modelling is a powerful and relatively recent set of computational techniques that allows simulations of concurrent and independent entities, namely, the agents 9 .Each agent is initialized at the beginning of a new simulation with a set of rules, encoded into its behaviours, and positioned in an environment that can be sensed by and react to the agent itself.The interactions among different agents and between the surrounding environment and the agents themselves can lead to the emergence of elaborated dynamics and provide a framework to model complex systems, such as the biological ones.MC methods, on the other hand, have been extensively used to simulate the interaction of radiation with matter 10 and can provide accurate estimates of, among others, dose depositions in biological structures.
Figure 1.Schematic representation of the implemented workflow.The ABM is used to generate the alveolar segment structure in homeostatic conditions.The data are then fed into the MC model which replicates the geometry and irradiates it.Information about the delivered dose to each cell are then loaded into the ABM which simulates the evolution of the system.If a multi-fractionation scheme is set up, the data exported from the ABM are used as an input to the MC.If not, the data can be directly analysed to determine whether RILF onsets and, if so, its severity.
In what follows we provide insights on the onset of the radiation-induced lung fibrosis (RILF) in a coupled 3D Agent-Based (ABM) -MC model of an alveolar segment.We updated our ABM of RILF 11,12 and implemented the alveolar segment geometry within the TOPAS-nBio MC extension 13 .As a result, our coupled model is now able to provide a mechanistic description of the depletion of the functioning distal airways of the lung as well as 3D spatial information about the delivered dose distributions.More specifically, this work details the implementation of the coupled model and provides a comparison between the ABM and the ABM-MC models' results.Moreover, the effects of different radiation qualities as well as the impact of the bystander threshold on the outcomes are shown.Finally, we demonstrate how the use of a multi-fractionation scheme affects the model's outcomes.Temporal fractionation techniques are widely adopted in the clinics [14][15][16] due to their ability to spare normal tissues without losing efficacy in tumour cells killing 17 .Our results, in agreement with previous experimental studies 18 , show a right shift towards higher doses for the same amount of damage when using 5 fractions with respect to a single fraction.

Characterization of the dose distributions in the alveolar segment
Before running the full ABM-MC model, different setups of the photon source for the irradiation of the alveolar segment were tested and the corresponding dose distributions compared.More specifically, doses were delivered using both an external beam with 4 coplanar fields and an isotropic source located at the center of the alveolar segment.The particle energy was set to Eg = 1 keV for all the experiments and the total number of histories (i.e.emitted particles per field) was set to 1 million, which ensured reasonable simulation times on a common laptop.Moreover, each experiment was performed 10 times and the results averaged.As can be seen from figures 2a and 2b, while the dose delivered with the external beam could be fitted to a gaussian distribution that peaks at around 50% of the maximum dose with thin (although short) tails, the dose distribution from the isotropic irradiation shows two distinct peaks (at 0-1% and 10% of the maximum dose) and a long right tail.Given the symmetry of the resulting dose distribution and its closer resemblance to the beam setups used in the clinics, the external beam was used in all the experiments performed in this study.We further analysed the effect of an increase in the number of histories in the dose distribution delivered by the external beam.Raising the number of particles from 4 millions (for 4 fields) to 30 millions (figure 2c) didn't alter the shape of the dose distribution substantially, but required a greater computational cost.Therefore, the lower value was set in the parameter file of the MC model and used in the following simulations.

Simulation outcomes: ABM-MC vs ABM
To characterise the coupled ABM-MC model, we compared the main simulation outcomes with those from our previous ABM model 11 .In particular, we evaluated the surviving fraction of the FSUs, the RSI, the early and the late ECM fluctuations in the ECM concentration as a function of the delivered dose.To ensure that the system reached the steady state (both in terms of cell number and substances concentration, see supplementary material), we simulated 1200 days.In the ABM-MC model the dose was delivered once at the beginning of the simulation and the parameters were set to the same values as those in the ABM model (notably: phagocytic fraction = 100%, phagocytic index = 1, apoptotic-tosenescent ratio = 0 and bystander threshold = 2).As can be seen from figure 3, both the early and the late components of the DECM (for the ABM-MC model) could be fitted by the equations ( 2) and (3) and are qualitatively in agreement with experimental data [18][19][20][21] .As expected, the DECM increases as a function of the dose and reaches a plateau around 20Gy for the early component and 30Gy for the late component.Due to the higher fraction of mesenchymal cells in the inflammatory phase (i.e.before the macrophages are able to clean the senescent cells), the DECMmax of the early phase is more than twice as much as that in the late phase and the plateau is reached at lower doses as the total amount of ECM that can be stored in the alveolar space is limited.Figure 3c shows the surviving fraction of the FSUs as a function of the dose.For doses equal to or smaller than 7.5Gy in the ABM-MC model the survival is close to 100%, which implies that no alveolus is fully depleted from the healthy AEC2, neither directly or indirectly.Moreover, fitting of the dataset with the equation (4) shows better agreement than that observed for the results from the ABM model alone as well as improved consistency at the lower doses.Finally, the late component DECM and the FSU surviving fraction were combined into the RSI, a surrogate measure of the RILF severity, presented in figure 3d.The dataset obtained from the ABM-MC model could be fitted using equation (3) and showed qualitative agreement with experimental findings 18 , although the RSI saturates at higher values than the FI.Overall, the outcomes from the ABM-MC model are in agreement with experimental results found in the literature.However, our previous model showed increased damage (i.e.lower surviving fractions and higher DECM) at the same doses, likely caused by narrower and more peaked dose distributions.

Temporal fractionation
To assess the ability of the alveolar segment model to repair the damaged tissue in between irradiations, we implemented a 5-fraction scheme using the ABM-MC model.The temporal fractionation involved the delivery of relatively small doses (from 0.5 to 8.5 Gy) 24 hours apart, thus allowing the M1 and M2 macrophages to gather (triggered by the secretion of MCP-1 from the damaged AEC2) and start removing the senescent cells.Moreover, the healthy AEC2 could replenish the depleted alveoli.After a 1-day simulation, the system was "frozen" and the structure of the alveolar segment loaded into TOPAS-nBio for the following irradiation.
The number of fractions and the doses were matched to those used in the experiments presented in the work by Zhou et al 18 , which was used to evaluate our previous ABM model.Figure 4 presents a comparison of the main outcomes from the ABM-MC model, using 1 and 5 fractions, as a function of the total delivered dose.The results were fitted using equations ( 3) and ( 4).As expected, a right shift in both the late component DECM and in the FSU survival was observed.Consequently, an isoeffect was measured at higher doses for the RSI when using 5 fractions with respect to 1 fraction.Similar RSI values for the two fractionation schemes were measured for doses equal to or lower than 10Gy and above 35Gy, due to the high variability in the FSU survival at low doses and the formulation of the RSI, that forces it to saturate at high doses.Of note, the FSU survival data for the single fraction deviated from the LQ fit at high doses (³ 30Gy).
The results from our ABM-MC model with 5 fractions could not only be fitted using the equations provided by previous experimental studies, but also matched the observations from Zhou et al. 18 regarding the right shift of the RSI (a surrogate of the severity of the RILF, as is the FI) as the number of fractions increased.

Impact of the bystander threshold
The bystander threshold dictates the minimum number of senescent AEC2 necessary to damage an healthy AEC2 that is located in their neighbourhood.It is, therefore, the main regulator for controlling the speed of the spread of the indirect damage.In fact, the probability that an healthy AEC2 will be indirectly damaged depends on the time spent in the neighbourhood of senescent cells and this, in turn,  22 ).In our previous work 11 , we determined the optimal values for the phagocytic fraction and phagocytic index, i.e. 100% and 1.Moreover, we showed how different values of the apoptotic-to-senescent ratio affected the survival of the FSUs.We assessed the best value for the bystander threshold as well (i.e. 2), but given that the results from the ABM-MC model exhibited decreased damage for the same doses with respect to our previous model, we evaluated the impact of a faster spread of the indirect damage (by setting the bystander threshold to 1) to determine if the effects could offset each other.As can be seen from figure 5, even though the ABM-MC model showed reduced damage with respect to the ABM one, it could not compensate the increased spread of the indirect damage.In fact, lowering the bystander threshold led to FSU survivals < 100% and DECM > 0 even at very low doses.As a consequence, the lowest RSI did not drop below 60 and none of the major outcomes could be fitted using equations ( 3) and (4).

Radiation qualities
The use of charged particles in the clinical setting has seen rapid spread in recent times due to the increased conformality of the dose deposition to the target volume with respect to the photons 23 .Moreover, experimental studies 24 reported differences in the dynamics of the normal tissue response between lung cancer patients treated with photons and protons.Contextually, we used TOPAS-nBio, coupled to our ABM model, to compare the effect of different radiation qualities on the alveolar segment.In particular, we simulated the irradiation of the lung structure with 2 million protons delivered with an external beam using 4 coplanar fields in one fraction.TOPAS-nBio can generate standard DNA damage files for each cell that is hit by one or more particles, thus providing punctual information about each and every interaction between the particle and the DNA.However, the process on a common laptop is relatively slow and the amount of data that is stored for each cell can be very large, thus making it prohibitive for thousands of cells as those in the alveolar segment.Therefore, our simulations didn't take into account the differences in the patterns generated by photons and protons when damaging the DNA of the AEC2 cells.Different simulation outcomes are thus solely the result of different dose distributions.We simulated 60MeV protons as previous studies have reported negligible differences in the DNA damage response between proton and photon irradiations in the plateau region for this energy 25 and very low linear energy transfer (LET, < 2.5 keV/µm) at short distances (< 5mm) 26 .This requirement ensured very small LET variations along the proton tracks (as the diameter of the alveolar segment is < 1mm) and allowed us to use the LQ parameters derived from the survival curve of the AEC2 irradiated with photons.Figure 6 shows the dose distribution for the protons expressed as a fraction of the maximum delivered dose.When compared with the photons histograms in figure 2, the protons distribution exhibits a more pronounced and narrower peak and the data points align well with a gaussian distribution.When comparing the main simulation outcomes from the different radiation qualities, a marked shift towards lower doses for the same effect (i.e.same DECM and FSU surviving fraction), shown in figure 7, was observed for the protons.This, in turn, resulted in higher RSI at the same dose, similar to what was found when comparing different fractionation schemes for the photons (see figure 4).
Given that the FSU survivals could be fitted by the LQ model from equation (4) for both the photons and the protons, we introduced a relative biological effectiveness 27 relative to the FSU survival (RBEFSU).The RBEFSU was defined as the ratio of the of the absorbed dose from the photons to the absorbed dose from the protons resulting in the same effect and provides a quantitative measure of the relative effectiveness of the different radiations in depleting the FSUs.The definition reads as follows: As can be seen from figure 7b, small variations were observed for the RBEFSU at the 3 isoeffects considered.While it is generally assumed an RBE = 1.1 for protons relative to photons for identical cell survivals 28 , with lower RBEs for higher reference doses, we measured increasing values for identical FSU survivals: 1.30, 1.31 and 1.34 for 50%, 37% and 10% survival, respectively.However, it is worth noting that for high doses the data points for the photons deviate from the LQ fit and therefore the RBEFSU at 10% survival might not be representative of the actual relative effectiveness of the radiation in depleting the FSUs.

Discussion
Despite being a crucial component in a large variety of cancer treatments, the efficacy of RT is intrinsically hindered by normal tissue toxicities 29 .Increasingly sophisticated NTCP models 7 have been implemented to estimate risk probabilities from patients' DVH, but the lack of mechanistic information constrains superficial representations of the underlying mechanisms and prevents patientsspecific parameters to be taken into account.
To address the aforementioned shortcomings, we have previously implemented an ABM that simulates the onset of RILF in an alveolar segment 11 .In this work, we outlined the development of a coupled ABM-MC model, where the alveolar segment structure was rebuilt using TOPAS-nBio 13 and linked to an updated version of our previous model, with a custom script in charge of handling the communication between the models.As discussed in our former work, the ability of the model to replicate experimental results was assessed via the FSU (i.e. the alveoli) survival, the ECM increase and the RSI as functions of the dose and the output datasets fitted using equations ( 3) and ( 4).More specifically, in the new implementation the ABM is used to generate the initial structure in healthy conditions and afterwards the cells' position, type and size are exported.Using this data as an input, the MC tool builds the structure in real-time and an irradiation is simulated based on the information provided via the parameter control system.Accordingly, realistic dose distributions can be simulated and dose depositions at cell scale registered.The generated data are then sent back into the ABM, where the absorbed doses are used to determine the fate of each cell using the LQ model.Subsequently, the ABM model is run and the onset of RILF simulated.Finally, to replicate temporal fractionation schemes the workflow can be executed multiple times by setting the number of fractions in the control script.We compared different setups of the photons' source by characterizing the resulting relative dose distributions in the alveolar segment.To maximize the dose homogeneity, the tested sources were either isotropic, located at the center of the simulated structure, or external, with 4 coplanar fields spaced 90 degrees apart.While the first setup resulted in a positive-skewed distribution with 2 peaks in the low-dose region and a long right tail, the one provided by the external beam was symmetric and peaked at 50% of the maximum dose.Moreover, figure 2c shows that an increase in the number of histories from 4 to 30 millions didn't alter the dose distribution significantly.In light of these results, the external beam was selected as the default photon source for the following simulations.In accordance with the results from our previous model, the outcomes of the ABM-MC approach (single fraction) qualitatively matched published experimental results.In particular, the DECM showed two distinct components, an early one at around 3 months from the irradiation and a late one as the model reached the steady state 30 .Both the datasets exhibited a sigmoidal response as the dose increased and could be fitted using equations ( 2) and (3).The surviving fraction of the FSUs was measured at the end of each simulation as the ratio of the number of alveoli with at least an healthy cell to the initial number of alveoli, expanding upon the idea introduced in the critical volume NTCP model 31,32 .The results were fitted using equation ( 4) and showed good agreement with the LQ model for doses < 20Gy.Finally, the previous metrics were combined into the RSI, which mimics the fibrosis index introduced by Zhou et al. 18 and provides a surrogate measure of the RILF severity, and could be fitted using equation (3).Overall, the ABM-MC model preserved the ability to simulate existing experimental results.Of note, despite showing consistency with the previous ABM model, it can be seen from figure 3 that all the curves were shifted towards higher doses for the same outcomes.We attribute these differences to the shape of the newly generated dose distributions which peak on the mean lung dose, but span a relatively large range with respect to those in the ABM.Moreover, in the old model the delivered dose was used to determine the initial damaged fraction which was the same for all the cells in an alveolus.Conversely, the ABM-MC simulates an actual dose distribution using the TOPAS-nBio module and each cell, given the absorbed dose, has a certain probability to be damaged, in accordance with the LQ model.Temporal fractionation schemes as a way to spare the normal tissue as sublethal damages are repaired and depleted areas repopulated are widely employed in the clinics [14][15][16] .Accordingly, we compared the results from 1-fraction with those from 5-fractions using the ABM-MC model.The ABM model was run for 1 day between fractions and the structure's data at the end of the simulation used to feed the MC model.After the 5 th fraction, the ABM model ran for 1200 days to match the simulation time from the 1-fraction simulation.As expected, the outcomes in Figure 4 show a shift towards higher doses for the alveolar segment irradiated with 5 fractions and could be fitted using equations ( 3) and (4).Notably, the RSI transition to higher doses from 1 to 5 fractions mimics the findings from 18 .The model is thus able to replicate normal tissue sparing as the macrophages clean the senescent cells and prevent excessive spread of the damage via bystander mechanisms, and the AEC2 repopulate the depleted alveoli.It is worth noting that the actual DNA damage is not taken into account by our model and therefore the repair of sublethal damages was not simulated.However, it was recently shown in 33 that residual DSBs in normal tissue cells exhibit low predictive power for cell survival and thus our assumptions remain valid.Finally, we noticed that delivering the dose in 5 fractions improved the ability of the LQ model to fit the FSUs survival with respect to the single fraction, especially at high doses.
The key role played by the bystander effect in the onset of RILF has been emphasized by multiple studies [34][35][36] and earlier works highlighted that, when radiation-induced lung toxicities are modelled mechanistically, inflammation-induced tissue damage should be taken into account 37 .In our previous work we detailed the implementation of the damage spread behavior for the senescent AEC2 agents, where the minimum number of senescent AEC2 needed to damage an healthy cell can be controlled via the bystander threshold parameter.Given that the ABM-MC model showed lower damage for the same dose with respect to our previous model, we investigated the effect of an increase in the speed of the damage spread on the model outcomes by lowering the bystander threshold.Figures 5a and 5b show that the reduced damage can't compensate for the increase in the bystander effect and, even at the lowest doses, the DECM didn't drop to 0 and the FSU survival was < 100%.This, in turn, resulted in a higher RSI with respect to that measured with the bystander threshold = 2. Notwithstanding this, the model proved to be robust, as it retained the ability to recover in-between fractions, as demonstrated by the right shift towards higher doses when a 5-fraction scheme was used.To further characterize the ABM-MC model, we simulated a monoenergetic 60MeV proton irradiation using the same setup for the source as that previously described for the photons.Remarkable differences with respect to the photons irradiation were found in the resulting dose distribution (see figure 6), which exhibited a narrower peak at 80% of the maximum delivered dose.As a result, despite the DNA damage was not taken into account, the model outcomes differed significantly from those obtained when the alveolar segment was irradiated with photons using a single fraction.Figure 7 shows increased and lower FSU survival for the same doses when protons were used, resulting in a higher RSI for doses > 10Gy and similarly to what Zhou et al. 38 observed when mice were irradiated using high LET particles.Furthermore, we compared the ability of the different radiation qualities to deplete the alveoli from the AEC2 by computing the RBE at different FSU survivals, a metric that accounts for long term effects as it was computed when the model reached the steady state and the alveolar segment was thus allowed to recover.We measured an average RBEFSU = 1.32, however, the FSU survival at high doses for the photons deviated substantially from the LQ fit and thus the value at 10% survival might be discarded.Given that for the organs with parallel architecture, such as the lungs, the FSU survival plays a key role when estimating the probability of radiation-induced toxicity, future studies could benefit from the RBEFSU as a tool for the comparison of different radiation qualities on the NTCP.Overall, these findings capture the unique feature of the model (and, more generally, the added value ABMs can bring to the field of radiation oncology) to be susceptible not only to variations in the values (i.e.different average doses, as equation-based models would be), but also to local changes (i.e.different dose distributions).In conclusion, this work could serve as a framework to investigate the effect of different fractionation schemes and dose distributions on the severity of RILF, while taking into account individualised parameters in the perspective of advancing precision medicine in RT treatments.Future developments might concentrate on taking into account the DNA damage and repair mechanisms in order to estimate the fractions of healthy, apoptotic and senescent cells more accurately.

The Agent-Based Model
Our 3D ABM of RILF in an alveolar segment was developed using the open-source platform BioDynaMo 39 and its implementation was detailed in our previous works 11,12 .The model replicates a small section of the distal airways in the human lungs, namely an alveolar segment, where three layers of alveoli are stacked in a cylindrical shape.Each layer includes 6 alveoli which are modelled as hollow spheres.The surface of each sphere, in turn, is lined by 3 main cell types: the mesenchymal cells (fibroblasts and myofibroblasts) in the outer layer, the macrophages (type 1 and 2, M1 and M2) in the inner layer and the epithelial cells (type 1 and 2) in the middle layer.Following a radiation-induced damage in the epithelial layer, type 2 alveolar epithelial cells (AEC2s) can either become apoptotic or turn into a damaged and eventually senescent state.By secreting a plethora of chemokines and cytokines 40 , the senescent cells can both damage the surrounding healthy AEC2 (where the minimum number of senescent cells needed to damage an healthy one is regulated by the bystander threshold parameter) and trigger an immune response led by M1 and M2.Provided by the capillaries (which are not simulated by our model), M1 and M2 (whose phagocytic fraction and phagocytic index are set with custom parameters) are gathered in the alveoli and remove the senescent cells 41,42 .At the same time, the healthy AEC2 increase their proliferation and differentiation rates in the attempt to replenish both their own and the AEC1 populations 43,44 .Ultimately, the disruption of the healthy epithelium and the secretion of growth factors from the type 2 macrophages concurrently stimulate the expansion of the mesenchymal compartment which in turn saturates the alveolus with extracellular matrix (ECM) components.The severity of the damage at later time points (e.g.months/years after the treatment) is measured via the RILF severity index (RSI), whose definition (proposed in our previous work 11 ) was inspired by the concept of functional subunits (FSU) in the critical volume model of NTCP 31,32 and by the fibrosis index (FI) introduced in the work of Zhou et al. 18 .The RSI reads as follows: Where the two factors describe the average increase in the concentration of the ECM across the whole simulation space and the decrease in the volume of surviving FSUs, represented by the alveoli.Other equations of note (which were described in more detail in our previous work and were used to fit the simulations' output in the Results section) are the following: -Increase in the average ECM concentration across the whole simulation space, early component (where ∆ !"# ,  and  $% represent the saturation value for the ECM increase, the steepness of the sigmoid and dose at 50% of the ∆ !"# , respectively) 19-21 -Increase in the average ECM concentration across the whole simulation space, late component, and RSI 18 -FSU survival probability, derived from the linear quadratic (LQ) 45 and critical volume NTCP models 31,32 (where  &'() is the total number of healthy AEC2 in an alveolus in homeostatic conditions) &-./,!"# = 1 −  3455,!"# = 1 − 5  3455,*655 7 "#$% 4)8 = 1 − 6 3455,*655 7 7 "#$% = 1 − 61 −  9:;9<; % 7 With respect to the previous version, the current model brings with it new features and noteworthy changes that will be detailed in what follows.Among them, modelling of the apoptotic AEC2s has been implemented to increase the accuracy of the simulations.Once an AEC2 has changed its state to apoptotic (either due to irradiation damage or aging), free movement is hampered and the time to removal from the simulation is drawn from a Poisson distribution 46 .At the beginning of a new simulation, each alveolus is initialized with the three cell compartments described before.Cells (i.e. the agents) are distributed at random positions on the alveolar surface and assigned type-specific behaviours.Among them, cell migration plays an important role and has been optimized in our latest model (see our previous work for additional details 11 ).Macrophages and mesenchymal cells travel along spherical arcs in random directions to patrol the alveolar space and maintain homeostatic ECM concentrations, respectively.AEC1 and AEC2, however, are capable of neighbourhood-informed migration and thus move preferentially towards depleted zones in order to repopulate them.Finally, damaged and senescent epithelial cells move in random directions at slower speed than the healthy ones, but might not be able to move at all at times.In our updated model, every cell movement is followed by a check on the cell's final position.If the cell doesn't happen to be located on the spherical surface it belongs to (mainly due to approximation errors when the spherical arc is computed), the cell is translated to the appropriate radial distance while keeping the polar and azimuthal angles fixed.
As in our previous model, the simulations run in a closed cubic space of side 2000 µm that encompasses the whole alveolar segment.A diffusion grid overlaps the simulation space and dissects it in smaller cubic voxels, whose number can be adjusted by the user.The simulated cells can both measure and change the concentration of the substances of the voxel they are located in.Once secreted, substances can both diffuse and decay and, in some cases, be depleted by other substances.Our model simulates, among others, such coupled substances, i.e. matrix metalloproteinases (MMP), tissue inhibitor of metalloproteinases (TIMP) and ECM (interaction mechanisms are outlined in our previous works 11,12 ).The system of reaction-diffusion partial differential equations (PDEs) for all the involved substances is then solved using a forward in time -central in space (FCTS) method with user-defined boundary conditions (BC).In its updated version, both the binding coefficients, the target and the BC for the substances in the model can be specified in the initialization phase.In particular, we set Neumann BC (as opposed to our previous model) with constant value, to ensure a net zero flux and mimic inter-compartmental communication.Moreover, we increased the duration of one time-step for the solution of the system of PDEs in the diffusion grid to 20 seconds.By doing so, we could shorten simulation times while still fulfilling both the stability (i.e.Courant-Friedrichs-Lewy) and positivity conditions.to characterize the inter-alveolar damage spread we implemented a new damage modality.More specifically, the model can now simulate heterogeneous damage patterns whereby the depletion of healthy AEC2 is limited to a portion of (user-defined) alveoli.The damage modality can be set using the parameter "damage_pattern" and the target alveoli specified with their unique index.

The Monte-Carlo Model
TOPAS-nBio is built on top of the Geant4-DNA 47 toolkit and extends the functionalities of the TOPAS platform 48 .It allows MC simulations at the microscopic and nanoscopic scale while simplifying the modelling of biological structures as well as consistent measurements of radiobiological quantities.To investigate the effect of realistic dose distributions on our model of alveolar segment, we reimplemented its structure from scratch using the TOPAS-nBio platform.Our extension consists of 3 new classes that are used to both build the geometry of the cylinder and score the dose delivered to its building blocks.When a new simulation is performed, an envelope for the alveolar segment is firstly built as a cubic box with side 900 µm.A for loop runs then through a vector containing the centers of the 18 alveoli (whose positions can be either hardcoded or provided in a separate file) and smaller containment boxes are generated.The modelling proceeds by instantiating a generic spherical structure with temporary radius for the cells and a second for loop runs through all the alveoli positions.Within its body, the alveoli are built by placing the cells on three spherical surfaces (as described in the previous section).More specifically, parameters for each cell (i.e.position in 3D, colour that indicates the cell type and size) are loaded from external files and, given the inherent repetitive nature of the geometry, the parameterization technique is used to allow for faster rendering and simulation times (see figure 8).As opposed to those in the ABM, cells in this model don't have behaviours, are static and differ only by colour and size.Hits within cells are processed using a custom scorer that provides both the total dose deposited within the cells as well as the average dose for each alveolus.By checking the cells' size, non-epithelial cells are filtered out and the computation of the dose is skipped.For the epithelial cells, however, the total energy deposited by the hitting particles is summed and converted to Gray (Gy) units.At the end of a simulation the total dose delivered to each cell is exported to a text file along with a unique identifier and the identifier of the alveolus the cell belongs to.Moreover, the average dose delivered to (the epithelial cells of) each alveolus is computed and exported together with the alveolus identifier to a second file.Given that our previous work identified the alveoli as the FSU of the lung, the latter can be useful to, for example, characterise the dose distribution on a bigger scale than the cellular one.Finally, simulations can be run using the parameter control system as described in the work by Schuemann et al. 13 .In the parameter file, the total world volume is defined as a cubic box with the same size as the simulation space of the ABM.The alveolar duct is then placed at its center and the envelope material set to vacuum, while we assumed that water is the sole component for the cells.With regard to the particle sources, we tested both an isotropic source located at the center of the alveolar duct and an external beam with one to four coplanar fields (see figure 9).Following the measurements of the average delivered dose per alveolus and outlined in the Results section, the homogeneity of the dose distribution was maximized by the external beam with four fields.The motion of the beam was implemented by using a step function with variable positions on the XZ plane and 90 degrees rotations about the Y-axis, while the cut-off shape was set to "rectangle".Energy and type of the particles to be included in the beam were also set via the parameter file as well as the angular and energetic spreads which were set to zero.Settings for our custom scorer were also adjusted via the parameter control system.

Coupling the models
The coupled model consists of the aforementioned AB and MC models each equipped with import/export functionalities that allow data exchange.In addition to this, overall control of the workflow is assigned to a bash script (see figure 10).The sequence starts with a 1000 steps long run of the ABM that does not require any external data to generate the alveolar segment structure in healthy conditions.We tuned our ABM by running 20 days Figure 8. Alveolar segment model built with TOPAS-nBio.The 18-alveoli structure is inscribed in a cubic box, while each alveolus has a spherical envelope.Cells have spherical shape and differ by size and colour.The outer space is empty, while cells are filled with water.The side of the outer box is 900 µm long, while the diameter of each alveolus measures 260 µm.
long simulations in homeostatic conditions and ensuring that both the number of cells and the average concentration of the substances remained constant and matched previously published results (more details on the tuning procedure as well as the parameters can be found in our previous work 11 ).To ensure the consistency of the results, 10 experiments are performed for each parameter set and the output is used to feed the next step.In particular, we implemented a BioDynaMo standalone operation 39 that can be activated using a parameter and is triggered at the end of each simulation.With it, the cells' parameters together with the position and the diameter of the alveoli and the concentration per voxel of all the simulated substances can be exported.The data files for the cells contain the cells' position, their radius, an RGB colour code that identifies the type and additional type-specific parameters, such as the number of phagocyted cells for the macrophages or the number of steps spent above the bystander damage threshold for the healthy AEC2.Moreover, the identifier of the alveolus each cell belongs to is used to generate a new file for each alveolus.Subsequently, the control script is used to set the total number of fractions as well as the steps for the ABM, the number of repeated experiments (where a new seed is selected for each run) and the target dose (through the number of histories of the MC model).The MC simulation is then run: TOPAS-nBio builds a new structure for the alveolar segment (using the data files mentioned above and the procedure described in the previous section), irradiates it and generates a summary of the dose deposited in each epithelial cell and the average dose per alveolus.It is worth noting that the duration of the irradiation is orders of magnitude shorter than the average cell's behavior time and therefore it assumed that cells don't change their position and/or state within it.The workflow proceeds with a new run of the ABM and custom parameters are used to enable reading of the external damage files.However, the concentration data are read in subsequent runs only as the duration of the preparatory simulation doesn't ensure steady state values for all the substances.Cells' data are read sequentially as well and used to rebuild the alveolar segment structure and assign the type-specific behaviours to the agents.For each epithelial cell the external damage file from the last step is checked and, if a match is found, the LQ model is used to determine the fate of the cell which will either remain healthy or become senescent or apoptotic.More specifically, a random number is drawn from a uniform distribution and compared against the AEC2 LQ survival curve value at the delivered dose (whose parameters were derived from previously published experimental data 50 , as detailed in our previous works 11,12 ).A new simulation is then performed (whose duration is set in the bash script) and, depending on the number of fractions, the output is used to feed either a longer run of the ABM or a new irradiation using the MC model.Preparatory simulations of the ABM used to generate the healthy structures, as well as those in between fractions, were run on a MacBook Pro laptop with a 2.3 GHz Quad-Core Intel Core i5

Figure 2 .
Figure 2. Dose distribution histograms for different setups in the alveolar segment model irradiated with photons.On the x-axis, doses are shown as percentages of the maximum delivered dose and the green curves are the best gaussian fits for the histograms data.Panels (a) and (c) correspond to 4-fields external beams with 4 million and 30 million histories (i.e.particles) respectively.In panel (b) the dose distribution resulting from an isotropic source (i.e. a point source located at the center of the alveolar duct) is reported.

Figure 3 .
Figure 3.Comparison of the major outcomes from the ABM model (black markers, red fit curves) and the ABM-MC model using photons in one fraction (red markers, green fit curves).Panels (a) and (b) show the increase in the average ECM concentration across the whole simulation space after 90 days (early component) and 1200 days (late component).In panel (c) the surviving fraction of the alveoli (FSU) in logarithmic scale at the end of the simulation is reported.The late DECM and the FSU survival are combined to provide the RSI, outlined in panel (d).

Figure 4 .
Figure 4. Comparison of the major outcomes from the ABM-MC model using photons in 1 fraction (black markers, red fit curves) and 5 fractions (red markers, green fit curves).Panels (a), (b) and (c) show the increase in the average ECM concentration (late component), the surviving fraction of the FSU and the RSI, respectively, at the end of the simulation.
the number of senescent neighbours exceeds the bystander threshold (as introduced in the work by McMahon et al.

Figure 5 .
Figure 5.Comparison of the major outcomes from the ABM-MC model using photons in 1 fraction (black markers, red fit curves) and 5 fractions (red markers, green fit curves) with bystander threshold = 2, 1 fraction (green markers) and 5 fractions (blue markers) with bystander threshold = 1.Panels (a), (b) and (c) show the increase in the average ECM concentration (late component), the surviving fraction of the FSU and the RSI, respectively, at the end of the simulation.

Figure 6 .
Figure 6.Dose distribution histogram for the alveolar segment model irradiated with protons.The dose was delivered using 4 fields from an external beam and 2 million histories in total.On the x-axis, doses are shown as percentages of the maximum delivered dose and the green curve is the best Gaussian fit for the histogram data.

Figure 7 .
Figure 7.Comparison of the major outcomes from the ABM-MC model using protons (black markers, red fit curves) and photons (red markers, green fit curves), both in 1 fraction.Panels (a), (b) and (c) show the increase in the average ECM concentration (late component), the surviving fraction of the FSU and the RSI, respectively, at the end of the simulation.Panel (c) provides also the RBE values for the surviving fraction of the FSUs at 50%, 37% and 10% survival.

Figure 9 .
Figure 9. Frontal and lateral irradiation of the alveolar duct in TOPAS-nBio with 1 keV gamma particles.

Figure 10 .
Figure10.Flowchart of the coupled ABM -MC model workflow.BioDynaMo runs the ABM for 1000 steps, generates the healthy structure and exports the alveoli positions and the cells data.The number of fractions, the number of steps for the ABM and the irradiation parameters are then set into the control script which triggers the first run of the MC model using TOPAS-nBio.The structure is loaded from the data exported during the previous step and the alveolar segment is irradiated.Consequently, the dose distribution is exported and the data used as input for the ABM to determine the cells' fate.A new run of the ABM is then performed and the segment structure together with the cells' state and their position are exported.If a new fraction is to be delivered, the data are fed into the MC model, otherwise the workflow terminates.