Optimization of two opposing neutron beams parameters in dynamical ( B / Gd ) neutron cancer therapy

We demonstrate how the therapeutic utility index and the ballistic index for dynamical neutron cancer therapy (NCT) with two opposing neutron beams form a nonlinear optimization problem. In this problem, the modulation frequencies ω and π of the beams and the relative time advance ε are the control variables. A Pareto optimal control vector ω* = (ω*, π*, ε*) for this problem is identified and reported for the first time. The utility index is shown to be remarkably periodically discontinuous in ε, even in the neighborhood of ε*.


Introduction
In both dynamical and/or stationary (B/Gd) neutron cancer therapy (NCT), there is always an unavoidable undesirable effect of irradiating some surrounding healthy tissues together with a targeted, B/Gd-loaded, see e.g.(Hosmane 2012), cancerous tumor, R. The ballistic index of neutrons, to be maximized in a NCT (in order to minimize this undesirable effect) is the ratio of the (B/ Gd) reaction rate inside R, to the total neutronic reaction rate inside the patient.This index was advanced by this author in (Haidar 2002(Haidar , 2018)), then tailored to the therapeutic setup of two opposing dynamical one-speed neutron beams in (Haidar 2017).Furthermore, in dynamical NCT, the neutron density distribution inside R is decomposable, (Haidar 2017(Haidar , 2018)), into three components: dissipative (associated with the useful dynamical reaction rate), periodic and stationary.Here the ratio of the dissipative-to-periodic neutron density wave, as a measure of the utility of time of the used modulated accelerators, was also developed by this author in (Haidar 2017(Haidar , 2018)), as another therapeutic index.
A standard technique conventionally used for boosting the ballistic index in stationary NCT is to irradiate the tumor several times with beams coming from different directions via rotating gantries (Smirnov andVorozhtsov 2016, Yoshida et al. 1976), or via rotating the patient himself.Nowadays, we already have operating compact superconducting medical cyclotrons that may rotate around the patient like e.g. the one in the Harper Hospital in Detroit USA.The entire structure of this cyclotron weighs less than 25 tons and this allowed the cyclotron to be mounted directly on a rotating gantry so that neutrons produced by any internal charged particle target can be directed at a supine patient from any desired direction.Medical accelerators are steadily becoming more compact and more flexibly temporally controllable (Yoshida et al. 1976), and it may not be unrealistic to conceive an NCT facility housing two accelerators, or two time-modulated neutron beams, that are directed by gantries towards a cancer patient.Alternatively, external-beam pulsing systems (Haidar 2002, Alvarez-Estrada andCalvo 2004), may also be employed for the same purpose.
It should be underlined, from the outset of this paper, that the diffusion model for neutron transport near boundaries of strongly absorbing R regions can only serve as a rudimentary approximation, see e.g.(Henry 1975), and a one-speed diffusion model makes the situation worse.It is also known that the multigroup transport operator is a much more difficult operator to handle analytically than the one-speed diffusion operator, particularly with time evolution applications, as in the present work.
In the therapeutic setup for dynamical (B/Gd) NCT by one, (Haidar 2018), or two, (Haidar 2017), opposing neutron beams, the state diffusion equation is linear in the neutron flux solution f(x,t) of the parabolic boundary value problem (BVP), (Haidar 2017), of one-speed neutron diffusion.The above pertaining two objective functions : therapeutic utility index ῆ(ω, ϖ, ε) and ballistic index ᾶ(ω, ϖ, ε) for a cancerous slab R of parameters {Σ a , D, l} turn out to be nonlinear functions of the components of the control vector ω = (ω, ϖ, ε).In the set of parameters, Σ a is the macroscopic one-speed neutron absorption cross section, D is the diffusion length for these neutrons, see e.g.(Henry 1975, Haidar 1985), a is the thickness of the cancerous slab and l is its extrapolated end, as illustrated in Figure 1.The closed analytical solvability of the pertaining BVP happens to lead to differentiable objective functions within the following nonlinear programming problem: . (1) The first nonlinear inequality constraint G 1 (ω) = ε < 2π/ ω means that the time advance should not exceed the period of the modulated reference neutron beam, otherwise if ε = k(2π/T 0 ) , k = 1, 2, 3, ... it becomes redundant.In this multi-objective (Pareto) optimization (M'silti and Tolla 1993) problem, the affine two-sided second inequality constraint, in which T 0 is the life time of thermal neutrons in R, is the same as the system of three affine inequality constraints It is well known that if it happens that the objective functions are conflicting, i.e. if ῆ(ω) cannot be improved, by varying ω, without degrading the value of ᾶ(ω), then the problem becomes entirely unsolvable.Otherwise, if ῆ(ω) and ᾶ(ω) are not conflicting, then the problem may accept Pareto optimal, see e.g.(M'silti and Tolla 1993, Miettinen 1998), solutions pending to satisfaction of the inequality constrains (2).Demonstrating that this what actually happens, is the purpose of this note.

Pareto optimization
To have a look at the posing problem from a different angle, we note that the results of (Haidar 2018) indicate that problem (1) can always be replaced (or cast into) the following analytically closed equivalent system of 10 nonlinear and linear equations: (3) with 6 control variables.Here σ = (σ 1 , σ 2 , σ 3 , σ 4 ) is a vector of slack variables, see e.g.(Sundaram 1996), necessary to transform the inequality G constraints to equality Q constraints.
Problem (3) is clearly overdefined and may, in general be unsolvable, unless when at least four of the above 10 equations happen to be inactive ( or insignificant).Luckily this turns out to take place for the posing therapeutic setup, as to be demonstrated in the result that follows.

Proposition 1.
A Pareto optimal control vector for the therapeutic optimization problem (1-2) is: Proof.As usual in Pareto optimization, see e.g.(Sundaram 1996, Miettinen 1998), the arguments we shall apply to the analysis of (1), are basically of graphical qualitative nature.
Step 1.Let us analyze qualitatively, first as a function of ϖ, the therapeutic utility index , (4) in which, , and all the symbols have their usual meaning, as in ( 70) of (Haidar 2018).
The basic modal ϖ -factor in the numerator of (4) is , in which is a bell-shaped convex function of ϖ, which is peaked at ϖ = 0 with a negative peak height of -(1/β n ) and peak width ~ (1/m 4 ) .A peak height that diminishes with increasing n.Hence most significant term in the double summation is M ^11 (ϖ, ε).
Consequently, as a sum (or difference, due to the sign pattern of cos(2n-1)π/4) of equicentral bell shapes, should be convex with a minimum at ϖ = 0 and stationarizing (with maxima) at ϖ  ± ∞.The numerator of (4), with conceived as a constant independent of ϖ, retains the convexity property of ϖ, retains the convexity property of As for the basic modal ϖ -factor in the denominator of (4), which is the term is another bell-shaped, flatter than , but concave function of ϖ, which is peaked at ϖ = 0 with a positive peak height of 1/β n .Then taking into consideration that i.e. sin[Ω mn (ϖ, ε)] ≈ 1, it can be argued that Ņ mn (ϖ, ε) remains concave ∀ϖ ∈ (-∞; ∞).Hence the most significant term in the double summation is Ņ 11 (ϖ, ε).Consequently, as a sum (or difference, due to cos(2n-1)π/4) of equicentral bell shapes, should be concave with a maximum at ϖ = 0 and stationarizing (with minima) at ϖ  ± ∞.The denominator of (4), with conceived as a constant independent of ϖ, retains the concavity property of Hence ῆ(ω, ϖ, ε), as a function of ϖ, being a convex-to-concave ratio of concentric functions should be convex with a minimum at ϖ = 0 and stationarizing (with maxima) at ϖ  ± ∞.
Second we consider (4), as a function of ω.Repeating essentially the same analysis on the basic modal ω -factor in the numerator of (4) : , leads to a conclusion that should be convex with a minimum at ω = 0 and stationarizing (maxima) at ω  ± ∞.The numerator of (4), with conceived as a constant independent of ω retains the convexity property of .
As for the basic modal ω -factor in the denominator of (4), which is the term is another bell-shaped, flatter than , but concave function of ω which is peaked at ω = 0 with a positive peak height of 1/β n .Then taking into consideration that , indicates that N mn (ω) remains concave ∀ ω ∈ [0, ∞].Also here the most significant term in the double summationis is N 11 (ω).
Then we observe that the denominator of (4), with conceived as a constant independent of ω, retains the concavity property of .Hence ῆ(ω, ϖ, ε), as a function of ω, being a convex-to-concave ratio of concentric functions should be convex with a minimum at ω = 0 and stationarizing (with maxima) as ω  ± ∞.
Third, we consider (4) as a function of ε.The basic modal ε -factor in the numerator of ( 4) is M ^mn (ϖ, ε) in which -e -βnε is concave with a maximum attained at ε  ∞.As a sum (or difference, due to the sign pattern of cos(2n-1) π/4) of such aligned concave functions, where M ^11 (ϖ, ε) as the most significant term, should be concave with maximum attained at ε  ∞.The numerator of (4), with conceived as a constant independent of ε, retains the concavity property of .Quite distinctively, the basic modal ε -factor in the denominator of (4), which is Ņ mn (ϖ, ε) happens to contain the term , ϖ  ± ∞, (5) which is almost periodic in ε.The double summation represents a half-range expansion of a certain periodic function F(ε) in sinusoidal Fourier series.The period of this F(ε) is 2π/ϖ.
The denominator of ( 4) is simply this periodic function added to a constant ( independent of ε) term: . Now ῆ(ω, ϖ, ε), as a function of ε, being a concave-to-periodic ratio of functions should be concave with a global maximum attained at ε  ∞, but periodically discontinuous, each half-period Step 2. Then we apply the same analysis for the ballistic index , ( 7) of (74), in (Haidar 2018), first as a function of ϖ.
The basic modal ϖ -factor in the numerator of ( 7) is , with in which is a bell-shaped convex function of ϖ which is peaked at ϖ = 0 with a negative peak height of -(1/β n ).A peak height that diminishes with increasing n, and the most significant term in the double summation being A ^11 (ϖ, ε).
Consequently, as a sum (or difference, due to the sign pattern (-1) n-1 ) of equicentral bell shapes, should be convex with a minimum at ϖ = 0 and stationarizing (maxima) at ϖ  ± ∞.The numerator of (4), with conceived as a constant independent of ϖ, retains the convexity property of .
Also the basic modal ϖ -factor in the denominator of (7), has the same convex behavior as A ^mn (ϖ, ε), while distinctively increasing with increasing n.Also convex are and , with a minimum at ϖ = 0 and stationarizing (maxima) at ϖ  ± ∞.
Then the denominator of (7), with conceived as a constant independent of ϖ, retains the convexity property of .
Hence ᾶ(ω, ϖ, ε), as a function of ϖ, being a convex-to-convex ratio of concentric functions should be a rather flattened convex (or quasi-convex) function with a minimum at ϖ = 0 and stationarizing (with maxima) at ϖ  ± ∞.Second, we consider (7), as a function of ω.Repeating essentially the same arguments to the basic modal ωfactors in the numerator and denominator of (7), which are respectively , and , leads to the conclusion that ᾶ(ω, ϖ, ε), as a function of ω,being a convex to convex ratio of concentric functions, should also be a rather flattened convex (or quasi-convex) function with a minimum at ω = 0 and stationarizing (with maxima) at ω  ± ∞.
Third, we consider (7) for ᾶ as a function of ε.The basic modal ε -factor in the numerator of ( 7) is A ^mn (ϖ, ε), in which -e -βnε is concave with a maximum attained at ε  ∞.As a sum (or difference, due to the sign pattern (-1) n- 1) of such aligned concave functions, where A ^11 (ϖ, ε) is the most significant term, should be concave with maximum attained at ε  ∞.The numerator of ( 7), with conceived as a constant independent of ε, retains the concavity property of .Quite distinctively, the basic modal ε -factor in the denominator of (7), which is B ^mn (ϖ, ε) contains the term -e -βnε which is concave with a maximum attained at ε  ∞.As a sum (or difference, due to (-1) n-1 ) of such aligned concave functions, where B ^11 (ϖ, ε) as the most significant term, and should be concave with maximum attained at ε  ∞.The denominator of (7), with conceived as a constant independent of ε, retains the concavity property of .Obviously ᾶ(ω, ϖ, ε), as a function of ε, being a concave to concave ratio of aligned functions, should also be a rather flattened concave (or quasi-concave) function with a minimum near ε = 0 and quickly attained (maximum) at ε > 0 and towards ε  ± ∞.
From the previous proof, it is clear how the sensitivities of ῆ and ᾶ to variations either in ϖ and/or ω are quite similar.The sensitivities of ῆ(ε) and ᾶ(ε) to variations in ε are, however, remarkably very different.

Remark 1.
In some cancer patients, transport of thermal neutrons by neutron guides or neutron optical fibers through the regions Λ and Π (with respective thicknesses l Λ and l Π and neutron macroscopic removal cross sections Σ Λ and Σ Π ), may turn out to be medically unfeasible.As an approximate substitute to solving the composite Λ ∪ R ∪ Π regional neutronics problem in this case, one can simply assume a planar attenuation, towards R, of the sources S(x,t) and Ş(x,t) respectively to and . (8) As a result, the entire analysis, reported in this note, holds true if we replace each å m and b ~m respectively with . (9)

Conclusion
In this note, we have demonstrated that the therapeutic indices in NCT, with two opposing neutron beams, are controllable in a Pareto optimization process.The existence of a Pareto optimal control vector is proved, for any given pair {S(t), Ş(t)} of beam source shape functions.
These facts pave the way towards launching the building of the first experimental benchmark for this new kind of NCT, on which, e.g. the impact of changes in the set {S(t), Ş(t)} on ω * , can further be investigated.

Figure 1 .
Figure 1.Sketch to illustrate the two opposing neutron beams in a dynamical (B/Gd) NCT setup; with x α and z α as the entrance points for fast neutrons to moderators and x β and z β are exit points for thermal neutrons from moderators.