Corresponding author: Nassar H. S. Haidar ( nhaidar@suffolk.edu ) Academic editor: Yury Korovin
© 2019 Nassar H. S. Haidar.
This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Citation:
Haidar NHS (2019) Optimization of two opposing neutron beams parameters in dynamical (B/Gd) neutron cancer therapy. Nuclear Energy and Technology 5(1): 17. https://doi.org/10.3897/nucet.5.32239

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 ε^{*}.
Accelerator Based Modulated Neutron Sources, OneSpeed Neutron Diffusion, Two Opposing Neutron Beams, Dynamical NCT, Pareto Optimization.
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/Gdloaded, see e.g. (
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 (
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. (
In the therapeutic setup for dynamical (B/Gd) NCT by one, (
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 multiobjective (Pareto) optimization (
G2(ω) = ϖ < 2π/ T_{0}.
G3(ω) = ω – ϖ ≤ 0
G4(ω) = – ω < 0. (2)
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. (
To have a look at the posing problem from a different angle, we note that the results of (
with 6 control variables. Here σ = (σ_{1}, σ_{2}, σ_{3}, σ_{4}) is a vector of slack variables, see e.g. (
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.
In addition to the notation of (
we shall also need: (β)^{+} = β + δβ , where δβ is an infinitesimal positive increment of β.
A Pareto optimal control vector for the therapeutic optimization problem (1–2) is: ω^{*} = (ω^{*} = ϖ^{*} = 2π/T_{0}, ε^{*} = (T_{0}/2)^{+}).
Proof. As usual in Pareto optimization, see e.g. (
Step 1. Let us analyze qualitatively, first as a function of ϖ, the therapeutic utility index
in which,
and all the symbols have their usual meaning, as in (70) of (
The basic modal ϖ – factor in the numerator of (4) is
in which
is a bellshaped 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 bellshaped, 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 convextoconcave 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 bellshaped, 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 convextoconcave 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 halfrange 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 concavetoperiodic ratio of functions should be concave with a global maximum attained at ε → ∞, but periodically discontinuous, each halfperiod π/ϖ of F (ε). So, if ε_{k} is the kth discontinuity of ῆ(ω, ϖ, ε), then
ε_{k} = k (π/ϖ), k = 1, 2, 3, ...
with
ῆ(ω, ϖ, ε_{k}) < ῆ(ω, ϖ, ε_{k+1}), ∀k (6)
Step 2. Then we apply the same analysis for the ballistic index
of (74), in (
The basic modal ϖ – factor in the numerator of (7) is
in which
is a bellshaped 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 convextoconvex ratio of concentric functions should be a rather flattened convex (or quasiconvex) 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 quasiconvex) 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,
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 quasiconcave) function with a minimum near ε = 0 and quickly attained (maximum) at ε > 0 and towards ε → ± ∞.
Step 3. The previous analysis indicates that maximization of ῆ(ϖ) and ᾶ(ϖ) occurs asymptotically as ϖ → ∞, but by the constraint 0 << ω < ϖ ≤ 2π/T_{0}, ϖ should not exceed 2π/T_{0}, i.e. ϖ^{*} = 2π/T_{0}, Pareto optimal. Also since both ῆ(ω) and ᾶ(ω) maximize asymptotically as ω → ∞ and ω << ϖ, then ω^{*} = ϖ^{*}, Pareto optimal.
The situation with ε is remarkably quite different. While ᾶ(ε) attains a stationarizing maximum quickly with increasing ε > 0, ῆ(ε) is concave with a global maximum attained at ε → ∞, but is periodically discontinuous at ε_{k} = k (π/ϖ), k = 1,2,3, ... . Moreover since ε < (2π/ω), then ε^{*} = ε_{1} ≈ π/ϖ^{*} = T_{0}/2. However, by structure of F (ε), the limit to the right of the ε_{k} discontinuity should correspond to a local maximum of ῆ(ε). For that reason ε^{*} = (T_{0}/2)^{+}, Pareto optimal. Here the proof ends.
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
As a result, the entire analysis, reported in this note, holds true if we replace each å_{m} and b~_{m} respectively with
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.
The author is grateful to two anonymous referees for their critical reading of an original version of this article and for a number of useful comments.