Computational study of the effect of hypoxia on cancer response to radiation treatment

We perform a computational study of the propagation of the oxygen concentration within a two-dimensional slice of a heterogeneous tumour region where the position and shape of the blood vessels are known. Exploiting the parameters space, we explore which effect is noticeable what concerns the formation of hypoxic zones. We use this information to anticipate a patient-specific radiation treatment with controlled response of the cancer growth.


Introduction
It has long been known that the lack of oxygen in cells (hypoxia) greatly can hinder the 2 effectiveness of a radiation treatment of cancer patients [2]. The concept of 3 hypoxic/normoxic fractions is insufficient for predicting radiation treatment 4 outcomes [3]. A good control of the tissue oxygenation is desired if one targets an 5 excellent radiation of localised cancer regions with a minimal damage on the 6 neighbouring tissue. This is the rationale for attempting to model, in detail, the tumour 7 oxygenation. It is worth noting that the literature on the modeling and simulation of 8 cancer radio-and chemo-therapies is large and keeps growing rather fast. At this stage, 9 we only mention [4,6,8] and references cited therein. 10 The amount of oxygen present locally in a tumour can be measured directly by 11 invasive methods (e.g. via sticking a needle). Pointwise indications can be obtained also 12 via imaging techniques, see e.g. [10] and references cited therein. In this context, we 13 support a non-invasive investigation of hypoxia by combining knowledge of the existing 14 vasculature and numerical estimations of the oxygen transport and metabolism. 15 From the modeling perspective, the key is understanding the transport of oxygen 16 through the blood vessels and through the surrounding heterogeneous environment. 17 Within this framework we focus the attention on the transport of oxygen outgoing from 18 the blood vessels travelling through both healthy and cancerous regions. The dynamics 19 appears to be simple: oxygen is supplied by blood vessels from which it is then diffused 20 in an environment where cells consume it. That is, as a first attempt, we can formulate 21 a macroscopic diffusion-reaction model for the diffusion and consumption of oxygen in 22 any given region when the arrangement of blood vessels and reasonable assumptions can 23 be stated on the structure of the effective diffusion coefficient of oxygen as well as on 24 the effective production rates by consumption.

25
Trusting previous works in this field (compare e.g. [9] and references cited therein), 26 our model equations are as listed in (1). We denote by Ω the spatial region which 27 October 17, 2020 1/9 includes the network of blood vessels and where the oxygen diffusion takes place during 28 a time interval (0, T f in ). Here T f in is the final observation time of the overall process.

29
Our modeling and computations are supposed to be valid until an a priori fixed time 30 instance T ≤ T f in (both T and T f in are supposed to be close to the moment when the 31 stationary state is reached). We denote by K = K(x, t) the mass concentration of 32 oxygen distributed in the region Ω \ V , where by V we denote the region occupied by the 33 blood vessels. Here x ∈ Ω \ V is the space variable while t ∈ (0, T ) is the time variable. 34 Our problem boils down to computing the oxygen concentration K(x, t) satisfying 35 the following model equations: In this context, n is the outer normal vector to ∂V , D = 2000µm 2 /s is an averaged 37 diffusion coefficient for oxygen in water, K 0 = 40mmHg 1 is the supply of oxygen by the 38 blood vessels 2 . Along the lines of [7], e.g., we assume for the production term in (1) the 39 Michaelis-Menten consumption structure: where C 0 = 15mmHg/s and K m = 1mmHg.

41
Interestingly, the model (1) is essentially a 3D-one; no correct reduction to a 2D 42 setting can be made simply because the network of blood vessels is too complex. It is 43 inherently a 3D object with no special structure (unless for particular cases). Hence, a 44 dimension reduction exercise from 3D to 2D cannot be done in general without 45 introducing too much uncontrolled errors in interpreting the geometry of the vasculature. 46 As far as we see, treating such a problem correctly is a challenging task that arises in tumours. This is a preparation step for handling a truly 3D setting.

57
(ii) Currently, we notice the lack of an efficient mesh generation process from slice 58 figures to mesh that could be imported into the computational platform FEniCS; 59 see [1] for the users manual. 60 1 This is only a typical value. We vary it in different ways as reported in the Results section. 2 "The vessel voxels represents the oxygen level on the surface of the vessel that is available for diffusion into the tissue. In our simulations, we use a default value of 40 mmHg (cf. [7]). However, it is important to note that the model allows for time and space variations in the choice of K 0 -for instance, a direct proportionality with the diameter of the vessel can be incorporated. 3 We point out a feature in the simulation which, with given probability, assigns an active vessel to have K 0 = 50mmHg instead of a lower K 0 which is the reference value set to all other active vessels. The probability for such an modification to occur is 1/3 as of now. We proceed this way mainly because the tumours use angiogenesis, which causes vessels to develop in a very unorganised manner, causing large differences in vessel length, trajectory and in tissue vascularisation, hence differences in the oxygen level would be expected.
October 17, 2020 2/9 The main question we ask is: (Q) Given the network as well as level of efficiency of 61 blood vessels in a tumour, can one detect the oxygen diffusion regions where hypoxia is 62 likely to be present? 63 We address this question via a simulation study of the model 1 when the spatial 64 distribution of blood vessels (occupying the region V ) is provided. Our results indicate 65 that the main challenge is to capture in detail the correct heterogeneous time and space 66 distribution of the oxygen concentration. Once available in real time, this information 67 can be used further to fine-tune a patient-specific radiation treatment.

75
By visual inspection of our available images on blood vessels, it seems that most of 76 the vessels are anyhow running perpendicular to the image surface and hence the effect 77 of vessels running 'along image surface' is supposed to be small. This allows us to 78 process good quality cross sections.

79
Generation of the computational mesh 80 We employ GMSH [5] to generate meshes that accurately capture the underlying vessel 81 structure from the 2D-slice images. We first obtain the vessels by thresholding the green 82 channel of the slice image so that all pixels with green channel having value greater

97
Firstly, we homogenise the Dirichlet boundary condition acting on the boundary of the region V occupied by the blood vessels, i.e. we set We rephrase the model equations as: where the production term in (3) is given by Based on (5), we observe that it makes sense to consider a parameter regime so that 99 K in − K 0 ≥ 0. In this context, we are going to choose a constant value for K 0 to 100 represent the contact with a "constant reservoir of oxygen" coming from the blood 101 vessels. However, the constant K 0 can vary between vessels as later indicated.

102
The corresponding weak formulation reads: Find u such that for all suitable test 103 functions ϕ the following integral identity holds: 104

105
We use the tailor-made mesh produced with the help of GMSH to discretise efficiently 106 our model equations. To handle (1) with corresponding initial and boundary conditions, 107 we use the open source computational platform FEniCS; see for details [1]. The 108 reference set of parameters is taken from [7,9]. Linearising and transforming to weak where ψ is test function, K n is the oxygen concentration evaluated at the previous time 111 step, while µ n = µ(K n ). To compute the correct approximate evolution of the overall 112 process, we need to ensure that (K − K n ) =: δK << K + K m as this was used to 113 linearise µ(K) from Eq. 2. We use ∆t = 20µs, which guarantees that the condition for 114 δK is satisfied. GMSH takes care in the automatic generation of the adaptive mesh, 115 which accounts for the precise distribution of the blood vessels.

116
The main output of our simulations consists of oxygen distribution profiles within 117 the tissue. As post-processing work, we evaluate the local oxygen density (we refer to it 118 in terms of histograms) as well as the expected irradiation response.

119
Oxygen distribution in tissue 120 Due to e.g. vasomotion, the set of active vessels, i.e., those delivering oxygen to tissue 121 as well as their oxygenation K 0 can vary. For a given slice in the tumour we model this 122 variation by randomly sampling a percentage (active vessel percentage takes values from 123 [15, 30, 50, 75, 100]) of vessels to be 'active', that is, they have a surface K 0 set to 124 non-zero value. Moreover, for any given active vessel there is a probability (1/3) that 125 the K 0 is set to high value (50 mmHg) to create variation in the vessel surface oxygen 126 concentrations.

127
We let the oxygen diffuse through the tissue while keeping the concentration at the 128 vessel surface constant. After the steady state is reached, i.e. (in practical terms) the

Response to irradiation 136
Based on the oxygen distribution at stationarity, say K ∞ (x), we can estimate the 137 required irradiation dose D r needed to ensure a given cell death percentage in the tissue. 138 We evaluate the dose (D99) required for cell death of 99%. The cell survival fraction (S) 139  around a point x depends on the local oxygenation K ∞ (x). After a dose D r 140 (compare [9] among others), we assume that it holds The α and β values are taken 0.03/Gy, and 0.003/Gy 2 respectively and OER m = 3, in 142 line with, [7,9,11]. To get the dose required for 99% cell death (D99) within tumour, we 143 solve numerically D r from the integral equation where ρ(x) is the cell number density (= constant), and N cells = Ω ρdA. The further components (see e.g. [4]). Furthermore, to bring the numerically estimated 160 information at the level of the patient (e.g. to suggest a personalised radiation 161 treatment (dosage, fractionation pattern, time schedule, etc.), we strongly believe that a 162 data-driven modeling needs to improve our predictions. We expect that perhaps MRI 163 scans of the patient's tumour can bring in such additional information.

164
There are several ways in which our results could be improved. We mention here 165 only three of the possible directions. However, these are as of now considered out of the 166 scope of this work. Some of these outlets will be considered in forthcoming publications. 167 • The surface concentrations could be sampled from a given distribution.

168
• We could use several slices of the tumour for averaging over the resulting oxygen 169 concentration histograms.

170
• The FEniCS code, approximating with FEM our model (1), can be extended in a 171 straightforward way to handle a 3D computational domain. However, we expect 172 the adaptation of the mesh (for instance via GMSH) to the 3D structure of the 173 blood vessels to be cumbersome.

174
• The visualisation of the hypoxic regions in 3D requires a special care. The hope is 175 to be able to detect visually the possibility of clogging in the oxygen transport 176 through Ω \ V .