{'fraction_non_alphanumeric': 0.03854293854293854, 'fraction_numerical': 0.014814814814814815, 'mean_word_length': 4.412114537444934, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Deterministic Partial Differential Equation Model for Dose Calculation in Electron Radiotherapy

Treatment with high energy ionizing radiation is one of the main methods in modern cancer therapy that is in clinical use. During the last decades, two main approaches to dose calculation were used, Monte Carlo simulations and semi-empirical models based on Fermi-Eyges theory. A third way to dose calculation has only recently attracted attention in the medical physics community. This approach is based on the deterministic kinetic equations of radiative transfer. Starting from these, we derive a macroscopic partial differential equation model for electron transport in tissue. This model involves an angular closure in the phase space. It is exact for the free-streaming and the isotropic regime. We solve it numerically by a newly developed HLLC scheme based on [BerCharDub], that exactly preserves key properties of the analytical solution on the discrete level. Several numerical results for test cases from the medical physics literature are presented.


Introduction
Mathematical methods play an increasing role in medicine, especially in radiation therapy. Several special journal issues have been devoted to cancer modeling and treatment, cf. [4,5,6,12] among others.
Together with surgery and chemotherapy, the use of ionizing radiation is one of the main tools in the therapy of cancer. The aim of radiation treatment is to deposit enough energy in cancer cells so that they are destroyed. On the other hand, healthy tissue around the cancer cells should be harmed as little as possible. Furthermore, some regions at risk, like the spinal chord, should receive a dose below a certain threshold.
Most dose calculation algorithms in clinical use rely on the Fermi-Eyges theory of radiation which is insufficient at inhomogenities, e.g. the lung. This work, on the other hand, starts with a Boltzmann transport model for the radiation which accurately describes all physical interactions.
Until recently, dose calculation using a Boltzmann transport equation has not attracted much attention in the medical physics community. This access is based on deterministic transport equations of radiative transfer. Similar to Monte Carlo simulations it relies on a rigorous model of the physical interactions in human tissue that can in principle be solved exactly. Monte Carlo simulations are widely used, but it has been argued that a grid-based Boltzmann solution should have the same computational complexity [9]. Electron and combined photon and electron radiation were studied in the context of inverse therapy planning, cf. [35,34] and most recently [36]. A consistent model of combined photon and electron radiation was developed [22] that includes the most important physical interactions. Furthermore, several neutral particle codes have been applied to the dose calculation problem, see [19] for a review and most recently [38].
In this paper, we want to study a macroscopic approximation to the mesoscopic transport equation. After the problem formulation in section 2, we derive the approximation of the macroscopic model in section 3. This approximation consists of a system of nonlinear hyperbolic partial differential equations, whose properties we briefly discuss. Due to the possibility of shock solutions, hyperbolic PDEs have to be solved with great care. In section 5, we introduce a scheme which is adapted to the problem at hand. Numerical results for tests from the medical physics literature are presented in section 6.

A deterministic model for dose calculation
A ray of high energy electrons that interacts with human tissue is subject to elastic scattering processes and inelastic ones. It is this latter process that leads to energy deposition in the tissue i.e. to absorbed dose.
To formulate a transport equation for electrons we study their fluence in phase space. Let ψ(r, ε, Ω) cos ΘdAdΩdεdt be the number of electrons at position r -r being a vector in 2D or 3D space-that move in time dt through area dA into the element of solid angle dΩ around Ω with an energy in the interval (ε, ε + dε). The angle between direction Ω and the outer normal of dA is denoted by Θ. The kinetic energy ε of the electrons is measured in units of m e c 2 , where m e is the electron mass and c is the speed of light.

Boltzmann transport equation
The transport equation can generally be formulated as [14] Ω · ∇ψ(r, ε, Ω) = ρ in (r) with σ in being the differential scattering cross section for inelastic scattering, and σ el the differential cross section for elastic scattering; σ tot in = S 2 σ in dΩ and σ tot el = S 2 σ el dΩ are the total cross sections for inelastic and elastic scattering, respectively; ρ in and ρ el are the densities of the respective scattering centers.
Explicit formulas for the cross sections that we used in this model can be found in section 2.3. They are based on the model developed in [22]. The energy integration is performed over (ε, ∞) since the electrons lose energy in every scattering event. Also, we consider only electron radiation. Equation (1) could also be used to model electrons which are generated by the interactions of photons with matter, as in [22]. In this case we would have an additional source term on the right hand side for the generated electrons.
Besides the transport equation one needs an equation for the absorbed dose. It was derived in [22] as an asymptotic limit of a model with a finite lower energy bound ε s > 0. The formula is exact if one chooses the lower energy limit ε s = 0, as we do here.
T being the duration of the irradiation of the patient and ρ the mass density of the irradiated tissue. If all quantities are calculated in SI units, equation (2) leads to SI units J/kg or Gray (Gy) for the dose. S is the stopping power related to the inelastic cross section. It is defined as

Continuous slowing-down approximation
Electron transport in tissue has very distinctive properties. The soft collision differential scattering cross sections have a pronounced maximum for small scattering angles and small energy loss.
This allows for a simplification of the scattering terms in the Boltzmann equation. The Fokker-Planck equation is the result of an asymptotic analysis for both small energy loss and small deflections. It has been rigorously derived in [31] and has been applied to the above Boltzmann model in [22]. However, some electrons will also experience hard collisions with large changes in direction and energy losses which have to be described by Boltzmann integral terms. Thus we only use an asymptotic analysis to describe energy loss, called continuous slowing-down approximation. This approximation has a greater domain of validity than the Fokker-Planck approximation. The Boltzmann equation in continuous slowing-down approximation (BCSD) is [25] Ω · ∇ψ(r, ε, Ω) = ρ in (r) A truncation in the energy space is introduced, that does not allow particles with arbitrary high energy, lim ε→∞ ψ(r, ε, Ω) = 0.
In the numerical simulations, we use a sufficiently large cutoff energy. Furthermore, we prescribe the ingoing radiation at the spatial boundary, ψ(r, ε, Ω) = ψ b (r, ε, Ω) for n · Ω < 0, where n is the unit outward normal vector.

Henyey-Greenstein Scattering Theory
The detailed interactions of electrons with atoms give rise to complicated explicit formulas for the scattering coefficients. Because of this, many studies use the simplified Henyey-Greenstein scattering kernel for elastic scattering [2], The parameter g, which can depend on r, is the average cosine of the scattering angle and is a measure for the anisotropy of the scattering. The case where g ≤ 1 matches an anisotropic scattering configuration.

Mott and Møller Scattering
A more realistic model for elastic and inelastic scattering of electrons in tissue has been developed in [22]. This model introduces material parameters (namely densities ρ e and ρ c , ionization energy ε B and effective atomic charge Z). The energy integration for inelastic scattering is cut-off at ε B . The model uses the Mott scattering formula for elastic scattering of an electron by an ion [29,26] with ϑ = arccos(Ω ′ · Ω), and ε is the outcoming electron energy in m e c 2 units. Here, α ≈ 1/137 is the fine structure constant, Z is the atomic number of the irradiated medium, r e is the classical electron radius. Z depends on r to account for heterogeneous media. To avoid an otherwise occurring singularity at ϑ = 0 a screening parameter can be introduced [39] that models the screening effect of the electrons of the atomic shell, denoted by a.
The inelastic scattering process is Møller scattering, where an electron impinges an atom that releases itself an additional electron e − + a → a + + 2 e − For this process, the electrons can be considered indistinguishable. The electron which has the higher energy after the collision is called primary electron, the other electron secondary. Due to kinematical reasons of the scattering processes the range of solid angles in Møller scattering is restricted. After the collision, the angle between the directions of the electrons is at most π/2. For an angle in [0, π/4], the electron with energy ε is the primary electron, for an angle in [π/4, π/2], it is the secondary electron. Therefore the Møller cross section can be written as where χ denotes the characteristic function of a set, is the Møller differential cross section of primary electrons and is the Møller differential cross section of secondary electrons. Here, In the simulations the model parameters ρ el , ρ in , ε B and Z are fitted to tabulated values taken from the database of the PENELOPE Monte Carlo code [32].

Partial Differential Equation Model
We will try to reduce the cost of solving system (1) by assuming a minimum entropy principle for the angle distribution of particles. This principle has been first proposed by Jaynes [24] as a method to select the most likely state of a thermodynamical system having only incomplete information. It has subsequently been developed in [28], [27], [1] and [15], among others, and has become the main concept of rational extended thermodynamics [30]. A full account and an exhaustive list of references on the historical development can be found in [21]. We define the first three moments in angle: where we note that ψ (0) is a scalar, ψ (1) is a vector and ψ (2) is a tensor. If we integrate the system (3) over Ω, we can derive the following equations, We have introduced the transport coefficients T el (r, ε) = πρ el (r)   These coefficients and the stopping power can be computed for both Henyey-Greenstein and Mott/Møller scattering. Explicit expressions can be found in [22,17].
The remaining problem is the computation of moment ψ (2) as a function of ψ (0) and ψ (1) . The Minimum Entropy M 1 closure for electrons [10] can be derived in the following way. To close the system we determine a distribution function ψ ME that minimizes the entropy of the electrons, under the constraint that it reproduces the lower order moments, By using this entropy, we have implicitly assumed that the electrons obey classical Maxwell-Boltzmann statistics. This is justified, since here quantum effects can be neglected.
Analogous to the calculations in [27] we can show that the entropy minimizer has the following form, ψ ME = a 0 exp(−Ω · a 1 ), where a 0 is a non-negative scalar, and a 1 is a three component real valued vector. This is a Maxwell-Boltzmann type distribution and a 0 , a 1 are (scaled) Lagrange multipliers enforcing the constraints. An important parameter is the anisotropy parameter α, whose norm is by construction less than or equal to one. If we compute the different moments of the distribution function given by (15) we obtain, In fact, these relations can be combined to give, or by taking the modulus, The relation (18) cannot be inverted explicitly by hand, i.e. we cannot express |a 1 | as a function of α in a closed form. However, this relation determines a unique solution which can in principle be computed. If we assume that we know a 1 , ψ (2) can be computed as where is a function of α by means of (18). For its efficient numerical evaluation, the Eddington factor has to be approximated. Several possibilities exist: • One could solve the closure relation (18) for |a 1 | e.g. by a Newton iteration in each step during the simulation.
• One could precompute a table that gives the Eddington factor χ as a function of α.
• One could approximate χ(α) by a suitable special function.
The second approach has been followed in [17]. It is advantageous only if the space in which one interpolates is low-dimensional. For more moments, this approach becomes more expensive, and the first approach appears to be more advantageous.
In some cases, an ansatz for χ can provide a good approximation. This is the approach we are following here. The Eddington factor χ can be approximated by a very simple rational function, This approximation is very accurate (the difference with exact curve is about 10 −15 ). The coefficients are given by

Properties of the System
In the literature, the system that has been thoroughly investigated (both analytically and numerically) is system (10) restricted to its conservative terms, without external sources, but with time-dependence. In the present work, we adapt a pseudo-time technique. We focus on the spatial discretization and use a standard discretization for the terms on the right-hand side. Thus we consider with the closure (18). The Eddington factor χ is shown in Figure 1. Furthermore, we show the system eigenvalues in two dimensions. In the isotropic regime (anisotropy parameter zero), they coincide with the P1 eigenvalues. On the other hand, in the case of free-streaming (|α| = 1), they coincide and have absolute value one. Thus the system (10) is hyperbolic and the speed of propagation is limited by one. Moreover the system is hyperbolic symmetrisable [15].
System (23) closed by the relation (19) has been analyzed thoroughly in [23]. There, solutions to Riemann problems are constructed and invariant regions are computed. Since the reconstruction (15) of the kinetic distribution ψ is always positive, it can be expected that system (19) -(23) must admit a positive solution ψ (0) and a limited flux α < 1. To our knowledge, however, there exists no proof of this fact. The invariant regions computed in [23] only cover a subset of all admissible values. For a related model [18], bounds were proved, but only in 1D and steady state. Nevertheless, we construct a scheme which preserves exactly the positivity of ψ (0) and the flux limitation, i.e. the convex set of the admissible states of the system (23) is [8] A = ψ (0) , ψ (1) : ψ (0) ≥ 0, |ψ (1) | ≤ ψ (0) .
In the absence of sources or boundaries, the total mass, momentum and energy are conserved. In addition, the minimum entropy system recovers the equilibrium diffusion regime as a relaxation limit for large absorption coefficients [13].
In a two-or three-dimensional geometry, we have in addition [8]: Let n be the unit normal vector to an interface; then the system exhibits two acoustic waves, with velocities λ L (n) and λ R (n), supplemented by a contact wave with velocity β (n). The quantity β · n satisfies the following inequality λ L (n) ≤ β (n) · n ≤ λ R (n). The Riemann invariants associated with the contact wave are {β , Π}. They are defined by the relations (24b)

Numerical Method
The properties of the continuous model should be reproduced by the numerical scheme. In particular the positivity and flux limitation constraints are fundamental. An HLL scheme [20] can be constructed [3,11,8], that satisfies the required properties. However such an approach cannot capture the contact discontinuity. To prevent this failure, an HLLC scheme [3] has been derived, that resolves the contact discontinuity and satisfies the physical constraints.
To complete this presentation of the numerical approximation, we mention that a suitable high order extension that preserves both the positivity and the flux limitation can be derived, relying on an appropriate limitation procedure.

An HLL scheme for the free transport M 1 angular moment system
In this section, we derive a Finite Volume method, issued from the HLL method [20] to discretize the free transport equation (23). Put in other words, we omit the source terms and we consider the one dimensional generic conservative system where U = ψ 0 ψ 1 and F stands for the flux of the M 1 system in the x space direction. We consider a structured mesh of size ∆x i , defined by the cells I i = x i−1/2 , x i+1/2 , where we have set x i+1/2 = x i + ∆x/2 , i ∈ Z, at time t n . As usual, we consider known a piecewise contant approximation U h (x,t n ), defined by U h (x,t n ) = U n i , x ∈ I i , ∀i ∈ Z. At initial time t = 0, we impose where U 0 is the initial data. This approximation evolves in time, involving a suitable approximate Riemann solver. In the HLL approach, the exact Riemann solver solution is substituted by a single approximate state (see Figure 2). Here b L and b R are relevant approximations of λ L and λ R , respectively. Let us introduce the proposed approximate solution: Moreover, the search of weak solutions leads to the Rankine-Hugoniot jump conditions At each interface x i+1/2 , we impose the above HLL approximate Riemann solver, assuming the CFL like condition (29) ensuring that the Riemann solvers do not interact in the case where b L,i+1/2 < 0 and b R,i−1/2 > 0 : We set U h (x,t + ∆t), at time t n + ∆t, the superposition of the non-interacting Riemann solutions. We define the updated approximation at time t n+1 by An easy computation gives where The robustness of the scheme, namely the positivity, the flux limitation, the total mass preservation, has been established for the HLL scheme (see [8] for further details). Finally, concerning the high order extension, we adopt a van Leer MUSCL technique [37], supplemented by a suitable slope limitation preserving these expected physical properties [7].

An accurate HLLC scheme
The HLL scheme has proved to be robust, however, its 2D extension fails when approximating contact waves. Several works [3,8] introduce a more accurate scheme, the HLLC scheme, based on a two state approximation, denoted by U * L and U * R . First, let us recall the relevant linearization that permits us to define an approximation with two intermediate states: on the one hand, the Rankine-Hugoniot conditions (27) are considered; on the other hand, they are supplemented by the continuity of the Riemann invariants accross the contact wave: where β x and Π are defined by the relation (24). The combination of both the Rankine-Hugoniot condition (27) and the relation (31) standing as the continuity of the Riemann invariants accross the contact wave, is sufficient to determine uniquely [3,8] the two approximate states U * L and U * R , together with their associated fluxesF L andF R . The proposed HLLC approximate solution can be written as Similar to the derivation of the HLL scheme, we integrate over a cell I i the juxtaposition of the non-interacting HLLC Riemann approximate solvers at each interface (projection step), in order to obtain the updated quantity This brief description of the HLLC scheme is now completed. It is able to capture exactly the contact wave, and satisfies the positivity, the flux limitation, and the total mass preservation.

Central Void
The first test case is taken from the medical physics literature [2]. We consider only elastic scattering, which is modeled by the Henyey-Greenstein kernel. Thus S = 0 and T in = 0. We compare the particle flux ψ (0) (x) obtained with the minimum entropy model (labeled M1) with a discrete ordinates solution of the transport equation (labeled SN) with sufficiently many angles (128). The method has been described in [16]. The test case consists of a one-dimensional geometry with three layers: optically thick, followed by optically thin followed again by optically thick. The layers have an equal depth of 40 mm. The scattering and absorption coefficients are σ s = 0.5 mm −1 , σ a = 0.005 mm −1 for the optically thick region, and σ s = 0.01 mm −1 , σ a = 0.0001 mm −1 for the optically thin region. Moreover, g = 0. Figure 3 shows the particle flux ψ (0) as a functon of space. Compared to the benchmark solution, the minimum entropy model slightly overestimates but nevertheless quite accurately describes the particle flux. In Figure 3 we also show the partice distribution function ψ(x, Ω), where Ω = (0, 0, µ) in 1D. The main difference is that for M 1 , the forward-peak of the incoming particles reaches further into the medium.

Two-dimensional Void-like Layer
Our second test case, again taken from [2], is a two-dimensional quadratic domain which contains a void-like layer, shown in gray in Figure 4(a). Again, we consider only elastic scattering modeled by the Henyey-Greenstein kernel.
We take σ s = 0.5 mm −1 and σ a = 0.005 mm −1 inside the square, and σ s = 0.01 mm −1 and σ a = 0.0001 mm −1 in the void-like ring. In both regions, g = 0. An isotropic source of particles is placed on the left boundary. In a 2D contour plot (Figure 4), the fluxes ψ 0 from the discrete ordinate method and from the minimum entropy method are virtually indistinguishable. The propagation into the medium, as well as the void-like layer are equally well resolved. A difference between the models only becomes apparent in a logarithmic plot of a cut through the center of the square at y = 50 mm. Figure 4 shows the particle flux along this line. The difference between both solutions is again of the order of one percent.

Electrons on Water Phantom
As a first test case that includes energy loss, we consider a 10 MeV electron beam impinging onto a slab of water. In Figure 5 we compare the results computed with our code to the dose computed by the state-of-the-art Monte Carlo code PENELOPE [32]. This code has been extensively validated against experimental results.
To obtain a good fit with the tabulated scattering data, we have fixed our model parameters for water as ε B = 16.0 eV, Z = 9.40, ρ el = 0.256 × 10 23 g/cm 3 , ρ in = 6.21 × 10 23 g/cm 3 . These parameters are directly inserted into the model (3), and subsequent derived models issued from (3). As boundary conditions, we have taken a very narrow Gaussian in energy, and a δ pulse in angle and computed the angular moments. PENELOPE was set up in a pseudo-1D setting with a large beam size perpendicular to the beam direction. In order to compare the different formulations of the models, both depth-dose curves in     dependence of ψ in the M 1 model. Both possibilities will be investigated further. However, we believe that this result can serve as a proof of concept of a PDE based modeling of dose computation.

CT Data
In our final test case, we compare our method with Monte Carlo results from PENELOPE using real patient CT data showing the hip bone. We took a two-dimensional slice of 6 × 6 cm from the three-dimensional CT data. A square region is split into 64×64 squares. In each of the squares, the material is described by its Hounsfield grey value G (x, y). The grey values can be translated into physical parameters as follows, ρ(x, y) = G (x, y) 1000 i.e. the densities ρ el and ρ in for water are multiplied by a specified factor. The region shows the hip bone and the density varies between 86% and 226% of the value of water. The boundary conditions were set up similar to the previous case, with three beams of width 2 cm, each consisting of 10 MeV electrons, impinging from the centers of three sides of the domain. Contour plots of the dose distribution are shown in Figure 7. There, we also show two cuts through the dose distribution. looking at the 2D dose distribution, the contour lines agree very well. Note that, although we have used 3 × 10 10 particles, there still is significant noise in the Monte Carlo results. The two cuts through the beam centers show that also quantitatively the independently computed dose distributions agree very well.  The computation time for the 3D Monte Carlo dose was 3 × 29 hours for 3 × 10 10 particles on a 3GHz Pentium 4 with 1 GB RAM. In 1D, the minimum entropy model took about 1 second, in 2D 4 seconds. Thus we expect a computation time of several seconds in a full 3D dose computation.
Again, this result shows that if our model is developed further, it may serve as an alternative to existing dose computation methods.