Deterministic treatment planning for ion beam therapy

Abstract

Techniques for deterministic treatment planning for ion beam therapy are disclosed. In an example method, a computing device receives information about a target volume and surrounding anatomy within a patient, each characterized by a depth and set of dimensions. The computing device defines a spatial grid including a first number of steps. The computing device computes, for each step of the spatial grid, a fluence spectrum and a radial distribution based on iterative solutions to the Lagrangian form of the Boltzmann transport equation. The computing device computes corrections to the fluence spectra based on nuclear interactions. The computing device then determines an absorbed dose distribution to the target volume and the surrounding anatomy and generates an ion beam therapy treatment plan based on the absorbed dose distribution.

Description

PATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT DETERMINISTIC TREATMENT PLANNING FOR ION BEAM THERAPY CROSS-REFERENCE TO RELATED APPLICATION

[0001] This international application claims priority to provisional application U.S. Ser. No.63 / 596,515 entitled “Deterministic Treatment Planning For Ion Beam Therapy” and filed on November 6, 2023, the entire disclosure of which is incorporated herein by reference for any purpose. GOVERNMENT SUPPORT CLAUSE

[0002] This invention was made with government support under CA225961 awarded by the National Institutes of Health. The government has certain rights in the invention. FIELD

[0003] The present application generally relates to ion beams in matter, and more particularly relates to techniques for deterministic treatment planning for ion beam therapy. BACKGROUND

[0004] Ion beam therapy can be used to treat cancer caused by an uncontrolled cell growth called a tumor. Charged ions are used to target and destroy tumor cells. When ions are accurately aimed at tumor cells they can damage the tumor cell DNA, thus inhibiting the cells’ ability to reproduce and grow. To accomplish this with efficacy, a treatment planning process is typically used to ensure that a therapeutic dose arrives at the tumor while limiting exposure to the surrounding tissue or organs.

[0005] Some early approaches to treatment planning for ion beam therapy involved pencil beam algorithms based on simple mathematical representations of ion beams and rough approximations of their properties. Such methods were guided by measured data and intuitive arguments and not by rigorous physical models. Moreover, treatment planning based on pencil beam algorithms is known to yield poor accuracy.

[0006] Some later approaches to treatment planning using ion beams involve the modeling of multiple scattering of ions as they pass through matter. For example, Molière’s theory of multiple scattering can be used to predict angular and radial dose distributions resulting from multiple Coulomb scattering. However, such methods have also thus farPATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT yielded inaccurate results because, for example, they did not account for energy straggling and nuclear reactions.

[0007] Some existing systems use Monte Carlo algorithms for all aspects of treatment planning. Monte Carlo algorithms can vary greatly in their performance. Some Monte Carlo algorithms may be accurate, employing rigorous physical models and computational techniques. However, such algorithms may be prohibitively slow for routine treatment planning. Some other Monte Carlo algorithms, sometime referred to as “fast Monte Carlo” algorithms, may perform fast enough to be suitable for ion beam treatment planning. However, fast Monte Carlo algorithms invariably rely on oversimplified physical models and less accurate numerical techniques. As a result, this approach increases systematic errors that are difficult to predict and quantify. In some cases, the resulting systematic errors can reach unacceptable levels that are unsuitable for ion beam treatment planning, certification, or regulatory compliance. Examples of oversimplifications may include, for instance, modeling transport without considering multiple scattering or energy loss straggling near the particle track end, where dose gradient remains strong. These and other oversimplifications may result in significant systematic error.

[0008] Additionally, statistical errors are inherent in all Monte Carlo algorithms. Any substantial reduction of statistical uncertainties is typically associated with a high computational cost. For example, reducing statistical uncertainties by half may require increasing calculation time by a factor of four. This computational cost is particularly high at low dose levels that normal tissues usually receive. This high computational cost may be unacceptable even for systems that rely on fast Monte Carlo algorithms.

[0009] Some existing systems for treatment planning for ion beam therapy utilize deep learning techniques. Such techniques may require a substantial amount of work and specialized training of clinical personnel and are thus not suitable for general deployment. Moreover, existing systems implementing deep learning techniques have achieved accuracy no better than existing fast Monte Carlo methods and are subject to the same statistical uncertainties. Therefore, these methods may not be useful for low dose levels. SUMMARY

[0010] Techniques for deterministic treatment planning for ion beam therapy are disclosed. The deterministic approach disclosed herein calculates fluence spectra for ion beams of a defined shape (e.g., Gaussian beams) by solving the Lagrangian form of thePATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT Boltzmann transport equation involving stepwise line integration along the particle path. The fluence spectra are calculated iteratively by approximating the fluence initially using the Vavilov distribution followed by a normal distribution for later steps. Angular fluence distributions are again calculated iteratively using the Molière distribution and converted to a radial distribution through a coordinate transformation. Contributions to fluence spectra due to nuclear interactions with hydrogen and heavier elements are also calculated using line integration to include fluence contributions from all nuclear reactions that occur along the path of an incident ion. Calculation of these contributions accounts for propagation of reaction products (including scattered ions or recoil ions) in the medium. The calculated fluence spectra are combined with the characteristics of the target material to produce dose distributions.

[0011] Treatment planning using purely deterministic methods improves on existing systems and provides advantageous technical effects, which can include eliminating statistical error and reducing systematic error, as well as being significantly faster. The elimination of statistical error through the elimination of the random quantities that characterize Monte Carlo methods make ion beam treatment safer. For example, the statistical error associated with Monte Carlo methods typically increases with decreasing dose. However, treatment plans must ensure that normal tissues and organs do not receive doses above allowed limits. Thus, developing a treatment plan that requires the use of Monte Carlo methods to plan doses to low-dose regions may not be possible to safely plan without unacceptably large safety margins that cannot utilize the full range of the maximum allowed dose. Relatively large statistical uncertainties in low-dose regions are unavoidable in Monte Carlo methods and are not always reported by the software. Such uncertainties pose the risk of a patient receiving a dose in excess of allowed levels due to a miscalculation. Deterministic methods do not have statistical errors and for this reason are more accurate in low-dose regions. The greater accuracy gained through the use of the techniques disclosed herein thus reduces the risk of overdose, allowing safer treatments that utilize more of the maximum allowed dose to be planned.

[0012] Moreover, the deterministic methods disclosed herein have significantly improved spatial resolution as compared with some existing Monte Carlo methods. For example, existing Monte Carlo methods may significantly underestimate dose in the vicinity of narrow Bragg peaks due to their relatively poor resolution which can cause undesired volume averaging effects. This may result in a treatment plan that underestimates thePATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT maximum dose to a critical volume, which can, in some cases, harm patients. The deterministic methods disclosed herein have a spatial resolution that may be more than a factor of ten higher than that found in some existing Monte Carlo-based systems, thus significantly reducing risk associated with underestimation of the peak dose.

[0013] The deterministic methods disclosed herein are also significantly faster than methods used by existing systems. Faster methods of calculation for treatment planning mean less waiting time for patients. Additionally, better treatment plan optimization can be attained when faster methods for dose calculation are available since plan optimization requires hundreds and thousands of dose calculations. Existing systems require simplifying assumptions to speed up the calculations, which can drive accuracy down. In contrast, with faster, deterministic methods, accuracy can be driven up by pushing algorithm parameters to values that ensure more accurate results or by implementing better physical models. Finally, faster, deterministic methods as disclosed herein can facilitate four-dimensional (4D) dose calculations and adaptive planning.4D dose calculation may account for anatomical changes during a treatment session, such as patient breathing. Similarly, adaptive planning refers to situations in which anatomical changes are encountered during the course of a treatment that may last for several weeks. If such changes are significant, the original treatment plan needs to be replaced by a new one. In both cases, additional dose calculations are necessary which may be prohibitively slow using existing systems. The methods of the present disclosure can be used to, for example, implement 4D treatment planning using currently available hardware.

[0014] The advantages pursuant to the faster speeds of the deterministic methods are not theoretical and have been prototyped on unspecialized hardware and comparatively slow scientific numerical frameworks. The methods disclosed herein have proven to be at least 60 times faster than some existing systems, exhibiting sub-second calculation times in concert with higher accuracies, as described above. Reimplementation of the prototyped methods using efficient coding and specialized hardware are expected to yield performance improvements up to hundreds of times faster still.

[0015] Another advantage of the deterministic methods disclosed herein is the intermediate calculation of fluence spectra en route to the calculation of the treatment beam distributions (e.g., distributions of ion energy or dose distributions). In existing Monte Carlo methods, the fluence spectra have much higher statistical uncertainties than the dose, leading to longer calculation times. Using deterministic methods described herein, these intermediatePATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT spectra can be calculated using interpolation methods at any number of points, whether closely or sparsely placed (i.e., an arbitrary spatial grid). These intermediate fluence spectra can be used for radiobiological modelling and for biological optimization of treatments, both growing trends in ion beam therapy. Many radiobiological models, for example, require the frequency or dose average linear energy transfer (LET). Using the deterministic methods of the present disclosure, these distributions can be calculated in a manner similar to the calculation of the absorbed dose, as described below. The fluence spectra can be used for advanced radiobiological models such as the Local Effect Model (LEM) or the Microdosimetric Kinetic Model (MKM), commonly used for planning carbon beam therapy.

[0016] In an example method for deterministic treatment planning for ion beam therapy, a computing device receives information about a target volume (e.g., a tumor) associated with a patient, the target volume located at a depth within the patient and characterized by a set of dimensions. The information may also include information about the patient’s organs that may be at risk of damage by the planned incident ions. The information may further include information about the patient’s particular anatomy in a volume surrounding the target volume, in which the planned incident ion beams may travel during treatment delivery, for calculation of doses to both the target volume and to surrounding organs at risk. The computing device determines an optimized number of ion beams needed for the treatment and characteristics of each beam (directions, energies, modulation etc.) based on all the input anatomical information. The computing device then defines a first spatial grid including a first number of discrete steps, each step having a step size. The computing device computes, for each step of the first spatial grid, a fluence spectrum based on a first iterative solution to the Lagrangian form of the Boltzmann transport equation and a radial distribution, which includes computing an angular distribution based on a second iterative solution to the Lagrangian form of the Boltzmann transport equation and converting the angular distribution to the radial distribution. The computing device determines an absorbed dose distribution based on the fluence spectra and the radial distributions for the steps of the first spatial grid and generating an ion beam therapy treatment plan based at least in part on the absorbed dose distribution.

[0017] In an example system including one or more processors and one or more computer-readable storage media storing instructions, the instructions, when executed by the one or more processors, can cause the one or more processors to perform operations including receiving information about a target volume associated with a patient and organs atPATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT risk, the target volume and organs at risk located at different depths within the patient and characterized by a set of dimensions. The instructions include determining the number of treatment ion beams, their directions and detailed characteristics of each beam, including ion energy distributions, beam collimation and modulation, and other beam parameters. Determination of optimized values for all these and other example beam characteristics for an individual patient can constitute a complex multidimensional optimization problem. Solution of the multidimensional optimization problem includes accounting for a particular patient’s anatomy, including the locations and dimensions of the target volume and nearby healthy organs at risk, among other factors. Thus, solutions of the multidimensional optimization problem typically involve multiple dose calculations. The instructions accordingly include the performance of a dose calculation in which a first spatial grid is defined, including a first number of discrete steps, each step having a step size. The instructions then include computing, for each step of the first spatial grid, a fluence spectrum based on a first iterative solution to the Lagrangian form of the Boltzmann transport equation and a radial distribution, in which determining the radial distribution includes determining an angular distribution based on a second iterative solution to the Lagrangian form of the Boltzmann transport equation and converting the angular distribution to the radial distribution. The instructions further include determining an absorbed dose distribution based on the fluence spectra and the radial distributions for the steps of the first spatial grid and generating an ion beam therapy treatment plan based at least in part on the absorbed dose distribution.

[0018] An example non-transitory computer-readable medium can storing instructions that, when executed by one or more processors, cause the one or more processors to perform operations including receiving information about a target volume associated with a patient and organs at risk, the target volume and such organs located at different depths within the patient and characterized by a set of dimensions. The instructions include determining the number of treatment ion beams, their directions and detailed characteristics of each beam, including ion energy distributions, beam collimation and modulation, and other beam parameters. Determination of optimized values for all these and other example beam characteristics for an individual patient can constitute a complex multidimensional optimization problem. Solution of the multidimensional optimization problem includes accounting for a particular patient’s anatomy, including the locations and dimensions of the target volume and nearby healthy organs at risk, among other factors. Thus, solutions of the multidimensional optimization problem typically involve multiple dose calculations. ThePATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT instructions accordingly include the performance of a dose calculation in which a first spatial grid is defined, including a first number of discrete steps, each step having a step size. The instructions then include computing, for each step of the first spatial grid, a fluence spectrum based on a first iterative solution to the Lagrangian form of the Boltzmann transport equation and a radial distribution, in which determining the radial distribution includes determining an angular distribution based on a second iterative solution to the Lagrangian form of the Boltzmann transport equation and converting the angular distribution to the radial distribution. The instructions further include determining an absorbed dose distribution based on the fluence spectra and the radial distributions for the steps of the first spatial grid and generating an ion beam therapy treatment plan based at least in part on the absorbed dose distribution.

[0019] In some embodiments, an apparatus is provided, which includes means for implementing part or all of the operations and / or methods disclosed herein.

[0020] In some embodiments, a computer program product is provided, which includes computer instructions that, when executed by a processor, implement part or all of the operations and / or methods disclosed herein.

[0021] The following detailed description, together with the accompanying drawings, will provide a better understanding of the nature and advantages of the claimed invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] The accompanying drawings, which are incorporated into and constitute a part of this specification, illustrate one or more certain examples and, together with the description of the example, serve to explain the principles and implementations of the certain examples.

[0023] FIG. 1 depicts an ion beam treatment system that may be used in connection with an example.

[0024] FIG.2 shows a detailed side view of an ion beam treatment system that may be used in connection with an example.

[0025] FIG.3 shows a process for that can be used for deterministic treatment planning for ion beam therapy, in accordance with some aspects.

[0026] FIG.4 shows a process for that can be used for deterministic treatment planning for ion beam therapy, in accordance with some aspects.PATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT

[0027] FIG. 5 shows a process for generating a treatment plan for ion beam therapy using the deterministic methods for dose distribution calculation described above, in accordance with some aspects.

[0028] FIG.6 shows a process for that can be used for administering a treatment plan generated using deterministic treatment planning for ion beam therapy, in accordance with some aspects.

[0029] FIG.7 shows a simplified block diagram of a computer system suitable for use in some examples. DETAILED DESCRIPTION

[0030] The following description of exemplary embodiments is presented for the purpose of illustration and description. It is not intended to be exhaustive or to limit the claimed embodiments to the precise form described, and persons skilled in the art will appreciate that many modifications and variations are possible. The embodiments (or examples) have been chosen and described in order to best explain their principles and practical applications to thereby enable others skilled in the art to best make and use various embodiments and with various modifications as are suited to the particular use contemplated. Radiation Therapy Systems

[0031] Ion beam therapy is a method for delivering a beam of protons or other ions (e.g., helium ions, carbon ions, etc.) to a target volume such as a patient’s tumor. Ion beams are generated outside the patient and are precisely directed to the tumor site. Ion beams can be configured to deposit the majority of their energy at a specific depth in tissue such that the tumor receives a high dose while the surrounding healthy tissues receive minimal radiation. Ion beams can damage the DNA of tumor cells and thereby cause the death of the malignant cells or prevent them from dividing.

[0032] FIGs. 1 and 2 depict an ion beam treatment system 100 that may be used in connection with an example. FIG.1 shows a perspective view of an ion beam treatment system 100. Ions are typically accelerated to the desired incident energy using a combination of components such as a linear accelerator, a cyclotron, or a synchrotron (not shown) and directed to a patient on a treatment couch 35. Some examples are presented in the context of proton beams, but it will be appreciated by those skilled in the art that the same principles apply to other heavy ion beams such as those made up of helium or carbon ions.PATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT

[0033] System 100 includes a support structure 10 supporting a rotatable gantry 20 with a treatment head 30. In some examples, system 100 includes a control unit (not shown) which includes control circuitry for controlling the different modes of operation of the ion beam treatment system 100. Ions derived from an ion source are accelerated externally to system 100 and guided to support structure 10 and into gantry 20 using an arrangement of vacuum tubes and electromagnets. Gantry 20 may include a series of bending magnets and focusing magnets and rotates around the patient to deliver a multi-direction dose to the target volume. Treatment head 30 may include steering magnets such as dipole magnets for precision control of the ion beam as the dose is delivered to the target volume.

[0034] FIG.2 shows a somewhat more detailed side view of ion beam treatment system 100. A patient P is shown lying on treatment couch 35. Ion beams generated as described above are emitted from a steering nozzle in treatment head 30 in a shaped beam 104, such as a Gaussian pencil beam. In a typical example, a patient plane 116, which is perpendicular to the page in FIG. 2, is positioned about one meter from treatment head 30. The rotatable axis of gantry 20 is located in patient plane 116, such that the distance between the exit window in treatment head 30 and isocenter 178 remains constant when gantry 20 is rotated. Isocenter 178 is a point located at the intersection between patient plane 116 and central axis 122 of beam 104. During treatment planning, an optimized position for the isocenter 178 within the tumor is determined and marked on a patient’s CT image. Patient P can be positioned on treatment couch 35 such that the marked isocenter precisely coincides with the physical isocenter 178 of the treatment ion beam treatment system 100. Multiple Scattering

[0035] A goal of treatment planning for ion beam therapy is to ensure that a therapeutic or prescribed dose arrives at the tumor while limiting exposure to the surrounding tissue or organs. Accurate models of the radiation dose to a tumor and surrounding tissue for a given beam energy and direction are among the building blocks of developing an effective treatment plan. At the same time, modeling the dose cannot take an excessively long time, in part because treatment planning is an iterative process that involves repeatedly modeling the dose. For example, radiation treatment planning may occur before treatment begins or, using the methods of the present disclosure, may occur during the course of treatment (including between treatment sessions) at any time before the treatment course is completed as when, for example, adaptive treatment planning methods are used.PATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT

[0036] FIG.3 shows a deterministic process 300 that can be used for determining the fluence spectra, angular distributions, radial distributions, and dose distributions used in treatment planning for ion beam therapy. Process 300 illustrates determination of these distributions for a particular ion beam of a given direction and incident kinetic energy. Extension of process 300 to a treatment plan including multiple ion beams from multiple directions is described more fully below with reference to process 600. In some examples, process 300 can be implemented in a computing device used for radiation treatment planning.

[0037] At block 302, the computing device receives information about a target volume and surrounding anatomy associated with a patient, the target volume and one or more elements of the surrounding anatomy each characterized by a depth within the patient and a set of dimensions. For example, the target volume can be a tumor to be treated, and the surrounding anatomy can include organs at risk or other anatomical information associated with a patient, such as information about the volume surrounding the target volume. For example, imaging obtained through techniques such as computed tomography (CT) or magnetic resonance imaging (MRI) may be used to delineate the target volume and organs at risk, and to determine their exact locations, shapes, and sizes. Additional imaging modalities, such as positron emission tomography (PET) or ultrasound, may be used to further refine the target volume characterization.

[0038] At block 303, the computing device defines characteristics of an ion beam to be directed toward the target volume. The characteristics can include ion energy distributions, beam collimation and modulation, and other beam parameters. In some embodiments, ion beam characteristics are defined as part of an iterative treatment planning process, as described below, and portions of process 300 can be repeated for each ion beam.

[0039] At block 304, the computing device defines a first spatial grid including a first number of discrete steps, each step having a step size, for a given initial kinetic energy of the ion beam. For example, a step size Δ^^^for numerical line integration along the ion path can be chosen so that the ion loses a fraction ^^ of its kinetic energy ^^^at the start of each step as it travels distance Δ^^^. For a given ^^ the step size is calculated as follows:where ^^^^^^is the range of an ion with an initial energy ^^. Pre-calculated or tabulated values of R(E) can be used, in some examples, to perform this calculation. In typical examples involving proton beams with initial energies in the range 40-220 MeV, an empirically derivedPATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT value of ^^ that results in optimized software performance is 0.05, but other values may be used in accordance with various design objectives.

[0040] Use of this method for determining step size Δ^^^may result in needlessly small steps in the vicinity of lower particle energies near the end of the ion beam track, which can negatively impact performance. A minimum step size Δ^^^^^can be used to prevent this. If step size Δ^^^becomes less than Δ^^^^^during definition of the first spatial grid, then Δ^^^can be set to Δ^^^^^. Thus, beyond a certain depth, all steps sizes may equal Δ^^^^^. A typical value of Δ^^^^^for an incident proton energy of 160 MeV is 0.1 cm.

[0041] Considerations when determining the step sizes Δ^^^include the fact that stepwise calculations may be simplified when particle energy loss per step is small, corresponding to a smaller step size, while overall calculations may be faster with larger step sizes. Other considerations are described below with reference to block 306. An optimized step size can be determined empirically by balancing these considerations. For example, an optimized step size can be determined through an iterative process by calculating fluence spectra using various step sizes and comparing the results with Monte Carlo simulations in view of both performance and accuracy.

[0042] In some examples, treatment planning can be based on representations of the patient as a voxelized phantom with a voxel size of 2-4 mm, with each voxel treated as a homogenous medium representing the average properties of all materials in the voxel. In some cases, the integration steps Δt can be chosen so that within each step the medium within the step is homogeneous. However, such a choice of the integration steps may be suboptimal for performance reasons or due to other constraints. In such cases, a step Δt may span more than one voxel, including multiple media. To account for this heterogeneity, equations involving range can be updated to account for heterogeneity. Additionally, multiplications by the step size Δt can be replaced by integration over a distance Δt. This integration does not add appreciably to the performance demands since in a voxelized phantom all the integrands are piecewise-constant functions.

[0043] At block 306, the computing device computes, for each step of the first spatial grid, a fluence spectrum based on a first iterative solution to the Lagrangian form of the Boltzmann transport equation, given by:PATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCTIn equation ^2^, t is the path length, Φ is the fluence, ^^ is the ion location, ^ Ω^⃗ is a unit vectorin the direction of ion travel, ^^ is the ion kinetic energy, ^^ is the total interaction cross section, and ^^^is the double differential scattering cross section.

[0044] Because equation ^2^ is integro-differential, it can be solved iteratively using line integration, integrating stepwise over the first spatial grid determined in block 306. For example, solutions to equation ^2^ can be used to iteratively calculate fluence spectra with:In equation (3), Φ^ା^^^^^ and Φ^^^^^ are fluence spectra at depths ^^^ା^ and ^^^ , where ^^^ା^ ^ ^^^and Φ ^^^ ∣ ᇱ^ ᇱ^ ^^ is the fluence spectrum at depth ^^^ା^ produced by an ion of energy ^^ thatoriginated at depth ^^ ᇱ^. If the range of an ion with energy ^^ is less than the step size ^^^ା^ െ ^^^,then Φ^^^ ∣ ^^ᇱ^ can be set to 0.

[0045] The path length that an ion has travelled when it reached a given depth can be assumed to be approximately equal to that depth. This simplifying assumption may be used for the determination of fluence spectra at points along the axis of the incident ion beam and is consistent with the narrow angular distribution of the ion beam. For example, for protons travelling in water the difference between the depth and the path length is less than 0.2%. This simplifying assumption is even more accurate for heavier ions such as helium and carbon ions.

[0046] Evaluation of the integral in equation (3) may involve selection of a conditionaldistribution Φ^^^^ ∣ ^^ᇱ^. In some examples, the Vavilov distribution can be used for calculationof the fluence spectrum after the first step, Φ^^^^^. Thereafter, for optimized performance, theconditional distribution Φ^^^^ ∣ ^^ᇱ^ can be approximated using a normal distribution. Thisapproximation is based on the asymptotic properties of the Vavilov distribution that tends to anormal distribution as the step size Δ^^ ൌ ^^^ା^ െ ^^^ increases. Therefore, the step size Δ^^determined in block 306 should be sufficiently large for the normal distribution to be an accurate approximation of the exact Vavilov formula at each subsequent step. In some other examples, the more accurate Vavilov distribution may be used for calculation of the fluence spectra for all steps or only for a few first steps, in accordance with treatment planning performance considerations.

[0047] The Vavilov distribution is a solution of the multiple scattering problem for charged particles, such as ions, that travel a distance t such that the energy losses are muchPATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT smaller than the initial energy of the particles. The solution is based on a relativistic form of the Rutherford formula for the scattering cross section given by:in which:and:In these equations, ^^^^ଶ is the electron rest energy (0.510999 MeV); ^^ ൌ ^^ / ^^ is the ratio ofion velocity to the speed of light; ^^^is the classical electron radius (2.817940 ∙10-13cm); ^^ is the density of the material (e.g., 1 g / cm3for water); ^^^is Avogadro's number(6.022141∙1023); ^^ is the number of electrons per molecule (e.g., ^^ ൌ 10 for water); and ^^ isthe molar mass of the material (e.g., 18.01528 g for water).

[0048] Using the Vavilov distribution, if ions begin with an initial energy ^^^and travela distance ^^, the fluence spectrum expressed in terms of the total energy lost ^^ ൌ ^^^ െ ^^^^^^ isgiven by:PATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT

[0049] In some examples, the normal distribution can be used for the conditionaldistribution Φ^^^^ ∣ ^^ᇱ^. For instance, the Vavilov distribution can be used for calculation ofthe fluence spectrum after the first step, Φ^^^^^. Thereafter, for optimized performance, theconditional distribution Φ^^^^ ∣ ^^ᇱ^ can be approximated using a normal distribution.Subsequent steps may use the normal distribution if it is a sufficiently accurate approximation of the Vavilov formula at each subsequent step. Alternatively, the Vavilov distribution may be used for calculation of the fluence spectra for all steps or only for a few first steps.

[0050] In some examples, a normal distribution with center ^‾^ and width ^^ଶ may beused for the conditional distribution Φ^^^^ ∣ ^^ᇱ^. The average energy ^‾^ can be determined usingthe continuous slowing down approximation. For an ion with energy ^^ᇱat the beginning of astep with size Δ^^, the average energy at the step end, ^‾^ is given by:In equation^13^, ^^^^^ᇱ^is the ion range as a function of ion energy and ^^ି^is the inversefunction. In this context, the “inverse function” refers to E as a function of R. In some examples, ^^^^^ᇱ^and the inverse function ^^ି^may be obtained from precalculated, tabulated ion ranges for various ion energies and materials. Because the precalculated, tabulated ion ranges only include discrete values, interpolation techniques that are known in the art can be used to obtain intermediate values of the continuous variables R and E. The corresponding normal distribution width can be given by:

[0051] Although the integral of equation (3) is calculated at discrete steps ^^^, ^^ଶ, … , ^^^along the first spatial grid that is optimized for the speed and accuracy of the method, fluencespectra may nevertheless be required at arbitrary depths ^^^ ^ ^^ ^ ^^^. Accordingly, aninterpolated fluence spectrum can be determined. For instance, equation^13^can be used to calculate the average ion energy ^‾^^at each step. Using the same equation, the average ionenergy ^‾^ at the desired depth ^^ is also calculated. The ^‾^^ at the two consecutive depths ^^^ and^^^ା^ nearest to ^^ on both sides ^^^^ ^ ^^ ^ ^^^ା^^ can be used to determine a weight:PATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCTand two energy shifts: Δ^^^ ൌ ^‾^^ െ ^‾^ ^16^The interpolated fluence spectrum at arbitrary depth z can then be given by: Φ^^^,^^^ ൌ ^^Φ^^^^, ^^ ^ Δ^^^^ ^ ^1 െ ^^^Φ^^^^ା^, ^^ െ Δ^^^ା^^ ^18^

[0052] In some examples, the computing device can calculate corrections to the fluence spectra determined in block 306 due to attenuation from nuclear interactions. For example, if an ion participates in a nuclear reaction, it is removed from the incident ion beam, since such an interaction typically changes its direction and separates it from the incident beam. Thus, an attenuation factor based on a total nuclear cross section and the step size can be used to account for such interactions. Calculation and application of the correction due to attenuation from nuclear interactions is described in more detail below in the section on Nuclear Interactions.

[0053] At block 308, the computing device computes, for each step of the first spatial grid, an angular distribution based on a second iterative solution to the Lagrangian form of the Boltzmann transport equation. The angular distribution may characterize the lateral spread of ions as the ions pass through matter due to multiple Coulomb scattering during each step. For example, as with the calculation of fluence spectra, the angular distribution can be calculated iteratively using:in which Φ^ା^^^Ω^⃗ ^ and Φ^൫^Ω^⃗ ᇱ൯ are angular distriubtions of fluence at depths ^^^ା^ and ^^^,where ^^^ା^ ^ ^^^ and Φ^൫^Ω^⃗ ᇱ ⋅ ^ Ω^⃗ ൯ is the angular distribution of fluence at depth ^^^ା^ producedby an ion that had direction ^ Ω^⃗ ᇱ when it originated at depth ^^^.

[0054] In some examples, equation ^19^ can be solved using the Molière distribution. The Molière distribution can be used to calculate the angular distribution of ions for a small step size Δt. For protons propagating in water, the Molière distribution can be expressed as:PATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT where:and:In these equations, Δ^^ is the step size, ^^ is the scattering angle, and ^^^is the Bessel function of the first kind of order zero. Equation (20) et seq. as shown are for protons propagating in water. For other ions and propagation media, the same distribution can be used with different values of certain interaction constants.

[0055] The dependence on step size Δt is included in distribution width ^^^ଶwhich is, for a chemical element j:For a compound such as water, this is given by:where ^^^is the weight fraction of the ^^-th element; ^^^is Avogadro's number; ^^ is the protonmass; ^^ is the classic ଶ^ al electron radius; and ^^ ൌ ^^ / ^^^^ . In some examples, the functionscan be precalculated for the ^^ range of interest.

[0056] In equation ^22^, parameter ^^ can be found by solving ^^ ൌ ^^ െ ln ^^ for agiven b. In some examples, precalculated values of ^^^^^^ may be used for b ranges of interest. Parameter ^^ can be calculated as follows:where parameter ^^^ଶfor a compound such as water, is as given above in equation ^24^ and ^^ is Euler's constant (approximately 0.5772156649). The parameter ^^^ଶcan be calculated for water (H2O) using:where:PATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCTin which ^^^is the atomic number of element ^^; ^^ is the fine structure constant (approximately 0.0072973525664); ^^்ிis the Thomas-Fermi constant (approximately 0.88534); and ^^ுி,^is the Hartee-Fock correction for element ^^.

[0057] The integral of equation ^19^ can be numerically evaluated by using the azimuthal symmetry of the scattering of ions due to the Coulomb force to express the Molière distribution as:where Θ is the angle between directions ^ Ω^⃗ and ^ Ω^⃗ ᇱ. The integral over ^^ᇱ can be evaluatednumerically using, for example, a 40-point Legendre quadrature. The integral over ^^ᇱcan be similarly evaluated numerically using, for example, a 20-point Chebyshev quadrature and the identity:

[0058] At block 310, the computing device converts the angular distribution to aradial distribution. For example, the radial distribution of Φ^^^, ^^^ can be derived from theangular distribution Φ^z,^Ω^⃗ ^ through certain geometric approximations. For instance, onestepwise geometric approximation includes using: ^‾^^ା^ ൌ ^‾^1 ^^2^tan^^‾^^ା^^ ^ tan^^‾^^^^^^^^ା^ െ ^^^^ ^33^PATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT at depth ^^^ା^, where ^^ is the radial distance from the beam axis and ^̅^ and ^̅^ denote average values. Using this step-wise geometric approximation, the angular distribution of equation ^31^ can be converted to a radial distribution using the conversion formula:

[0059] As with the fluence calculated at arbitrary depths ^^^ ^ ^^ ^ ^^^ above in equation^18^, the radial distribution can be similarly approximated at arbitrary z. For the radialdistribution between Φ^^^^ ,^^^ and Φ^^^^ା^, ^^^, the interpolation formula is given by:^^ െ ^^^^ ൌ ^^^^ା^ െ ^^^^35^Φ^^^, ^^^ ൌ ^1 െ ^^^Φ^^^^ ,^^^ ^ ^^Φ^^^^ା^, ^^^ ^36^

[0060] In some examples, for calculation of absorbed dose using a physical model of a realistic ion beam, the incident ion beam may be modeled as a Gaussian pencil beam with width ^^ described in two dimensions asIn polar coordinates, this can be written:Integral^38^can be numerically evaluated over ^^ᇱusing, for example, an 80-point Legendrequadrature. Integral^38^can likewise be transformed using the identity from equation^32^and numerically evaluated over ^^ᇱusing, for example, an 80-point Chebyshev quadrature. In cases where the incident beam is better modeled by a beam shape other than a Gaussian distribution, the same calculations can be used, substituting an expression that characterizes the incident fluence of the non-Gaussian beam for equation (37).

[0061] At block 312, the computing device determines whether iteration over the first spatial grid is complete (i.e., whether all steps of the line integral of equation (3) are complete). If the particular step is the last step, then execution of the loop over blocks 308 to 312 can terminate. Otherwise, the process returns to block 306 to compute the fluence spectrum and angular distribution for the next step.PATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT

[0062] At block 314, the computing device determines an absorbed dose distribution to the target volume and the surrounding anatomy based on the fluence spectra and the radial distributions for the steps of the first spatial grid according to the depth and dimensions of the target volume and the surrounding anatomy. For example, the Gaussian pencil beam representation of equation ^38^ may be used to calculate the dose to the target volume at a particular location ^^. The absorbed dose can be calculated using:where ^^ is the medium density and ^^ is the ionic stopping power for a material at point ^^. Absorbed dose is typically reported either as dose-to-water or dose-to-muscle. (Some care providers may have guidelines that specify using one or the other.) In those cases, the ionic stopping power and density either for water or for muscle are used in equation (39)^32^. In some examples, the ionic stopping power of water or other media, such as muscle, may be obtained from precalculated or tabulated data.

[0063] Similar calculations can be performed to obtain the frequency or dose average linear energy transfer (LET), useful distributions for radiobiological modelling and for biological optimization of treatments. These distributions are given by:where ^^^^^,^^^ is proton LET for the material at point ^^; ^^ி^^^^ is the frequency average LET; and ^^^^^^^ is the dose average LET.

[0064] At block 316, the computing device generates an ion beam therapy treatment plan based at least in part on the absorbed dose distribution for the incident ion beam. In a typical implementation the calculations of process 300 are combined with contributions from ion beams from multiple directions and / or incident kinetic energies (which can be calculated in the manner described above) to constitute the ion beam therapy treatment plan. For example, the depth and set of dimensions of the target volume may be used to establish a target volume geometry. The isocenter 178 can be configured such that it is located at an optimized location within the target volume and further such that the axis of all incident ion beams pass throughPATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT the isocenter 178. Each incident ion beam enters the patient through the skin. For each step along the path along the beam axis beginning from the skin to the depth where all ions stop, a corresponding fluence distribution with a two-dimensional extent can be determined using the techniques of the present application. For example, the solution of the Boltzmann equation ina medium for the fluence at coordinates (^^, ^^) as shown in, for example, equation (38), can beused in concert with the absorbed dose at ^^ given by equation^39^to calculate the absorbed dose at all points where the beam deposits significant dose, including the target volume and organs at risk located near the beam path. For non-Gaussian or arbitrary beam profiles, equation (38) can be modified to account for the non-Gaussian beam profile.

[0065] The treatment plan generated in block 316 includes irradiating both the target volume as well as organs at risk located near the ion beam path. As described above, the step size can be adjusted according to performance requirements, and the interpolation formulae described with respect to block 306 can be used to determine the absorbed dose at arbitrary locations near the beam path. For instance, in some examples, the target volume may be a tumor. Irradiation of the tumor using the ion beam may be desired to weaken or eliminate the tumor. As described below with reference to process 600, multiple beams targeting the tumor from different directions can be used sequentially or simultaneously. The techniques of the present disclosure may be used to determine the absorbed dose at various surface and interior locations of a particular tumor geometry as well as the surrounding regions including organs at risk. In some cases, calculation of the absorbed dose along each step of the first spatial grid may be sufficient. However, additional resolution may be desired in the vicinity of certain parts of the tumor. For example, in cases involving heterogeneous tumors, tumors close to vital organs, or moving tumors, higher resolution may be desired. In those cases, the interpolation formulae may be used. Moreover, in addition to estimating the dose to the tumor, estimation of the dose to the surrounding tissue is necessary to avoid damage to healthy tissue. The same methods and interpolation formulae as those used for tumors may be used for precision calculation of dose to these areas. Nuclear Interactions

[0066] Process 300 models fluence spectra and absorbed dose distribution based on multiple Coulomb scattering. The accuracy of dose distribution modeling can be further improved by accounting for additional physical processes. For example, the accuracy of the fluence spectra can be improved by accounting for nuclear interactions, which can attenuate the ion beam or produce additional scattering effects.PATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT

[0067] For instance, deterministic methods can be used to calculate corrections to the fluence spectra determined in block 308 due to attenuation from nuclear interactions. For each step of the first spatial grid, an attenuation factor based on a total nuclear cross section and the step size can be determined. The attenuation factor may account for various nuclear processes, including elastic interactions of ions with hydrogen atoms; elastic interactions of ions with atoms heavier than hydrogen, such as carbon, nitrogen, oxygen, etc.; and inelastic interactions of ions with atoms heavier than hydrogen. In some cases, certain nuclear attenuation effects can be ignored if their effect is sufficiently negligible. For instance, inelastic interactions of protons with hydrogen atoms are negligible in some cases. The cross sections for these reactions arerespectively. For some applications, precalculated or tabulated values of these cross sections may be suitable. The total nuclear cross section is given by:

[0068] Ions that undergo nuclear interactions may experience a change of energy or direction of travel, or they may be absorbed. Thus, for each step of the first spatial grid, a corrected fluence spectrum based on the fluence spectrum determined at block 308 of process 300 and the attenuation factor. Such ions can be removed from the ion beam step-by-step by applying to the fluence Φ^^^^ the attenuation factorwhere Δ^^ is the step length and is the total nuclear cross section determined above.

[0069] FIG. 4 shows a process 400 for using deterministic methods to calculate corrections to the fluence spectra determined in block 308 due to nuclear scattering effects. In some examples, process 400 can be implemented in a computing device used for radiation treatment planning. The elements of process 400 may be performed following the determination of the fluence spectra as described above with reference to FIG.3. In a typical example, process 400 includes steps for correcting and refining the fluence spectra of process 300 to account for nuclear scattering.

[0070] At block 402, the computing device determines a second spatial grid, including a second number of discrete steps, each step including a step size. Because nuclear interactions are less frequent than Coulomb scattering, a coarser spatial grid can be used for calculating the nuclear scattering corrections as compared with the first spatial grid determined at block 306 of process 300.

[0071] In one example, the step size for the second spatial grid for nuclear interactions modelling may begin with a relatively large step Δtmaxat shallow depths, which can bePATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT gradually reduced to a Δtmin towards the track end. For instance, the step size can be reduced linearly or using a fractional energy loss method comparable to the method described above with reference to block 306. In a typical example, for a proton with an incident energy of 160 MeV, a Δtmax of 2 cm and a Δtmin of 0.4 cm may yield performant results in implementing the techniques of process 400 described herein.

[0072] At block 404, the computing device determines, for each step of the second spatial grid a first contribution to the fluence spectrum from elastic interactions of ions with hydrogen atoms. For example, a correction to the fluence at a depth ^^ can be estimated by calculating contributions to scattering across a spectrum of energy bins constituting a fluence spectrum histogram, similar to the Monte Carlo surface tally method, a computational technique that uses statistical sampling to estimate the particle fluence on a surface.

[0073] In an example calculation involving collisions between protons and hydrogenatoms, contributions to fluence spectra at a depth ζ due to incident protons ^^^ ൌ 1^ and recoilprotons ^^^ ൌ 2^ can be found by considering scattered incident protons and recoil protons, forwhich the distance ^^^is less than the continuous slowing down approximation range R(^^^) of an ion with energy ^^ , whe ^ି௧ ^ re ^^^ൌఓ^,^ೌ್. If the cosine of the scattering angle in the center of mass frame is ^^ୡ୫, then: ^^^ାఓ^^lab,^ ൌ ^ଶ ^43^with corresponding kinetic energies:A contribution ofห may be added to each energy bin associated withenergy ^^^∗for scattered ions that contribute to the fluence, as described above.

[0074] In this example, the correction at arbitrary depth ^^ including integrating over ^^^^to account for the angular distribution of scattered protons is given by:PATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT spectrum at depth ζ produced by each of the two scatteredis the total cross section for elastic ^^ → ^^ scattering at point ^^;^^^∗(^^, ^^ୡ୫^ is the energy of a proton with initial energy ^^^^^^ୡ୫^ after it travels distance ^^୧;^^^^^^, ^^ୡ୫^ accounts for attenuation due to nuclear interactions of scattered protons as theytravel distance ^^୧; and ^^^is set to zero, if ^^୧< 0 or ^^୧> ^^^^^^^, otherwise it is calculated as ^^^ൌ⋅ ^^^^. In some examples, the accuracy of this formula can be improved byintegrating ^^ேover distance ^^^instead of multiplication. In some examples, the contributions to the fluence spectramay be represented as a histogram. In that case,^^൫^^ െ ^^∗^ ^^^, ^^^^^൯ represents the contribution to an energy bin to which ^^^∗belongs. While this example applies to protons colliding with hydrogen atoms, a similar procedure can be used to determine corrections to fluence spectra for heavier ions colliding with protons, using instead the mass of the particular ion instead of the hydrogen mass.

[0075] Equation Error! Reference source not found. is a line integral calculated step-wise using steps Δt in analogy to Monte Carlo methods, in which a particle is traced along its path. Methods described herein, however, involve deterministic rather than random steps and are further optimized for best performance. The use of approximations such as zeroing the attenuation due to nuclear interactions outside of the given range can improve computational efficiency while still maintaining adequate accuracy. Equation Error! Reference source not found. can be evaluated as a line integral by following the incident ion along its path. This is performed for a discrete set of ion energies accurately representing fluence Φ൫^^^,^^ᇱ൯, in which the contribution for each energy is weighted accordingly. The integral over ^^^^can be numerically evaluated using, for example, an 80-point Legendre quadrature.

[0076] At block 406, the computing device determines, for each step of the second spatial grid, a second contribution to the fluence spectrum from elastic interactions of ions with atoms heavier than hydrogen. Different ions interacting with different atoms may undergo many different reaction types. For simplicity of description, cases where the ions are protons are considered. Calculations for other ions can follow the same procedure, with the difference being in the number and types of particles produced in the interaction and in probability distributions of their initial energies and directions of travel. In an example calculation involving collisions between protons and heavier atoms, the operations of block 406 are performed in a manner similar to the operations of block 404, except that the proton-to-heavyPATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT atom scattering involves only one scattered particle (the proton) instead of two, and the scattering of the primary incident proton is highly anisotropic. In this case, the displacement of the recoil atom, A, can be disregarded as negligible.

[0077] For incident ions heavier than protons, the displacement of the recoil atom can be, in some cases, significant. In that case, an algorithm similar to that for protons colliding with hydrogen atoms can be used, with additional computations to account for the number and types of particles produced in the interaction and the probability distributions of the ions’ initial energies. These differences can be fully accounted for by using the appropriate distributions when the integration of equation (46) is performed.

[0078] Continuing the example above, the example calculation for incident protons willbe given. The differential cross section for targets with an atomic mass number ^^ ^ 62 can begiven by: ^.^ଷ ^ ^.^^^ ^.ଷଷ^^ exp െ14.4^^^^ ^ 1.4^^ exp^െ10^^^ ^46^where ^^ is the invariant momentum transfer in GeV / c2. The parameter s is calculated using the sequence: ^^^^^ ൌ^47^^^௧^௧ ൌ ^^^ ^ ^^ ^ ^^ଶ ^48^in which ^^^^^and ^^^^are proton momenta in the laboratory and center of mass frames, respectively, in GeV / c; ^^^,^^ଶare the rest energies of the proton and an atom ^^ (e.g., for an incident proton,ൌ 0.938247 GeV, and ^^ଶൌ 14.999^^^for oxygen); and ^^ is proton kinetic energy in the laboratory frame.

[0079] For incident ions heavier than protons, similar computations can be used, with the substitution of the mass of the ion and the appropriate differential cross section. Moreover, in some reactions more than one particle can emerge from the collision. In such cases the procedure can be extended to account for the multiple particles.PATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT

[0080] The contribution to the fluence of scattered ions after elastic scattering with heavy atoms can be calculated using:This integral over ^^ is a line integral that can be calculated stepwise using steps Δt as determined in block 402. The integral over µcm can be calculated using a quadrature optimized for a differential cross section of the particular ion type. For example, for incident protons, a 40-point Laguerre quadrature can be used. To perform this integration, the cosine of the scattering angle ^^′^^^and the kinetic energy of a proton after scattering ^^′^^^^^^, for a given ^^^^can be calculated. These calculations can be performed using relativistic kinematics, following the steps described below.

[0081] For example, the total energy and relativistic momenta of the incident ion in the laboratory and center of mass frames, ^^^^^and ^^^^, are given by equations (47)-(50). In examples where the incident ion travels parallel to the z-axis,PATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT^^′^^^^^^ ൌ ^^ ′^^^^^^,^^^^^^ െ ^^1 ^59^

[0082] At block 408, the computing device determines whether iteration over the second spatial grid is complete. If the particular step is the last step, then execution of process 400 can terminate. Otherwise, the calculation described above continues at block 402. As the elements of process 400 may be performed following the determination of the fluence spectra as described above with reference to FIG. 3, block 402 may be performed following the performance of the next iteration of the steps of process 300.

[0083] Process 400 including corrections for nuclear scattering to the calculated fluence spectra omits corrections due to inelastic scattering of incident ions with atoms. For example, the cross sections for inelastic scattering between protons and hydrogen is zero for proton energies below 300 MeV. For heavier target atoms, after an inelastic collision between an incident ion and the heavier atom, the nucleus may emit photons, neutrons, or low energy protons. In some implementations, the contributions by these effects to the fluence spectra may be omitted as the contribution to the total dose may be negligible (e.g., on the order of ∼0.1%). Where this is the case, omitting corrections due to inelastic scattering of incident ions with atoms thus helps to optimize performance without a noticeable loss of accuracy. In cases where contributions from inelastic nuclear interactions are not negligible, computing contributions from inelastic nuclear interactions can be implemented using an algorithm similar to that for elastic nuclear interactions described above. Ion Beam Treatment

[0084] The methods of calculating dose distribution described above with reference to FIGs.3 and 4 can used as part of a treatment planning process. FIG.5 shows a process 500 for generating a treatment plan for ion beam therapy using the deterministic methods for dose distribution calculation described above. In some examples, process 500 can be implemented in a computing device used for radiation treatment planning. Process 500 may be performed using a computing device, for example, by a health care provider generating a treatment plan for ion beam therapy for a patient with a tumor. Process 500 includes selection of an ion to use for irradiation, determination of a number of ion beams and their incident ion energies, directions, and so on. Then, the techniques of the present disclosure can be used to determinePATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT the absorbed dose from the beams at various surface and interior locations of a particular tumor geometry, as well as the surrounding regions including organs at risk.

[0085] At block 502, a target volume is defined. The target volume is typically a tumor on or within the body of a human patient, experimental animal, training phantoms, and so on. The target volume may be located at a particular depth within the patient. Its location and shape can be determined from a computed tomography (CT) image, sometimes referred to as a planning CT image. For example, a qualified radiation oncology specialist may draw tumor boundaries on the CT image. For more accurate delineation of the target volume, additional imaging modalities can be used, such as magnetic resonance imaging (MRI) or positron emission tomography (PET).

[0086] For precise positioning of a patient for treatment, an optimized location of the isocenter (the point toward which all beams will be directed) is identified (e.g., marked on CT image), typically within the tumor. The patient will be positioned for treatment so that the marked isocenter precisely coincides with the location of the physical isocenter 178 of the ion therapy system 100. For example, the patient may be imaged in the treatment position using x- ray or CT imaging techniques and then necessary adjustments can be made.

[0087] Incident ion beams will travel through the patient’s body before they reach the tumor and after they exit it. Thus, in addition to identification of the target volume (e.g., the tumor) the shapes and locations of all organs at risk and the complete anatomy of the volume surrounding the target volume is also determined. This additional anatomical information can be added, for example, to the planning CT image. The planning CT image can be used to determine an optimized ion beam plan.

[0088] The irradiation of the tumor may proceed with a series of ion beams applied with varying directions and initial energies. Additionally, the ion itself may be selected as part of treatment planning. Examples of ion beams include proton beams, helium ion beams, carbon ion beams, or beams consisting of other heavy ions. Different ions may have varying penetration depths or Bragg peaks, relative biological effectiveness, scattering characteristics, availability, cost, and so on. Selection of the ion may be based on factors including the specific clinical scenario, desired treatment precision, biological effects, and available resources. Applications of the methods disclosed herein to heavier ions, as opposed to protons, may involve replacement of the proton mass and charge in certain equations by those of carbon or helium, or the like. In examples of an ion beam apparatus capable of providing different types of ion beams, the ion beam selected should match the type of ion assumed during developmentPATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT of the treatment plan. In some examples, an ion beam apparatus may be selected based on the desired type of ion beam.

[0089] At block 504, a number of ion beams and one or more characteristics of each ion beam on the ion beam therapy treatment plan are selected. The one or more characteristics may include beam energies, directions, modulation, and so on. As described above with reference to processes 300 and 400, the incident or initial ion beam energy is a parameter used for stepwise calculation of the fluence spectra as a function of depth (i.e., steps) as well as the angular and radial distributions for a particular ion beam. These distributions are then used to determine a predicted absorbed dose. Thus, the initial ion beam energies, directions, modulations, etc. are selected in accordance with the desired dose according to the patient’s anatomy and specifically the position and dimensions of the target volume and organs at risk.

[0090] For example, configured ion beams may be the first in a series of ion beam applications planned for a course of treatment. For instance, it may be desired to apply the ion beam to the tumor from a diversity of incident directions sequentially or simultaneously, which can spread the dose over a larger volume while still maintaining a high dose to the target volume, since all incident beams meet at the isocenter 178 within the target volume (e.g., tumor). The use of multiple incident ion beam directions can also be chosen so to avoid directly hitting critical organs. Accordingly, treatment plans may include incident beams from a diversity of directions. The direction selected at block 504 may thus be one of a number of planned applications of the ion beam to the target volume. For instance, ion energy and beam direction may be iteratively selected by a human operator with computer assistance or automatically by executing optimization software to obtain the desired dose to the target volume while keeping exposure of surrounding healthy tissues within acceptable limits.

[0091] At block 506, a dose distribution to the target volume and the surrounding region is computed. For example, process 300 can be used to calculate the dose distribution based on the fluence spectra for a given incident ion beam kinetic energy and direction. The nuclear corrections of process 400 can be applied to improve the accuracy of the computed dose for the ion beam. If multiple incident beams are planned (as is usually the case), the calculated dose distribution can be determined, e.g., as a sum of dose distributions from the different beams. The dose distribution may be displayed on a suitable graphical user interface (GUI) for inspection by health care providers engaged in treatment planning. For instance, the dose distribution may be represented as a two- or three-dimensional graph or as a region on aPATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT three-dimensional display. The dose distribution may similarly be represented numerically or using any other suitable representation.

[0092] At block 508, the dose distribution is compared to a dose prescription. For example, a health care provider such as a radiation oncologist may determine an appropriate dose to damage or destroy the tumor according to certain known criteria. The prescription dose can accordingly be used as a reference for iteratively adjusting the parameters such as number of beams, incident ion energies, directions, and the duration of time that the beam(s) remain on. In some examples, the dose may be cumulative. In that case, a dose from each direction and / or energy may be summed or otherwise combined to attain the total dose. In addition, the radiation oncologist, or other provider, will consider dose that may be delivered to surrounding healthy or at-risk organs to avoid complications stemming from incident ion beam radiation.

[0093] At block 510, the dose distribution to the surrounding region is compared to predefined limits. For example, published limits on safe doses to various regions, organs, and tissues may be used as the predefined limits. A computing device may be used to highlight or otherwise warn concerning instances where the calculated dose distribution exceeds the predefined limits. In general, efforts should be made to made to minimize the dose to the anatomy and organs surrounding the target volume.

[0094] At block 512, if the dose distribution is incommensurate with the dose prescription or the dose distribution to the surrounding region exceeds the predefined limits, the process continues at block 504. For example, the dose distribution to the target volume may locally or regionally too high, too low, or inconsistent with the prescription. Likewise, the dose distribution to the surrounding region may exceed the predefined limits, which can result in damage to healthy tissues or organs. In these cases, the process 500 is continued at block 504, iteratively determining dose distributions at various incident energies, directions, and durations (among other parameters) until a suitable combination of parameters for the particular segment of the treatment plan is determined. Treatment planning can continue until an acceptable or optimized set of incident energies, directions, and other parameters are determined, and at block 514, the treatment plan can be finalized, which can include, for instance, writing a set of control parameters for an ion beam apparatus to an appropriately-formatted data file. The control parameters can be used by a computing device to administer the ion beam therapy.

[0095] Process 500 can be used to plan treatment before treatment begins or at a later date (at any point before the treatment is completed). In some embodiments, planning can occur before or between treatments, e.g., in 4D planning or adaptive planning. For instance, ifPATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT prior to executing a treatment plan, it is observed that a tumor moves outside of certain predefined limits or safety-related margins (e.g., due to breathing), a 4D CT image can be used to perform 4D planning to account for motion of the tumor and update the treatment plan accordingly. Treatment is often conducted over several sessions (e.g., over a period of days or weeks), and process 500 can be used to plan each session.

[0096] FIG.6 shows a process 600 for that can be used for administering a treatment plan generated using deterministic treatment planning for ion beam therapy. In some examples, process 600 can be implemented in a computing device used for radiation treatment planning. Process 600 includes steps for calculation and refinement of the fluence spectra as described in, for example, processes 300 and 400.

[0097] At block 602, the patient is positioned relative to an ion beam apparatus 100, so that the location of the isocenter in a target volume marked on a planning CT image, precisely corresponds to the location of the physical isocenter 168 of the ion treatment device 100. In some implementations, real-time x-ray or CT imaging can be used to further facilitate accurate positioning.

[0098] For instance, the target volume may be a tumor. Accurate positioning of the patient with respect to the ion beam apparatus 100 ensures that the tumor aligns accurately with the treatment plan’s predetermined specifications. Irradiation of the tumor with high precision and accuracy are highly relevant to successful treatment and minimization of collateral damage to surrounding healthy organs and other anatomy. In some cases, specialized imaging tools or markers may be used to refine the patient positioning.

[0099] At block 604, a treatment plan file is loaded. The treatment plan file is a computer file that includes detailed instructions for delivery of the treatment plan as determined in process 600 by the ion beam therapy device 100. These instructions can be executed automatically by a computing device that controls operation of the ion beam therapy device 100 to deliver each ion beam, sequentially or simultaneously. Treatment delivery may be initiated by a human operator, such as a radiation oncology technician, after the patient is positioned for treatment and the position is verified using, for example, by x-ray imaging, as described in process 500 above.

[0100] At block 606, the ion beam is administered to the patient based on the ion beam treatment plan determined in process 500. Administration of the ion beam is guided by the ion beam therapy treatment plan, which includes parameters such as the direction, content, energy, intensity, duration, and modulation of the ion beam. In some examples, the ion beam energy orPATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT intensity can be adjusted in real-time using beam modulation techniques, following instructions contained in the treatment plan. Feedback mechanisms may be employed to ensure that all parameters of the delivered beam are the same as those specified in the treatment plan. For example, integrated imaging systems within the ion beam apparatus can provide real-time feedback to ensuring that the ion beam remains aligned with the target volume.

[0101] Some examples may include modules for 4D ion beam administration which can account for temporal changes and movements, such as respiratory motion during beam delivery, to ensure that the target volume receives the prescribed dose even as it moves in real- time. Likewise, some examples can include adaptive planning modules. Adaptive planning may allow for adjustments to the treatment plan during the administration of course of a treatment based on data received, for example, during daily or weekly checks, such as anatomical changes or tumor shrinkage. Adaptive planning can thus be used to adjust the ion beam treatment plan during the course of administration. Computer System Implementation

[0102] Processes described herein can be implemented in computer systems of various designs. FIG.7 shows a simplified block diagram of a computer system 700 suitable for use in some examples of the present disclosure. Computer system 700 includes a number of different subsystems interconnected via a system bus 775. The core subsystems include an input / output (I / O) controller 771, a system memory 772 (e.g., DRAM, SRAM, PROM, and / or other computer-readable media), and a central processor 773. Central processor 773, which can be implemented using one or more programmable integrated circuits (including single-core and / or multi-core microprocessors) controls operations of computer system 700 by executing program code that can be stored (at least temporarily) in system memory 772. Accordingly, central processor 773 can communicate with each subsystem and can control the execution of instructions from system memory 772 or storage device(s) 779, as well as the exchange of information between subsystems. Similarly, any of the data mentioned herein can be delivered from one component to another component and can be output to (or input from) the user. In some examples, central processor 773 may be coupled to one or more coprocessors, such as one or more graphics processing units (not shown) that are designed for high-throughput parallel processing.

[0103] I / O controller 771 allows other components to be communicatively coupled to central processor 773, and central processor 773 can receive input from other componentsPATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT and / or send output to other components via I / O controller 771. Accordingly, additional subsystems such as printer 774; user input device(s) 778 (e.g., keyboard, mouse, etc.); storage device(s) 779 (e.g., various computer-readable media such as hard disk drives or other fixed storage devices, removable disks, removable solid-state memory devices such as USB thumb drives, etc.); monitor 776, which is coupled to display adapter 782; and the like may be communicably coupled to central processor 773. Peripherals and I / O devices, which may couple to I / O controller 771, can be connected to the computer system using various interconnect standards known in the art, such as serial port 777. Wireless local-area connectivity (e.g., via Bluetooth or Wi-Fi or the like) may also be supported.

[0104] In some examples, network interface 781 may be provided to enable communication between computer system 700 and other computer systems, e.g., via Ethernet, Wi-Fi, or the like. Network interface 781 may support connection to a local area network and / or to a wide-area network such as the internet. Thus, for example, processes 300, 400, 500, 600 and other processes described herein can be implemented in one instance of computer system 700, which can communicate treatment plans to another instance of computer system 700 local to ion beam treatment system 100.

[0105] In some examples, computer system 700 is implemented as a single computer apparatus with some or all of the subsystems described above. In some examples, a single instance of computer system 700 can include multiple instances of the same components or subsystems, e.g., connected together by an internal interface. In some examples, two or more instances of computer system 700 (which can be configured alike or differently as desired) can communicate over a network. In such examples, one instance can be considered a client and another instance a server.

[0106] Various features described herein, e.g., methods, apparatus, computer-readable media and the like, can be realized using any combination of dedicated components and / or programmable processors and / or other programmable devices. The various processes described herein can be implemented on the same processor or different processors in any combination. Where components are described as being configured to perform certain operations, such configuration can be accomplished, e.g., by designing or connecting electronic circuits to perform the operation, by programming programmable electronic circuits (such as microprocessors) to perform the operation, or any combination thereof. Further, while the examples described above may make reference to specific hardware and software components, those skilled in the art will appreciate that different combinations ofPATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT hardware and / or software components may also be used and that particular operations described as being implemented in hardware might also be implemented in software or vice versa.

[0107] Any of the software-implemented components or functions described in this application may be realized in the form of software code to be executed by a processor; such code may be created using any suitable computer language such as, for example, Java, C++ or Perl using, for example, conventional or object-oriented techniques. The software code may be stored as a series of instructions or commands on a computer readable storage medium; suitable media include random access memory (RAM), a read only memory (ROM), a magnetic medium such as a hard-drive or a floppy disk, or an optical medium such as a compact disk (CD) or DVD (digital versatile disk), flash memory, and the like. The computer readable medium may also be a combination of multiple such media. Computer readable storage media encoded with the program code may be packaged with a compatible device or provided separately from other devices (e.g., as a separately packaged computer readable storage medium or via an internet download operation that results in the program code being stored on a computer readable storage medium of the device that downloaded it). Any such computer readable storage medium may reside on or within a single computer product (e.g., a hard drive, a CD, or an entire computer system) and may be present on or within different computer products within a system or network. Such programs may also be encoded and transmitted using carrier signals adapted for transmission via wired, optical, and / or wireless networks conforming to a variety of protocols, including the Internet. (It is noted that “storage” of programs or data is distinct from propagation of programs or data using transitory media such as carrier waves.)

[0108] Any of the methods described herein may be totally or partially performed with a computer system including one or more processors, which can be configured to perform the steps or operations. Thus, examples of the present disclosure can include computer systems configured to perform the steps or operations of any of the methods described herein, potentially with different components performing different steps or operations (or different groups of steps or operations).

[0109] Embodiments may be implemented by using a computer program product, comprising computer program / instructions which, when executed by a processor, cause the processor to perform any of the methods described in the disclosure.PATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT

[0110] The foregoing description of some examples has been presented only for the purpose of illustration and description and is not intended to be exhaustive or to limit the disclosure to the precise forms disclosed. Numerous modifications and adaptations thereof will be apparent to those skilled in the art without departing from the spirit and scope of the disclosure.

[0111] Reference herein to an example or implementation means that a particular feature, structure, operation, or other characteristic described in connection with the example may be included in at least one implementation. The disclosure is not restricted to the particular examples or implementations described as such. The appearance of the phrases “in one example,” “in an example,” “in one implementation,” or “in an implementation,” or variations of the same in various places in the specification does not necessarily refer to the same example or implementation. Any particular feature, structure, operation, or other characteristic described in this specification in relation to one example or implementation may be combined with other features, structures, operations, or other characteristics described in respect of any other example or implementation.

[0112] Use herein of the word “or” is intended to cover inclusive and exclusive OR conditions. In other words, A or B or C includes any or all of the following alternative combinations as appropriate for a particular usage: A alone; B alone; C alone; A and B only; A and C only; B and C only; and A and B and C. EXAMPLES

[0113] These illustrative examples are mentioned not to limit or define the scope of this disclosure, but rather to provide examples to aid understanding thereof. Illustrative examples are discussed above in the Detailed Description, which provides further description. Advantages offered by various examples may be further understood by examining this specification.

[0114] As used below, any reference to a series of examples is to be understood as a reference to each of those examples disjunctively (e.g., "Examples 1-4" is to be understood as "Examples 1, 2, 3, or 4").

[0115] Example 1 is a method for ion beam treatment planning, the method comprising: receiving information about a target volume and surrounding anatomy associated with a patient, the target volume and one or more elements of the surrounding anatomy each characterized by a depth within the patient and a set of dimensions; receiving parameters forPATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT an ion beam, the parameters including an initial ion beam energy, beam intensity, beam modulation and beam duration; defining a first spatial grid comprising a first plurality of discrete steps, each step having a step size; computing, for each step of the first spatial grid: a fluence spectrum based on a first iterative solution to the Lagrangian form of the Boltzmann transport equation; and a radial distribution, wherein computing the radial distribution includes: computing an angular distribution based on a second iterative solution to the Lagrangian form of the Boltzmann transport equation; and converting the angular distribution to the radial distribution; determining an absorbed dose distribution to the target volume and the surrounding anatomy based on the fluence spectra and the radial distributions for the steps of the first spatial grid according to the depth and dimensions of the target volume and the surrounding anatomy; and generating an ion beam therapy treatment plan based at least in part on the absorbed dose distribution to the target volume and the surrounding anatomy.

[0116] Example 2 is the method of example(s) 1, wherein the ion is a proton, a helium ion, or a carbon ion.

[0117] Example 3 is the method of one or more of the preceding examples, wherein the target volume is a tumor.

[0118] Example 4 is the method of one or more of the preceding examples, further comprising: positioning the patient relative to an ion beam apparatus, wherein the patient is positioned to place the physical isocenter of the ion beam apparatus within the target volume; configuring the ion beam apparatus based on the ion beam therapy treatment plan, including: determining a number of ion beams; and for each ion beam, determining one or more characteristics based on the ion beam therapy treatment plan; and administering the ion beam to the patient based on the ion beam therapy treatment plan.

[0119] Example 5 is the method of one or more of the preceding examples, wherein the one or more characteristics include at least one of ion beam energy, ion beam collimation, ion beam modulation, or ion beam duration.

[0120] Example 6 is the method of one or more of the preceding examples, wherein the first iterative solution to the Lagrangian form of the Boltzmann transport equation for a particular step comprises integrating a distribution that characterizes a slowing of an ion due to multiple Coulomb scattering during the particular step.

[0121] Example 7 is the method of one or more of the preceding examples, wherein the distribution is the Vavilov distribution.PATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT

[0122] Example 8 is the method of one or more of the preceding examples, wherein the distribution is a normal distribution.

[0123] Example 9 is the method of one or more of the preceding examples, wherein the second iterative solution to the Lagrangian form of the Boltzmann transport equation for a particular step comprises integrating a distribution that characterizes the lateral spread of ions as the ions pass through matter due to multiple Coulomb scattering during the particular step.

[0124] Example 10 is the method of one or more of the preceding examples, wherein the distribution is the Molière distribution.

[0125] Example 11 is the method of one or more of the preceding examples, further comprising determining, for each step, a spatial distribution of beam fluence based on the radial distribution.

[0126] Example 12 is the method of one or more of the preceding examples, wherein the spatial distribution of beam fluence is a Gaussian pencil beam.

[0127] Example 13 is the method of one or more of the preceding examples, further comprising: determining, for each step of the first spatial grid: an attenuation factor based on a total nuclear cross section and the step size; and a corrected fluence spectrum based on the fluence spectrum and the attenuation factor.

[0128] Example 14 is the method of one or more of the preceding examples, further comprising: determining a second spatial grid, comprising a second plurality of discrete steps, each step comprising a step size; and determining, for each step of the second spatial grid: a first contribution to the fluence spectrum from elastic interactions of ions with hydrogen atoms; and a second contribution to the fluence spectrum from elastic interactions of ions with atoms heavier than hydrogen, wherein each ordered step of the second spatial grid is larger than a corresponding ordered step of the first spatial grid.

[0129] Example 15 is the method of one or more of the preceding examples, wherein determining the absorbed dose distribution based on the fluence spectra and the radial distributions comprises integrating, over all energies, the product of fluence spectra and an ionic stopping power of a material representative of the target volume and the surrounding anatomy, and a density of the material.

[0130] Example 16 is the method of one or more of the preceding examples, wherein determining the first spatial grid comprises: determining a fractional energy loss per step; determining a minimum step size; and determining each step size using the fractional energyPATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT loss per step and an ionic range in a material representative of the target volume and the surrounding anatomy.

[0131] Example 17 is the method of one or more of the preceding examples, further comprising: determining the fluence spectrum between first consecutive steps of the first spatial grid using a first interpolation technique; and determining the radial distribution between second consecutive steps of the first spatial grid using a second interpolation technique.

[0132] Example 18 is the method of one or more of the preceding examples, further comprising: determining a voxelized representation of the target volume and the surrounding anatomy; determining that at least one step size is greater than a size of at least one pair of voxels, wherein the at least one pair of voxels comprise different materials; and determining corrected fluence spectra and corrected radial distributions based on the different materials.

[0133] Example 19 is a system comprising: one or more processors; and one or more computer-readable storage media storing instructions which, when executed by the one or more processors, cause the one or more processors to perform operations including: receiving information about a target volume and surrounding anatomy, the target volume and one or more elements of the surrounding anatomy each characterized by a depth within a patient and a set of dimensions; receiving parameters for an ion beam, the parameters including an initial ion beam energy, beam intensity, beam modulation and beam duration; defining a first spatial grid comprising a first plurality of discrete steps, each step having a step size; computing, for each step of the first spatial grid: a fluence spectrum based on a first iterative solution to the Lagrangian form of the Boltzmann transport equation; and a radial distribution, wherein determining the radial distribution includes: determining an angular distribution based on a second iterative solution to the Lagrangian form of the Boltzmann transport equation; and converting the angular distribution to the radial distribution; determining an absorbed dose distribution to the target volume and the surrounding anatomy based on the fluence spectra and the radial distributions for the steps of the first spatial grid according to the depth and dimensions of the target volume and the surrounding anatomy; generating an ion beam therapy treatment plan based at least in part on the absorbed dose distribution; determining a number of ion beams; and for each ion beam, determining one or more characteristics based on the ion beam therapy treatment plan.

[0134] Example 20 is a non-transitory computer-readable medium storing instructions that, when executed by one or more processors, cause the one or more processors to performPATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT operations including: receiving information about a target volume and surrounding anatomy, the target volume and one or more elements of the surrounding anatomy each characterized by a depth within a patient and a set of dimensions; receiving parameters for an ion beam, the parameters including an initial ion beam energy, beam intensity, beam modulation and beam duration; defining a first spatial grid comprising a first plurality of discrete steps, each step having a step size; computing, for each step of the first spatial grid: a fluence spectrum based on a first iterative solution to the Lagrangian form of the Boltzmann transport equation; and a radial distribution, wherein determining the radial distribution includes: determining an angular distribution based on a second iterative solution to the Lagrangian form of the Boltzmann transport equation; and converting the angular distribution to the radial distribution; determining an absorbed dose distribution to the target volume and the surrounding anatomy based on the fluence spectra and the radial distributions for the steps of the first spatial grid according to the depth and dimensions of the target volume and the surrounding anatomy; generating an ion beam therapy treatment plan based at least in part on the absorbed dose distribution; determining a number of ion beams; and for each ion beam, determining one or more characteristics based on the ion beam therapy treatment plan.

[0135] Example 21 is an apparatus, comprising: means for implementing the operations of the method of any of example(s)s 1-18.

[0136] Example 22 is a computer program product comprising computer instructions that, when executed by a processor, implement the operations of the method of any of example(s)s 1-18.

Claims

PATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT CLAIMS That which is claimed is:

1. A method for ion beam treatment planning, the method comprising: receiving information about a target volume and surrounding anatomy associated with a patient, the target volume and one or more elements of the surrounding anatomy each characterized by a depth within the patient and a set of dimensions; receiving parameters for an ion beam, the parameters including an initial ion beam energy, beam intensity, beam modulation and beam duration; defining a first spatial grid comprising a first plurality of discrete steps, each step having a step size; computing, for each step of the first spatial grid: a fluence spectrum based on a first iterative solution to the Lagrangian form of the Boltzmann transport equation; and a radial distribution, wherein computing the radial distribution includes: computing an angular distribution based on a second iterative solution to the Lagrangian form of the Boltzmann transport equation; and converting the angular distribution to the radial distribution; determining an absorbed dose distribution to the target volume and the surrounding anatomy based on the fluence spectra and the radial distributions for the steps of the first spatial grid according to the depth and dimensions of the target volume and the surrounding anatomy; and generating an ion beam therapy treatment plan based at least in part on the absorbed dose distribution to the target volume and the surrounding anatomy.

2. The method of claim 1, wherein the ion is a proton, a helium ion, or a carbon ion.

3. The method of claim 1, wherein the target volume is a tumor.

4. The method of claim 1, further comprising: positioning the patient relative to an ion beam apparatus, wherein the patient is positioned to place the physical isocenter of the ion beam apparatus within the target volume;PATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT configuring the ion beam apparatus based on the ion beam therapy treatment plan, including: determining a number of ion beams; and for each ion beam, determining one or more characteristics based on the ion beam therapy treatment plan; and administering the ion beam to the patient based on the ion beam therapy treatment plan.

5. The method of claim 4, wherein the one or more characteristics include at least one of ion beam energy, ion beam collimation, ion beam modulation, or ion beam duration.

6. The method of claim 1, wherein the first iterative solution to the Lagrangian form of the Boltzmann transport equation for a particular step comprises integrating a distribution that characterizes a slowing of an ion due to multiple Coulomb scattering during the particular step.

7. The method of claim 6, wherein the distribution is the Vavilov distribution.

8. The method of claim 6, wherein the distribution is a normal distribution.

9. The method of claim 1, wherein the second iterative solution to the Lagrangian form of the Boltzmann transport equation for a particular step comprises integrating a distribution that characterizes the lateral spread of ions as the ions pass through matter due to multiple Coulomb scattering during the particular step.

10. The method of claim 9, wherein the distribution is the Molière distribution.

11. The method of claim 1, further comprising determining, for each step, a spatial distribution of beam fluence based on the radial distribution.

12. The method of claim 11, wherein the spatial distribution of beam fluence is a Gaussian pencil beam.

13. The method of claim 1, further comprising:PATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT determining, for each step of the first spatial grid: an attenuation factor based on a total nuclear cross section and the step size; and a corrected fluence spectrum based on the fluence spectrum and the attenuation factor.

14. The method of claim 1, further comprising: determining a second spatial grid, comprising a second plurality of discrete steps, each step comprising a step size; and determining, for each step of the second spatial grid: a first contribution to the fluence spectrum from elastic interactions of ions with hydrogen atoms; and a second contribution to the fluence spectrum from elastic interactions of ions with atoms heavier than hydrogen, wherein each ordered step of the second spatial grid is larger than a corresponding ordered step of the first spatial grid.

15. The method of claim 1, wherein determining the absorbed dose distribution based on the fluence spectra and the radial distributions comprises integrating, over all energies, the product of fluence spectra and an ionic stopping power of a material representative of the target volume and the surrounding anatomy, and a density of the material.

16. The method of claim 1, wherein determining the first spatial grid comprises: determining a fractional energy loss per step; determining a minimum step size; and determining each step size using the fractional energy loss per step and an ionic range in a material representative of the target volume and the surrounding anatomy.

17. The method of claim 1, further comprising: determining the fluence spectrum between first consecutive steps of the first spatial grid using a first interpolation technique; and determining the radial distribution between second consecutive steps of the first spatial grid using a second interpolation technique.PATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT 18. The method of claim 1, further comprising: determining a voxelized representation of the target volume and the surrounding anatomy; determining that at least one step size is greater than a size of at least one pair of voxels, wherein the at least one pair of voxels comprise different materials; and determining corrected fluence spectra and corrected radial distributions based on the different materials.

19. A system comprising: one or more processors; and one or more computer-readable storage media storing instructions which, when executed by the one or more processors, cause the one or more processors to perform operations including: receiving information about a target volume and surrounding anatomy, the target volume and one or more elements of the surrounding anatomy each characterized by a depth within a patient and a set of dimensions; receiving parameters for an ion beam, the parameters including an initial ion beam energy, beam intensity, beam modulation and beam duration; defining a first spatial grid comprising a first plurality of discrete steps, each step having a step size; computing, for each step of the first spatial grid: a fluence spectrum based on a first iterative solution to the Lagrangian form of the Boltzmann transport equation; and a radial distribution, wherein determining the radial distribution includes: determining an angular distribution based on a second iterative solution to the Lagrangian form of the Boltzmann transport equation; and converting the angular distribution to the radial distribution; determining an absorbed dose distribution to the target volume and the surrounding anatomy based on the fluence spectra and the radial distributions for the steps of the first spatial grid according to the depth and dimensions of the target volume and the surrounding anatomy; generating an ion beam therapy treatment plan based at least in part on the absorbed dose distribution;PATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT determining a number of ion beams; and for each ion beam, determining one or more characteristics based on the ion beam therapy treatment plan.

20. A non-transitory computer-readable medium storing instructions that, when executed by one or more processors, cause the one or more processors to perform operations including: receiving information about a target volume and surrounding anatomy, the target volume and one or more elements of the surrounding anatomy each characterized by a depth within a patient and a set of dimensions; receiving parameters for an ion beam, the parameters including an initial ion beam energy, beam intensity, beam modulation and beam duration; defining a first spatial grid comprising a first plurality of discrete steps, each step having a step size; computing, for each step of the first spatial grid: a fluence spectrum based on a first iterative solution to the Lagrangian form of the Boltzmann transport equation; and a radial distribution, wherein determining the radial distribution includes: determining an angular distribution based on a second iterative solution to the Lagrangian form of the Boltzmann transport equation; and converting the angular distribution to the radial distribution; determining an absorbed dose distribution to the target volume and the surrounding anatomy based on the fluence spectra and the radial distributions for the steps of the first spatial grid according to the depth and dimensions of the target volume and the surrounding anatomy; generating an ion beam therapy treatment plan based at least in part on the absorbed dose distribution; determining a number of ion beams; and for each ion beam, determining one or more characteristics based on the ion beam therapy treatment plan.

21. An apparatus, comprising: means for implementing the operations of the method of any of claims 1-18.PATENT Attorney Docket No.: 090723-1472827 Client Reference No.: MDA23-045PCT 22. A computer program product comprising computer instructions that, when executed by a processor, implement the operations of the method of any of claims 1-18.

Images