Monte Carlo calculation optimization method and neutron capture therapy system

Abstract

A method for optimizing Monte Carlo calculations and a neutron capture therapy system, the method including at least the steps of: simulating a particle transport process using a transport calculation method; calculating particle dose expectation values ​​using a dose counting method; and calculating an average dose expectation value for a voxel grid based on the dose expectation values ​​of each particle. By optimizing the simulation of the particle transport process or the dose counting method, the Monte Carlo calculations are optimized, thereby reducing the time required for developing a treatment plan.

Description

[Technical Field]

[0001] The present invention relates to radiation therapy systems and methods, and more particularly to Monte Carlo optimization methods and neutron capture therapy systems. [Background technology]

[0002] With the development of atomic science, radiation therapy using cobalt-60, linear accelerators, electron beams, etc. has become one of the main methods of cancer treatment. However, conventional photon and electron beam therapy kill tumor cells due to the physical limitations of the radiation itself, but also damages many normal tissues in the beam path. Furthermore, because tumor cells have different degrees of sensitivity to radiation, conventional radiation therapy often has poor therapeutic effects on malignant tumors that have a certain degree of radiation resistance (e.g., glioblastoma multiforme and malignant melanoma).

[0003] To reduce radiation damage to surrounding normal tissues, the concept of targeted chemotherapy has been applied to radiation therapy. Furthermore, for tumors with high radiation resistance, radiation sources with high relative biological effectiveness (RBE), such as proton therapy, heavy ion therapy, and neutron capture therapy, are currently being actively developed. Among these, neutron capture therapy combines the above two concepts. For example, boron neutron capture therapy (BNCT) offers a better cancer treatment option than conventional radiation by combining the specific accumulation of boron-containing drugs in tumor cells with precise beam control.

[0004] In boron neutron capture therapy, boron ( 10 B) Taking advantage of the large capture cross section of the contained drug for thermal neutrons, 10 B(n,α) 7 Li neutron capture and fission reactions4 He and 7 Two types of heavy charged particles of Li are generated, and the total range of the two types of particles is equivalent to the size of a single cell, so that the radiation damage to the living body is limited to the cellular level. When boron-containing drugs are selectively accumulated in tumor cells, by combining them with an appropriate neutron radiation source, the goal of locally killing tumor cells can be achieved without causing excessive damage to normal tissues.

[0005] In order to ensure that radiation particles kill as many cancer cells as possible while minimizing damage to normal cells, a CT or PET scan is usually performed before treating the patient, and information on the tissue material of the human body is obtained from the scan results. A computational model is then constructed based on the material information and the radiation source to simulate the transport process of radiation particles within the human body, and finally the dose distribution of the radiation particles within the human body is obtained. The optimal dose distribution plan for the patient is then selected as the patient's treatment plan.

[0006] Currently, dose calculation modules in treatment planning systems primarily use the Monte Carlo method to simulate radiation particles. While conventional radiation therapy requires simulation of the motion processes of photons and electrons, BNCT treatment requires simulation of the motion processes of neutrons and photons. The Monte Carlo method offers high accuracy because it can accurately model the random walk of radiation particles within the human body and perform realistic simulations. However, its slow convergence speed and long calculation times limit its use. Therefore, the Monte Carlo dose calculation process needs to be optimized and accelerated so that physicists can obtain dose distribution results more quickly. Summary of the Invention

[0007] To overcome the drawbacks of the prior art, a first aspect of the present invention provides a method for optimizing Monte Carlo calculations, including at least the steps of: simulating a particle transport process using a transport calculation method; calculating particle dose expectation values ​​using a dose counting method; and calculating an average dose expectation value for a voxel grid based on the dose expectation values ​​of each particle. By optimizing the simulation of the particle transport process or the dose counting method, the Monte Carlo calculations are optimized, thereby reducing the time required for developing a treatment plan.

[0008] The transport calculation method includes at least the steps of acquiring particle information, calculating a macroscopic cross-sectional area value based on the particle information, acquiring a next transport point or collision point of the particle based on the macroscopic cross-sectional area value, and determining a survival state of the particle. Further, the step of acquiring the next transport point or collision point of the particle based on the macroscopic cross-sectional area value includes sampling the next transport point or collision point of the particle and updating the particle information.

[0009] Furthermore, the step of determining the viability of the particle includes at least the steps of determining whether to recalculate the macroscopic cross-sectional area value if the particle is viable, and stopping the simulation of particle transport if the particle is non-viable.

[0010] The particle information includes at least position information, direction information, and energy information of the particles.

[0011] Furthermore, the material of the voxel grid in which the particle is located is determined based on the position information and direction information of the particle, and a macroscopic cross-sectional area value is calculated based on the energy information of the particle and the material of the voxel grid in which the particle is located.

[0012] Furthermore, the step of determining whether to recalculate the macroscopic cross-sectional area value includes at least the steps of: determining whether the energy of the particle has changed; determining whether the material of the voxel grid where the particle is located has changed; if the energy of the particle or the material of the voxel grid where the particle is located has changed, recalculating the macroscopic cross-sectional area value based on the energy of the particle and the material of the voxel grid where the particle is located; and if the energy of the particle and the material of the voxel grid have not changed, obtaining the next transport point or collision point of the particle based on the macroscopic cross-sectional area value. If the material of the voxel grid where the particle is located has not changed and the energy has not changed, the calculation result of the macroscopic interface in the previous step is used as is, and the next transport point or collision point of the particle is directly obtained based on the macroscopic cross-sectional area value, thereby reducing the number of calculations of the macroscopic cross-sectional area value and thereby reducing the calculation time.

[0013] The dose counting method includes at least the steps of: calculating the track length of each particle in each voxel grid; calculating the dose expectation value of the particles in each voxel grid based on the track length of each particle and the corresponding macroscopic cross-sectional area value; and calculating the average dose expectation value of the voxel grid.

[0014] Furthermore, the step of calculating the mean dose expectation value of the voxel grid comprises calculating using Equation (1):

number

[0015] Furthermore, when one particle passes through the voxel grid twice, the dose expectation value of the particle is counted as two passes, and the total number of particles N remains unchanged. By omitting the merging of dose expectation values ​​when the same particle passes through the same voxel grid twice and directly calculating them as two independent particles, the total number of particles N remains unchanged, ensuring the accuracy of the dose calculation result while reducing the dose calculation time, thereby reducing the calculation time. Furthermore, the dose counting method further includes at least a step of calculating a dose error of the voxel grid.

[0016] Furthermore, the step of calculating the dose error of the voxel grid includes calculating using Equation (2):

number

[0017] Furthermore, if one particle passes through the voxel grid twice, the expected dose value of the particle is counted as two times, and the total number of particles N does not change.

[0018] Furthermore, the particles include neutrons and photons.

[0019] Another aspect of the present invention provides a neutron capture therapy system including an image acquisition module for acquiring medical images of an irradiated body, and a treatment planning module for constructing a three-dimensional voxel model based on the medical images and formulating a treatment plan, wherein the treatment planning module includes at least a transport calculation unit for simulating a particle transport process, and a dose counting unit for collecting statistics of dose expectation values ​​of the particles and calculating an average dose expectation value of a voxel grid based on the dose expectation value of each of the particles.

[0020] Furthermore, the neutron capture therapy system further includes a dose calculation unit for calculating tissue dose values ​​based on average dose expectation values ​​of the voxel grid, and the treatment planning module is further used to develop a treatment plan based on the tissue dose values.

[0021] The neutron capture therapy system provided by the present invention optimizes Monte Carlo calculations by optimizing the simulation of particle transport processes or the dose counting method, thereby shortening the time required for treatment planning. [Brief explanation of the drawings]

[0022] [Figure 1] 1 is a flowchart of a Monte Carlo calculation according to the present invention. [Figure 2] 1 is a flowchart of a simulation of a particle transport process according to the present invention. [Figure 3] 1 is a flow chart of the statistics of dose expectation values ​​of particles of the present invention. [Figure 4] 1 is a flowchart of specific statistics of dose expectation values ​​of particles according to the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0023] DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS Hereinafter, embodiments of the present invention will be described in more detail with reference to the drawings so that those skilled in the art can implement the present invention by referring to the text of the specification.

[0024] Radiation therapy is a common means of cancer treatment today. In recent years, neutron capture therapy has been increasingly applied as an effective means of cancer treatment today, among which boron neutron capture therapy is the most common. Neutrons for boron neutron capture therapy can be supplied by nuclear reactors or accelerators. Boron neutron capture therapy (BNCT) is a method of ionizing radiation using boron ( 10 B) Utilizing the large capture cross section of the contained drug for thermal neutrons, 10 B(n,α) 7 Li neutron capture and fission reactions 4 He and7 The boron neutron capture reaction generates two types of heavy charged particles: Li and Li. 10 B(n,α) 7 The Li neutron capture nuclear reaction equation is shown. These two types of heavy charged particles have an average energy of approximately 2.33 MeV, a high linear energy transfer (LET), and a short range. The linear energy transfer and range of the alpha particle are 150 keV / μm and 8 μm, respectively. 7 The Li-heavy-loaded particles have a range of 175 keV / μm and 5 μm, and the combined range of these two particles is roughly equivalent to the size of a single cell, limiting radiation damage to the living body at the cellular level. When boron-containing drugs are selectively accumulated in tumor cells, combining them with an appropriate neutron radiation source can achieve the goal of precisely killing tumor cells without excessive damage to normal tissues.

[0025] Whether the neutron source for boron neutron capture therapy is derived from a nuclear reactor or from the nuclear reaction of charged particles with target material, a mixed radiation field is generated. That is, the beam contains neutrons and photons ranging from low to high energy. In boron neutron capture therapy for deep-seated tumors, the greater the content of other radiation sources (excluding epithermal neutrons), the greater the proportion of non-selective dose deposition in normal tissues. Therefore, it is necessary to minimize these sources of radiation that cause unnecessary dose deposition. To further understand the neutron dose distribution within the human body, in addition to the quality factor of the beam in air, an embodiment of the present invention uses a human head tissue phantom to calculate the dose distribution, and the quality factor of the beam in the phantom is used as a reference for neutron beam design.

[0026] The International Atomic Energy Agency (IAEA) has proposed five beam quality factors in air for neutron sources used in clinical boron neutron capture therapy. These five proposals can be used to compare the advantages and disadvantages of different neutron sources and can also serve as a reference for selecting neutron production pathways and designing beam shapers. The five proposals are as follows: Epithermal neutron flux>1×10 9 n / cm 2 s Fast neutron contamination<2×10 -13 Gy-cm 2 / n Photon contamination<2×10 -13 Gy-cm 2 / n Thermal to epithermal neutron flux ratio<0.05 Epithermal neutron current to flux ratio>0.7 NOTE: The energy range of epithermal neutrons is 0.5 eV to 40 keV, while the energy range of thermal neutrons is less than 0.5 eV and the energy range of fast neutrons is greater than 40 keV.

[0027] In neutron capture therapy, treatment planning is extremely important in order to maximize cancer cell killing while reducing radiation damage to normal tissues. However, the Monte Carlo method currently used for treatment planning has a slow convergence speed and long calculation times, limiting its use. Therefore, the present invention proposes an optimization method for Monte Carlo calculations that can increase calculation speed. This will be explained below with reference to specific examples.

[0028] As shown in FIG. 1, a first aspect of the present invention provides a Monte Carlo optimization method for increasing calculation speed, which includes at least the steps of: simulating particle transport processes using a transport calculation method; calculating particle dose expectation values ​​using a dose counting method; and calculating an average dose expectation value for a voxel grid based on the dose expectation values ​​of each particle. Optimizing the particle transport process simulation or dose counting method optimizes the dose calculation method, thereby shortening the time required for treatment planning. The principle of Monte Carlo dose calculation is to simulate the random motion of each radiation particle in the human body or a medium and then calculate the particle's dose contribution during motion. When simulating the motion of a single particle, a conventional Monte Carlo algorithm first randomly samples the initial attributes of the source particle, calculates the macroscopic cross-section value of the particle in the medium, and then calculates the collision probability and calculates the next transport or collision point based on the macroscopic cross-section value. If the particle disappears after the collision (is absorbed or blocked), the calculation stops, i.e., the simulation of the particle's motion stops. On the other hand, if the particle is still alive, the above procedure is repeated. To ensure the accuracy of the dose calculation, the voxel model used for dose calculation is characterized by a very large number of grids (e.g., on the order of one million) and small size (e.g., 1 mm). Therefore, there is a high probability that the same material will be filled in the same grid. If the particle energy is constant, the calculated macroscopic cross-sectional area value will not change when it moves to an adjacent grid of the same material.

[0029] As shown in FIG. 2 , in the Monte Carlo optimization method provided by the present invention, the transport calculation method includes the steps of obtaining particle information; calculating a macroscopic cross section value based on the particle information; obtaining a next transport or collision point of the particle based on the macroscopic cross section value; The method includes at least a step of determining the survival state of the particles. The steps of acquiring particle information are as follows. Each radiation particle moves randomly within the human body or medium, and to simulate the particle transport process, particle information must be obtained, and the particle information includes at least position, direction, and energy. In an alternative embodiment, since both neutrons and photons provide radiation to the human body during treatment, the particles whose transport process is simulated include neutrons and photons when formulating a treatment plan.

[0030] The steps for calculating the macroscopic cross section value based on the particle information are as follows. The macroscopic cross-sectional area value of a particle is related to the energy of the particle and the tissue material of the voxel grid where the particle is located. After obtaining particle information, the material of the voxel grid where the particle is located is determined based on the position information and direction information of the particle, and the macroscopic cross-sectional area value is calculated based on the energy information of the particle and the material of the voxel grid where the particle is located.

[0031] The steps for obtaining the next transport or collision point of the particle based on the macroscopic cross section value are as follows: The collision probability is calculated based on the macroscopic cross section value, and the particle's next transport or collision point is sampled. When the particle's next transport or collision point is sampled, the particle information needs to be updated because the particle's position, direction, or energy has changed.

[0032] The steps for determining the viability of a particle are as follows: During the transport process, particles may disappear due to collisions, or may move outside the range of the simulation space, such as outside the 3D voxel model. In such cases, there is no need to continue simulating the particles, and transport simulation is performed only on particles that remain alive within the 3D voxel model.

[0033] After sampling the particle's next collision or transport point based on the macroscopic cross section value, the particle's survival state needs to be determined.

[0034] As shown in FIG. 3, in one alternative embodiment, if the particle is viable, it is determined whether to recalculate the macroscopic cross-section value, and if the particle is non-viable, the simulation of particle transport is stopped.

[0035] More specifically, determining whether to recalculate the macroscopic cross-sectional area value includes at least determining whether the energy of the particle has changed, determining whether the material of the voxel grid where the particle is located has changed, and if the energy of the particle or the material of the voxel grid where the particle is located has changed, recalculating the macroscopic cross-sectional area value based on the energy of the particle and the material of the voxel grid where the particle is located, and if the energy of the particle and the material of the voxel grid have not changed, obtaining the next transport point or collision point of the particle based on the macroscopic cross-sectional area value. If the material of the voxel grid where the particle is located has not changed and the energy has not changed, the calculation result of the macroscopic cross-sectional area in the previous step is used as is, and the next transport point or collision point of the particle is directly obtained based on the macroscopic cross-sectional area value, thereby reducing the number of calculations of the macroscopic cross-sectional area value and thereby reducing the calculation time.

[0036] As shown in FIG. 4, in the optimization method of Monte Carlo calculation provided by the present invention, the dose counting method is: The method includes at least the steps of: calculating the track length of each particle in each voxel grid; calculating the expected dose of the particle in each voxel grid based on the track length of each particle and the corresponding macroscopic cross-sectional area value; and calculating the average expected dose of the voxel grid.

[0037] In Monte Carlo dose statistics, doses are often calculated by estimating track length. Each time a particle moves within a grid, it leaves a track. Based on the track length and macroscopic cross-section, the collision probability of the particle within the grid can be calculated. Multiplying this by the average energy deposition of each particle collision yields the energy deposition of the track. Adding and normalizing the energy deposition of all particles within a grid yields the expected dose for that grid. Traditional dose counting methods first simulate each particle's track, record the track length and the grid number it passes through using two arrays, and then calculate the expected dose for each grid based on the track length and macroscopic cross-section. If a particle passes through a grid multiple times, i.e., if the same number exists in the grid number array, the expected dose values ​​corresponding to the same number are merged, and finally, the expected dose value for that particle is added to the total dose. Because the number of grids a particle passes through is enormous, merging the expected dose values ​​for the same grid takes a long time. The dose counting method proposed in this invention saves a lot of time in dose calculation by eliminating the process of merging doses when the same particle passes through the same voxel grid twice. The steps to calculate the average dose expectation value of the voxel grid are as follows: Calculate using formula (1).

number

[0038] When a particle passes through a voxel grid twice, the particle's expected dose is counted twice, and the total number of particles, N, remains unchanged. By omitting the merging of expected doses when the same particle passes through the same voxel grid twice and calculating them directly as two independent particles, the total number of particles, N, remains unchanged, ensuring the accuracy of the dose calculation results while reducing the dose calculation time, thereby reducing computation time. Suppose the i-th particle passes through the grid twice and the remaining expected doses are t1 and t2, then ti = t1 + t2. If the doses of the same particle are calculated as two independent particles without merging, the result of ti remains unchanged. This means that the accuracy of the dose calculation is not affected, but the dose calculation time is effectively reduced.

[0039] In another alternative embodiment, the dose calculation method further comprises at least the step of calculating a dose error for a voxel grid.

[0040] Calculating the dose error of the voxel grid includes calculating using equation (2).

number

[0041] When one particle passes through the voxel grid twice, the expected dose value of the particle is counted as two times, and the total number of particles N does not change.

[0042] Another aspect of the present invention provides a neutron capture therapy system including an image acquisition module for acquiring medical images of an irradiated body, and a treatment planning module for constructing a three-dimensional voxel model based on the medical images and formulating a treatment plan, wherein the treatment planning module includes at least a transport calculation unit for simulating a particle transport process, and a dose counting unit for collecting statistics of particle dose expectation values ​​and calculating an average dose expectation value of a voxel grid based on the dose expectation value of each particle.

[0043] In another alternative embodiment, the neutron capture therapy system further includes a dose calculation unit for calculating tissue dose values ​​based on average dose expectation values ​​of the voxel grid, and the treatment planning module is further used to develop a treatment plan based on the tissue dose values.

[0044] The neutron capture therapy system provided by the present invention optimizes Monte Carlo calculations by optimizing the simulation of particle transport processes or the dose counting method, thereby shortening the time required for treatment planning.

[0045] Although specific exemplary embodiments of the present invention have been described above to facilitate understanding of the present invention by those skilled in the art, the present invention is obviously not limited to the scope of the specific embodiments, and as long as various changes are made within the spirit and scope of the present invention as defined and determined by the appended claims, these changes will be obvious to those skilled in the art, and all of these changes will fall within the scope of the protection claimed by the present invention.

Claims

1. simulating a particle transport process using a transport calculation method; calculating particle dose expectation values ​​using a dose counting method; and calculating an average expected dose value for the voxel grid based on the expected dose value for each particle.

2. The transportation calculation method includes: acquiring particle information; calculating a macroscopic cross section value based on the particle information; obtaining a next transport or collision point of the particle based on the macroscopic cross section value; 2. The method for optimizing Monte Carlo calculations according to claim 1, further comprising at least a step of determining the survival state of particles.

3. The step of obtaining a next transport point or collision point of the particle based on the macroscopic cross section value includes:

3. The method for optimizing Monte Carlo calculations according to claim 2, further comprising the step of sampling the next transport or collision point of the particle and updating the particle information.

4. The step of determining the viability of the particle comprises: If the particle is viable, determining whether to recalculate the macroscopic cross section value; 3. The method for optimizing Monte Carlo calculations according to claim 2, further comprising at least the step of: stopping the simulation of particle transport if the particle is in a non-viable state.

5. 5. The method for optimizing Monte Carlo calculations according to claim 4, wherein the particle information includes at least position information, direction information, and energy information of the particles.

6. 6. The method for optimizing Monte Carlo calculations according to claim 5, characterized in that the material of the voxel grid in which the particle is located is determined based on the position information and direction information of the particle, and a macroscopic cross-sectional area value is calculated based on the energy information of the particle and the material of the voxel grid in which the particle is located.

7. The step of determining whether to recalculate the macroscopic cross-sectional area value comprises: determining whether the energy of the particle has changed; determining whether the material of the voxel grid in which the particle is located has changed; if the energy of the particle or the material of the voxel grid in which it is located changes, recalculating the macroscopic cross-sectional area value based on the energy of the particle and the material of the voxel grid in which it is located; and obtaining the next transport or collision point of the particle based on the macroscopic cross section value if the energy of the particle and the material of the voxel grid are unchanged.

8. The dose counting method includes: calculating track lengths for each particle within each voxel grid; calculating an expected dose for particles within each voxel grid based on the track length and corresponding macroscopic cross section value for each particle; 2. The method of claim 1, further comprising at least the step of: calculating an average dose expectation value for a voxel grid.

9. said step of calculating an average dose expectation value for a voxel grid comprising: Calculating using equation (1), [Equation 1] where i is the particle number and t ｉ 9. The method of optimizing Monte Carlo calculations according to claim 8, characterized in that: ∑ i = ...

10. said step of calculating an average dose expectation value for a voxel grid comprising:

10. The method for optimizing Monte Carlo calculations according to claim 9, further comprising the step of: when one particle passes through the voxel grid twice, the dose expectation value of the particle is counted as two passes; and the total number of particles N remains unchanged.

11. The dose counting method includes: further comprising calculating a dose error of the voxel grid, including calculation using Equation (2); [Equation 2] where i is the particle number and t ｉ 9. The method of optimizing Monte Carlo calculations according to claim 8, characterized in that: ∑ i = ...

12. The step of calculating dose error for a voxel grid comprises:

12. The method for optimizing Monte Carlo calculations according to claim 11, wherein when one particle passes through the voxel grid twice, the dose expectation value of the particle is counted as two passes, and the total number of particles N remains unchanged.

13. 2. The method of optimizing Monte Carlo calculations according to claim 1, wherein said particles include neutrons and photons.

14. an image acquisition module for acquiring medical images of the subject; a treatment planning module for constructing a three-dimensional voxel model based on medical images and formulating a treatment plan; A neutron capture therapy system, characterized in that the treatment planning module includes at least a transport calculation unit for simulating a particle transport process, and a dose counting unit for collecting statistics of dose expectation values ​​of the particles and calculating an average dose expectation value of a voxel grid based on the dose expectation value of each of the particles.

15. 15. The neutron capture therapy system of claim 14, further comprising a dose calculation unit for calculating tissue dose values ​​based on average dose expectations of a voxel grid, wherein the treatment planning module is further used to develop a treatment plan based on the tissue dose values.

Images