Method and apparatus for simulating charged particle transport in an external magnetic field, and its application in radiotherapy planning

Abstract

A method for simulating the transport of charged particles in an external magnetic field B is provided, in which the unit direction vector of the charged particle, according to Newton's second law u′(t) = F(u(t), η(t)), includes a force term F, which depends on the unit direction vector u(t) and further nonlinearly depends on the particle's motion variable η(t). According to this method, Newton's second law is solved using an implicit finite difference method, where, over the duration of a time step, the force term is approximated by interpolated values ​​of the unit direction vector and the motion variable, where the interpolation of each quantity is a combination of the values ​​of that quantity at the endpoints t and t+Δt of the time step. This method can be used to calculate the dose when a radiation delivery system 300 delivers a treatment plan P to a patient 310.

Description

Technical Field

[0001] This disclosure relates to the field of numerical simulation of charged particle transport. It addresses a problem in the field of radiotherapy: generating radiotherapy plans for delivery in the presence of an external magnetic field. Background Technology

[0002] Treatment planning systems are expected to find a treatment plan that describes how to deliver a prescribed dose of radiation to a tumor while minimizing radiation to surrounding at-risk tissues and organs. To this end, treatment planning systems can calculate the dose a patient will receive when a linear accelerator operates according to a provisional treatment plan. This is particularly relevant in so-called magnetic resonance-guided radiotherapy (MRgRT), where a magnetic field is used to perform real-time magnetic resonance imaging (MRI) of the patient. It is well known that external magnetic fields also affect the dose delivered to patient tissues (dose distribution shift, electron return effect, electron flow effect). Furthermore, the use of magnetic forces for local dose adjustment is envisioned; for example, a suitable local magnetic field can deflect charged particles around difficult-to-locate dangerous organs. Therefore, it is desirable to provide a dose calculation method that also considers magnetic effects.

[0003] Monte Carlo dosimetry calculations are well known to provide accurate approximations and enable high-level detail in material analysis and precise descriptions of particle interactions. Such dosimetry calculation methods are made possible through Monte Carlo-based particle transport simulations. Implementation of this dosimetry method can include particle transport simulations as a component. An early method for simulating electron transport in the presence of electric and magnetic fields was proposed in AF Bielajew's "Electron Transport in E and B Fields" (see: TM Jenkins et al., Monte Carlo Transport of Electrons and Photon, Ettore Majorana International Scientific Series, Vol. 38, Plenum Press, 1988). Here, the Lorentz force is integrated by evaluating the change in direction based on the particle's orientation and energy prior to the step (1-P method). Bielajew's method is used in the EGSnrc and GPUMCD software.

[0004] A newer Monte Carlo-based method for simulating particle transport, incorporating magnetic effects, is disclosed in VNMalkov et al.'s paper, "Charged particle transport in magnetic fields in EGSnrc," Med. Phys., vol. 43, no. 7, pp. 4447-4458 (2016). This method employs a three-point integral method (3-P method) that incorporates multiple scattering when evaluating orientation changes.

[0005] Furthermore, the PENELOPE code system uses analytical solutions for circular orbits. See F. Salvat et al., PENELOPE-2018: A Code System for Monte Carlo Simulations of Electron and Photon Transport, OECD Nuclear Energy Agency, NEA / MBDAV / R(2019)1(2019).

[0006] One quality factor for comparing the performance of particle transport simulation methods is the maximum step size they support under given conditions. Referring to equations 19.32–19.37 in Bielajew’s work above, these equations introduce various dimensionless parameters that must be kept small for convergence to be expected. Generally, the required smallness of the parameters can be achieved by reducing the step size, but this will result in a corresponding increase in computational complexity for simulating particle paths of a given length. For the parameter δ, which represents the change in direction (in radians) per step, it has been shown that 3-P integrals perform satisfactorily up to δ = 0.2, while 1-P integrals are generally limited to even smaller values ​​of δ ≤ 0.02. Both Malkov (2016) and PENELOPE recommend limiting δ to ≤ 0.02. Summary of the Invention

[0007] The purpose of this disclosure is to enable methods and apparatus to accurately and efficiently simulate the transport of charged particles in an external magnetic field. In particular, the methods and apparatus should be able to operate with acceptable large step sizes. Another objective is to provide such methods and apparatus suitable for supporting radiation therapy dose calculations. A further objective is to solve the simulation of charged particle transport governed by nonlinear equations of motion. Yet another objective of this disclosure is to provide a method for calculating the dose distribution in an irradiated medium based on a treatment plan including machine-oriented instructions.

[0008] At least some of these objectives are achieved by the invention as defined in the independent claims. The dependent claims relate to advantageous embodiments of the invention.

[0009] In a first aspect of the invention, a method is provided for simulating the transport of charged particles in an external magnetic field B according to a model in which the unit direction vector of the charged particles is determined by Newton's second law.

[0010] u′(t)=F(u(t),η(t)),

[0011] This includes a force term F, which depends on the unit direction vector u(t) and further nonlinearly on the particle's motion variable η(t). In this method, Newton's second law is solved using an implicit finite difference method, where the force term is approximated over a time step duration Δt by interpolated values ​​of the unit direction vector and the motion variable. The interpolation of each quantity is a combination of its values ​​at the endpoints of time steps t and t+Δt.

[0012] A key advantage of this method is its ability to operate over a wide range of step sizes. More precisely, numerical evaluations show that it does not suffer any significant loss of accuracy even when the direction changes by as much as δ = 0.3 radians per time step. Furthermore, the method proposed in the first aspect can be readily integrated with other methods—such as those used in the applicant's RayStation. TM The software application combines the random hinge method used in the PENELOPE code system. It can also be applied to each individual point in the 3-P method, potentially enabling larger strides.

[0013] In some embodiments, the interpolation used to approximate the force term is a linear combination of endpoint values. Specifically, the interpolation can be an average value, thereby:

[0014]

[0015] Each of these options provides a simple and robust approximation of Newton's second law.

[0016] In some embodiments, the simulation method further includes simulating particle interactions, such as multiple scattering between particles and a medium. This interaction is simulated as discrete events occurring between two time steps using a finite difference method.

[0017] The charged particles in the simulated transport can be any charged particle or atomic nucleus, possibly including positrons, pions, and muons. In some embodiments, the charged particles are electrons, protons, helium ions, or carbon ions.

[0018] The method according to the first aspect can be integrated into a dose distribution calculation method that takes a treatment plan as its input, by which charged particles irradiate a medium in the presence of a magnetic field. The method includes simulating particle transport involving interactions with the medium, wherein the dose distribution is calculated as a result of the simulated interactions.

[0019] In a second aspect of the invention, a treatment planning system configured to perform the methods described above is provided. The invention further relates to a computer program comprising instructions for causing a computer, or particularly a treatment planning system, to perform the methods described above. The computer program may be stored or distributed on a data carrier. As used herein, "data carrier" can be a temporary data carrier, such as modulated electromagnetic waves or light waves, or a non-transient data carrier. Non-transient data carriers include volatile and non-volatile memories, such as permanent and non-permanent storage media of magnetic, optical, or solid-state types. Still within the scope of "data carrier," such a memory may be fixedly mounted or portable.

[0020] This disclosure relates to "Newton's second law" in a general sense and its "force term". In some embodiments, the differential equation to be solved is not the classical form of Newton's second law with a velocity v(t) in units of 1 m / s, but rather Newton's second law formulated for a unit direction vector.

[0021]

[0022] The unit for this quantity can be 1 / s. The classical form of Newton's second law can be transformed into a unit direction vector formula by substitution of variables; see Bielajew (1988). In this formula, although the right-hand side of Newton's second law has a different unit than that of classical forces, it is still referred to as the force term.

[0023] Generally, unless otherwise expressly defined herein, all terms used in the claims are to be interpreted according to their ordinary meaning in the art. All references to “a / an / element, device, component, measure, step, etc.” are, unless otherwise expressly stated, publicly interpreted as referring to at least one instance of an element, device, component, measure, step, etc. Unless expressly stated otherwise, the steps of any method disclosed herein need not be performed in the exact order described. Attached Figure Description

[0024] Aspects and embodiments will now be described by way of example with reference to the accompanying drawings, in which:

[0025] Figure 1 This is a flowchart of a method for calculating dose distribution according to an embodiment;

[0026] Figure 2 This is a block diagram of the treatment plan system;

[0027] Figure 3 Details of the radiation delivery system are shown;

[0028] Figure 4 A multi-leaf collimator suitable for shaping radiation beams is shown;

[0029] Figure 5 and Figure 6 It is two projections of the electron trajectory simulated based on a simplified model and using one-point (1-P) and two-point (2-P) integration methods with equal step δ = 0.2 compared to the small-step 1-P integral that is understood as the basic reality;

[0030] Figure 7 This is the differential dose distribution, which shows the dose calculated by integrating 1-P with δ = 0.3. Figure 9 The extent to which the basic true dose distribution shown (calculated by small-step 1-P integral) is represented;

[0031] Figure 8 This is another differential dose distribution, which shows the extent to which the dose calculated by 2-P integral with δ = 0.3 according to an embodiment of the invention is the same as the fundamental true dose distribution; and

[0032] Figure 9 It shows Figure 7 and Figure 8 The basic true dose distribution referenced. Detailed Implementation

[0033] Various aspects of this disclosure will now be described more fully below with reference to the accompanying drawings, in which certain embodiments of the invention are illustrated. However, these aspects may be embodied in many different forms and should not be construed as limiting; rather, these embodiments are provided by way of example so that this disclosure will be thorough and complete, and will fully convey the scope of all aspects of the invention to those skilled in the art. Throughout the specification, similar numerals refer to similar elements.

[0034] One intended application of the charged particle transport simulation method according to the present invention is a dose calculation method 100, which can be adapted for MRgRT or magnetically enhanced dose delivery. Another application of the charged particle delivery simulation method is to calculate the trajectory of electrons in a magnetic field that has been generated for macroscopic deflection purposes (e.g., around dangerous organs). For this second purpose, charged particles entering the human body with energies in the range of 150-200 MeV may be deflected by a magnetic flux density of approximately 5-10 T.

[0035] Method 100 can be used to calculate the dose distribution in a medium to be irradiated with charged or neutral particles according to a treatment plan. For example, the medium may be exposed to external photon radiation, and charged particles are released when photons interact with the medium. When irradiation is performed in the presence of an external magnetic field B, the Lorentz force acts on each charged particle. The charged particles being transported can be any charged particle or atomic nucleus. For example, the charged particles can be electrons, protons, helium ions, or carbon ions.

[0036] The calculated dose distribution can be expressed as a dose as a function of position, wherein the medium is divided into dose grids, such as three-dimensional grids of voxels with a size of 1 to 5 mm, preferably less than 2 mm.

[0037] exist Figure 1 The method 100, visualized as a flowchart, can be executed on a general-purpose computer, a computer suite, or a dedicated processing system. Specifically, the method 100 can be implemented in a radiotherapy planning system.

[0038] Figure 2 An exemplary treatment planning system 200 is depicted in block diagram form. The treatment planning system 200 includes an interface 210, a memory 220, and processing circuitry 230. Data communication between these components may be via an internal data bus, or, in a distributed implementation, via a data network. In such a distributed implementation, components of the treatment planning system 200 may be provided as a networked (“cloud”) resource. The interface 210 is configured to communicate with external parties such as other processors or operators. Using the interface 210, the treatment planning system 200 can obtain a treatment plan for which dose distributions are to be calculated, output the calculated dose distributions for operator viewing, receive configuration inputs from the operator, and so on. Needless to say, if a treatment plan has already been generated within the treatment planning system 200, it can be transmitted internally without going through the interface 210. The memory 220 may store configuration settings, data describing the radiation delivery system 300, descriptions of the patient and their prescribed doses, and program code 221 executable by the processing circuitry 230. The processing circuitry 230 may include one or more processors. For example, it includes a central processing unit (CPU), a graphics processing unit (GPU), and / or a hardware accelerator.

[0039] Within method 100, the first step 110 includes obtaining a treatment plan, which includes a treatment plan suitable for controlling the radiation delivery system 300 (see [link to method 100]). Figure 3The instructions in the treatment plan can be machine-oriented (or machine-level) instructions that define the behavior of actuators in the radiation delivery system 300—including actuators associated with a radiation source, radiation source gantry, or a patient bed capable of rotating. The instructions in the treatment plan can typically be executed by these actuators or by control circuitry in the radiation delivery system 300 without substantial alteration. Typically, the instructions in the treatment plan are not expressed in terms of clinical effect or delivered dose. To calculate the dose delivered to the medium when the instructions are executed, method 100 can be relied upon. A possible use of method 100 is to calculate the dose distribution that would result if a provisional treatment plan were executed, for example, as validation of the treatment plan. If the calculated dose distribution is unsatisfactory, it may be necessary to modify the treatment plan or regenerate it with adjusted settings.

[0040] In step 114, a particle transport simulation is performed. The particle transport simulation can be a computer-implemented Monte Carlo simulation. At a general level, this simulation includes sampling incident particles (i.e., neutral or charged particles impacting the medium to be irradiated), simulating the interaction between the medium and the particles, and simulating further interactions between the medium and further particles created by said interaction. The incident particles—including their position, orientation, and / or energy—are sampled according to the treatment plan. In some embodiments, the incident particles are not sampled directly from the treatment plan, but from a flux field calculated based on the treatment plan in the aforementioned optional step 112. (The concept of a flux field is...) Figure 4 (As illustrated in the description.) In some embodiments, particle transport simulation 114 may include the following sub-steps:

[0041] - The particles sampled according to the treatment plan will be transported to the boundary of the dose grid.

[0042] - Sample the distance for the next interaction and transport the particle there.

[0043] - Sample the type of interaction that will occur based on the total cross-section. This interaction can be multiple scattering. The total effect of multiple scattering over the distance the particle travels can be simulated as discrete interaction events (condensed history).

[0044] - Initiate any new particles created during an interaction (catastrophic event) in the simulation model, and sample the particle's energy and orientation based on the differential cross-section after the interaction. This initiation may include instantiating a data structure representing the new particle.

[0045] - Place newly launched particles on the stack for later simulation.

[0046] - Continue transporting the particle until it loses all its energy or leaves the dose grid.

[0047] In the next step 116, the dose distribution resulting from the simulated interactions is calculated. More precisely, the energy loss of electrons as they move through voxels can be deposited in the voxels, and the density of the material in the voxels is converted into a dose. In some embodiments, particular care is taken at voxel boundaries. For example, bookkeeping associated with energy deposition may be simpler (or inconsistent) if the step size is adjusted to avoid occasional voxel boundary crossings. More details can be found in Malkov's (2016) paper.

[0048] The technique for simulating the transport of charged particles according to the present invention is used in step 114, and more precisely, when the particle to be transported is charged. This technique may be characterized by a two-point (2-P) integral method, which considers energy loss and magnetic field bending (or deflection) when calculating magnetic field deflection in a single simulation step. Similar to some existing methods, multiple scattering and other interactions are treated separately from magnetic field bending. However, since the two-point integral is more accurate than the single-step method, the magnetic field bending size can be larger without any significant loss of accuracy. As can be seen from Bielajew's work above, the described method is effective if the step is short enough such that (i) energy loss is small, (ii) the magnetic field is approximately constant, and (iii) directional change is minimal. If we assume small energy loss and a constant magnetic field, the allowable step size for the angular change δ radians per step is...

[0049]

[0050] Where q is the particle's charge, c is the speed of light, and the magnetic field has been decomposed into... in

[0051]

[0052] Numerical evaluations of the model have been performed for step sizes corresponding to δ = 0.01, 0.02, 0.05, 0.1, 0.15, 0.2, and 0.3, and satisfactory accuracy has been achieved up to at least δ = 0.2. For step size s, the simulation step will have a duration given by the following formula.

[0053]

[0054] Where β(t) is the particle velocity at time t, and its unit is the speed of light c.

[0055] According to the embodiment, Newton's second law applies to the unit direction vector u(t) of a charged particle: u′(t) = F(u(t), η(t)).

[0056] The force term F(u(t), η(t)) is solved using an implicit finite difference method, where it is approximated by constant values ​​over the duration Δt of the time step. This force term depends on the unit direction vector u(t) and the motion variable η(t). More precisely, it is approximated by interpolated values ​​of the unit direction vector and the endpoint values ​​of the motion variable (t and t+Δt are the endpoints of the time step), as follows:

[0057]

[0058] Where Δu represents magnetic curvature (deflection) such that u(t+Δt) = u(t) + Δu. In this expression, p(x,y) represents an interpolation function that linearly or nonlinearly combines the independent variables x and y. In one embodiment, the interpolation function is a linear combination of parameters:

[0059]

[0060] Where p1 and p2 are arbitrary positive constants. Here, p1 ≠ p2 or p1 = p2. In the latter case, the interpolation function maps the independent variables x and y to their arithmetic mean:

[0061]

[0062] For nonlinear F, this approximation differs significantly from the approximation used in the classical Heun method, namely:

[0063]

[0064] The motion variable η(t) can be total energy, kinetic energy, linear momentum, or velocity. If the force term equals...

[0065]

[0066] The total energy E(t) can then be identified as the motion variable η(t). According to the implicit finite difference method in this embodiment, Newton's second law is approximated as follows:

[0067]

[0068] The expression can be simplified by introducing symbols.

[0069] u(t)=u || +u ⊥ ,

[0070] u(t+Δt)=u || +(u ⊥ +Δu)N′,

[0071]

[0072] Wherein, component u || Parallel to the magnetic field Therefore, we can conclude that:

[0073]

[0074] The cross product is rewritten in matrix form as in

[0075]

[0076] (The symbol u is used to represent the unit direction vector component and its column matrix representation.) The result of multiple cancellations and groupings of the coefficients is...

[0077]

[0078] get:

[0079]

[0080] Where I is the identity matrix. Finally, inserting N′=1 (particle energy conservation under Lorentz force) and re-identifying the cross product simplifies this to

[0081]

[0082] Or, equivalently,

[0083]

[0084] The above work has shown that the new position of the electron can be calculated based on the deflection Δu.

[0085]

[0086] At this stage, the following standardizations may be applied optionally:

[0087]

[0088] This limits the propagation of floating-point inaccuracies. Recall that the parallel and perpendicular components refer to u(t).

[0089] The total energy value E(t+Δt) at the later endpoint (appearing in α) can be calculated, for example, by evaluating the stopping power function of the distance s traveled by the particle during the time step. In this disclosure, the stopping power S(T) at the kinetic energy T is defined as the average energy loss per unit path length and can be modeled as...

[0090]

[0091] Where N is the number of scattering targets per unit volume, and This is the differential cross section interacting with the energy loss W. This stopping power can be divided into collision stopping power and radiation stopping power, with collision stopping power dominating among the other contributions for particle energies up to 10 MeV. See MJ Berger, “Electron stopping powers for transport calculations,” in: TM Jenkins et al. (eds.), Monte Carlo Transport of Electrons and Photon, Ettore Majorana International Scientific Series, Vol. 38, Plenum Press, 1988. In embodiments of the invention, the energy loss can be modeled as a confined collision stopping power.

[0092] Figure 3 An example of a radiation delivery system (or radiotherapy machine) 300 is schematically shown, configured to operate according to a radiotherapy plan P, which may have been determined by a treatment planning system. As mentioned initially, the treatment plan P is generated with the aim of delivering a prescribed radiation dose to at least one target while delivering radiation not exceeding a threshold amount to any nearby organs of danger. The radiation delivery system 300 is configured to perform therapy with respect to a treatment volume 310 within the patient's body by controlling at least one radiation source 305 of the system 300 according to the treatment plan P. The radiation source 305 may be, for example, a linear accelerator. It may be integrated into a bench 315.

[0093] The radiation delivery system 300 may further include means for generating a magnetic field with a defined orientation. For example, the magnetic field may be parallel to the direction of the radiation source 305. In another example, the magnetic field is oriented along the patient's longitudinal axis (foot to head).

[0094] Figure 4An exemplary multi-leaf collimator (MLC) 400 is shown, which can be used as a blocking device in a radiation delivery system 300 to control radiation from at least one radiation source 305 therein. The MLC 400 has a set of collimator blades 400-1, 400-2, ..., 400-(n-1), 400-n, arranged in opposite pairs. Typically, each blade has a width between 5 mm and 10 mm. The MLC 400 is arranged in the path of the radiation field and, by selectively moving the collimator blades 400-1, 400-2, ..., 400-(n-1), 400-n in the direction of travel T, forms a shape that can change the exposure flux region 410 of the field on the beam plane 420. Blade adjustment can be performed continuously or in discrete steps. The flux transmitted through the openings in the MLC 400 is uniform at a given time point. However, the modulation flux distribution in the beam plane 420 can be generated by the temporal superposition of multiple fluxes with openings of different shapes in the MLC 400.

[0095] In the dynamic multi-leaf collimation (DMLC) therapy to which this invention applies, the collimator blades 400-1, 400-2, ..., 400-(n-1), 400-n move during irradiation. In the segmented multi-leaf collimation (or so-called step-ejection) therapy to which this invention also applies, the collimator blades 400-1, 400-2, ..., 400-(n-1), 400-n are stationary during irradiation. Here, the radiation beam is turned off when the collimator blades 400-1, 400-2, ..., 400(n-1), 400-n are repositioned.

[0096] Intensity-modulated radiotherapy (IMRT) and volume-modulated arc therapy (VMAT) are two important applications of DMLC treatment. IMRT treatment consists of a series of continuously delivered static beams. Irradiation is stopped when a gantry 315 holding at least one radiation source 305 is rotated to a position for the next beam. VMAT treatment delivery is performed via a gantry 315 that rotates continuously during irradiation. Therefore, the treatment energy flux here is delivered over an arc of up to 360 degrees. In VMAT, the radiotherapy plan can consist of multiple potentially overlapping arcs. IMRT and VMAT treatments are typically delivered through a set of identical treatment sections, for example, thirty treatments per day over a six-week period.

[0097] Ion-based radiotherapy involves irradiating the patient with protons or heavier ions. This is typically done in the form of pencil beam scanning (PBS), which involves delivering radiation to the patient in the form of beams from different beam directions. Each beam comprises multiple points, up to several thousand, each defined by its location, orientation, weight, and energy. The beam orientation depends on the position and angle of the gantry, as well as the position and orientation of the recliner. Before entering the patient, the orientation of each point can be adjusted using magnets so that the beam will cover a larger area of ​​the patient. The point weight indicates the number of particles included in the point, and the energy determines how far the particles will travel within the treatment volume 310. Thus, the orientation and energy determine where the Bragg peak of that point will be located within the patient's body; that is, where most of its energy will be deposited.

[0098] Figure 5 This is an xz-plane view of the electron trajectory simulated according to a simplified model, where only energy loss due to water-based blocking power and angular deflection caused by the Lorentz force from a 1.5T magnetic field in the y-direction are considered. The electron initiates with a kinetic energy of 2MeV, where... The 1-P integration method with δ = 0.2 was used, and according to an embodiment of the invention, the 2-P integration method with the same δ = 0.2 was used to make the electron step. The results are plotted using the symbols × and +, respectively. The solid line plots represent the calculation results using the 1-P method with δ = 0.001. Since the 1-P method is known to perform well in the "simple" case with very small step sizes (i.e., when particles are transported in very small increments), these results can be considered representative benchmarks or basic realities. Figure 5 As can be seen, the trajectory simulated using the 2-P method differs much less from the actual trajectory. Therefore, it can be assumed that the 2-P method can operate with a much larger stride δ than the traditional 1-P method without causing any significant loss of accuracy.

[0099] Figure 6 This is a yz-plane view of the same simulated electron trajectory. This seems to confirm... Figure 5 A good impression of the performance of the inventive 2-P method.

[0100] Figures 7 to 9Simulations were conducted using a model consisting of a 20×20×20 cubic centimeter water column containing an air cylinder. The air cylinder had a radius of 2 cm, located 8 cm above the origin (x) and 3 cm to the side (z), with the longitudinal direction of the cylinder corresponding to the magnetic field direction. The idea behind the model (phantom) was to evaluate how different simulations compare around the air cavity. Due to abrupt shifts in the medium density, differences in dose calculations were expected near the air cavity. Simulations were performed using the air cavity model with a 5×5 cm field having incident energies of 0.5, 2, and 6 MeV and magnetic flux densities of 1.5 and 0.35 T respectively. The statistical uncertainty of the simulations was set to 0.1% to minimize random noise, ensuring that only systemic differences were observed. The simulations were compared by visually examining the differences between simulations with δ = 0.01 and simulations with δ = 0.3.

[0101] The comparison results between 6MeV and 1.5T are as follows: Figure 7 and Figure 8 As shown. Compared to the basically true results ( Figure 9 Compared to ), the 1-P integral with a large δ ( Figure 7 ) seems to introduce system differences, while the 2-P integral with equal δ ( Figure 8 The 1-P method has the most pronounced drawbacks in high gradient regions, at the dose accumulation point of the beam entering the imagined patient, within and around the air cavity, and at the beam edge. In high gradient regions, the difference is on the order of 1% or less, which may be acceptable. However, the 2-P integral produces almost no systematic difference, only some seemingly random differences in the air cavity, which may be attributed to random noise and some overestimation only within the air region. One might conclude that the 2-P integral method can use a larger δ value than the 1-P integral method without affecting the accuracy of the simulation. This can significantly reduce the simulation computation time compared to the conservative choice of δ = 0.2 for the 1-P integral method, which is the recommended value in the EGSnrc implementation in Malkov (2016).

[0102] The foregoing has primarily described various aspects of this disclosure with reference to several embodiments. However, as will be readily understood by those skilled in the art, other embodiments are equally possible within the scope of the invention as defined by the appended claims, in addition to those disclosed above.

Claims

1. A computer-based method for supporting radiotherapy dose calculation, the method being configured to simulate an external magnetic field based on a model. B The transport of charged particles in the model, where the unit direction vector of the charged particle is determined by Newton's second law. Including force terms F The force term F Depends on the unit direction vector Furthermore, it depends non-linearly on the motion variables of the particle. The method includes: Newton's second law is solved using an implicit finite difference method, in which the force term is approximated over the duration of a time step by interpolated values ​​of the unit direction vector and the motion variables, where the interpolation of each quantity is the endpoint of that quantity at that time step. t , t +Δ t The combination of values ​​at the location.

2. The method according to claim 1, wherein, At least one of the interpolations is a linear combination of endpoint values.

3. The method according to claim 2, wherein, At least one of the interpolations is the average of the endpoint values.

4. The method according to any one of claims 1 to 3, wherein, The motion variables of the particle It is one or more of the following: total energy, kinetic energy, linear momentum, velocity.

5. The method according to claim 4, wherein, The motion variable at the later endpoint is calculated by evaluating the stopping power function of the distance traveled by the particle. The value of .

6. The method according to claim 4, wherein, The motion variable is the total energy. And the force term is equal to ， in, q It is the charge of the particle. c It is the speed of light, and It is a unit vector that makes .

7. The method according to claim 5, wherein, The motion variable is the total energy. And the force term is equal to ， in, q It is the charge of the particle. c It is the speed of light, and It is a unit vector that makes .

8. The method according to claim 6 or 7, wherein, The finite difference method evolves the unit direction vector as follows: ， in, It is the unit direction vector component perpendicular to the magnetic field, and 。 9. The method according to any one of claims 1 to 3, further comprising: Simulate discrete interaction events that occur between two time steps of the finite difference method.

10. The method according to claim 9, wherein, The simulated interactions include multiple scattering.

11. The method according to any one of claims 1 to 3, wherein, The charged particles are electrons, protons, helium ions, or carbon ions.

12. The method according to any one of claims 1 to 3, wherein the method is performed as part of a particle transport simulation in the irradiated medium according to a radiotherapy plan.

13. A method (100) for calculating the dose distribution in an irradiated medium according to a treatment plan, the method comprising: Obtain a treatment plan (110) that includes instructions suitable for controlling the radiation delivery system (300); Perform particle transport simulation (114), including: - Sample the incident particles according to the treatment plan. - Simulate the interaction between the medium and the particles, and - Simulate further interactions between the medium and other particles created by the aforementioned interactions. Calculate the dose distribution resulting from the simulated interaction (116). The particle transport simulation includes performing the method according to any one of the preceding claims.

14. The method according to claim 13, wherein, The particle transport simulation is a Monte Carlo simulation.

15. The method according to claim 13 or 14, further comprising: The flux field is calculated based on the treatment plan (112). The incident particles are sampled based on the flux field.

16. A treatment planning system (200) includes a memory (220) and a processing circuit (230) configured to perform the method of any one of claims 1 to 15.

17. A computer program product comprising instructions that, when executed by a computer, cause the computer to perform the method according to any one of claims 1 to 15.

Images