A spatial measure-valued model for radiation-induced DNA damage kinetics and repair under protracted irradiation condition

In the present work, we develop a general spatial stochastic model to describe the formation and repair of radiation-induced DNA damage. The model is described mathematically as a measure-valued particle-based stochastic system and extends in several directions the model developed in Cordoni et al. (Phys Rev E 103:012412, 2021; Int J Radiat Biol 1–16, 2022a; Radiat Res 197:218–232, 2022b). In this new spatial formulation, radiation-induced DNA damage in the cell nucleus can undergo different pathways to either repair or lead to cell inactivation. The main novelty of the work is to rigorously define a spatial model that considers the pairwise interaction of lesions and continuous protracted irradiation. The former is relevant from a biological point of view as clustered lesions are less likely to be repaired, leading to cell inactivation. The latter instead describes the effects of a continuous radiation field on biological tissue. We prove the existence and uniqueness of a solution to the above stochastic systems, characterizing its probabilistic properties. We further couple the model describing the biological system to a set of reaction–diffusion equations with random discontinuity that model the chemical environment. At last, we study the large system limit of the process. The developed model can be applied to different contexts, with radiotherapy and space radioprotection being the most relevant. Further, the biochemical system derived can play a crucial role in understanding an extremely promising novel radiotherapy treatment modality, named in the community FLASH radiotherapy, whose mechanism is today largely unknown.


Introduction
Radiotherapy is, today, a widely used treatment against cancer [Thariat et al., 2013].Conventional radiotherapy is based on X-rays, i.e. photons, but in the last decades constantly increasing attention has been devoted to advanced radiotherapy treatment with ions [Durante and Paganetti, 2016].Ion beams have many essential features making them preferable compared to photons, related mostly to the extremely localized energy released in tissues which can lead to a superior biological effect than X-rays.The effect of radiation on biological tissue has been studied by the community over the last decades, and DNA is believed to be the most sensitive target to radiation so DNA damage is the most relevant biological vehicle that leads to cell killing induced by radiation, [Durante and Loeffler, 2010].Despite the potential superiority of hadrons in theory, additional research is crucial to incorporate this treatment modality into clinical practice fully.One of the primary obstacles to the widespread use of hadrons is in fact accurately estimating the biological effect caused by radiation, a crucial aspect to account for in order to prescribe the best possible treatment.Mathematical models have thus been developed over the years to understand and accurately predict the biological effect of ions on biological tissue, [Bellinzona et al., 2021, Hawkins, 1994, Hawkins and Inaniwa, 2013, Kellerer and Rossi, 1974, Herr et al., 2015, Pfuhl et al., 2020, Cordoni et al., 2021], focusing on the DNA damage Double Strand Breaks (DSB).Such mathematical approaches focus on developing models that describe the formation, evolution, and interaction of DSB, with the final goal of predicting the probability that a certain cell survives radiation.
To date, very few models in the context of radiotherapy have a robust mathematical and probabilistic background even if the community widely acknowledges stochastic effects play a major role in the biological effect of radiation.In fact, despite the early development of stochastic models for the description of the kinetic repair of radiation-induced DNA damages, [Sachs et al., 1990, Albright, 1989], the radiobiological community soon drifted to developing deterministic models of damage repair assuming Poisson fluctuations of the number of damages around the average values, [Bellinzona et al., 2021].This type of modelization is strictly linked to a linear-quadratic description of the relation between the logarithm of the cell-survival probability and the absorbed dose, a physical quantity that describes the energy deposited by the particles over the mass of the biological tissue traversed by the particles.Although such models provide a fast way to assess the cell survival fraction, which is a key aspect for a use of such models in clinical applications in which the run time of a model is extremely relevant, in recent years, the need began to be felt for more robust modeling from a purely probabilistic point of view.From a mathematical point of view, the Generalized Stochastic Microdosimetric Model (GSM 2 ) recently introduced in [Cordoni et al., 2021, Cordoni et al., 2022b, Cordoni et al., 2022a], appears to be a general mathematical model, that includes several relevant stochastic effects emerging in the creation, repair, and kinetics of radiation-induced DNA-damages, [Cordoni et al., 2022a].GSM 2 considers two types of DNA lesions X and Y, representing respectively lesions that can be repaired and lesions that lead to cell inactivation.In the current context, the specific exact meaning of sub-lessons is left unspecified.This is because there are mainly two different ways that cells can be affected by radiation.One is the creation of DNA Double-Strand Breaks (DSB) from two Single-Strand Breaks (SSB), and the other is the formation of chromosome abnormalities from pairs of chromosome breaks, [Kellerer and Rossi, 1974].Both of these mechanisms are important in understanding how cells respond to radiation and can be described by the model developed in this work.
The present paper aims at extending GSM 2 to include a spatial description allowing for reaction rates that depend on the spatial position, lesion distance, and density.In fact, a true spatial distribution of DNA damage inside the cell nucleus is today almost completely missing in existing models.At the same, it is widely known in the community that the spatial distribution of DNA damages strongly affects the probability that a cell repairs the induced damages, so spatial stochastic plays a major role in the modelization of the repair of radiation-induced DNA damages.We will thus model the spatial distribution of DNA damages as a general measure-valued stochastic particle-based system, characterizing existence and uniqueness as well as some relevant martingale properties that, as standard, will play a crucial role in the derivation of the large system limit.Stochastic particle-based systems have been long studied in the mathematical community, [Bansaye and Méléard, 2015, Popovic et al., 2011, Pfaffelhuber and Popovic, 2015].Recently, a lot of attention has been devoted to studying the spatial non-local stochastic particlebased system, [Bansaye and Tran, 2010, Fontbona and Méléard, 2015, Champagnat and Méléard, 2007, Ayala et al., 2022], where a measure-valued stochastic process describes the population.The population can interact according to a specific rate leading to either the creation or removal of individuals.Mathematically, these systems are described by Stochastic Differential Equations (SDE) driven by Lévy-type noises that besides a diffusive component include jump operators in the form of Poisson random measure, that account for the creation and removal of individuals from the population, [Bansaye and Méléard, 2015].
Most results focus on birth-and-death spatial processes, meaning that at each time, at most, a single individual can be born or die.In this setting, pairwise interactions, involving either the creation or removal of more than one individual, are not allowed.Such interactions are relevant in many biological and chemical applications so a general mathematical theory that extends and generalizes the birth and death process is greatly desirable.Recently, few papers appeared that include pairwise reactions, [Isaacson et al., 2022, Lim et al., 2020, Popovic and Veber, 2020], but none of these deal with the existence and uniqueness or regularity of results for such systems, where instead the focus is mostly on the large-population limit.
The developed model includes some key features that make the mathematical treatment of the spatial model non-trivial.First, given the application considered, where clusters of DNA lesions are more difficult to repair by the cell and has been recognized as one of the main factors that lead to cell inactivation in radiobiology, [Kellerer and Rossi, 1974], we will include pairwise interaction and second-order rates, meaning that a couple of lesions can interact to create an unrepairable lesion that inactivates the cell.It is worth stressing that, already, many existing radiobiological models include parameters to account for the interaction of damages, [Hawkins, 1994, Sato and Furusawa, 2012, Hawkins and Inaniwa, 2013, Bellinzona et al., 2021, Cordoni et al., 2021, Cordoni et al., 2022a].Still, none have a true mathematical spatial formulation and often rely on fixed domains to limit pairwise interaction within a certain distance neglecting nonetheless any true spatiality inside a fixed domain.The latter approach can be restrictive and may lead to overfitting with the inclusion of unnecessary parameters.One of the proposed model's main strengths is that it considers a true spatial distribution of lesions, allowing for true pairwise interaction that can depend on the distance between lesions, which is a novel and important aspect of the model.As mentioned, the existence and uniqueness results for second-order systems are rare and yet a general theory is missing, so the derived results represent a novelty both from a radiobiological as well as a mathematical perspective.
Another key aspect of the studied model is that we explicitly consider the case of protracted irradiation, that is we consider the situation in which a continuous radiation field induces a random number of lesions in the cell.Such a situation is non-trivial from a purely mathematical perspective as the generation of a random number of damages must be considered.Nonetheless, it is extremely relevant to include protracted irradiation in a biological model since it allows us to better estimate the kinetics repair of radiation-induced damage with benefits both in radioprotection and clinical application.Existing radiobiological models account to a certain extent for the protracted irradiation case, [Inaniwa et al., 2013, Manganaro et al., 2017], and a similar rate has already been considered in [Cordoni et al., 2021] in a non-spatial setting.In generalizing the setting to include a spatial description, we need to describe a spatial energy deposition pattern within the cell nucleus.We will build such a theory properly generalizing some existing approaches.Nonetheless, a robust theory to account for the spatial formation of radiation-induced DNA lesions is missing and future efforts will be made to derive such a theory.In fact, another relevant aspect of the studied model is the computation of the spatial distribution of radiation-induced DNA damages.Such a description is relevant both in the protracted case and in the instantaneous irradiation case as it describes the initial damage distribution.In particular, the initial damage distributions, X 0 , and Y 0 can be computed using different methods.One possible approach is the use of Monte Carlo (MC) track structure codes, [Nikjoo et al., 2006], to simulate the passage of charged particles in biological tissue and their energy release and to thus estimate the DNA damage distribution caused by radiation.MC track structure codes have been shown to be effective in accurately characterizing DNA damage formation, [Goodhead, 1994, Ottolenghi et al., 1995, Cucinotta et al., 2000, Chatzipapas et al., 2022, Kyriakou et al., 2022, Zhu et al., 2020, Thibaut et al., 2023], however, once the initial damage distribution is computed, in order to assess the cell survival probability, these models typically neglect the spatial distribution of damages and focus on average values described by Ordinary Differential Equations (ODE).The model developed in this research is unique in that it is able to fully exploit the accuracy of the spatial distribution of damages as predicted by MC track-structure codes.Further, since MC track structure codes simulate all the energy released by a particle along its path, which is referred to in the community as track, as well as all secondary energy releases associated with the original particle, the computational time is extremely demanding.To shorten the computational time, a threshold on energy release can be applied.so that all events that release lower than a certain energy are neglected and incorporated into the deposition that has originated it.Such an approach is called condensed history MC, [Agostinelli et al., 2003], and it provides accurate results of energy deposition at a lower computational time compared to MC track structure codes.
An alternative approach to MC track-structure codes would be to develop an analytical model for DNA damage formation and distribution.Such a model would be less accurate, but less computationally expensive.Currently, the Local Effect Model (LEM), [Friedrich et al., 2012], and the Microdosimetric Kinetic Model (MKM), [Hawkins, 1994, Kase et al., 2007], which are the only models used in Treatment Planning Systems (TPS), take into account the spatial distribution of the absorbed dose without MC codes.These models are based on the Amorphous Track (AT) model, [Kase et al., 2007], which parametrizes the dose distribution around a track of a particle.However, to eventually assess the cellsurvival probability, both models make extensive use of fixed domain so that a true spatial distribution of damages is again neglected.It is worth further stressing that, although the AT model can be used to compute the imparted dose in a fast fashion, it is based on several assumptions, such as the so-called track-segment condition (assumes tracks do not lose energy when traversing the cell nucleus) and uniform radiation fields (cylindrical geometry is often assumed for the cell nucleus and tracks are perpendicular to the cell nucleus).Being the former approach based on track-structure codes beyond the scope of the present work, we will focus on the latter one.Nonetheless, future research will focus on developing a comprehensive analytical model for DNA damage formation that accurately describes energy and spatial stochasticity.
The developed spatial DNA-damage model is expected to play a relevant role in the modelization of a novel radiotherapeutic technique, named in the community FLASH radiotherapy [Favaudon et al., 2014].Up to 2014, most of the mechanisms happening in the interaction of radiation with biological tissue were believed to be known, and therefore in the last two decades models focused on specific applied aspects and general and robust mathematical theories were strongly believed to be unnecessary given the overall understanding of the problem at hand.Starting in 2014 a series of ground-breaking papers [Favaudon et al., 2014, Montay-Gruel et al., 2017, Vozenin et al., 2019] finally showed that an increase in the rate of delivery of ions at Ultra-High Dose Rates (UHDR), namely a high amount of energy released in a small fraction of time, spares healthy tissue and yet maintains the same effect on the tumor.This effect, which was completely unexpected, represents the final goal of any radiotherapeutic treatment [Esplen et al., 2020, Griffin et al., 2020].All available models brutally failed to predict such peculiar effects and hundreds of publications have appeared recently trying to understand the mechanism at the very core of the FLASH effect [Labarbe et al., 2020, Abolfath et al., 2020a, Abolfath et al., 2020b, Liew et al., 2021, Petersson et al., 2020]; many physical explanations and mathematical models have been proposed in the last five years, but up to date, no model is believed to be capable neither of predicting nor understanding the origin of the FLASH effect [Weber et al., 2022].Two facts are today believed to be at the very core of the FLASH effect: (i) this effect has its origin in a spatial interaction of ions, that involves besides the physics of ions and biology, also chemistry, and (ii) nonlocal effects of ions while traversing different cells and how their spatial interaction affects the overall chemical and biological environment.The model developed in the present research could, when coupled with an adequate description of the chemical environment affected by radiation, help unravel the mechanism behind the FLASH effect.In order to do that, we couple the spatial model with a reaction-diffusion equation that describes the evolution of the chemical system.In the following treatment, we will not specify a particular chemical description.This is because there are several possible choices, and the choice depends on the specific application.The chemical stage can be broadly divided into two stages: (i) the heterogenous chemical stage and (ii) the homogeneous chemical stage.The former is characterized by a heterogeneous spatial distribution and occurs immediately after particles hit the cell nucleus, i.e. between 10 −12 seconds and 10 −6 seconds.The latter is characterized by a homogeneous distribution and occurs after the heterogeneous stage, i.e. between 10 −6 seconds to 10 0 seconds.Mathematically, this means that while the homogeneous stage can be described by ODEs that characterize the time evolution of the concentration of chemicals within the domain, [Labarbe et al., 2020, Abolfath et al., 2020a], the heterogeneous stage requires advanced mathematical tools since the system is highly nonlinear and the reactions occur locally.Therefore, most of the literature is devoted to the development of simulation codes, [Clifford et al., 1986, Pimblott et al., 1991, Boscolo et al., 2020, Ramos-Méndez et al., 2020].To date, there is no general mathematical formulation via local reaction-diffusion PDEs that exists in literature, even though it could provide an accurate representation of the system and fast computational time.Also, general results of well-posedness that cover relevant non-local chemical systems are not available in the literature due to the highly complex mathematical formulation needed.For these reasons, a deep mathematical study of such a system is left for future research.Relevant systems to study could include a spatial non-local version of the homogeneous chemical systems presented in [Labarbe et al., 2020, Abolfath et al., 2020a] or an analytical formulation in terms of highly dimensional non-local reaction-diffusion PDEs of [Boscolo et al., 2020, Ramos-Méndez et al., 2020].In this paper, we will instead limit ourselves to considering a general reaction-diffusion PDE with coefficients satisfying rather general assumptions that could in principle include several relevant examples.In particular, almost any chemical description includes bimolecular reactions, meaning that the resulting PDE has quadratic terms.A significant effort has been made in the literature to study reaction-diffusion PDEs under the most general assumptions on the coefficients in order to include as many examples as possible, [Pierre, 2010].In this direction, a mass control assumption has been typically seen as a general condition that allows obtaining the existence and uniqueness of the equations in many cases dropping the standard global Lipschitz condition.For the sake of simplicity, we will consider the system introduced in [Fellner et al., 2020], which allows for quadratic growth of the coefficients, adding random discontinuity due to the effect of radiation.Therefore, we prove the existence and uniqueness of a stochastic particlebased system coupled with a reaction-diffusion system with random jumps.We will show later that this system allows generalizing some relevant systems.A future effort will be devoted to the study of general local reaction-diffusion systems similar to the ones studied in [Isaacson et al., 2022].
At last, we will also characterize the large-system behavior.Such a limiting system can be obtained with standard arguments proving the tightness of the measure and identifying the listing process, which can be proved to admit a unique solution.Although the techniques are standard, the result is new in literature since no stochastic system allowing for pairwise interaction and creation of random numbers of particles has never been studied before.In fact, The resulting governing equation can be useful to study the behavior of the system at high doses, where the number of damages within the cell increases arbitrarily.It is worth stressing that, the high-dose case is recognized to be non-trivial and most of the existing models fail to predict the behavior of the system at high doses.For this reason, often, correction terms are included in the model to better math experimental data, [Bellinzona et al., 2021].
The main contributions of the present paper are: (i) to provide a general mathematical description of a spatial model governing the formation and kinetics of radiation induced-damages; (ii) to study the well-posedness of a measure-valued stochastic particle system with pairwise interaction and random creation of damages; (iii) to propose a multi-scale model to couple biology and chemistry that could possibly describe the FLASH effect; (iv) to study the large-system limit of the system with pairwise interaction and protracted irradiation.

The microdosimetric master equation
The main goal of the present section is thus to introduce the classic setting for the GSM 2 [Cordoni et al., 2021, Cordoni et al., 2022b].GSM 2 models the time-evolution of the probability distribution of the number of lethal and sub-lethal lesions denoted by (X(t); Y(t)), where X and Y are two N−valued random variables counting the number of the lethal and sub-lethal lesion, respectively.In the following, we will consider a standard complete filtered probability space " Ω; F; (F t ) t≥0 ; P " satisfying usual assumptions, namely right-continuity and saturation by P-null sets.
We thus assume that a sub-lethal lesion X can undergo three different pathways: (i) at rate r a sublethal lesion is repaired, (ii) at rate a a sub-lethal lesion is left unrepaired by the cell and thus it becomes a lethal lesion and (iii) at rate b two sub-lethal lesion form a cluster that cannot be repaired by the cell and thus become a lethal lesion.Any lethal lesion leads to cell inactivation.These three pathways can be summarized as follows Denoting by p(t; y ; x) = P ((Y(t); X(t)) = (y ; x)) ; the probability to have at time t exactly x sub-lethal lesion and y lethal lesions, following [Cordoni et al., 2021], we can obtain the microdosimetric master equation (MME) x − 1)bp(t; y ; x)] ; p(0; y ; x) = p 0 (y ; x) ; (2) where we have denoted the creation operator as In [Cordoni et al., 2022b] it is derived a closed-form solution for the survival probability as predicted by the MME (2), defined as the probability of having no lethal lesions Y. Further, GSM 2 is closely connected with one of the most used radiobiological models to predict the survival probability of cell nuclei when exposed to ionizing radiation, that is the Microdosiemtric Kinetic Model (MKM) [Hawkins, 1994].The main equations of the MKM describe the time-evolution of the average value ȳ , resp.x, of the number of lethal, resp.sub-lethal, lesions, and are given by (3) The model further assumes that ȳ is the average of a Poisson random variable so that by describing the average values we have complete knowledge of all the moments.
To obtain a suitable analytical solution to the equations (3), it is often assumed that (a + r)x >> 2bx, so that above equation is reduced to This highlights why in the high dose case the MKM must be corrected including additional terms.In fact, even if it is typically true that (a + r) >> 2b, at sufficiently high doses, the number of lesions x increases so that (a + r)x does not dominate anymore 2bx and therefore the omission of the term in equation ( 4) becomes non-negligible.
Further, it has been shown in [Cordoni et al., 2021], that the average of the MME coincides with the MKM equations (3) under a suitable mean-field assumption, that us which in turn coincides exactly with the requirement that X follows a Poisson distribution.It has thus been shown in [Cordoni et al., 2022a] that the GMS 2 is able to give a more general description of many stochastic effects relevant in the formation and repair of radiation-induced DNA lesions that play a crucial role in estimating the surviving probability of a cell nucleus.
It can be further shown, [Weinan et al., 2021, Bansaye andMéléard, 2015], that equation (2) describes the time evolution for the probability density function associated with the following stochastic differential equation Above in equation ( 5), N Y (ds; dz) and N X (ds; dz) are two independent Poisson point measure with intensity ds dz on R + × R + , see, e.g.[Applebaum, 2009].The main of the present work will be to provide a spatial description of the SDE (5) so that X and Y are replaced by random measures.

On the initial distribution
In order to later generalize the initial damage distribution, we introduce in the current section the distribution introduced in [Cordoni et al., 2021, Cordoni et al., 2022b, Cordoni et al., 2022a].For a detailed treatment, we refer the interested reader to the mentioned papers or to [Bellinzona et al., 2021].
Among the most powerful approaches to describe the formation of DNA lesions is using microdosimetry [Zaider et al., 1996].Microdosimetry is the branch of physics that investigates the energy deposition in domains comparable to cell nuclei, that is of the order of some microns.At that scale, energy deposition is purely stochastic, so the main object used in microdosimetry are random variable and their corresponding distributions.Over the years many models have been developed based on microdosimetric principles, [Kellerer andRossi, 1974, Zaider et al., 1996], and both the MKM and GSM 2 assess the formation of DNA lesions using microdosimetry, [Hawkins, 1994, Bellinzona et al., 2021, Cordoni et al., 2021, Cordoni et al., 2022b].
The main microdosimetric quantity of interest from the point of view of radiobiological models is the specific energy z, [Zaider et al., 1996].The specific energy z is the ratio between energy imparted by a finite number of energy depositions " over the mass m of the matter that has received the radiation, that is z = " m : The stochastic nature of " implies that also z is inherently stochastic.The single-event distribution f 1 (z) of energy deposition on a domain, [Zaider et al., 1996], is the probability density distribution describing the energy deposition due to a single event, typically a particle traversing the domain.Such distribution is associated with a random variable Z that describes the specific energy imparted on a certain domain of mass m.The average values of the random variable Z, referred to in the literature as fluence-average specific energy, that is the mean specific energy deposition, is typically denoted in literature as z F .By additivity property, the specific energy distribution resulting from tracks can be computed convolving times the single event distribution, [Zaider et al., 1996].Therefore, the distribution f of the imparted energy z is computed iteratively as We denote by p e ( |D; z F ) a discrete probability density distribution denoting the probability of registering events.Typically such distribution is assumed to be dependent on the total dose absorbed by the mass and the fluence of the incident particles.The standard assumption is that, since events are in a microdosimetric framework assumed to be independent, the distribution p e is a Poisson distribution of average D zF so that we have Therefore, microdosimetry postulates that the actual energy deposition on a certain domain can be obtained via the multi-event specific energy distribution At last, given a certain specific energy deposition z by events, the induced number of lethal and sub-lethal lesions is again a random variable, with a discrete probability density function denoted by p.In general the average number of lethal, resp.sub-lethal, lesions is assumed to be a function of z, namely »(z), resp.-(z).Again, by independence on the number of created lesions, such distribution is assumed to be a Poisson distribution.Overall, the probability of inducing x sub-lethal and y lethal lesions can be computed as, [Cordoni et al., 2021], or assuming Poissonian distributions for suitable functions »(z) and -(z).These quantities summarize the free-radical reactions that result in a lesion.It is a function of the type of ionizing particle, details of the track structure, radical diffusion, and reaction rates, the point in the cell cycle, and the chemical environment of the cell.In the following, we will explicitly model these functions so that they depend on chemical concentration.The classical assumption, which has been also considered in [Cordoni et al., 2021], is to assume such functions to be linear in z.Notable enough, it has been shown in [Cordoni et al., 2022b] that, also assuming a Poissonian distribution on both p e and p, the resulting discrete probability density function ( 8) is not a Poisson distribution; as a matter of a fact, it has been shown to be a microdosimetric extension of the so-called Neyman distribution, [Neyman, 1939], which is a well-known distribution in radiobiological modeling to treat the number of radiation-induced DNA damages.To have a better grasp on the distribution (7), it can be described by a stochastic chain of interconnected events: (i) given a certain dose D and fluence average specific energy z F , a given random number of events is registered in a cell nucleus; then (ii) such events deposits a certain random specific energy z = z 1 + • • • + z .At last, (iii) the specific energy deposited z induces a random number of lethal and sub-lethal lesions y and x.
3 The spatial radiobiological model The current section aims at generalizing the radiobiological model as introduced in Section 2 to consider a spatial measure-valued process.Consider a closed bounded regular enough domain Q ⊂ R d , d ≥ 1, which should represent a cell nucleus.We assume that Q has a smooth boundary @Q, and denote by n(q) the outward normal direction to the boundary @Q at the point q.
We consider two possible types of DNA damage, S = {X; Y}, where X denotes sub-lethal lesions and Y are lethal lesions.We assume sub-lethal and lethal lesions can undergo three different pathways, a, b, and r , as introduced in Section 2.
We consider thus a process that lives in the state space encoding the i-th lesion position q i and type s i .For a metric space E, we define by M F (E) the space of finite measure over E, endowed with the weak topology; given a regular enough function Also, we denote by M(E) the space of point measure over E, defined as In general, in the following, we will often consider either E = P or E = Q; if no confusion is possible we will omit the subscript in the scalar product.
Fix a finite time horizon T < ∞, for t ∈ [0; T ], we define the concentration measure of lesion at time t, as the total number of lesions at time t.We further denote by X (t) and Y (t) the marginal distributions Analogously to the previous notation, N X (t), resp.N Y (t), denote the total number of lesions of type X, resp.Y, at time t.
Besides lesion concentration we will often use a vector listing all lesions in the system; thus, given a system state (t), we denote by H( (t)) := " `q1;X (t); X ´; : : : ; " q N X (t);X (t); X " ; `q1;Y (t); Y ´; : : : ; " q N Y (t);Y (t); Y " ; 0; : : : the position and type of all lesions in the system at time t.It is worth stressing that, since lesions of the same type are indistinguishable, the chosen ordering is arbitrary and there is no ambiguity in H( (t)).

The model
Each lesion i, characterized by its position and lesion type P i = (q i ; s i ), can move and undergo three different pathways.Such rates can be described by the system and can be characterized as follows: (i) -Repair each lesion in the class of sub-lethal lesions X can repair at a rate r : Q × R → R + ; r `q; ˙Γr q ; ¸´; that depends on the spatial position of the i−th lesion and on the concentration of the system.A sub-lethal lesion that repairs disappear from the system.
The repair rate r is associated to a Poisson point measure The index i ∈ N 0 gives the sampled lesion in X to repair.The corresponding intensity measure associated with N r is We denote with Ñr the compensated Poisson measure defined as Ñr (ds; di; d") := N(ds; di; d") − -r (ds; di; d") : (ii) -Death each lesion in the class of sub-lethal lesions X can die at a rate a : Q × R → R + ; a `q; ˙Γa q ; ¸´; that depends on the spatial position of the i−th lesion and on the concentration of the system.A sub-lethal lesion that dies generates a lethal lesion Y at a new position q according to the probability distribution m a (q), q ∈ Q.
The death rate a is associated to a Poisson point measure The index i ∈ N 0 gives the sampled lesion in X to die and become a lethal lesion in Y in position q sampled from m a (q|q 1 ), q ∈ Q.The corresponding intensity measure associated with N a is We denote with Ña the compensated Poisson measure defined as (iii) -Pairwise interaction two lesions in the class of sub-lethal lesions X can interact at a rate that depends on the spatial position of the (i 1 ; i 2 )−th lesions and on the concentration of the system.Two sub-lethal lesions that interact can either (i) die with probability p, generating a lethal lesion Y at a new position q according to the distribution m b (q), q ∈ Q, or (ii) repair with probability 1 − p and disappear from the system.The probability p depends also on the positions of the sampled lesions, namely The pairwise interaction rate b is associated to two Poisson point measure The index i = (i 1 ; i 2 ) ∈ N 0 × N 0 gives the sampled lesions in X to either become a lethal lesion in Y in position q sampled from m(q|q 1 ; q 2 ), q ∈ Q, or repair and be removed from the system.The corresponding intensity measures associated with N b;p and N b;1−p are We (iv) -Spatial diffusion each lesion of type X and Y moves around the domain Q with diffusion term and drift term In the following, we will also denote with S + `Rd ´the space of symmetric non-negative d × d matrices.
To describe lesion motion we introduce a countable collection of standard independent Brownian motion `W n;X (t) ´n∈N and `W n;Y (t) ´n∈N on R d .Brownian motion is assumed to reflect with normal derivative at the boundary of the domain Q.In particular, denote by T n and T n+1 two successive jump times of the process , and assume that at time T n we have ), the number of lesions remains constant so that the process is solely subject to the diffusive component.Thus, for any t ∈ [T n ; T n+1 ) each lesion evolves according to the following SDE with reflection at the boundaries where we denoted by dW (t) the integration in the sense of Itô.
In the following, we will consider a filtered and complete probability space ing standard assumptions, namely right-continuity and saturation by P-null sets.In particular, (F t ) t∈R+ is the filtration generated by the processes defined in (i) − (ii) − (iii) − (iv ) as well as a M × M−valued initial distribution = ( X 0 ; Y 0 ).Remark 3.1.Notice that, compared to the original interaction rates as introduced in [Cordoni et al., 2021], we included in the present version of the model a further possible pathway, namely This is done since, as noted in early versions of advanced radiobiological models [Sachs et al., 1992], pairwise interaction of damages, can result also in correct repairs; such a process is called for instance in [Sachs et al., 1992] as complete exchange.△ Through the paper we will assume the following hypothesis to hold: Hypothesis 3.2.1. Jump components: (2.i) the repair rate r is uniformly bounded over compact subsets, that is for N ≥ 0 it exists r such that sup (2.ii) the death rate a satisfies a linear growth condition, that is, there exists a positive constant ā such that, for all q ∈ Q it holds (2.iii) the pairwise interaction rate b satisfies a linear growth condition, that is, there exists a positive constant b such that, for all q 1 and q 2 ∈ Q, it holds (2.iv) the pairwise interaction death p is a probability, that is, for all q 1 and q 2 ∈ Q, it holds p(q 1 ; q 2 ) ∈ [0; 1] ;

Diffusive components:
(2.i) there exist positive constants L X and L Y such that, for any q 1 , q 2 ∈ Q, it holds 3. Kernel components: is continuous and uniformly bounded, that is, there exists a constant Γh such that sup q∈Q sup ( q;s)∈Q×S Γ h q (q; s) < Γh < ∞ ; (3.ii) the sampling measure m(q|q 1 ; q 2 ) is a probability measure, that is In the following we consider the next class of cylindrical test functions: for F ∈ C 2 b (R × R), that is, F is bounded with continuous and bounded second order derivative, and for f X ; f Y ∈ C 2 0 (Q), that is f X and f Y are continuous with bounded second order derivative in the domain variable Q, satisfying In the following, we will denote by @ @x , resp.@ @x , the derivative with respect to the first argument, resp.the second argument, of the function F .Also, ∇, resp.∆, resp.T r , resp.Hess, denotes the gradient with respect to the space variable q, resp.the Laplacian operator with respect to the space variable q, resp.the trace operator, resp.the Hessian matrix.Cylindrical functions (13) are a standard class generating the set of bounded and measurable functions from M × M into R. Remark 3.3.(i) A natural choice for the kernel Γ would be to assume that only nearby mass affects the overall rate; in such a case we have, for q ∈ Q, Γ q (q; s) := ½ {|q− q|<›} (q; s) : Therefore, only lesions that are at most distant › from the position q where the reaction happens to participate in the reaction.
(ii) Regarding the pairwise interaction rate b, it is natural to assume that b depends only on the separation distance between two lesions, that is, it exists a function b : R → R + ; b (q) ; such that b (q 1 ; q 2 ) = b (q 2 ; q 1 ) = b (|q 1 − q 2 |) : Further, it is natural to assume that the closer two lesions are, the more likely to interact they are.Therefore, relevant choices for the rate b are for instance a step interaction rate b (q) := b½ {|q|<"} ; Whereas the former rate ( 14) models the case where only lesion closer than " can interact, the latter rate (15) considers that the rate of interaction decreases exponentially as the lesions are more distant from each other.At last, as noted in [Kellerer and Rossi, 1974], enhanced short-range interaction can be modelled using b (q) := b1 for suitable constants, so that the rate of interaction declines fast but still has a fat tail at larger distances.Similarly, it is reasonable to assume that also p depends on the distance between the interacting lesions.
(iii) Possible choices for the sampling measure m(q|q 1 ; q 2 ) would be to assume that, whenever two sublethal lesions at (q 1 ; q 2 ) reacts, a lethal lesion is created randomly on the segment connecting q 1 and q 2 .For instance, the following probability distributions can be considered, △ With the previous notation, we can thus introduce the following weak representation for the spatial radiation-induced DNA lesion repair model, given f X and f Y ∈ C 2 0 (R), we have ×N b;p (ds; di 1 ; di 2 ; dq; d" 1 ; d" 2 ) : (17) Definition 3.3.1.We say that (t) = ( X (t); Y (t)) as defined in equations 9-10 is a spatial radiationinduced DNA damage repair model if = ( (t)) t∈R+ is (F t ) t∈R+ −adapted and for any f X and f Y ∈ C 2 0 (Q) equation ( 17) holds P−a.s.
The above process is characterized by the following infinitesimal generator where LF (f X ;f Y ) ( ) is the infinitesimal generator of the reaction terms, whereas L d F (f X ;f Y ) ( ) is the infinitesimal generator of the diffusive part of equation ( 17).
In particular, we have It further holds that Regarding the infinitesimal generator of the reaction terms it holds where we have denoted by Q2 := Q 2 \ {(q 1 ; q 2 ) : q 1 = q 2 } :

Stepwise construction of the process
In the present Section, we provide a step-wise construction of the process.Such construction, besides being relevant from a theoretical point of view, is particularly important in implementing a simulation algorithm for the process defined in the previous section.Notice that, using assumptions 3.2, we have that the rate at which a, b and r happen is bounded in the uniform norm.Thus, between the occurrence times of the jump components, each lesion moves according to the diffusive generator D X and D y .

Well-posedness and martingale properties
We can in the present Section the existence and uniqueness of solutions to the above-introduced model.
Theorem 3.4.Let X 0 and Y 0 two independent random measures with finite p−th moment, p ≥ 1, that is it holds Then, under Hypothesis 3.2, for any T > 0, there exists a pathwise unique strong solution to the system (17) in D ([0; T ]; M × M).Also, it holds In particular, the process in Definition 3.3.1 is well-defined on R + .
Proof.Since the jump times are isolated, the construction of X (t) and Y (t) can be done pathwise inductively along the successive jump times.In particular, denote by T m and T m+1 two successive jump times of the process , and assume that at time T m we have N X (T m ), resp.N Y (T m ), lesions of type X, resp.Y.As noted above, for t ∈ [T m ; T m+1 ) the number of lesions remains constant so that the process is solely subject to the diffusive component as described in equation ( 12).Using Hypothesis 3.2, conditionally on F Tm , equation ( 12) can be seen as a purely diffusive SDE with globally Lipschitz coefficients on R d×N X (Tm) × R d×N Y (Tm) , so that the process `XiX (t); Y iY (t) ´iX=1;:::;N X (Tn); iY=1;:::;N Y (Tn) ; admits a unique strong solution for t ∈ [T m ; T m+1 ).
Define then, for n ≥ 0, fi X n := inf{t ≥ 0 : ˙1; X (t) ¸≥ n} ; and set for short f X n := t ∧ fi X n .We can construct a solution algorithmically in [0; T ).For t ≥ 0, noticing that the number of lesions in X can only decrease, we have, Regarding Y, we have that, Taking the expectation in equation ( 25) and using estimate (24) together with we have that, for some C > 0 that can take possibly different values, ¸p« : (26) From Gronwall lemma it thus follows that it exists C > 0 depending on p and T but independent of n, such that E sup In fact, if that was not the case, we can find T 0 < ∞ such that This would in turn yields which contradicts equation ( 27).A similar argument holds for X.Using thus Fatou's lemma we can let n → ∞ as proving thus (23).
At last, since the above claim holds also for p = 1, we have that so that the process can be constructed step by step between consecutive jumps and the sequence of jump times (T m ) m∈N goes to infinity and the process in well-defined.The proof is thus complete.
Remark 3.5.Notice that if, instead of conditions (2:ii) − (2:iii) in Hypothesis 3.2 we require the weaker conditions Theorem 3.4 would follow analogously with the only difference that existence and uniqueness can be proved only up to a sufficiently small finite horizon time T 0 < ∞ rather than on the whole real line R + .
In particular, equation ( 28) does not hold.To see that, consider the truncated death rate and pairwise interaction rate a n (q; v ) := a(q; v )½ {v ≤n} ; b n (q; v ) := b(q; v )½ {v ≤n} : By the boundedness of the rates a n and b n , we have that Theorem 3.4 is valid and existence and uniqueness hold true up to a stopping time fi n .We further clearly have that fi n ≤ fi n+1 , so that, if fi n → ∞ as n → ∞, we have the existence and uniqueness for any time horizon T , whereas if on the contrary we have that fi n → T 0 , we have an explosion of the solution in finite time. △ The next result state a martingale property for the spatial GSM 2 introduced in previous sections.
Theorem 3.6.Assume that Hypothesis 3.2 holds true and that X 0 and Y 0 are two random measures independent with finite p−th moment, p ≥ 2. Then: (i) is a Markov process with infinitesimal generator L defined by (18); Then, the process is a cádlág martingale starting at 0; are cádlág L 2 −martingale starting at 0 with predictable quadratic variation given by Proof.(i) to show that is a Markov process is standard using the fact Poisson point processes have independent increments.Then, for any function f X and f Y ∈ C 2 0 (Q), we have the representation given in equation ( 17).By compensation, we can reformulate equation ( 17) as where MX and MY are local-martingales accounting for the noises W , N r , N a and N b .A straightforward computation shows that for F ∈ C 2 b (R × R), dividing equation ( 34) by t, taking the limit as t ↓ 0 and taking the expectation we finally have that LF (f X ;f Y ) ( ) has the expression as given in equation ( 18).
(ii) using condition (30) we have can infer that (31) is integrable and well-defined.Using point (i) we can finally conclude that (31) is a cádlág martingale.

On the initial distribution
As clear by the description given in Section 2.1, the initial distribution considered lacks any spatial distribution on the dose and on the formation of the lesion in the cell nucleus.To compute the initial lesion distribution 0 , we need to generalize the treatment given in Section 2.1 to include in equation ( 7) a spatial description.A possible mathematical formulation of the initial lesion computation would be the following.For a better understanding, we will provide a step-wise construction of such distribution: (i) given a certain dose D and fluence average specific energy z F , a random number of events in a cell nucleus is sampled from a distribution p e .A typical assumption would be, due to the independence of events, to assume p e a Poisson distribution with average D zF ; (ii) the events are distributed randomly over the cell nucleus.Under an isotropic and uniform random field, the distribution can be assumed to be uniform over the domain, or in a more general setting, the distribution can be sampled from an a priori calculated distribution of tracks using a condensed history MC code, [Agostinelli et al., 2003].A similar distribution has been for instance calculated in [Missiaggia et al., 2021, Missiaggia et al., 2022]; (iii) for any event i, i = 1; : : : ; , a certain specific energy z i is sampled according to the single-event microdosimetric specific energy distribution f 1 (z); (iv) for any event i, i = 1; : : : ; , with specific energy deposition z i , the number ‰ X i and ‰ Y i of sub-lethal and lethal lesion respectively is sampled from a distribution p.A typical assumption would be to assume such distribution the product of two independent Poisson distributions of average »(z i ) and -(z i ) respectively, for some suitable functions » and -; i the number of sub-lethal and lethal lesion respectively.Thus sample the positions (q X ; q X ) ∈ Q |(‰ X ;‰ Y )| , according to a distribution " ‰ ( q X ; q Y ˛‰X ; ‰ Y ).A reasonable choice for such distribution would be to distribute z i spatially around the track using the Amorphous Track model, [Kase et al., 2007], which is a parametrization of the radial dose distribution around a particle track.In particular, denoting by AT i (q) the normalized radial dose distribution representing the probability of depositing a certain dose in a domain, for any Q 1 ⊂ Q, the relative dose absorbed in Q 1 is thus given by Then, the probability density distribution describing the probability of creating a lesion in Q 1 is given by Remark 3.7.As mentioned in the introduction, a further choice would be to use track structure code to simulate the spatial distribution of lesions within a cell nucleus.Several papers have been published in the literature showing how this can be achieved, [Chatzipapas et al., 2022, Kyriakou et al., 2022, Zhu et al., 2020, Thibaut et al., 2023].Nonetheless, none of these papers then asses the biological effect of given radiation using a true spatial biological model, so the accuracy in the description of the geometry of the biological target is lost.

The protracted irradiation
In the present Section, we assume a further rate besides the interaction rates described in Section 3.Such a rate accounts for the formation of a new random number of both lethal and sub-lethal lesions due to protracted irradiation.We thus have the following system of possible pathways where ‰ X and ‰ Y are two suitable (possibly correlated) N−valued random variables.The last pathway, namely the dose-rate ḋ is reasonable to be assumed strictly positive up to a certain time horizon T ir r , representing the irradiation period.We thus have the following: (v) -protracted irradiation at a certain dose rate ḋ a random number ‰ X and ‰ Y sub-lethal and lethal lesions, respectively, are created in Q.This process can be described by a spatial compound random measure We assume that the random measure " admits a probability measure of the form with p(‰ X ; ‰ Y ) a discrete probability distribution on N 2 representing the probability of inducing (‰ X ; ‰ Y ) sub-lethal and lethal lesions and " ‰ a spatial distribution representing the probability of creating (‰ X ; ‰ Y ) sub-lethal and lethal lesions at positions (q X ; q Y ) ∈ Q |(‰ X ;‰ Y )| .We will further denote for short the marginal distributions by The protracted irradiation rate ḋ is associated to a Poisson point measure The corresponding intensity measure associated to N ḋ is We denote with Ñb the compensated Poisson measure defined as Ñḋ (ds; d‰ X ; d‰ Y ; dq; d" 1 ; d" 2 ) := N ḋ(d s; d‰ X ; d‰ Y ; dq; d" 1 ; d" 2 ) − -ḋ(d s; d‰ X ; d‰ Y ; dq; d" 1 ; d" 2 ) : Remark 4.1.The protracted irradiation can be interpreted as an improved description of the initial distribution.In particular, if ḋ is sufficiently large and T irr is sufficiently small, only the creation of new damages due to the protracted irradiation happens before any of the other pathways can happen.This is typically the case in the clinical scenario, where the dose rate usually dominates the biological interaction rates; such a situation is referred to as conventional dose rate.In this case, it is reasonable to assume that the initial distribution of lesions X 0 and Y 0 in the instantaneous irradiation, that is ḋ = 0, coincides with the distribution of lesions under protracted irradiation at time T 0 .For this reason, a typical distribution " can be obtained following the description provided in Section 2.1 in the particular case of a single event hitting the domain, that is = 1.It is further worth stressing that there are some relevant situations where an explicit treatment of the effect of protracted irradiation can play a relevant role: (i) a split dose irradiation treatment, where the treatment is split into several treatments with a smaller dose to favorite normal tissue recovery between treatments, (ii) space radioprotection, characterized by extremely low dose rates exposure over a long period and (iii) FLASH radiotherapy.Both (i) and (ii) are situations where is fundamental to model the entire spatial distribution of radiation-induced damages over a relatively long time period so that the inclusion of a specific protracted irradiation rate is necessary.For what instead concern (iii), it will be explicitly treated in Section 4.1.△ We will assume the following to hold.
Hypothesis 4.2.4. protracted dose rate components: (4.i) the protracted irradiation rate ḋ is positive and finite; (4.iii) the random measure p admits finite p−moments, that is, for p ≥ 1, it holds Therefore, the resulting process is characterized by the process defined in equation ( 17) with the addition of the random measure N ḋ.In particular, denote for short by L X B and L Y B the process introduced in equation ( 17), then the process under the effect of protracted irradiation is characterized by the following weak representation, given f X and f Y ∈ C 2 0 (R), we have We thus augment the probability space with the processes defined in (v ).We thus have the following definition.
We thus have the following well-posedness result.
Theorem 4.3.Let X 0 and Y 0 two random measures with finite p−th moment, p ≥ 1, that is it holds Then, under Hypothesis 3.2-4.2,for any T > 0 it exists a unique strong solution in D ([0; T ]; M × M) to the system (17).Also, it holds Proof.The proof proceeds with similar arguments as in the proof of Theorem 3.4, taking also into account the protracted irradiation term.
For t ≥ 0, we have, Notice that, for any positive integer x and y it holds (42) so that, using equation ( 42) into equation ( 41), we have that, for some C > 0 that can take possibly different values, Gronwall lemma implies that there exists C > 0 depending on p and T but independent of n, such that Similar arguments can be used to prove that The proof thus follows in a straightforward manner proceeding as in the proof of Theorem 3.4.
The above process is characterized by the following infinitesimal generator where ), h ∈ {r; a; b; ḋ}, are the infinitesimal generators introduced in equations ( 20)-( 21), whereas L ḋ is defined as We thus have the martingale properties and representation corresponding to Theorem 3.6.
Theorem 4.4.Assume that Hypothesis 3.2-4.2holds true and that X 0 and Y 0 are two random measures independent with finite p−th moment, p ≥ 2. Then: (i) is a Markov process with infinitesimal generator L defined by (46); Then, the process is a cádlág martingale starting at 0; are cádlág L 2 −martingale starting at 0 with predictable quadratic variation given by Proof.The proof is analogous to the proof of Theorem 3.6.

The bio-chemical system under protracted irradiation
As mentioned above, before focusing on a macroscopic limit of the spatial DNA repair model, we will consider a different setting relevant to the considered application.As commented in Section 2.1, the functions » andusually include information regarding the chemical environment and radical formation.In the following treatment, we make this assumption explicit, assuming that the formation of new damages depends on the chemical environment described by a set of reaction-diffusion equations.It is worth stressing that in general, the energy deposition of the particle also affects the chemical environment so the above-mentioned reaction-diffusion equation also includes a term dependent on the energy deposition.As the chemical evolves on a much faster time scale, the concentration of chemical species will be described by a set of parabolic PDE, with a random discontinuous inhomogeneous term due to the effect of radiation.We will not consider a specific model for the chemical environment, but on the contrary, we will assume a general mass control hypothesis that includes many possible systems proposed in the literature.Future investigation will be specifically devoted to analyzing and implementing the highly dimensional chemical system, including the homogeneous chemical stage also the heterogeneous one.
Assume a set of L chemical species, then, for i = 1; : : : ; L, the concentration of the i−th species  i evolves according to 8 > < > : We consider the following.
(v') -protracted irradiation at a certain dose rate ḋ a random number ‰ X and ‰ Y sub-lethal and lethal lesions, respectively, are created in Q.This process is described by a random measure that depends on chemical concentration We assume that the random measure " admits a decomposition of the form " ‰ ( q X ; q Y ˛‰X ; ‰ Y )p `‰X ; ‰ Y ˛(q; t) ´; with p(‰ X ; ‰ Y ) a discrete probability distribution on N 2 .
(vi) -chemical environment for all i = 1; : : : ; L, the random discontinuous inhomogeneous term G is defined as

:
Remark 4.5.The random function Z i represents the energy deposited by an event and can be computed as described in Section 2.1.Regarding instead the dependence of p on the chemical environment, several possible choices can be made.What is currently known is that the actual number of damage created by certain radiation depends on the chemical environment.Therefore, considering the description of the initial damage distribution given in Section 2.1, a meaningful choice would be to assume that given a certain specific energy deposition, the average number of induced lethal and sub-lethal lesions described in (iv ) depends on the chemical environment, that is we have -(z i ; ) and »(z i ;  i ).△ We will assume the following to hold.
Hypothesis 4.6.5. chemical environment components: (5.i) for all i = 1; : : : ; L, the random initial condition  0;i is bounded and non-negative P−a.s., that is  0;i (q; !) ∈ L ∞ (Q) ; and  0;i (q; !) ≥ 0 ; for a.e.q ∈ Q ; P − a.s. ; and has finite p−th moment, p ≥ 1, that is (5.ii) there exist constants C 0 and C 1 such that (5.iii) for all i = 1; : : : ; L, f i is locally Lipschitz and f i () ≥ 0 ; for all  = 0 ; (5.iv) for all i = 1; : : : ; L, there exist " > 0 and K > 0 such that (5.v) for all i = 1; : : : ; L, the random function Z i is bounded and non-negative P−a.s., that is Z i (q; !) ∈ L ∞ (Q) ; and Z i (q; !) ≥ 0 ; for a.e.q ∈ Q ; P − a.s. ; and has finite p−th moment, p ≥ 1, that is As mentioned in the Introduction, we will not focus on specific examples of chemical environments.Nonetheless, since almost any chemical model contains second-order reaction rates, we are forced to consider more general assumptions than the standard global Lipschitz condition.For this reason, we rather consider a mass control condition in Hypothesis 4.6(ii).Such assumptions can be immediately seen to hold in a reaction-diffusion description of the systems introduced in [Abolfath et al., 2020a, Labarbe et al., 2020].△ Therefore, the resulting process is characterized by the process defined in equation ( 38) with the addition of the chemical system as defined in equation ( 52).We thus have the following representation Augment then the filtration with the processes defined in (v ′ ) − (v i).
We thus have the following well-posedness result.
Theorem 4.8.Let X 0 and Y 0 two random measures with finite p−th moment, p ≥ 1, that is it holds Then, under Hypothesis 3.2-4.2-4.6, for any T > 0 it exists a unique strong solution of the system (17 . Also, it holds Proof.The main step of the proof follows alike Theorems 3.4-4.3,noticing that between consecutive jump times T m and T m+1 , under Hypothesis 4.6 using [Fellner et al., 2020, Theorem 1.1], equation ( 52) admits a unique strong solution in C . In fact, since Z is bounded and non-negative P−a.s.we can infer that for any jump time T m it holds  i (q; T m ) := lim t↑Tm  i (q; t) + Z(q) ;  i (q; T m ; !) ∈ L ∞ (Q) ; and  i (q; T m ; !) ≥ 0 ; for a.e.q ∈ Q ; P − a.s.: The rest of the proof follows Theorems 3.4-4.3.
Remark 4.9.It is worth noticing that, using [Fellner et al., 2020, Theorem 1.1], we can infer that, between consecutive jump times T m and T m+1 , it holds . Nonetheless, since we cannot require Z to be smooth, we cannot conclude that  ∈ D " (0; T ); `C2 0 (Q) ´L" as at the jump times T m , (T m ; q) can only be shown to be bounded.△

The large population limit
In the following we assume the model coefficients depend on a parameter K and study the behavior of the measure-valued solution K = ` X;K ; Y;K ´studied in previous Sections as K → ∞.
We thus consider a sequence of the initial measure so that 1 K " X;K 0 ; Y;K 0 " → `uX 0 ; u Y 0 ´and we assume that the rates of the model as introduced in Sections 3-4 are rescaled as follows b K (q 1 ; q 2 ; v ) := 1 K b " q 1 ; q 2 ; v K " ; ḋK := K ḋ : We thus aim at characterizing the limit as K → ∞ of the rescaled measure As above we also define for short the marginal distributions Analogous arguments to the one derived in Sections 3-4 show that ` X;K (t); Y;K (t) ´t≥0 is a Markov process with infinitesimal generator of the form and (1 − p(q 1 ; q 2 )) b(q 1 ; q 2 ; ) We thus can infer from Theorem 4.4 the following martingale property for the rescaled system.
Lemma 5.1.Consider K ≥ 1 fixed, assume that Hypothesis 3.2-4.2holds true and that X 0 and Y 0 are two random measures independent with finite p−th moment, p ≥ 2.Then, the processes M X and M Y defined for f are cádlág L 2 −martingales starting at 0 with predictable quadratic variation given by Proof.The proof is analogous to the proof of Theorem 3.6.
We assume the following.
Hypothesis 5.2.6. rescaled system: (6.i) the initial measure u K 0 = (u X;K 0 ; u Y;K 0 ) converges in law for the weak topology on M (6.ii) all the parameters r K , a K and b K are continuous on the corresponding space, that is either and they are Lipschitz continuous w.r.t. the last variable, that is there exist positive constants L r , L a and L b such that, for all q ; q 1 ; q 2 ∈ Q and v Next is the main result of the present section.
Then for all T > 0, the sequence `uK ´K∈N converges in law in D ([0; T ]; M F (Q)) to a deterministic continuous measure-valued function in C ([0; T ]; M F (Q)), which is the unique weak solution of the following non-linear integro-differential equation, The steps of the proof of Theorem 5.3 are standard in the literature, [Isaacson et al., 2022, Bansaye andMéléard, 2015], or also [Isaacson et al., 2022] for the case of a bimolecular reaction.Nonetheless, general results that include the case treated in the present paper are not available in the literature since [Isaacson et al., 2022] cannot include zero-th order reactions whereas other results do not include pairwise interaction.Therefore, the study of the large-population limit of a spatial model with pairwise interaction and random generation of lesions has never been considered in the literature, so existing results cannot be directly applied to the present case.As customary in literature, for the sake of readability the proof of Theorem 5.3 will be divided into four steps.

Step 1: uniqueness of solution
The first result concerns the proof of the uniqueness of the limiting process (63).Proof.Arguments similar to the proof in Theorem 4.3 imply that neglecting negative terms and using Grownall's lemma, the following estimate holds Denote by T X and T Y the semigroups generated by the operators L X d and L Y d .We thus have for any bounded and measurable functions f such that f We thus have ˛(b(q 1 ; q 2 ; u 2 ) − b(q 1 ; q 2 ; u 1 )) T X (t − s)(f X (q 1 ) + f X (q 2 )) ˛≤ C sup 1; u 1 (t) − u 2 (t) : (67) We therefore can infer, using estimates (67) in equation ( 66), that ˙f X ; u 1 (s) − u 2 (s) ¸ds : Using thus Gronwall's lemma, we can finally conclude that for all t ≤ T , sup f X ∞≤1 ˙f X ; u 1 (s) − u 2 (s) ¸= 0 ; from which the uniqueness follows.

5.2
Step 2: propagation of moments ˙1; u Y;K (t) ¸3 ≤ C : Taking the limit fi X n and fi Y n as K → ∞, Fatou's lemma implies equation (68).Estimate (69) thus follows from the fact that f X and f Y are bounded.

Step 3: tightness
In the following, we will denote by Λ K the law of the process u K = (u X;K ; u Y;K ).We then have the following.
Theorem 5.6.Assume Hypothesis 3.2-4.2-5.2 hold, then the sequence of laws `ΛK ´K∈N on D ([0; T ]; M F ) is tight when endowed with the vague topology.
Proof.We equip M F with the vague topology; to prove the tightness of the sequence of laws Λ K , using [Roelly-Coppoletta, 1986], it is enough to show that the sequence of laws of the processes ˙f X ; u X;K ¸and ˙f Y ; u Y;K ¸are tight in D ([0; T ]; R) for any function f X and f Y .As standard, in order to accomplish this, we employ the Aldous criterion [Aldous, 1978] and Rebolledo criterion [Joffe and Métivier, 1986].

5.4
Step 4: identification of the limit Theorem 5.7.Assume Hypothesis 3.2-4.2-5.2 hold, denote by Λ the limiting law of the sequence of laws `ΛK ´K∈N .Denote by u the process with law Λ.Then u is a.s.strongly continuous and it solves equation (63).
Proof.The fact that u is a.s.strongly continuous follows from the fact that ˛˙f X ; u X;K (t) ¸− ˙f X ; u X;K (t − ) ¸˛≤ 1 K : Consider t ≤ T , f X , f Y and u = (u X ; u Y ) ∈ D ([0; T ]; M F × M F ); define the functionals q; u) + a(q; u)] f X (q)u X (s)(dq)ds+ + Z t 0 Z Q2 b(q 1 ; q 2 ; u)(f X (q 1 ) + f X (q 2 ))u X (s)(dq 1 )u X (s)(dq 2 )ds+ A " ‰ X ( q X ˛‰X )dq X dp X (‰ X )ds ; p(q 1 ; q 2 )b(q 1 ; q 2 ; u)f Y (q)m b (q|q 1 ; q 2 )d qu X (s)(dq 1 )u X (s)(dq 2 )ds+ We aim to show that, for any t ≤ T , Using Lemma 5.1 we have that M X (t) = Ψ X t;f X (u which concludes the proof. Theorem 5.8.Assume Hypothesis 3.2-4.2-5.2 hold, then the sequence of laws `ΛK ´K∈N on D ([0; T ]; M F ) is tight when endowed with the weak topology.
Proof.Using [Méléard and Roelly, 1993], since the limiting process is continuous, repeating the above calculations with f X = f Y = 1, we obtain that the sequences `˙1; u X;K ¸´K and `˙1; u Y;K ¸´K converge to ˙1; u X ¸and to ˙1; u Y ¸in D ([0; T ]; R).

Conclusions
In the present paper, we introduced a general stochastic model that describes the formation and repair of radiation-induced DNA damage.The derived model generalizes the previously studied model, [Cordoni et al., 2021, Cordoni et al., 2022b, Cordoni et al., 2022a], including a spatial description, allowing for pairwise interaction of cluster damages that might depend on the distances between damages, as well as a general description of the effect of radiation on biological tissue under a broad range of irradiation condition.We studied the mathematical well-posedness of the system as well as we characterize the large system behavior.Further, we believe that the derived model could play a role in describing an effect of recent interest in radiobiology, called the FLASH effect.In particular, at extremely high rates of particle delivery, it has been empirically seen that the unwanted effects of radiation on healthy tissue decrease while the killing effect on the tumor is maintained.Although numerous studies on the topic, the physical and biological mechanism at the very core of this effect is today largely unknown.It is nonetheless believed that spatial interactions of particles can play a major role in the origin of this effect.Therefore, the model introduced in the present research can have the generality to provide a stochastic description of the effect that spatial interdependence between particles can have on biological tissue.

Lemma 5. 5 .
Assume that Hypothesis 3.2-4.2-5.2 holds, then for any T > 0, there exists a constant C := C(T ) > 0 which depends on T such that the following estimates holdsup K E sup t∈[0;T ] ˙1; u X;K (t) ¸3 ≤ C ; sup K E sup t∈[0;T ] ˙1; u Y;K (t) ¸3 ≤ C : ˙f X ; u X;K (t) ¸3 ≤ C ; sup K E sup t∈[0;T ] ˙f Y ; u Y;K (t) ¸3 ≤ C : (69)Proof.The proof follows from computations similar to what is done in the proof of Theorem 4.3 obtaining a similar estimate as in equation (27).That is we have,E sup t∈[0;T ∧fi X n ] ˙1; u X;K (t) ¸3 ≤ C ; E sup t∈[0;T ∧fi Y n ]