Spatial optimization for radiation therapy of brain tumours

Glioblastomas are the most common primary brain tumours. They are known for their highly aggressive growth and invasion, leading to short survival times. Treatments for glioblastomas commonly involve a combination of surgical intervention, chemotherapy, and external beam radiation therapy (XRT). Previous works have not only successfully modelled the natural growth of glioblastomas in vivo, but also show potential for the prediction of response to radiation prior to treatment. This suggests that the efficacy of XRT can be optimized before treatment in order to yield longer survival times. However, while current efforts focus on optimal scheduling of radiotherapy treatment, they do not include a similarly sophisticated spatial optimization. In an effort to improve XRT, we present a method for the spatial optimization of radiation profiles. We expand upon previous results in the general problem and examine the more physically reasonable cases of 1-step and 2-step radiation profiles during the first and second XRT fractions. The results show that by including spatial optimization in XRT, while retaining a constant prescribed total dose amount, we are able to increase the total cell kill from the clinically-applied uniform case.


Introduction
Glioblastomas are the most aggressive, and unfortunately most common, form of primary brain tumour [1][2][3][4][5]. They are characterized by rapid growth and invasiveness, yielding survival times that seldom exceed a year [6]. Because of this, treatments for glioblastomas are swift and aggressive, usually involving a combination of surgical intervention, chemotherapy, and external beam radiation therapy (XRT). Furthermore, the tendency for recurrence of glioblastomas after surgery makes postoperative chemotherapy and XRT a crucial part of effective treatments. Although current treatment plans do often extend survival time, they are far from perfect and leave much room for improvement. However, while these efforts focus on optimal scheduling of radiotherapy, they do not include a similarly sophisticated spatial optimization.
Non-uniform dose distributions are not a new concept in radiation oncology. Phenotypic variations across the volume of a tumour can result in differing levels of radiation effectiveness. In particular, hypoxic regions cause a difference in cell radio-sensitivity, making radiation less a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 effective in those areas. This leads to a technique called "dose-painting" in which different regions are prescribed a different dose to combat the reduced effect (see [7,8] for an overview). Dose painting divides the tumour into a small, discrete number of regions allowing a different dose to be prescribed to each area. Additionally, we may need to apply non-uniform radiation to areas with different cell density, irrespective of radio-sensitivity.
Several previous works have addressed the dependence of the optimal beam shape on the density profile of the tumour, (see [9][10][11][12] for example). A possible criterion for optimality (used in our calculation) is minimizing the total number (N) of tumour cells remaining after application of XRT; another is to minimize the tumour control probability (TCP), e −N (see below, Eq (7)). Naturally either minimization must satisfy a number of physical constraints. The constraint on the total radiation dose (used in the paper) leads to a particular shape of the beam, which was also obtained in previous studies [9,11,12]. By contrast, a constraint on the 'mean dose,' weighted by the local cell density, leads to a uniform beam profile [10].
In this paper, we study the following question: given a maximum allowable total dose to administer, what dose profile results in the maximum global TCP. As an equivalent metric to the TCP, we instead minimize the total number of surviving tumour cells after radiotherapy. We then go further, making our conclusions clinically reasonable by discretizing this optimal profile and considering multiple radiation fractions. While we do not deal with tumour heterogeneity in particular, we do include an alternative death mechanism which can account for different radio-sensitivies resulting from tumour heterogeneity such as hypoxia. Other works have examined optimization of radiation therapy incorporating normal tissue complication probability (NTCP) in addition to TCP [8,13]. The NTCP is an important measure to consider as it quantifies the chance of problems arising in nearby organs at risk. While we do not include NTCP in this work, an interesting topic for future study would be to incorporate TCP and NTCP in a single mathematical optimization.
Naturally, when discussing optimization in a given system, mathematical modelling is an invaluable tool. In the context of brain tumours, mathematical models have been widely used by many different investigators (for example, [14][15][16][17][18]). While originally developed for investigations of brain tumours, we anticipate that these models and the following approach and results can be easily generalized to study other types of tumours. In the following, we build upon a host of previous works (see [19] for a review) which show that the natural growth of glioblastomas can be well-described by two governing parameters, ρ and D n , which describe growth and invasion processes respectively. The two mathematical models most commonly used to describe growth of tumours are the so-called exponential and logistic growth laws. The more general logistic growth naturally provides a better description of tumour proliferation and stabilization. However, during the relatively short time frame before the growth of the tumour begins to plateau, exponential growth is a reasonable approximation. Including a linear (Laplacian) diffusion term, commonly used to model tumour invasion, we arrive at the well-established equation for tumour growth @nðx; tÞ @t ¼ D n r 2 nðx; tÞ þ rnðx; tÞ 1 À anðx; tÞ n max Here, nðx; tÞ is the tumour cell density at positionx ¼ ðx 1 ; � � � ; x d Þ, and r 2 ¼ P d a¼1 @ 2 =@x 2 a is the Laplacian operator. We have introduced d as the number of dimensions; d = 3 in three dimensions, while for certain computations we focus on the two dimensional case of d = 2. Exponential growth corresponds to a = 0, while a = 1 leads to logistic growth to a maximal density n max .
Prior to treatment, two MRIs are commonly performed: one diagnostic and one to aid in treatment planning. Using the measures for velocity of growth and tumour abnormality obtained from these two scans alone, the tumour-specific values for ρ and D n can be estimated [18].
We begin by presenting the general mathematical model, followed by studies of a number of special cases. First, the optimal continuous profile is derived for the cases of one and two radiative fractions with both exponential and logistic cytotoxic action. The insights from these cases are subsequently extended to show that the optimal tumour cell density is uniform for any radiation procedure. From there, we move to optimization of non-continuous radiation profiles, and derive the optimal scenario for the cases of one-step and two-step profiles for single or multiple fractions, constrained individually or together. Finally, the non-continuous optimization is extended to include a logistic cytotoxic action for one and two step radiation profiles.

Model derivation
Consider a tumour that has evolved according to Eq (1) for some time, such that its cell density profile is given by nðx; t 0 Þ. Now we intend to apply XRT to this tumour. To do so, we introduce a function f ðx; tÞ which we call the cytotoxic profile. The action of most therapeutic interventions is to remove a fraction of existing cells, which is incorporated in our model by adding a term À gf ðx; tÞnðx; tÞð1 À bnðx; tÞ=n max Þ to the right hand side of Eq (1). For radiation, the parameter γ is a measure of the radiation rate, and can be written as g ¼ a _ D (1/day) where α is the linear coefficient in the well-known Linear-Quadratic model for radiation efficacy, and _ D is the dose rate applied to the tumour. For simplicity, we do not include the quadratic term of the Linear-Quadratic model into the model as it does not affect the qualitative shape of the radiation profile, which is the focus of this study. This is clear mathematically since the inclusion of the quadratic component simply perturbs the value of γ in Eq (2) below, which acts as a scaling constant. The parameter b, much like a in the growth term, simply differentiates between exponential (b = 0) and logistic (b = 1) cell killing. It is also known that saturated tumours with a high cell-density and low proliferation are less affected by radiation than low cell-density tumours. A logistic death term is mathematically able to capture this behaviour as a higher cell density reduces the magnitude of the final term. This addition modifies the governing equation to Normal XRT occurs in a series of bursts called fractions. We therefore impose that f ðx; tÞ is an on-again-off-again function in time such that f ðx; tÞ ¼ 0 for all times except during the scheduled fractions. This greatly simplifies the solution to the nonlinear partial differential equation, as during the time interval Δt over which each fraction is applied, the processes of cell division and spreading have little effect on the cell density. The natural growth of tumours proceeds on scales of days, months, or years, while radiation fractions occur over a scale of minutes. Thus, during each fraction, the first two terms on the right hand side of Eq (2) can be mathematically neglected, leading to @n @t Spatial optimization for radiation therapy Integrating this (now ordinary) differential equation in the interval t 2 [t 0 , t 0 + Δt], leads to in the exponential case (b = 0), and in the logistic case (b = 1). We focus on the first and second fractions of XRT and further impose a simple upper bound for f ðx; tÞ to adhere to patient safety standards. We write this constraint as 0 � f ðx; tÞ � C for some C. The final constraint on f ðx; tÞ limits the total dose received by the patient. This constraint is mathematically represented by where the time integral is over the entire treatment length. This constrains the total radiation dose over all fractions (where we expect the inequalities to be saturated to achieve maximal effectiveness).
Our goal is to determine the function f ðx; tÞ that minimizes the total number of tumour cells after the final fraction, N(T), obtained by integrating the tumour cell density as To contrast, Brahme and Agren [9] instead optimized the TCP, which they defined as (using our notation) TCP = e −N(T) where N(T) is the number of cells surviving the treatment. Regarding nðx; t 0 þ DtÞ as the mean of a probability density for cells following XRT, TCP is the probability that all cells are exterminated. Clearly, one can see that minimizing N(T) is equivalent to maximizing e −N(T) . Furthermore, for the simplest case, we indeed achieve the same result.

Continuous profile optimization
We begin by analysing continuous cytotoxic profiles for one fraction of radiation. In the case of exponential death, the result matches that of previous studies, and from it we can derive a precise shape of the profile. In the case of logistic death, the resulting profile is more complex, but remains mathematically consistent with the previous case. A similar analysis of the second fraction of radiation is significantly more challenging as the effect of normal tumour progression must be taken into account. While it is possible, the results are far less interesting since the optimal profile will leave a uniform (or flat-topped) cell density (proof of this in supporting information). The second fraction is mathematically examined in the supporting information for the case of exponential death. The second fraction for the logistic death case.
Optimal profile with one fraction of exponential death. To proceed in the case of exponential cytotoxic action, we impose the first constraint in Eq (6) with a Lagrange multiplier λ, which requires extremizing where we have set γ = Δt = 1 for convenience. Solving the resulting Euler-Lagrange equation for f ðx; t 0 Þ, we find the optimal profile with λ chosen such that Eq (6) is satisfied. Using this cytotoxic profile would give us the optimal cell kill from the radiation fraction. Not surprisingly, following the simplifications leading to Eq (9), the above result is independent of the parameters ρ and D n . The appearance of the logarithm is merely a consequence of the killing effect of therapy being proportional to the number of existing cells, with the cytotoxic profile appearing in the exponent of Eq (4). This result has been previously derived by many others ( [9][10][11] for example). An interesting consequence of this optimal profile is that it leaves the resulting cell density as a uniform distribution. Stavreva et al [10] applied an extremization of TCP subject to a constraint on the mean dose defined as D v / R nðx; t 0 Þf ðx; t 0 Þd d x, i.e. weighting the dose according to the local cell density. It is clear that if this constraint is used in Eq (8) in place of the unweighted net dose, the resulting beam profile will be uniform.
While a useful starting point, the cytotoxic profile in Eq (9) is not guaranteed to satisfy the constraints of 0 � f ðx; tÞ � C. In particular, it leads to unphysical negative values when nðx; t 0 Þ < l. To better understand the limitations of this result, and how to overcome them, let us consider the simplest case of a Gaussian profile arising from radially symmetric growth of a single cell in exponential growth: assuming that the tumour begins with a single oncogenic transformation or single-cell metastasis, modelled by a delta-function cell density, exponential growth for a time t 0 leads to the cell density profile where r is the radial distance from the initial cell (tumour center). The width of the Gaussian profile is s ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 2D n t 0 p , while n 0 ¼ e rt 0 =ð2ps 2 Þ d=2 is the cell density at its center. Eq (9) now leads to a parabolic cytotoxic profile. Cutting off the negative portions of the parabola leads to the (semi-circular) profile We have introduced the parameters f m = ln(n 0 /λ) and r m ¼ s ffi ffi ffi ffi ffi ffi ffi 2f m p to indicate the maximum value of the cytotoxic profile, and the radius over which it is applied, respectively. The total radiation dose in this fraction is then given by where S d is the d-dimensional solid angle, with S 3 = 4π and S 2 = 2π. Using f m ¼ r 2 m =ð2s 2 Þ from Eq (11), we conclude that the optimal radius of the semi-circular beam is given by while its maximal intensity equals If the above value of f m exceeds the maximum allowed intensity of C, we should instead use with r m (F) now constrained by the maximal allowed density according to Eq (12). In either case, the cell density profile after application of XRT will attain a flat-top profile, as A flat post-XRT profile is also predicted for any initial tumour cell density nðx; t 0 Þ, although the volume over which the beam is applied will be different. In the best outcome, the density profile will be below the (single-cell) threshold for future growth. If not, additional fractions need to be applied. Optimal profile with one fraction of logistic death. Mathematically, the switch from an exponential to a logistic cell death does not change much. We now simply use Eq (5) instead of Eq (4) in our optimization. Again imposing the constraint of Eq (6) using a Lagrange multiplier and setting γ = Δt = 1 for convenience, we arrive at the optimal profile 2l þ n max þ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi n max ðn max þ 4lÞ Like the previous single-fraction case, this result is independent of ρ and D n . It also is not guaranteed to satisfy the bounds of 0 � f ðx; tÞ � C. As a quick check on the validity of this result, one can let n max ! 1 and see that Eq (9) is recovered. Note that the optimal profile again has a uniform or flat-toped form. This can be seen by substituting Eq (17) into Eq (5).

Discrete profile optimization
Computational methods exist for determining how to administer a heterogeneous dose, primarily sub-volume boosting and dose painting by numbers [7,8]. Unfortunately, these methods are limited by the mathematical optimization techniques as well as the physical reasonability of their results. Specifically, current technology is not capable of producing beams with a high degree of precision for clinical use such as those derived above in the Continuous Profile Optimization. The standard procedure is to coalesce several beams on the location of the tumour, creating an area of high radiative strength. As such, a reasonably practical non-uniform beam profile is a step function. In the following, we consider cases of both one-step and two-step radiation profiles. We emphasize that our aim is to explore the feasibility of spatial optimization and to gain qualitative insights, and thus make mathematical simplifications throughout to reflect this. The first simplification is to consider radially symmetric profiles in two dimensions. In the following, we assume exponential growth (a = 0) leading to a Gaussian profile, as in Eq (10). A summary of the parameter values used in calculations can be seen in Table 1.
D1: One-step radial profile. The simplest step-function case of XRT involves a uniform beam of radius r 1 and strength f 1 , applied for a duration Δt at time t 0 , i.e.
The goal is minimize Approximating the PDE as an ODE as before, the tumour cell density distribution immediately after the fraction is obtained as Integrating this result gives the total number of cells as The constraint on the total beam flux can be imposed through a Lagrange multiplier λ to create an augmented N, Extremizing with respect to r 1 , f 1 , and λ leads to the following system of equations: After eliminating f 1 and λ from the above equations, we arrive at the following implicit expression for r 1 : Values of r 1 that satisfy Eq (25) can be used to find a corresponding f 1 , together specifying the optimal N(r 1 , f 1 ). We assume that the duration of radiation is approximately 10 minutes, or Δt = 0.007 (days), and use the linear model for the radiation g ¼ a _ D, where α is the radiobiological parameter and _ D is the dose rate. For α = 0.08 (1/Gy) (taken from an average of the values obtained in [17]) and the standard dose rate of 5 (Gy per radiation time), we obtain γ = 60 (1/day). Using the parameter values C = 2.5 (corresponding to the maximum allowed dose rate of 5 (Gy per radiation time)), F 0 = 25 (mm 2 ), and a range of (n 0 , σ) pairs chosen such that N(t = t 0 ) = 10 7 (cells) is constant, we can solve for r 1 and calculate the corresponding f 1 . The location of the optimal radius r 1 is plotted in S1 Fig. As expected, larger values of σ require radiation over a wider radius. However, the increase of r 1 with σ is sub-linear, and quite well fitted by r 1 / σ 1/2 . This is precisely the scaling behavior predicted in Eq (13) for the semi-circular beam shape in d = 2. The scaling of the optimal radius with tumour size thus appears to be robust, irrespective of the beam shape. This procedure is done for 5 different (n 0 , σ) pairs in Table 2 and the profiles can be seen in S3 Fig (Supporting Information). Note that for the σ = 1 and σ = 2 cases, the optimal r 1 falls below ffi ffi ffi ffi ffi 10 p (mm). However, due to our constraints of F 0 ¼ 25 ¼ f 1 r 2 1 and 0 � f(r, t) � C = 2.5, we have a lower bound of r 1 � ffi ffi ffi ffi ffi 10 p (mm). Also note that the calculated values in Table 2, and in the following data tables, are truncated decimals, but the values of N are calculated using more precise solutions. Therefore, inserting the values for r 1 and f 1 found from Table 2 into Eq (20) will not necessarily exactly reproduce the listed values of N. Optimal one-step profile for one fraction of radiation with exponential growth and death.n 0 and σ are the magnitude and standard deviation of the initial Gaussian tumour cell density, which starts with 10 7 total cells. D2: Two-step radial profile. Now we consider the more elaborate example of a 2-step cytotoxic profile, but still applied in only one fraction. Introducing two new variables, r 2 and f 2 , the radiation profile is: where f 2 applies a different dosage to the outer region of the tumour. Adding the second radial arc modifies the constraint on dosage to We need to minimize the total cell number, N, with respect to the four-parameter set (r 1 , r 2 , f 1 , f 2 ). Using the Lagrange multiplier method as above results in a system of 5 coupled transcendental equations. Unfortunately, this set of equations has many local extrema, making the result heavily dependant on the initial guess used in the computational solver. This makes identification of the global extreme difficult. To avoid this, we use a Monte Carlo method to identify the optimal shape of the radiation beam. First, a random radiation beam is generated by selecting a point from the parameter space (r 1 , r 2 , f 1 , f 2 ). The selection of this point is made by randomly assigning values to 3 of the parameters within their acceptable ranges, then calculating the fourth according to Eq (27) such that the dose constraint is satisfied. This generated parameter set defines the candidate radiation beam. The total cell number resulting from each beam is calculated and compared to the previous minimum. This is done many times (*10 9 in our simulations) such that N converges to a global minimum. The point in the parameter space which generates this minimum is then inserted back into the equations for the optimal profile to check that they are indeed satisfied. Such a point parameterizes the optimal profile (in subsequent optimizations the parameter space will become more complex, however this method will remain the same). The numerical values resulting from this method are summarized in Table 3, and the resulting cell density profiles are shown in S4 Fig (Supporting Information). As a check on the numerical optimization, we do find that the the addition of the second radial arc leads to better treatment outcomes. Once the constraint on maximal value of f 1 = 2.5 is no longer operative (for σ > 2 mm), the resulting cell density profiles are close approximations to the flat-top profiles expected from a semi-circular beam. The optimal two step radial profile thus represents a crude approximation to the semi-circular beam. Optimal two-step profile for one fraction of radiation with exponential growth and death. σ is the standard deviation of the initial Gaussian tumour cell density, which starts with 10 7 total cells. r 1 , r 2 , f 1 and f 2 are the two radii and strengths of the steps in the cytotoxic radiation profile. N(t + dt) is the final tumour cell number at the end of radiation. D3: Two fractions individually constrained. We next consider application of a second fraction of radiation, with dosage separately constrained on each fraction. Describing the second fraction requires introducing four new variables to fully parametrize f(r, t). We modify the previous notation by adding a second index to each of the radii and strengths of f(r, t), with the first index indicating the step-number and the second corresponding to the fraction number. Thus, r 1 , r 2 , f 1 , and f 2 from before become r 11 , r 21 , f 11 , and f 21 respectively, with corresponding r 12 , r 22 , f 12 , and f 22 for the second fraction. The interval between the two fractions will be labelled τ, such that the second fraction takes place over the interval [t 0 + Δt + τ, t 0 + 2Δt + τ]. Note that we assume the lengths of the two fractions to be the same. Defining t � 0 2 ½t 0 ; t 0 þ Dt� and t � 1 2 ½t 0 þ Dt þ t; t 0 þ 2Dt þ t�, we can write f(r, t) during the separate fractions as In this section, we consider the fractions constrained individually, such that Since our constraint on f ðr; t � 0 Þ is the same as before, the optimal f ðr; t � 0 Þ does not change from that in Table 2. We can use the same optimization procedure for the second fraction, but with a modified starting cell density profile. In order to deal with profiles with simple analytic expression, we further assume that the time interval τ between the fractions is small enough to neglect spatial migrations described by the diffusion term (mathematically, this can be expressed as the condition D n τ � σ 2 , which can be derived from a simple scale analysis). If so, the density profile simply grows exponentially, by a factor e ρτ without changing its spatial form, and immediately before the second fraction is given by nðr; t 0 þ dt þ tÞ ¼ n 0 e rt e À r 2 2s 2 e À gf 11 Dt 0 � r � r 11 n 0 e rt e À r 2 2s 2 e À gf 21 Dt r 11 < r � r 21 n 0 e rt e À r 2 2s 2 The density profile immediately after application of the second fraction is then given by: For each set of radii, r 11 , r 21 , r 12 , and r 22 , the cell density is a piecewise continuous function. The total number N can then be obtained as before by integration, as an explicit analytic expression. For the same initial combinations of (n 0 , σ) as above, we search for the radii (r 21 , r 22 ) that optimize the second fraction using our Monte Carlo method (the optimal radii (r 11 , r 12 ) are naturally the same as obtained previously in Table 3). The results of this minimization are summarized in Table 4.
Note that if diffusion is ignored, the final cell density after a second fraction of radiation is given by nðxÞ ¼ n 0 ðxÞ exp frt À g½f 1 ðxÞ þ f 2 ðxÞ�Dtg. From this expression it follows that in this limit, the order in which the two fractions are applied is not important. Indeed the same conclusion applies to any number of fractions, each separately optimized. Thus it is not necessary to use the XRT profile that is optimal at the time of its application, as long as there are planned future fractions that boost the overall amount of radiation at each point to the optimal value. This freedom provides an additional tool for therapeutic planning.
D4: Two fractions with overall constraint. As a final example within this class, with nondiffusive exponential growth between the two fractions, we consider the case when the overall dose is constrained, i.e. as The optimization procedure can be carried out as before through a Monte Carlo search. However, each step of the search now involves an exploration in the space of 8 variables (as opposed to separate searches in a space of 4 variables), subject to one constraint. Note also that (except for the replacement of 2F 0 for F 0 ) this is exactly the search that would be performed for a single fraction with a 4-step profile. The optimal values of these parameters are given in Table 5.
Observe that for each σ, the optimal N with the single overall constraint is either better, or the same as, in the individually constrained case. From the perspective of optimization, this is not unexpected as the latter also explores the subset of the space available to the former. In view of this, the surprising result may appear to be that for σ = 1 and 2 (see Table 6) the minimum occurs in the separately constrained space. However, the structure of the general optimization problem is such that for D n = 0, the same result should hold with one or more constraints (supporting infomation). We are thus unable to conclude if the unequal partition of flux between the two fractions (in cases of σ = 3, 4, 5) is correct or a computational artifact (we indeed find many solutions close to the optimum, so finding the true optimum requires considerable computation and precision). Optimal two-step profile for two separately constrained fractions of radiation with exponential growth and death. σ is the standard deviation of the initial Gaussian tumour cell density, which starts with 10 7 total cells. r 11 , r 21 , r 12 , r 22 , f 11 , f 21 , f 12   Optimal two-step profile for two mutually constrained fractions of radiation with exponential growth and death. σ is the standard deviation of the initial Gaussian tumour cell density, which starts with 10 7 total cells. r 11 , r 21 , r 12 , r 22 , f 11 , f 21 , f 12  D5: Tumour densities from logistic growth. The Gaussian density profiles employed to model the initial cell densities above are appropriate only for small undeveloped tumours. For contrast, in this section we consider density profiles resulting from the logistic growth (a = 1 in Eq (1)). Such profiles are obtained by evolving the partial differential equation starting with a localized initial condition. While there is no exact analytic form for the resulting solution, as an analytical approximant, we use the form where a 1 , b 1 , and c 1 are fit parameters. Logistic growth causes the tumour to form a flat-top profile as the density in the centre approaches n max . Eq (31), which is known as a Fermi Function, also describes a flat-top form, which is why we use it for fitting. A numerical implementation of Eq (1) was used to generate the fit parameters. For this procedure and more information, see supporting information. The values of the fitted parameters are given in Table 7. We consider the same radial step-functions as previous, and the optimization procedures are carried out in the same manner as before. The results for application of a single 2-step fraction are reported in Table 8. If this treatment is followed by a second, separately constrained 2-step fraction, the resulting parameters for the second application are given in Table 9.
D6: Tumour densities from logistic growth and logistic death. The previously examined cases of exponential cytotoxic action do not incorporate the phenomenon of radiation having a larger effect on faster-proliferating cells than on slower-proliferating cells due to the cell's position in the cell cycle and factors such as oxygen concentration [22]. If we instead consider logistic action (b = 1), then some aspects of such tendency are reproduced. This qualitative Distribution of radiation flux over two radiation fractions based on standard deviation of initial tumour cell density.
The '1st Fraction' number represents the amount of dose (out of 50) that is allocated to the 1st radiation fraction and the '2nd fraction' the amount allocated to the 2nd radiation fraction.
https://doi.org/10.1371/journal.pone.0217354.t006 Spatial optimization for radiation therapy behavior is important to incorporate because the vast majority of tumours exhibit some form of treatment resistance. This resistance is conferred to a tumour through heterogeneities of various traits across its volume, such as cell-type, phenotypic expression, and stem-ness. Of particular interest to radiation therapy, tumour hypoxia is a prominent feature that leads to increased radioresistance. As we will see, inclusion of logistic cytotoxic action reproduces this behavior and dramatically changes the results. We note that this is a purely phenomenological effect and not meant to describe the underlying biology at play. We make the same simplifications as before leading to the piecewise equivalent to Eq (5), The results of the Monte Carlo optimization are summarized in Table 10. Optimal one-step profile for one fraction of radiation with logistic growth and exponential death. σ is the standard deviation of the initial Gaussian tumour cell density, which starts with 10 7 total cells. r 1 , r 2 , f 1 , and f 2 are the radii and strengths of the steps in the cytotoxic radiation profiles. N(t + dt) is the final tumour cell number at the end of radiation.
https://doi.org/10.1371/journal.pone.0217354.t008 Optimal one-step profile for one fraction of radiation with logistic growth and exponential death. σ is the standard deviation of the initial Gaussian tumour cell density, which starts with 10 7 total cells. r 11 , r 21 , r 12 , r 22   Optimal two-step profile for one fraction of radiation with logistic growth and death. σ is the standard deviation of the initial Gaussian tumour cell density, which starts with 10 7 total cells. r 1 , r 2 , f 1 , and f 2 are the radii and strengths of the steps in the cytotoxic radiation profiles. N(t + dt) is the final tumour cell number at the end of radiation. https://doi.org/10.1371/journal.pone.0217354.t010 Spatial optimization for radiation therapy Mathematically, the results here differ dramatically from the exponential case in which the optimal profile focused the radiative strength in the center of the tumour where the cell-density was larger. With the incorporation of logistic cell death, the effectiveness of radiation in the center of the tumour diminishes. The optimal profile must now balance this with still attacking areas with high cell densities. We see that in the case with the smallest initial spread of cells, the optimal cytotoxic profile achieves this by focusing strength in a ring through the middle of the tumour. For the other two cases, the density was far more spread out initially, leading to an optimal profile which focused on cells at the outside of the tumour. This was because the effect of the logistic death was stronger than the decreased cell density.
In the case of exponential death, we were able to continue our discussion and optimize the second fraction as we did the first. In the logistic death case, the computations in the optimization method become significantly more challenging. For that reason, we do not pursue the second fraction of radiation here. Table 11 shows a summary of the discrete optimization cases undertaken in this work. It outlines the method, constraints, and details of each of the six cases.

Conclusion and discussion
In this paper we pose the question of how to spatially shape a sequence of XRT treatments to best eliminate tumour cells. To answer this question, we need to know (i) how the tumour grows in time; (ii) how it responds to treatment; and (iii) what constraints apply to radiation dosage. Answers to all questions need to be expressed in mathematical terms, which necessitates simplifications and approximations. We have relied on assumptions and mathematical models commonly used in the literature, and hope our general results are insensitive to choice of model.
The most important result follows from the assumption that the effect of XRT is to destroy a fraction of existing tumour cells, proportional to its strength. This assumption is embodied by the term À gf ðx; tÞnð1 À bn=n max Þ in Eq (2). It immediately leads to the conclusion that the optimal beam profile should depend non-uniformly on the tumour density at the time of its application. Shaping the beam to such a precise form is likely impossible, and may further violate constraints on the strength of the beam. As such, optimal beam profiles can be sought within certain constraints on the shape of the beam. In particular, we consider XRT profiles obtained by superposing several radially symmetric beams. Clinically, the normal procedure is to irradiate the entire radius of the visible tumour with a uniform beam. As can be seen from the various tables above, a uniform f(r, t) will not lead to optimal cell kill. Using a spatially optimized f(r, t) can reduce the total number of cells within a tumour significantly from what is currently done clinically. The metric of TCP is not the only quantity that can be used to measure radiation efficacy. Others such as normal tissue complication probability (NTCP), conformity index, homogeneity index, and survival time have been used by other researchers [23][24][25][26]. We focused our analysis here on TCP as it is intuitively simple and mathematically straightforward. We stress that this is an exploratory study meant to discern qualitative behaviour, and note that a more complete study should include analysis and comparison of other methods of radiation treatment planning. More work in this direction would be interesting and similar mathematical tools could theoretically be used.
For more than one fraction of XRT, it is necessary to account for the growth of the tumour in between the two treatments. A simple commonly used model is logistic growth, depending on three parameters ρ, D n and n max . In the various studies above, we have mostly neglected the change in the shape of the tumour between two fractions. However, this needs to be included in a more comprehensive study. Given all our results, we can propose a procedure for spatial optimization of radiation as follows: • Image tumour twice. The change in shape of the tumour can then be used to model its growth mode; e.g. in order to deduce the parameters ρ, D n , and n max . The most important ingredient for shaping the optimal beam is the cell density nðxÞ.
• Determine dose prescription and treatment schedule for the tumour. This provides constraints on individual and total allowed dosage, {F n } and the time intervals between fractions, {τ n }.
• Determine the physical limitations of the radiation apparatus; e.g. how many radial steps is the apparatus capable of accurately producing and superposing.
• Optimize the first radiation fraction using the above optimization procedure.
• Use the growth model deduced from the initial two images to model the growth of the tumour cells between fractions. The deduced cell density profile prior to each fraction can then be used as input to optimize the shape of the beam in that fraction. The assumptions on evolution of the profile after a fraction and in between fractions are likely to introduce errors. Ideally, more imaging can be done to obtain more precise cell density profiles at intermediate times.
This exploratory work examines spatial variations in the shape of radiation beams, and proposes improvements to radiation therapy by shaping beams according to the pre-treatment cell-density. The advantage of this method is the opportunity to reduce the total tumour cell number below that obtained in uniform radiation. In addition to the initial tumour density profile, the method relies only on few tumour-specific parameters in order to carry out the optimization. In its current implementation, the method does not take into account other factors in radiation planning such as nearby organs at risk, scheduling, or other scoring indices used to compare dose distributions. Such factors could potentially be included in a more elaborate model, a direction that we hope to pursue in future. A major limitation of the model is the exclusion of heterogeneity across the tumour volume as tumours usually vary in cell type, stemness, phenotypic expression, and perhaps most importantly, hypoxia. Although we added a logistic cytotoxic action term to account for different radiation efficacy across the area, this is relatively simplistic and could be improved with a more thorough specification of tumour heterogeneity. Another direction that was not examined here is the interaction of this spatiallyvarying radiation treatment with other therapies such as chemotherapy. We hope that this exploratory work sheds light on the spatial variation of radiation beams and inspires further work in the area. Cell densities resulting from the optimal radiation profiles from one 1-step fraction of exponential growth and death for the initial density deviations of (from top left to bottom right) σ = 1, σ = 2, σ = 3, σ = 4, and σ = 5. The red dots are those below the detectable threshold in imaging, of 5 cells/mm 2 , and the blue dots those above. Note the different y-axis scales to highlight the resulting shape. (EPS) S4 Fig. Cell densities from two one-step fractions of exponential growth and death. Cell densities resulting from the optimal radiation profiles from one 2-step fraction of exponential growth and death for the initial density deviations of (from top left to bottom right) σ = 1, σ = 2, σ = 3, σ = 4, and σ = 5. The red dots are those below the detectable threshold in imaging, of 5 cells/mm 2 , and the blue dots those above. Note the different y-axis scales to highlight the resulting shape. (EPS) S5 Fig. Cell densities from one two-step fraction of exponential growth and death. Cell densities resulting from the optimal radiation profiles from two 2-step fraction of exponential growth and death for the initial density deviations of (from top left to bottom right) σ = 1, σ = 2, σ = 3, σ = 4, and σ = 5. The red dots are those below the detectable threshold in imaging, of 5 cells/mm 2 , and the blue dots those above. Note the different y-axis scales to highlight the resulting shape. (EPS) S6 Fig. Cell densities from two two-step fractions of exponential growth and death. Cell densities resulting from the optimal radiation profiles from two 2-step fraction of exponential growth and death for the initial density deviations of (from top left to bottom right) σ = 1, σ = 2, σ = 3, σ = 4, and σ = 5. The red dots are those below the detectable threshold in imaging, of 5 cells/mm 2 , and the blue dots those above. Note the different y-axis scales to highlight the resulting shape. (EPS) S7 Fig. Cell densities from one two-step fractions of logistic growth and exponential death. Cell densities resulting from the optimal radiation profiles from one 2-step fraction of logistic growth and exponential death for the initial density deviations of (from left to right) σ = 1, σ = 3, and σ = 5. The red dots are those below the detectable threshold in imaging, of 5 cells/mm 2 , and the blue dots those above. Note the different y-axis scales to highlight the resulting shape. (EPS)

S8 Fig. Cell densities from two two-step fractions of logistic growth and exponential death.
Cell densities resulting from the optimal radiation profiles from two 2-step fractions of logistic growth and exponential death of (from left to right) σ = 1, σ = 3, and σ = 5. The red dots are those below the detectable threshold in imaging, of 5 cells/mm 2 , and the blue dots those above. Note the different y-axis scales to highlight the resulting shape. (EPS) S9 Fig. Cell densities from one two-step fractions of logistic growth and death. Cell densities resulting from the optimal radiation profiles from one 2-step fraction of logistic growth and death of (from left to right) σ = 1, σ = 3, and σ = 5. The red dots are those below the detectable threshold in imaging, of 5 cells/mm 2 , and the blue dots those above. Note the different yaxis scales to highlight the resulting shape. (EPS)