Physica Medica xxx (2016) xxx–xxx

Contents lists available at ScienceDirect

Physica Medica
journal homepage: http://www.physicamedica.com

Original paper

Development of Monte Carlo based real-time treatment planning systemwith fast calculation algorithm for boron neutron capture therapy
Kenta Takada a,⇑, Hiroaki Kumada a, Peng Hong Liem b, Hideyuki Sakurai a, Takeji Sakae aabFaculty of Medicine, University of Tsukuba, 1-1-1, Tennodai, Tsukuba, Ibaraki 305-8575, Japan
NAIS Co., Inc., 416, Muramatsu, Tokai-mura, Naka-gun, Ibaraki 319-1112, Japana r t i c l e

i n f oArticle history:
Received 20 July 2016
Received in Revised form 3 November 2016
Accepted 4 November 2016
Available online xxxx
Keywords:
BNCT
PHITS
Fast calculation
Real-time dose evaluationa b s t r a c t
Purpose: We simulated the effect of patient displacement on organ doses in boron neutron capture
therapy (BNCT). In addition, we developed a faster calculation algorithm (NCT high-speed) to simulate
irradiation more efficiently.
Methods: We simulated dose evaluation for the standard irradiation position (reference position) using ahead phantom. Cases were assumed where the patient body is shifted in lateral directions compared to
the reference position, as well as in the direction away from the irradiation aperture.
For three groups of neutron (thermal, epithermal, and fast), flux distribution using NCT high-speed with avoxelized homogeneous phantom was calculated. The three groups of neutron fluxes were calculated forthe same conditions with Monte Carlo code. These calculated results were compared.
Results: In the evaluations of body movements, there were no significant differences even with shifting
up to 9 mm in the lateral directions. However, the dose decreased by about 10% with shifts of 9 mm in a
direction away from the irradiation aperture.
When comparing both calculations in the phantom surface up to 3 cm, the maximum differences
between the fluxes calculated by NCT high-speed with those calculated by Monte Carlo code for thermal
neutrons and epithermal neutrons were 10% and 18%, respectively. The time required for NCT high-speed
code was about 1/10th compared to Monte Carlo calculation.
Conclusions: In the evaluation, the longitudinal displacement has a considerable effect on the organ
doses.
Conclusions: We also achieved faster calculation of depth distribution of thermal neutron flux using
NCT high-speed calculation code.
Ó 2016 Associazione Italiana di Fisica Medica. Published by Elsevier Ltd. All rights reserved.

1. Introduction
Boron neutron capture therapy (BNCT) has attracted attention
as a cancer therapy method that can be implemented at the cellular level. To perform BNCT, a treatment planning system is used forirradiation simulations in a manner similar to other radiotherapies
such as X-ray therapy and particle therapy. A team including
researchers from the University of Tsukuba is now developing a
new treatment planning system for BNCT (with the tentative nameof Tsukuba Plan) using a Monte Carlo algorithm [1,2].
BNCT has characteristics that differ from X-ray or particle beam
therapies, for example, the fact that a rotational gantry is notmounted. For this reason, since the irradiation position is aligned⇑ Corresponding author.
E-mail address: k-takada@md.tsukuba.ac.jp (K. Takada).

to the fixed irradiation aperture, there may be cases where
maintaining a fixed patient body position will be difficult. In addition, because BNCT irradiation times are relatively long comparedwith those of other radiation therapy modalities, the extended
immobilization can be uncomfortable for patients. Thus, patient displacement during irradiation can occur, and affect the organ dose(e.g., intra-fractional variation). In regards to the patient displacement during BNCT treatment, Wielopolski et al. evaluated the effectof the patient body shifting (0.5, 1.0, 2.0 cm) in the case of BNCT
treatment for brain tumor [3]. An approach to raise the immobilization precision of the patients in BNCT has also been studied [4].
This research aims to establish more efficient BNCT treatment
with a focus on the changes in the dose for a body displacementof less than 1 cm. Depending on the irradiation site and immobilization accuracy, body displacements less than 1 cm are contemplated as intra-fractional variation for radiation therapy [5,6]. If

http://dx.doi.org/10.1016/j.ejmp.2016.11.007
1120-1797/Ó 2016 Associazione Italiana di Fisica Medica. Published by Elsevier Ltd. All rights reserved.

Please cite this article in press as: Takada K et al. Development of Monte Carlo based real-time treatment planning system with fast calculation algorithm
for boron neutron capture therapy. Phys. Med. (2016), http://dx.doi.org/10.1016/j.ejmp.2016.11.007

2

K. Takada et al. / Physica Medica xxx (2016) xxx–xxx

the dose to organs at risk (OARs) most often or target dose changes
owing to intra-fractional variation in the treatment, a framework
for evaluating those changes becomes necessary. In particular,
real-time evaluation for OARs is necessary to improve BNCT treatment because the tolerance dose for OARs determines the irradiation time.
However, real-time evaluation is very difficult using conventional Monte Carlo calculations because of their time-consuming
nature. Therefore, we also investigated the higher-speed calculation algorithm by using neutron diffusion theory.
2. Materials and methods
2.1. Evaluation of organ dose uncertainties owing to patient body
movement
This simulation was carried out using a human body head phantom (PB-1, Kyoto Kagaku Co., Ltd) and the Tsukuba Plan [1,7]. The
Tsukuba Plan can simulate dose distribution using alternative radiation source information to enable irradiation simulations of various BNCT irradiation fields. In this evaluation, we assumed that thetherapy was carried out using the Japan Research Reactor No. 4(JRR-4) epithermal neutron beam [8].
For cases where the patient’s body could easily move, we performed an irradiation simulation for the left maxillary sinus withpatients in the seated position. The target was assumed to be atumor dominating the left maxillary sinus. For a situation where
the irradiation method involved the phantom forehead touchinga fixed aperture, the target center was set so that the conditions
for the beam center to pass through are those of the standard position (reference position). We evaluated the maximum dose to theskin, the limiting factor of the irradiation time, and the minimum
dose to the tumor from the dose volume histogram for each body
position. Table 1 shows the parameters required for BNCT
Table 1
Various parameters used in the dose evaluation.

⁄⁄⁄

Parameters

Tumor

Boron concentration (ppm)
CBE⁄
RBE⁄⁄

72
24
3.8
1.35
Hydrogen dose: 2.5
Nitrogen dose: 2.5
Gamma dose: 1.0

Normal tissue

Skin
28.8
2.5

treatment such as boron concentration, compound biological effectiveness, and relative biological effectiveness from published literature regarding the JRR-4 beam [9]. In the reference position, theirradiation time until the skin tolerance dose becomes 15 Gyweighted (Gyw) was determined. As shown in Fig. 1, patient body
motion away from center of the beam was simulated at steps of
1.5 mm up to a distance of 9.0 mm to the left and right of the reference position (lateral displacements). Furthermore, body
motions were also evaluated up to a distance of 9.0 mm from theaperture, with steps of 1 mm (longitudinal displacements). For aposition shifted from the reference position, a simulation was performed wherein irradiation was performed for the same time as
the irradiation time used for the reference position.
The particle and heavy ion transport code system (PHITS) was
selected as the Monte Carlo calculation code [10], and this code
was executed using a high performance, parallel computer cluster.
The dose distribution with three hundred million particle histories
was calculated using the PHITS code for each position. For each
body position, the calculations were performed five times with different initial random number seeds.
2.2. Development of a higher speed BNCT calculation algorithm, andneutron flux distribution calculation
Dose evaluation in real-time using Monte Carlo calculations ofchanges to dose distributions owing to patient body movement
during BNCT irradiation is extremely difficult because it is verytime consuming. Therefore, we developed a higher-speed BNCT
calculation algorithm (NCT high-speed). The basic concept of thiscode is based on neutron diffusion theory [11]. In the algorithm,
we divided the continuous neutron energies into multi-groups
and calculated the spatial distribution of the neutrons by solving
the fixed-source neutron diffusion equations. In this research, we
adapted the eight energy groups used for neutron beam calculations. To develop the code, and to enable an adaptation of such calculations to BNCT dose calculations for a human body, readings ofthe human body voxel model were created by averaging pixel dataof medical images.
The fixed source multi-group neutron diffusion problem can be
written as r  Dg ðrÞr/g ðrÞ þ Rt;g ðrÞ/g ðrÞ ¼

G
X

Rs;g0 !g ðrÞ/g0 ðrÞ þ Q g ðrÞ

g 0 ¼1

CBE: compound biological effectiveness.
RBE: relative biological effectiveness.

g ¼ 1; 2; . . . ; G:

ð1ÞFig. 1. Simulation geometry for the patient body motion. Dd: the distance from the final collimator to the surface of the head phantom.

Please cite this article in press as: Takada K et al. Development of Monte Carlo based real-time treatment planning system with fast calculation algorithm
for boron neutron capture therapy. Phys. Med. (2016), http://dx.doi.org/10.1016/j.ejmp.2016.11.007

K. Takada et al. / Physica Medica xxx (2016) xxx–xxx

where £ðrÞ and Q ðrÞ are the space (r) dependent group neutron flux
and the external group neutron source, respectively, andRt ðrÞ; Rs ðrÞ; and DðrÞ denote the space dependent total, scatteringcross sections, and the spatial diffusion constant, respectively.
These group constants are prepared so that Eq. (1) produces solutions close to the ones of the higher order neutron transport
method. The left and right hand sides of Eq. (1) represent neutron
loss and production terms, respectively. The first and second termsof the left hand side of Eq. (1) denote neutron loss from leakage andtotal removal from energy group g, respectively, while the first termat the right hand side denotes neutron production from scattering
into energy group g.
The external group space-dependent neutron source is represented by Qg(r).
In the NCT high-speed method, we evaluated the depth distributions of group neutron flux within a 20  20  20 cm cubic(8000 cm3) homogenous phantom (0.5 cm cubic voxel model) consisting of soft tissue defined by the International Commission on
Radiation Units and Measurements (density: 1.03 g/cm3). For theneutron beam used in the evaluation, we used the abovementioned JRR-4 epithermal neutron beam. The evaluation geometry is shown in Fig. 2. The neutron beam used in BNCT is a wide
range continuous energy spectrum [12]. In the NCT high-speed
method, we handled and calculated the energy structure as eightgroup structures, as shown in Table 2. We performed evaluationsof depth distributions within the phantom, and compared themwith the values calculated using PHITS, with regards to the neutron
fluxes of three energy groups: fast neutrons, epithermal neutrons,
and thermal neutrons. For the PHITS calculations, we implemented

3

one hundred and fifty million particle histories for the same voxel
model using a parallel computing system of 160 cores. This number of histories was determined by conditions that adequately
reduce the statistical uncertainties (<3%) of the Monte Carlo calculation results up to a depth of 15 cm within the phantom. The calculation mesh size was set as 0.5 cm3.
3. Results
3.1. Calculation results of patient body movement in BNCTThe dose evaluations using PHITS related to body movement in
the left and right directions are shown in Fig. 3. The evaluations
related to body movement in a direction away from the irradiation
aperture are shown in Fig. 4. When compared with the reference
position dose, and considering the error owing to body movementin the lateral direction when shifting up to 9 mm from the reference position, the maximum skin dose and the minimum target
dose were not significantly different. In contrast, the results related
to body movement in a direction away from the aperture, whenseparated by 9 mm from the reference position, gave the maximum skin dose and the minimum target dose of 10.4% and 10.1%
below the reference position dose, respectively.
3.2. Calculation result by NCT high-speed calculation algorithm for
BNCTRegarding the depth distribution of the three groups of neutron
flux within a homogeneous phantom, the calculation result of theNCT high-speed method compared with the PHITS method isshown in Fig. 5. Both results were normalized by the values atthe phantom surface.
The neutron fluxes for the three groups calculated using theNCT high-speed method showed the same trend as the depth distribution calculated with the PHITS method. However, large differences, in particular of the fast neutron flux, were seen in regions
deeper than 3 cm in the phantom. To evaluate the dose for OARs,
the thermal neutron and epithermal neutron flux are important.
In shallow regions up to 3 cm from the phantom surface, the maximum differences between the fluxes calculated by the NCT highspeed method with those calculated by PHITS for thermal neutron
and epithermal neutrons were 10% and 18%, respectively. In addition, the NCT high-speed method required only 2.3 min of computing time, while the PHITS method required 26.9 min.
4. Discussion
4.1. The impact of body movement on BNCT

Fig. 2. Simulation geometry for comparison of the non-Monte Carlo calculation
algorithm and the Monte Carlo calculations. The phantom consisted of homogeneous soft tissue.

Table 2
Neutron energy structure used in the high-speed calculation algorithm.
Group

Energy range (eV)

1
2
3
4
5
6
7
8

8.21E+05
6.74E+04
1.23E+03
1.76E+01
6.83E01
2.77E01
8.53E02
1.00E05

1.00E+07
8.21E+05
6.74E+04
1.23E+03
1.76E+01
6.83E01
2.77E01
8.53E02

Based on the evaluation results performed for BNCT in theseated position, we found that there was little effect on organ
doses even with a position shift of several millimeters in either
the left and right directions. However, there was a considerable
effect on the doses if the position shifts in the direction moving
away from the irradiation aperture. This is thought to be causedby beam divergence and the distance from the source to patient
position.
While BNCT irradiation facilities can be divided into reactorbased and accelerator-based facilities, in either case various materials (moderators, shield) are distributed in the upper stream of theirradiation aperture [13,14]. There is a mutual interaction (mostly
neutron scattering) between these materials and the neutrons,
producing a non-zero angular distribution of components in theneutron beam; thus, the beam straightness is slightly distorted,
especially in the case of accelerator-based BNCT.

Please cite this article in press as: Takada K et al. Development of Monte Carlo based real-time treatment planning system with fast calculation algorithm
for boron neutron capture therapy. Phys. Med. (2016), http://dx.doi.org/10.1016/j.ejmp.2016.11.007

4

K. Takada et al. / Physica Medica xxx (2016) xxx–xxx

Fig. 3. Calculation results of dose variation for skin and target due to body motion in the lateral direction. Error bars show the statistical uncertainties by calculations with 5
times of Monte Carlo calculations. Bar graphs: weighted dose. Symbols: dose deviations in the case of lateral displacements compared to reference position.

Fig. 4. Calculation results of dose variation for skin and target due to body motion of the direction away from the irradiation aperture. Error bars show the statisticaluncertainties by calculations with 5 times of Monte Carlo calculations. Bar graphs: weighted dose. Symbols: dose deviations in the case of longitudinal displacements
compared to reference position.

In our evaluation for longitudinal displacements, the relationship of Dd with the dose deviation for organ doses did not obey
the inverse square law. The inverse square law is obeyed the total
distance from the source to the patient with patient dose. The distance from the source to the final collimator depends on the BNCT
facilities. Thus, the relationship between the patient longitudinal
displacements with the dose also depends on the facilities.
4.2. The fast calculation algorithm for BNCTIn this evaluation, we adopted a calculation algorithm based on
the multi-group neutron diffusion theory: the NCT high-speed
method for BNCT evaluations. The BNCT dose calculation is commonly conducted by using a stochastic (Monte Carlo) method;

however, the NCT high-speed method is a deterministic one, andis not based on the Monte Carlo method. When adapting to realworld BNCT treatment planning, there is a need for the NCT
high-speed method to perform particle transport calculations
using a complex, realistic human head model. Fig. 6 shows the spatial distribution of the thermal neutron flux calculated with theNCT high-speed method with a realistic voxelized head phantom.
This figure shows that the NCT high-speed method is able to calculate the three-dimensional spatial distribution of thermal neutron
flux, even in complex shapes constructed from medical imaging
data.
However, the calculation accuracy of the NCT high-speed
method is inferior to the Monte Carlo method. Particularly, large
differences were found in the deeper region of the phantom. In

Please cite this article in press as: Takada K et al. Development of Monte Carlo based real-time treatment planning system with fast calculation algorithm
for boron neutron capture therapy. Phys. Med. (2016), http://dx.doi.org/10.1016/j.ejmp.2016.11.007

K. Takada et al. / Physica Medica xxx (2016) xxx–xxx

5Fig. 5. Comparison of the results of the calculated depth distributions for the three groups of neutron flux in a homogeneous phantom. Lines indicate the results of the NCThigh speed method (non-Monte Carlo calculation), and symbols give the results of the Monte Carlo calculations.

Fig. 6. Calculation results of the thermal neutron distributions using the NCT high-speed calculation algorithm in a complex human-shaped phantom model.

BNCT, the irradiation time is almost determined by the dose limitations of OARs such as skin and mucosa. If a patient position shift
has occurred, the dose calculation should be performed again tothese non-deep-sited organs with faster calculation. Therefore, itis thought that the NCT high-speed calculation can adapt to theshallow (e.g. superficial) regions under the present conditions.
In this research, calculation results for neutron fluxes were only
represented, of course, the absorbed dose should be calculatedusing an energy-dependent coefficient factor from neutron flux.
4.3. Challenges to developing a real-time dose evaluation system for
BNCT
Our final goal is to develop BNCT dose evaluations in real-time,
however, a number of issues remain.
For example, the depth distributions of the thermal neutron flux
within the homogeneous phantom, calculated with the NCT highspeed method, are not identical to the corresponding distributions
calculated by the PHITS method. Thus, there is a need to enhance

the calculation precision of the NCT high-speed method. One
method to do this is by using more neutron energy groups thanthe current eight group structure. Because the calculation time will
also increase when the energy structure is more detailed, there is aneed to find a good balance between the calculation precision andthe calculation time.
In addition, there is a need for a method of monitoring patientbody movements during treatment, as well as a method thatavoids patient body movement as much as possible. When such
movements are unavoidable, the direction and the degree of movement should be measured and fed back into the system for dose reevaluations. With the collection of this type of data, we expect thata real-time BNCT dose evaluation system can be realized.

5. Conclusion
In this research, we used the Monte Carlo calculation code
PHITS to evaluate the effects of patient position shift on the skin

Please cite this article in press as: Takada K et al. Development of Monte Carlo based real-time treatment planning system with fast calculation algorithm
for boron neutron capture therapy. Phys. Med. (2016), http://dx.doi.org/10.1016/j.ejmp.2016.11.007

6

K. Takada et al. / Physica Medica xxx (2016) xxx–xxx

and target doses during BNCT treatment. The results showed thatthe effect of lateral position shift is small, but that the organ doseschanges considerably with distance changes in the direction awayfrom the irradiation aperture.
In addition, we introduced the principles and the present development status concerning a high-speed calculation algorithm, theNCT high-speed method, which can achieve relatively high-speed
BNCT calculations as compared to Monte Carlo calculations. In
the future, we plan to implement further refinements to enhance
the possibility of realizing real-time BNCT dose evaluations.
Conflict of interest
The authors report no conflicts of interest.
Acknowledgements
Numerical calculations were performed in part using the Cluster of Many-core Architecture processor (COMA) at the Center for
Computational Sciences, University of Tsukuba.
This work was supported by JSPS KAKENHI Grant Numbers
JP26860394, JP16H05236.
References[1] Kumada H, Matsumura A, Sakurai H, Sakae T, Yoshioka M, Kobayashi H, et al.
Project for the development of the linac based NCT facility in University ofTsukuba. Appl Radiat Isot 2014;88:211–5.
[2] Kumada H, Takada K, Yamanashi K, Sakae T, Matsumura A, Sakurai H.
Verification of nuclear data for the Tsukuba plan, a newly developed treatment
planning system for boron neutron capture therapy. Appl Radiat Isot
2015;106:111–5.

[3] Wielopolski L, Capala J, Pendzick NE, Chanana AD. Patient positioning in static
beams for boron neutron capture therapy of malignant glioma. Radiat Med
2000;18:381–7.
[4] Kiger W, Albritton JR, Lu X-Q, Palmer MR. Development and application of an
unconstrained technique for patient positioning in fixed radiation beams. Appl
Radiat Isot 2004;61:765–9.
[5] Cacicedo J, Perez JF, Ortiz de R, del Hoyo O, Casquero F, Gómez-Iturriaga A,
et al. A prospective analysis of inter- and intrafractional errors to calculate CTV
to PTV margins in head and neck patients. Clin Transl Oncol 2015;17:113–20.
[6] Pang PP, Hendry J, Cheah SL, Soong YL, Fong KW, Wee TS, et al. An assessmentof the magnitude of intra-fraction movement of head-and-neck IMRT cases
and its implication on the action-level of the imaging protocol. RadiotherOncol 2014;112:437–41.
[7] Takada K, Kumada H, Isobe T, Terunuma T, Kamizawa S, Sakurai H, et al.
Whole-body dose evaluation with an adaptive treatment planning system forboron neutron capture therapy. Radiat Prot Dosimetry 2015;167:584–90.
[8] Kumada H, Yamamoto K, Matsumura A, Yamamoto T, Nakagawa Y, Nakai K,
et al. Verification of the computational dosimetry system in JAERI (JCDS) forboron neutron capture therapy. Phys Med Biol 2004;49:3353–65.
[9] Horiguchi H, Nakamura T, Kumada H, Yanagie H, Suzuki M, Sagawa H.
Investigation of irradiation conditions for recurrent breast cancer in JRR-4.
Appl Radiat Isot 2011;69:1882–4.
[10] Sato T, Niita K, Matsuda N, Hashimoto S, Iwamoto Y, Noda S, et al. Particle andheavy ion transport code system, PHITS, version 2.52. J Nucl Sci Technol
2013;50:913–23.
[11] Liem PH. Development and verification of BATAN’S Standard, two-dimensional
multigroup neutron diffusion codes (BATAN-2DIFF). Atom Indonesia
1994;20:1–19.
[12] Auterinen I, Serén T, Anttila K, Kosunen A, Savolainen S. Measurement of freebeam neutron spectra at eight BNCT facilities worldwide. Appl Radiat Isot
2004;61:1021–6.
[13] Sakurai Y, Kobayashi T. The medical-irradiation characteristics for neutron
capture therapy at the heavy water neutron irradiation facility of Kyoto
University Research Reactor. Med Phys 2002;29:2328–37.
[14] Capoulat ME, Minsky DM, Kreiner AJ. Computational assessment of deepseated tumor treatment capability of the 9Be(d, n)10B reaction for acceleratorbased boron neutron capture therapy (AB-BNCT). Phys Med 2014;30:133–46.

Please cite this article in press as: Takada K et al. Development of Monte Carlo based real-time treatment planning system with fast calculation algorithm
for boron neutron capture therapy. Phys. Med. (2016), http://dx.doi.org/10.1016/j.ejmp.2016.11.007