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

Nassar H. S. Haidar

doi:10.3897/nucet.5.32239

Abstract

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 ε^*.

Keywords

Accelerator Based Modulated Neutron Sources, One-Speed Neutron Diffusion, Two Opposing Neutron Beams, Dynamical NCT, Pareto Optimization.

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, 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, 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, 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 and Vorozhtsov 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 and Calvo 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)

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.

The first nonlinear inequality constraint G₁(ω) = ε < 2π/ ω means that the time advance should not exceed the period of the modulated reference neutron beam, otherwise if ε = k (2π/T₀) , 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₀ is the life time of thermal neutrons in R, is the same as the system of three affine inequality constraints

G2(ω) = ϖ < 2π/ T₀.

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. (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 σ = (σ₁, σ₂, σ₃, σ₄) 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.

In addition to the notation of (Haidar 2018, Haidar 2017), namely: f (x,t) = vg (x,t), with v as the neutron speed and g (x,t) is the pertaining neutron density, åm = ℵ_R_;Λa_m, b~m = ℵ_R_;Πb_m, m = 0,1,2,3, ... with am = am (S (t),ω) and bm = bm (Ş (t),ϖ) are Fourier coefficients for the modulated two opposing neutron beams, with respective P and Ϸ periods and respective beam source periodic shapes S (t) and Ş (t), acting on an Λ ∪ R ∪ Π composite, as illustrated in Fig. 1, and ℵ_R_;Λ and ℵ_R_;Π are coupling factors (Haidar 2017, 2018) between the R region and the respective adjacent Λ and Π tumor-free regions,

we shall also need: (β)⁺ = β + δβ , where δβ is an infinitesimal positive increment of β.

Proposition 1

A Pareto optimal control vector for the therapeutic optimization problem (1–2) is: ω^* = (ω^* = ϖ^* = 2π/T₀, ε^* = (T₀/2)⁺).

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⁴) . A peak height that diminishes with increasing n. Hence most significant term in the double summation

is M^₁₁ (ϖ, ε).

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 Ņ₁₁(ϖ, ε). 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₁₁(ω).

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^₁₁ (ϖ, ε) 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 π/ϖ of F (ε). So, if ε_k is the k-th 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

, (7)

of (74), in (Haidar 2018), first as a function of ϖ.

The basic modal ϖ – factor in the numerator of (7) is

,

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^₁₁(ϖ, ε).

Consequently, as a sum (or difference, due to the sign pattern (–1)ⁿ^–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)ⁿ^–1) of such aligned concave functions, where A^₁₁(ϖ, ε) 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)ⁿ^–1) of such aligned concave functions, where B^₁₁(ϖ, ε) 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 ε → ± ∞.

Step 3. The previous analysis indicates that maximization of ῆ(ϖ) and ᾶ(ϖ) occurs asymptotically as ϖ → ∞, but by the constraint 0 << ω < ϖ ≤ 2π/T₀, ϖ should not exceed 2π/T₀, i.e. ϖ^* = 2π/T₀, 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 ε^* = ε₁ ≈ π/ϖ^* = T₀/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₀/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

. (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.

Acknowledgements

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.

References

Alvarez-Estrada RF, Calvo ML (2004) Neutron fibers and possible applications to NCT. Applied Radiation and Isotopes 61: 841–844. https://doi.org/10.1016/j.apradiso.2004.05.039

Haidar NHS (1985) Neutron fission rate measurements across a poisoning discontinuity in water. Journal of Nuclear Science and Technology 22(6): 484–490. https://doi.org/10.1080/18811248.1985.9735684

Haidar NHS (2002) Composite surface integral transforms in neutron therapy of cancer. International Journal of Theoretical Physics 41(1): 109–129. https://doi.org/10.1023/A:1013277701565

Haidar NHS (2017) On dynamical (B/Gd) neutron cancer therapy by accelerator based two opposing neutron beams. Journal of Nuclear Medicine and Radiation Therapy 8(6): 345. https://doi.org/10.4172/2155-9619.1000345

Haidar NHS (2018) On dynamical (B/Gd) neutron cancer therapy: accelerator-based single neutron beam. Pacific Journal of Applied Mathematics 9(1): 9–26.

Henry AF (1975) Nuclear-reactor analysis. The MIT Press, Cambridge, Massachusetts.

Hosmane N (2012) Boron and Gadolinium neutron capture for cancer treatment. World Scientific, Singapore. https://doi.org/10.1142/8056

Miettinen K (1998) Nonlinear multiobjective optimization. Springer, Berlin. https://doi.org/10.1007/978-1-4615-5563-6

M’silti A, Tolla P (1993) An interactive multiobjective nonlinear programing procedure. European Journal of Operational Research 64: 1115–125. https://doi.org/10.1016/0377-2217(93)90012-C

Smirnov V, Vorozhtsov S (2016) Modern compact accelerators of cyclotron type for medical applications. Physics of Particles and Nuclei 47(5): 863–883. https://doi.org/10.1134/S1063779616050051

Sundaram RK (1996) A first course in optimization theory. Cambridge University Press, Cambridge. https://doi.org/10.1017/CBO9780511804526

Yoshida Y, Nakayama H, Nakai K, Yamazaki T (1976) Beam-pulsing system for the IMS cyclotron. Nuclear Instruments and Methods in Physics Research A 138(4): 579–588. https://doi.org/10.1016/0029-554X(76)90003-3