A physicochemical model of X-ray induced photodynamic therapy (X-PDT) with an emphasis on tissue oxygen concentration and oxygenation

X-PDT is one of the novel cancer treatment approaches that uses high penetration X-ray radiation to activate photosensitizers (PSs) placed in deep seated tumors. After PS activation, some reactive oxygen species (ROS) like singlet oxygen (1O2) are produced that are very toxic for adjacent cells. Efficiency of X-PDT depends on 1O2 quantum yield as well as X-ray mortality rate. Despite many studies have been modeled X-PDT, little is known about the investigation of tissue oxygen content in treatment outcome. In the present study, we predicted X-PDT efficiency through a feedback of physiological parameters of tumor microenvironment includes tissue oxygen and oxygenation properties. The introduced physicochemical model of X-PDT estimates 1O2 production in a vascularized and non-vascularized tumor under different tissue oxygen levels to predict cell death probability in tumor and adjacent normal tissue. The results emphasized the importance of molecular oxygen and the presence of a vascular network in predicting X-PDT efficiency.

Cancer is the second cause of death in the world after cardiovascular disease 1 .Among the various methods proposed for cancer treatment, radiotherapy (RT) utilizes ionizing radiation to fatally damage tumor cells.High tissue penetration of X-rays induce DNA breaks, resulting efficient cancer cell destruction 2 .Simultaneously, ionizing radiation induce side effects for surrounding normal tissues 3 .So, it is a challenging problem to balance between inhibiting tumor growth and reducing side effects to normal tissue.To achieve this goal, the concept of therapeutic efficiency has been defined as the ratio of cell death in tumor to healthy tissue.There are different factors affect the efficiency of RT such as tissue morphological characteristics (includes tissue oxygenation, Cell division speed, …), nature of radiation (dose, linear energy transfer, …) and presence of substances that affect the cell response to radiation (like radiosesitizer and radioresistant agents) 4 .Among these, radiosensitizers are chemical agents increase lethal effects of radiation 5 and significantly lowering the radiation exposure necessary for tumor regression 6 .
In contrast to RT, photodynamic therapy (PDT) is a relatively new cancer treatment approach utilizes intermediate energies of electromagnetic spectrum (optical radiation) to produce cytotoxic ROS 7 .PDT interactions need three components of oxygen, PS and light with an appropriate wavelength to damage cells 8 .In the absence of each component, PDT interactions will not initiate.Thus, PDT is recognized as a safe, low toxicity and highly selective method of cancer treatment 7 .However, PDT is restricted to accessible tissues due to limited penetration depth of light.
X-PDT, on the other hands, is a novel approach can be used for treatment of cancer 9 .X-PDT utilizes nanocomposites include some nanoscintillators to locally convert X-rays to optical luminescence and activate conjugated PSs.In presence of oxygen, the activated PS, generate ROS to damage cancer cells 8,10 .Therefore, there are two important factors that affect the efficiency of X-PDT: first, scintillator light yield transforming X-ray into light, and then efficiency of the emitted photons being absorbed by the PS.
Moreover, metallic nanoscintillators increase local dose by selectively scattering and absorbing X-rays causing localized damage to DNA 5,11 .As a result, X-PDT combines RT and PDT to improve therapeutic outcome and reduces radiation damage to normal tissues 12 .X-PDT attacks both DNA and membrane of cancerous cells causing more effective damage than conventional RT or PDT alone 13 .With the growth of nanotechnology, the idea of using nanoscintillators for activating adjacent PSs introduced by Chen et al. in 2006 14 .Since then further efforts have been made by different researchers as reviewed in 9 .However, there are still some challenges for the clinical use of X-PDT 9 and mathematical modeling will be able to response some of them.
Morgan et al. in 2009 15 made the first attempts to estimate the amount of produced ROS by X-PDT.Although many uncertainties existed in their simulations, it was a first attempt and raised valid questions about the practicality of X-PDT as an effective treatment option 9 .
Bulin et al. in 2015 16 proposed an alternative model using GEANT4-based Monte Carlo simulations to numerically estimate the spatial energy distribution in a macroscopic volume of water loaded with nanoscintillators.The results showed at physiological NP concentrations, biological effects of radiation are primarily mediated by water, the most abundant molecule in biological tissues and absorption of energy by NPs is mainly driven by inelastic electron scattering and is nearly independent of the nanoparticle X-ray stopping power 9 .Following these results, Klein et al. in 2019 17 proposed a simplified electron cross section model to describe luminescence yield of NPs in a mix environment including water and nanoscintillators.They estimated that luminescence yield of NPs can be approximated as a function of the radiation dose, NP concentration, scintillator light yield and the electron cross-sections for tissue and the NP material.This model provides an upper bound for the actual number of scintillation photons emitted and is a better predictor of X-PDT efficacy.However, electron cross-section model as well as before theoretical models of X-PDT are not able to fully explain the X-PDT efficacy reported during pre-clinical tests 9 .The current models most predict the number of 1 O 2 per cell over a wide excitation energy range, which is below the required 10 7 1 O 2 PDT threshold for cell killing.This is while experimental studies demonstrated that X-PDT is more than just a PDT derivative but it is essentially a PDT and RT combination 12 .Moreover, the current models of X-PDT most focused on the interaction of ionizing radiation with nanoscintillators and generated light with PSs to estimate 1 O 2 concentration and predict produced 1 O 2 based on some physical parameters such as radiation dosage, nanoscintillator concentration, scintillator light yield and tumor volume [15][16][17][18] .Indeed, none of them investigate effects of tumor microenvironment on the treatment efficiency while it is known that some chemical parameters like molecular oxygen has a significant impact on the success or failure of cancer treatments 19,20 .
In the present study, we attempted to model cancerous cells treated with X-PDT through a feedback of tissue oxygen content.For this purpose, we first developed the mathematical model of tumor growth and then stablished a link between treatment parameters and physiological conditions of tumor environment.Due to the key role of molecular oxygen in RT effectiveness 21 and 1 O 2 production in PDT 15 , we do take care of oxygen and its dynamics as the main nutrient in tumor environment.We have developed a physicochemical model of X-PDT to incorporate molecular oxygen as one of the main factors determining the amount of produced ROS and also effectiveness of X-PDT 18 .In this regards, we considered a five compartment model (Fig. 1) including chemical transport, cancer cells dynamics, angiogenesis, blood flow and treatment models.

Results
Figure 2 shows initial results of the coupled model after the growth of 9 initial tumor cells includes (a) the tumor stages at the end of avascular growth phase after 41.6 days (b) inhibitor distribution produced randomly around the parent vessel as well as sprouts initiated from the center of each semi-circle (c) distribution of vessel radiuses and (d) oxygen production rates by the new microvessel structure after 141.6 days.According to the Fig. 2(a), tumor includes 3221 cells that 2304 quiescent cells placed at the internal layers and 917 proliferative cells at the external layer.Production of VEGF by the hypoxic tumor cells initiates angiogenesis process with sprouting 8 sprouts along the parent vessel and production of inhibitor around each new sprout (Fig. 2b).The activated As shown in Fig. 2(c), the newly formed vessels have smaller diameters than the older ones (across x axis) and also vessel diameters vary depending on their positions along the parent vessel (across y axis) depending on the total applied stimulations S total .Different branching patterns among 8 activated sprouts are also related to the various amount of S total applied to each sprout tip.Contribution of new vessels in providing molecular oxygen is shown in Fig. 2(d) and illustrates that lower radiuses of vessel segments as well as more segment connections provide more oxygen to the tumor.Figure 3 shows distribution of oxygen in tumor microenvironment after 141.6 days.In a vascularized tumor of Fig. 3(a), a mass of oxygen and other nutrients transferred by the new vessels in response to the VEGF produced by hypoxic tumor cells.On the other hand, in a non-vascularized tumor of Fig. 3(b), that angiogenesis is not initiated (B O2 = 0), the lack of oxygen inside the tumor formed a layer of quiescent cells.The maximum level of oxygen in this situation corresponds to the minimum level of oxygen in the vascularized tumor of Fig. 3(a).Now we investigate X-PDT efficiency for the two conditions of vascularized and non-vascularized tissue of Fig. 3. Figure 4 shows different cells of a vascularized tumor as well as normal tissue killed by 1 O 2 during X-PDT with dose 4 Gy and nanocomplex concentration of 4 mg/ml.We assume that half of drug concentration remains in the adjacent normal tissue.As illustrated, cells in regions with higher oxygen concentrations targeted more effectively by 1 O 2 both in tumor and normal tissue.According to the model considerations, new microvascular network initiated from the parent vessel on the left border of the domain, as a result cells in the left hand of the grid killed mostly by 1 O 2 .Based on these results, in this vascularized tumor, about 22.73% of tumor cells and 8.39% % of normal cells damaged due to 1 O 2 toxicity.Among these, 2.54% of proliferative tumor cells (Fig. 4a) and 20.19% of quiescent tumor cells (Fig. 4b) have been damaged.Moreover, there is an isolated area inside the tumor where no cells have been killed and the treatment efficiency has been zero.Indeed, oxygen distribution determines the pattern of cell killing and hypoxia in this region acts like a defensive wall that protects tumor cells against 1 O 2 toxicity.This isolated area is not necessarily placed in the center of the tumor but is placed in the region with minimum oxygen distribution.As well, the lack of oxygen in a non-vascularized tumor does not provide necessary 1 O 2 threshold for cell killing.As a result, no tumor cells are destroyed by this X-PDT component in a non-vascularized tumor in the same condition.
Killing probability of 1 O 2 in a vascularized and non-vascularized tumor under different tissue oxygen concentrations investigated in Fig. 5.As shown, killing probabilities depend quantitatively and qualitatively on the tissue oxygen level as well as tissue oxygenation properties.With increasing background oxygen concentration, probability of cell killing also increases both in vascularized and non-vascularized tumor, until the level of oxygen concentration is too high (C O2 = 10 CO 2ref ) that there is no difference between a vascularized and nonvascularized tumor in killing probability.
Tissue oxygenation properties also determine killing probabilities for tumor and adjacent normal tissue.As shown in Fig. 5, for a constant oxygen background concentration, killing probability differs depending on the tissue oxygenation conditions (Vascularized or non-vascularized).In a vascularized tissue, killing probability is distributed among different parts of the tumor (and normal tissue), while in a non-vascularized tissue different layers of tumor have almost similar killing probabilities.Moreover, in a non-vascularized tissue 1 O 2 production is not sufficient to effectively destroy cancer cells except when tissue oxygen level becomes too high (Fig. 5b,  d, f).It should be noted that in high oxygen levels (C O2 = 10 CO 2ref ), there is not high differences between the presence of a vascular network or not (Fig. 5e, f).
According to the results of Fig. 5, almost 50% of tumor cells will be killed by 1 O 2 toxicity in sufficient oxygen concentration (C O2 = 10 CO 2ref ) in a vascularized tumor, but this is not all X-PDT efficiency because RT component will also target the cells, simultaneously.Figure 6 shows tumor and normal cells killed by RT in the presence of 4 mg/ml nanoscintillator concentration and radiation dose 4 Gy in a normal oxygen level (CO 2ref ).According to the results, 37.26% of tumor cells killed by RT + RS that includes 3.59% proliferative cells and 33.67% quiescent cells.An interesting result is that the isolated area of tumor cells in Fig. 4(b) that 1 O 2 failed to destroy them, now  Figure 7 shows the efficiency of RT, RT + RS, PDT and X-PDT under different regimes of radiation doses and/ or drug concentrations in a vascularized tumor with a constant oxygen concentration of CO 2ref .According to the results, RT alone (0-4 Gy) did not induce significant cell death, while RT in presence of radiosensitizers (RT + RS) increase cell death in drug concentrations of more than 2 mg/ml.PDT results estimate the amount of cell death happened due to 1 O 2 production during X-PDT and X-PDT results summarize the total effects of RT + RS, 1 O 2 toxicity and synergy effects of them, simultaneously.As shown, X-PDT is more effective than RT + RS and PDT alone at all doses and concentrations.As a comparison, in radiation dose of 4 Gy and drug concentration of 4 mg/ml, the cell viability dropped to 84.11% in RT only, 70.27% in PDT, 62.74% in RT + RS and 33.5 in X-PDT.
According to the Fig. 6(c), almost all of normal cells (98%) targeted by RT + RS in presence of 2 mg/ml nano-scintillator and radiation dose 4 Gy.Although most of surrounding normal cells are damaged, they are able to reconstruct themselves by cell repair mechanisms.Moreover, we assumed a strict situation that half of drug concentration remains in normal tissue under irradiation.However, treatment conditions include drug concentrations remain in healthy tissues should be optimization to reduce side effects.In pharmaceutical experiments, many efforts have been made to design ideal drug delivery systems so that drug loading in tumor cells with higher pH becomes more than normal low-pH tissues 22,23 .Figure 8 shows a moderately sensitive normal tissue and tumor cells before, just after and 4.16 days after X-PDT.Because DNA repair pathways of cancerous cells has been deregulated 24 , they are unable to repair themselves (Fig. 8) as a result, in most cases, fractionated treatments lead to complete elimination of cancerous cells.In the current situation, although some of tumor cells (35.08%) still remains after the treatment (Fig. 8e), it is possible to omit remaining cells by applying a low dose RT alternately.
Figure 9 compares viability of tumor and normal cells under different X-PDT regimes.In comparison to the control group that receive RT only (C = 0 & D = 0 ), cell viability decreases as NP concentration increases.Minimum cell viability of tumor equals to about 35.08% that is related to the NP concentration of 4 mg/ml and radiation dose 4 Gy (D = 4 & C = 4).Although increasing drug concentration compared to increasing radiation dose makes more reduction in tumor cell viability, this trend is opposite in normal cells.In other words, normal cells are more sensitive to changing in radiation doses compared to changes in NP concentration.Applying dose 4 Gy in all drug concentrations decreases normal cell viability below 20%, means that it is very toxic for them while the maximum reduction for tumor cells achieved in this dose and NP concentration 4 mg/ml.Therefore, it is a decision making problem to compromise between killing tumor cells and unwanted normal tissue side effects.
For better illustration, Tables 1 and 2 present percentages of tumor and normal cells killed by 1 O 2 in a vascularized tissue under different oxygen concentrations and vascular network densities, respectively.Consulting the literature, we considered two reference values for tissue oxygen concentration (CO 2ref ) and vessel branching age (ψ ref ) and compared the related results with the results in 0.25, 0.5, 2 and 4 times variations.According to the results, cell killing rate increases as tissue oxygen levels and microvascular network densities increase.Background oxygen concentration models inherent tissue differences in oxygen levels and we intended to predict why treatment efficiency differs among different cancers like colon compared to lung.Moreover, long-term effects of microvascular network densities appears in tissue oxygen levels.According to the simulations, when the tissue oxygen concentration increases to 4 CO2 ref , 46.47% of tumor cells killed due to 1 O 2 toxicity that is 2 times larger than killing rate in CO2 ref condition.However, the rate of normal cell killing in this situation is more than 3 times compared to CO2 ref situation.As a result, a high tissue oxygen level although increases tumor destruction, it simultaneously affect healthy tissue destruction more effectively.Thus, optimizing treatment conditions,

Discussion
In this paper, we considered available oxygen concentration in the vicinity of X-PDT drugs to evaluate production rate of singlet oxygen in tumor microenvironment.We also estimated cell killing probability in a vascularized and non-vascularized tissue for different levels of oxygen concentrations.Combination effects of RT in targeting the isolated cells escaping from PDT, emphasized that X-PDT is more than just a PDT derivative but it is essentially a PDT and RT combination 12 .Introducing physicochemical model of X-PDT with considering tissue oxygen concentration estimated killing rates in different tumor layers (proliferative and quiescent) and illustrated why some of the cells escape from 1 O 2 toxicity.In this paper, we improved the existing models of X-PDT to investigate the effects of tissue oxygenation to find out the cause of different X-PDT outcomes in various tissues from the point of their oxygen content, although it did not greatly improve the number of produced 1 O 2 .Considering variable vessel diameters, cause to produce a more realistic vascular network in providing molecular oxygen.The results showed that death rates of the cells as well as mortality patterns depend on the tissue oxygen concentration and oxygenation properties such as micro vessel density.Hope that such studies have great importance to clinical applications of X-PDT in the near futures.

Mathematical modeling Overview
Development of cancer begins with a mutation in certain key genes of a normal cell.One of this mutation results is cell's runaway from normal homeostatic mechanisms causing irregular cell proliferation.After a while, over successive divisions, a cluster of about 10 6 cancer cells forms that compete for space and nutrients.These new invading cells change the concentration of underlying molecules result in hypoxia at internal layers of tumor.Hypoxic cells secrete tumor angiogenic factors (TAFs) causing endothelial cells (EC), lining in nearby vessel walls, to form new blood vessels from the existing ones (angiogenesis) 25 .After triggering angiogenesis, blood flow penetrates into the network and provides oxygen and nutrients to the hypoxic cells.Angiogenesis process is a key feature of malignant tumors enables them to invade surrounding tissues (metastasis) and leads to challenges in their treatment.In the following, first we simulate growth of some initial tumor cells along with the consumption of chemicals include oxygen, glucose and hydrogen production to obtain appropriate concentration of molecules at the time of treatment.To model vascular tumor growth, we coupled our previous model 26 to the nutrient diffusion model 27 and estimate available tissue oxygen.Then, the edited model of X-PDT applied to estimate viability of cancerous cells through a feedback of tumor micro-environment.Details of each model explained in the following.

Chemical transport
We considered the key diffusible chemicals of oxygen, glucose, hydrogen ions and VEGF.Each tumor cell is assumed to consume nutrients (glucose and oxygen) and excrete metabolic waste (hydrogen).Tumor cells also secrete VEGF under hypoxic conditions.Differential equation describes the concentration of such molecules is 28 : where i = O 2 , G, H, denotes molecules of oxygen, glucose and hydrogen ion respectively.D i are diffusion constants and functions A i and B i describe the consumption rate of molecules by the tumor cells and production rate of them by the vasculature, respectively.Also: where W i is the uptake rate of molecule i by cancer cells, ϕ denotes the density of tumor cells, q i is the vessel permeability to molecule i, R is microvessel radius, C b,i and C w,i are intravascular and wall concentrations of molecule i.Terms T P and E P,k are Boolean indicators represent a tumor cell locating at node P and a capillary segment connecting node P with its neighboring nodes, respectively.The summation is taken over point P and its adjacent nodes (k = N, S, W, E).
For the simplicity, it is assumed that the intravascular concentrations of oxygen, glucose and hydrogen ions are preserved throughout the vasculature, so C b,i has a constant value.Finally, the interstitial VEGF transport equation is defined as 27 : where D v is diffusion coefficient of VEGF, η P and θ P are Boolean indicators of the existence of hypoxic tumor cells and ECs on a grid point P, respectively, ε v is take up rate of VEGF by ECs and v is its natural decay rate.

Cell dynamics
This model consists of three different cell types includes tumor cells, normal tissue and vasculature.The internal metabolic activities of normal and tumor cells cover glycolysis, production of acid and oxidative phosphorylation 28 .Normal cells only proliferate if density of neighbors falls below a certain value ( ~ 80%) but tumor cells can proliferate until they fill all the available spaces 29 .We adopted the previously developed models of tumor metabolism to incorporate core metabolic activities of tumor and normal cells 27,29 .Using Michaelis-Menten kinetics, oxygen consumption W o is determined by: where C O2 (x,y) is oxygen concentration,V O2 is maximum oxygen consumption and K O2 is the half maximum oxygen concentration.Glucose consumption is also driven by: where A 0 = 29 5 V O2 is the baseline production rate of ATP, and C G and K G are defined as glucose concentration, and half maximum concentration, respectively.The coefficient P G is a multiplier representing the Warburg effect (i.e.altered glucose metabolism seen in many tumor cells) where P G = 1 corresponds to normal glucose consump- tion, and P G > 1 denotes more glucose consumption than needed to meet normal ATP demand.Defining the target ATP production rate as W A A 0 , where: And where parameter K H accounts for proton buffering in the tumor microenvironment.
In the following, a set of cellular automaton (CA) rules are implemented into each cell to update its status in each time step (Fig. 10).One tumor cell per grid point is permitted either be a proliferating, quiescent, dead or treated ones.We developed some predefined rules 27 to incorporate progression of cells in cell cycle and model the effects of radiation cell cycle sensitivity 21 .Cell cycle is consisted of four stages: Mitosis (M), Gap1 (G1), Synthesis (S), and Gap2 (G2).In each stage, the cell undergoes a different process.During Mitosis, the cell divides (1)

Angiogenesis
The first step of angiogenesis is nearby vessel sprouting in response to factors secreted from hypoxic tumor cells (like VEGF).These sprouts then progress in extracellular matrix (ECM).Thus, two steps can be considered in angiogenesis modeling include sprouting and progression.To model sprouting, we use our previous model of adaptive sprout spacing (ASS) to estimate the number, position and time of activated sprouts along the parent vessel 26,32 .According this model, there are two conditions for sprouts to form: 1) If the concentration of activator around each point of the parent blood vessel become greater than or equal to a chosen trigger value, and 2) the inhibitor concentration at that point be less than or equal to a chosen inhibitor threshold value.Thus, If A and I denote the angiogenesis activator and inhibitor concentrations, respectively, the required conditions for sprouting will be formulated as: where the subscripts specify the location on the grid and the superscripts specify the time steps.
Assuming that the activator is produced by hypoxic tumor cells, simply diffuses and decays.Inhibitor is also produced by each formed sprout and diffuses into the tissue around of it to inhibit formation of new sprouts with a random distance and then decays.Thus, the partial differential equations describing activator and inhibitor concentrations are as follows: where D a and D i are activator and inhibitor diffusion coefficients, a and i are their decay rates, respectively, all taken as constants 33,34 .
Activated sprouts then progress in ECM in three ways of random motility, chemotaxis and haptotaxis 35 .Three key variables considered for progression of activated sprouts include vascular endothelial growth factor (VEGF), fibronectin and ECs.The non-dimensional equation governing the distribution of ECs is defined as 35 : where e, C v and f denote non-dimensional endothelial cell density, VEGF and fibronectin concentrations, respec- tively.D e , x and ρ are non-dimensional EC diffusion, chemotaxis, and haptotaxis coefficients, respectively.Using standard finite-difference methods, discretized form of the partial differential equations result in coefficients of the five-point finite-difference stencil to generate the probabilities of movement of an individual cell in response to its local milieu.The equation describing the interaction of sprouts with ECM fibers is: where β and γ are positive constants expressing fibronectin production and proteolytic activity of EC tip, respec- tively.As the activated sprouts progress in ECM, they branch and anastomosis to form a novel capillary structure.Then blood flows in the hollow capillaries to supply tumor for chemicals.Moreover, the interaction of blood with capillary network remodel their structure via wall shear stress and upregulate microvessel radius.

Blood perfusion
Vascular network formed by angiogenesis model in previous section is a hollow network just determines the path of the endothelial cells.In order to estimate the amount of chemicals transferred by new vessels, blood flow analysis is required.
For non-Newtonian fluids, apparent viscosity μ app is defined as the slope of the rheological curve at a specific shear rate.The well-known Poiseuille relation is applied to determine blood flow Q, while blood viscosity μ app (R, H D ) is defined as a function of capillary radius and hematocrit.Assuming blood flow as an incompressible flow, then mass conservation equation is applied for all nodes of the network as: When blood flows through a flexible vessel, resistance stresses will be generated depend on the vessel radius.On the other hand, when new vessels are forming, their radius can be determined depends on the entered shear stresses as: where Δt is the time step.We developed the chemical transport Eq. ( 3) to incorporate the effect of vessel diameter on the microvessel structure and amount of oxygen transferred by the novel network.Then, according to the amount of available oxygen and in the presence of specified concentrations of drug and radiation dose, cell viability will be computed.Therefore, Eq. ( 14) is used to update the new vessel radius as R new = R old + R .The total stimulus on a considered segment corresponds to the sum of each individual stimulus, S wss , S p , and S m representing the wall shear stress, intravascular pressure, and metabolic stimulus, respectively as follows 36 : We used our previous algorithm of vessel adaptation 26 that employs an auxiliary function to predict network remodeling at specified future time steps.The updated network with new capillary radiuses, branches and flow parameters includes updated shear stresses, pressure differences and flow rates substituted with previous ones at each time step.Finally, production rates of molecules at new vascular network updated according to the new parameters.

Treatment model
Previous studies have experimentally demonstrated that X-PDT is essentially a PDT and RT combination 12 .The two modalities target different cellular components (cell membrane by PDT and DNA by RT), leading to enhanced therapy effects.The efficacy of X-PDT is therefore due to synergy effects of RT and PDT 37 .Current models have focused only on the singlet oxygen concentration to predict treatment efficiency of X-PDT and none of them considered RT component.In the following, we model the contribution of X-PDT components and investigate the role of tissue oxygen concentration and oxygenation properties in cell viabilities.

PDT component
Absorption of a photon by PS in its ground state, promotes it to a singlet excited state S1 (Fig. 11).PS then lose its energy by emitting fluorescence and return to the ground state or transfer to a triplet state T1 that has a longer lifetime.This process known as intersystem crossing (ISC) is an essential feature of a good PS.In the triplet state, PS may return to its ground state by emitting phosphorescence or transfer its energy to molecular oxygen (Type II) or transfer an electron to the biomolecules (Type I), which in subsequent reactions produces ROS 38 .Table 3. Parameter definitions and units of X-PDT process described in Fig. 11.  3.

Variable
For modelling PDT component, we focus on the production of 1 O 2 as the dominant molecule during type II of PDT reactions.Because the dynamic processes of PDT are known to be very fast ( ~ μs or less) 38 , the simplified equation of 1 O 2 production in the time scale of a few seconds to hours, incorporates the mechanism developed by Georgakoudi et al. 39 and the corresponding X-PDT energy transfer mechanism proposed by Klein et al. 17

as:
Where N scint is density of scintillation photons emitted by a dilute suspension of nanoscintillator and C O2 (x,y) is oxygen concentrationin at the point (x,y).K p is decay rate of triplet PS to the ground state and K ot is the bimolecular rate of triplet PS quenching by 3 O 2 .We assume that all photons are converted to singlet oxygen molecules.N scint describes energy transfer in X-PDT and depends on some physical parameters like radiation dosage d, nanoscintillator concentration C sc , scintillator light yield Y sc as 17 : where μ/ρ is the electron cross-sections for tissue and the nanoscintillator material can be obtained from the ESTAR database, which is maintained by the National Institute of Standards and Technologies (NIST).
Finally, the corresponding number of 1 O 2 production per cells ( N 1O2 ) computed by multiplying the cell volume as: Available molecular oxygen plays an important role in treatment efficiency so that quantum yield of 1 O 2 is affected by the amount of oxygen consumption during X-PDT 19 .The physicochemical model of X-PDT can predict treatment failure at hypoxic tumor layers.However, for any (x,y) of the network if C O2 (x,y) > > K p /K ot , treatment efficiency depends only on the physical parameters of Eq. ( 17).As well, if C O2 (x,y) < < K p /K ot , production rate of 1 O 2 is concerned by oxygen content of tissue, leading to unpredictable treatment effect 40,41 .Therefore, the correspondence between the oxygen content of the tissue and the selected photosensitizer can be effective in treatment efficiency.In this regards, the oxygen content of the tissue for a specific photosensitizer in different states of presence or absence of a vascular network, as well as different tissue oxygen levels, has been investigated in the results.

Estimation of compensation coefficient
Proposed physicochemical model of X-PDT by considering oxygen concentration may not greatly improve the number of produced 1 O 2 , but it could illustrate dependence of 1 O 2 yield to the available molecular oxygen.Therefore, in order to model cellular death caused by 1 O 2 , we used invers engineering to compute the required coefficient for compensating the lack of 1 O 2 number from the lethal threshold value 4 × 10 7 reported by Neider et al. 42 .We assume this compensation coefficient η c represents all the factors that is predicted in experimental observations but have not yet been considered in modelling.Some of these factors include non-optical forms of energy transfer between nanoscintillator and conjugated photosensitizers 9 , catalyzing the radiation-induced formation of radicals and ROS by NPs 43 , enhancing radiation dose within nanometers of NP surfaces 44 .The role of all these factors applied by product a compensation coefficient η c in Eq. (18).

RT component
As mentioned before, each nanocomposite used in X-PDT contains a nanoscintillator and some conjugated PSs.While nanoscintillator is a light source for PS, it has also some radiation sensitizing properties 11 .In this part, we consider the effects of radiation at existence of nanoscintillators.Using nanoscintillators (radiosensitizers) can have many different mechanisms of action on effective killing of cancerous cells and improvement of treatment efficiency.Cardilin et al. 6 proposed a model based on LQ theory and lump all of radiosensitizer processes together as having the net action of linear stimulating the radiation-induced mass transfer to make the model generally applicable.According to this model, the proportion of cancer cells surviving an irradiation dose d and plasma concentration of the radiosensitizer C equals to: where α (Gy −1 ) and β (Gy −2 ) are cell-specific radiosensitivity parameters and b is a pharmacodynamic parameter associated with the radiosensitizing effects.
We call this part of X-PDT as "RT + RS component" and couple this linear stimulatory function with the tumor growth model as shown in Fig. 10 in order to estimate this component of X-PDT at the time of treatment.We also consider cell cycle sensitivity under irradiation 21  www.nature.com/scientificreports/α/β = 3.1 45 to estimate treatment side effects.Moreover, we assume that normal cells have the ability to repair damages caused by radiation while the repair pathways in tumor cells is demolished due to mutation 46 .

Initial and boundary conditions
Simulations were carried out on a 200 × 200 grid, which is a discretization of the unit square [0, 1] × [0, 1], with a space step of h = 0.005, representing a tissue with dimensions of 5 mm × 5 mm.Consulting the literature a grid size of Δx = Δy = 25 μm was set which corresponds to the approximate size of a tumor cell 47 .We initiallly assume 9 proliferative tumor cells with random stages and appropriate ages placed at the center of the domain.The rest of the lattice space is filled by normal cells with a density of 80%.A parent vessel is lining on the left border of the domain.Initial pH value of the tumor microenvironment is set to 7.4.Initial concentrations of 5.2 × 10 −6 mol/L and 5 × 10 −3 mol/L are set for oxygen and glucose, respectively 27 .Also, the nodal pressure difference of blood in preexisting vessels is set to ΔP = 100 mmHg with the radius of 14 µm.All capillary segments are assumed to be of initial radius and length of 2 µm and 25 µm, respectively.No flux boundary conditions are adopted for the interstitial diffusion equations.Other required parameter values are given in Table 4.

Figure 1 .
Figure 1.Overview of the model includes (a) Block diagram of the five compartment model and (b) Schematic of simulation space.

Figure 2 .
Figure 2. A snapshot of initial results of the coupled model includes: (a) tumor cell stages at the end of avascular growth phase (b) inhibitor distribution produced randomly around the parent vessel (c) vessel radius (d) Contribution of the produced capillary network in oxygen supply.

Figure 3 .
Figure 3. Distribution of oxygen in tumor and adjacent normal tissue, in a (a) vascularized tumor and (b) nonvascularized tumor.

Figure 4 .
Figure 4. Proliferative (a) and quiescent (b) tumor cells as well as normal tissue (c) killed by 1 O 2 during X-PDT in a vascularized tumor.The white spots in (a) & (b) are tumor cells and in (c) are normal cells.The black spots are background area in all.

Figure 5 .
Figure 5. Probability of cell killing by 1 O 2 during X-PDT in a vascularized and non-vascularized tissue as well as different tissue oxygen levels under dose 4 Gy and nanocomplex concentration of 4 mg/ml.

Figure 6 .
Figure 6.Different cells killed by RT + RS component of X-PDT includes (a) Proliferative and (b) Quiescent tumor cells as well as (c) Normal tissue.The white spots in (a) & (b) are tumor cells and in (c) are normal cells.The black spots are background area in all.

Figure 7 .
Figure 7.Comparison of mean tumor cell viability under different treatments of RT, RT + RS, PDT and X-PDT in a vascularized tumor with a constant tissue oxygen concentration of CO 2ref .

Figure 8 .
Figure 8.The results of X-PDT on a moderately sensitive normal tissue (a-c) as well as tumor cells (d-f) before (a & d), just after (b & e) and 4.16 days after X-PDT (c & f).The yellow color in the first row indicates normal cells and in the second row indicates tumor cells.The blue color in both rows shows the background area.

Figure 9 .
Figure 9.Comparison of tumor and normal cell viability under different X-PDT regimes with this assumption that half of these concentrations remains in adjacent normal tissue.

Figure 10 .
Figure 10.Model overview includes coupling of chemical diffusion model, cellular automaton rules and cancer cell dynamics.It is also shown how to get feedback from the model to apply the treatment at the desired time.

1 O 2
formation in type II of PDT reactions are focused for modeling PDT component of X-PDT.

1 ϕ
K otThe ratio of the monomolecular decay rate of the triplet state photosensitizer to the biomolecular rate of the triplet photosensitizer quenching by 3 O 2 11.9 × 10 −6 (PpIx) mol L −Density

Table 1 .
1ercentages of different tumor cells and normal tissue killed by1O 2 during X-PDT under various tissue oxygen concentrations in dose 4 Gy and NP concentration of 4 mg/ml.

Table 2 .
Percentages of different tumor cells and normal tissue killed by 1 O 2 during X-PDT under various vascular network densities in dose 4 Gy and NP concentration of 4 mg/ml.