Advertisement

Uncertainty quantification analysis and optimization for proton therapy beam lines

      Highlights

      • Uncertainties in proton therapy beam delivery.
      • Uncertainty Quantification methods in beam dynamics.
      • Analysis of errors in transport lines and their effects on the beam properties.
      • Optimized beam optics solutions using UQ methods.

      Abstract

      Since many years proton therapy is an effective treatment solution against deep-seated tumors. A precise quantification of sources of uncertainty in each proton therapy aspect (e.g. accelerator, beam lines, patient positioning, treatment planning) is of profound importance to increase the accuracy of the dose delivered to the patient. Together with Monte Carlo techniques, a new research field called Uncertainty Quantification (UQ) has been recently introduced to verify the robustness of the treatment planning. In this work we present the first application of UQ as a method to identify typical errors in the transport lines of a cyclotron-based proton therapy facility and analyze their impact on the properties of the therapeutic beams. We also demonstrate the potential of UQ methods in developing optimized beam optics solutions for high-dimensional problems. Sensitivity analysis and surrogate models offer a fast way to exclude unimportant parameters frcomplex optimization problems such as the design of a superconducting gantry performed at Paul Scherrer Institute in Switzerland.

      Keywords

      1. Introduction

      The conformal dose distribution that proton therapy can achieve is influenced by different uncertainties in both treatment planning and beam delivery. An underestimation of such uncertainties leads to more severe consequences (in terms of treatment quality) in proton therapy than in conventional photon therapy [
      • Paganetti H.
      Proton therapy physics.
      ].
      Imaging artifacts, tissue heterogeneity, patient positioning, intrafractional anatomical changes and organ motion are the main sources of uncertainty that compromise the quality of a treatment plan. These factors primarily impact the range of the proton beam with a reduced coverage of the tumor and higher dose to healthy tissues [
      • Paganetti H.
      Range uncertainties in proton therapy and the role of Monte Carlo simulations.
      ]. In beam delivery, uncertainties are mainly related to the instability of accelerator or transport lines, misplacement of magnets and errors in magnetic field strength. In consequence the expected beam properties (i.e. energy, intensity, position and size) at the isocenter are altered and the treatment quality is hence reduced. For example, in case of active spot scanning technique, the misplacement of a single pencil beam leads to a sensible deterioration of the dose homogeneity over the tumor as shown in Fig. 1.
      Figure thumbnail gr1
      Fig. 1Dose variation due to the misplacement of a single proton pencil beam (green dashed curve) with respect to the regular beam pattern (red curves). The violet dashed line represents the desired flat dose over the tumor. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
      A precise quantification of possible sources of uncertainty in each proton therapy aspect is of extreme importance to increase the robustness of the dose delivered to the patient [
      • Paganetti H.
      Proton therapy physics.
      ]. Several approaches can be adopted to mitigate the presence of uncertainties. In particular, after its initial use in the design phase, the development of a precise beam dynamics model allows a better understanding of the accelerator and beam transport system performances. Quality Assurance (QA) programs guarantee that the delivery of the dose fulfill rigorous medical protocols and recommendations [

      ICRU Bethesda. Clinical Proton Dosimetry Part I: beam production, beam delivery and measurement of absorbed dose. Technical Report ICRU-59, 1998.

      ]. Additionally, since several years, Monte Carlo codes are used to verify the robustness of the treatment plans [
      • Mairani A.
      • Böhlen T.T.
      • Schiavi A.
      • Tessonnier T.
      • Molinelli S.
      • Brons S.
      • Battistoni G.
      • Parodi K.
      • Patera V.
      A Monte Carlo-based treatment planning tool for proton therapy.
      ]. More recently a new research field called Uncertainty Quantification (UQ) has been also introduced to understand and manage uncertainties in the treatment planning [
      • Hartman C.M.H.
      Sensitivity and uncertainty analysis in proton therapy treatment plans using Polynomial Chaos Methods.
      ,
      • van der Voort S.R.
      Application of polynomial chaos in proton therapy. Dose distributions, treatment parameters, robustness recipes and treatment planning.
      ].
      In this work we profit from the general flexibility of UQ methods and apply them, for the first time, in the context of beam dynamics for proton therapy transport lines. A short introduction on UQ theory is presented in Section 2. We discuss a first application of UQ in Section 3 with the identification of typical errors in transport lines of a cyclotron-based proton therapy facility and the analysis of their impact on the properties of the therapeutic beams. Additionally, in Section 4, we demonstrate the potential of UQ in developing optimized beam optics solutions for high-dimensional problems such as the design of a superconducting gantry carried out at Paul Scherrer Institute in Switzerland.

      2. Uncertainty Quantification theory

      The outcomes of a numerical model are influenced by uncertainties in the formulation of the physical models, simplified assumptions, unknown effects and interactions among variables that are not properly considered.
      A well-established way to quantify the effect of uncertainties is provided by Monte Carlo (MC) sampling, where simulations of the original model are repeated several times with a random variation of the input parameters. These methods have a wide range of applications but tend to be rather slow especially in case of complex models.
      Recently a research field called Uncertainty Quantification (UQ) has been used together with MC techniques. Part of the so-called predictive science, UQ estimates the impact of uncertainty propagation through a numerical model or experiment using probability, statistics and applied mathematics. It has a quite broad range of applications like weather and climate modeling, nuclear reactor design, studies of biological systems [
      • Smith Ralph C.
      Uncertainty quantification. Theory, implementation, and applications.
      ] and recently also accelerator [
      • Adelmann A.
      On nonintrusive uncertainty quantification and surrogate model constructionin particle accelerator modeling.
      ] and medical physics [
      • Hartman C.M.H.
      Sensitivity and uncertainty analysis in proton therapy treatment plans using Polynomial Chaos Methods.
      ,
      • van der Voort S.R.
      Application of polynomial chaos in proton therapy. Dose distributions, treatment parameters, robustness recipes and treatment planning.
      ]. The UQ methods are divided in two categories: intrusive and non-intrusive. Intrusive methods require a reformulation, in stochastic way, of the original deterministic equations that describe a model. On the contrary, non-intrusive methods use existing codes considered as black-box. As for MC sampling, non-intrusive UQ requires repetitive runs of the original model for different input parameters.
      In this work we only focus on non-intrusive UQ methods, where the original numerical model (called high-fidelity model) is described as a set of Polynomial Chaos Expansion (PCE), the so-called surrogate model. The PCE surrogate model evaluates the response of a Quantity of Interest (QoI) to uncertainties in the input parameters. Following the derivation in [
      • Adelmann A.
      On nonintrusive uncertainty quantification and surrogate model constructionin particle accelerator modeling.
      ], a QoI u is described as a series expansion of basis vectors Ψi , function of input parameters ξ with corresponding coefficients αi as:
      u(ξ)=i=0αiΨi(ξ),
      (1)


      where ξ=(x1,x2,,xd) is the vector of input parameters xj (with j = 1, , d) that represent possible sources of uncertainty. Due to practical limitations, the PCE series in Eq. (1) is truncated up to a certain order p.

      2.1 Training points

      To build the surrogate model, repetitive runs of the high-fidelity model are performed at specific set ξ of input parameters referred as training points. Each input parameter xj varies in a defined domain between a minimum and maximum value. Within this variation range, the training points are selected from probability distributions, e.g. Uniform, Gaussian, Gamma etc. The number of training points depends on the order p of the surrogate model to be built and the number d of input parameters as:
      Ntrain=(p+1)d.
      (2)


      Depending on the complexity of the high-fidelity model, the Ntrain runs represent the most expensive part of the UQ analysis in terms of time and computation resources.

      2.2 PCE basis

      Being truncated at an order p, the PCE surrogate model contains only a limited number of basis vectors represented by a sequence of orthogonal polynomials. According to the Askey scheme [
      • Schoutens W.
      The askey scheme of orthogonal polynomials.
      ], several families of polynomials can be used as basis vectors depending on the type of probability distribution adopted for the input parameters as summarised in Table 1.
      Table 1Families of polynomials for different distributions of input parameters.
      DistributionPolynomialSupport
      GaussHermite[−,]
      UniformLegendre[−1,1]
      ExponentialLaguerre[0,]
      BetaJacobi[−1,1]
      In this work the PCE surrogate models are constructed from a uniform distribution of input parameters with Legendre polynomials as basis vectors. To construct the basis vectors (Ψi(ξ) in Eq. (1)) of a PCE at order p with d input parameters, it is necessary to introduce the multi-indices M = (m1,,md) of degree
      |M|=k=1dmkp.
      (3)


      Each multivariate basis vector Ψi(ξ) is the product of univariate polynomials as:
      Ψi(ξ)=k=1dψmkk(xk),
      (4)


      where ψmkk(xk) is, in our case, the univariate Legendre polynomial of degree mk associated to the input xk.

      2.3 Coefficients

      The number of coefficients to be evaluated depends on the order p of the PCE and on the number d of input parameters as:
      Ncoeff=(d+p)!d!p!.
      (5)


      Due to the orthogonality of the polynomial basis, the coefficients αi can be computed from the projection equation:
      αi=u(ξ),Ψi(ξ)Ψi(ξ),Ψi(ξ),
      (6)


      where the i-th coefficient is obtained by the orthogonal projection of u onto the i-th basis polynomial. The integral at nominator in Eq. (6) is evaluated numerically by quadrature rules [
      • Smith Ralph C.
      Uncertainty quantification. Theory, implementation, and applications.
      ]. At the denominator, the orthonormality of the univariate Legendre polynomials propagates to the multivariate ones:
      Ψi(ξ),Ψj(ξ)=δij,
      (7)


      where δij denotes the Kronecker delta.

      2.4 Statistical information

      The orthogonality of the basis vectors provides a complete statistical analysis of the surrogate model. For each QoI, mean E[·] and variance Var[·] are directly given from the coefficients αi as:
      E[u(ξ)]=α0
      (8)


      Var[u(ξ)]=i=1Ncoeffαi2
      (9)


      The confidence interval to 95% is calculated from Eq. (9) and t-Student distribution as:
      Conf.95%=±t·Var[u(ξ)]d,
      (10)


      where the t-Student parameter is adjusted with the number d of input parameters [
      • Taylor J.R.
      An introduction to error analysis.
      ].

      2.5 Sensitivity analysis

      Together with the surrogate model, UQ also provides a sensitivity analysis that identifies the contribution of each input parameter to the variability of a certain QoI. This analysis allows verifying the model robustness with respect to different inputs or determining whether the model can be simplified by fixing insensitive parameters [
      • Smith Ralph C.
      Uncertainty quantification. Theory, implementation, and applications.
      ].
      The sensitivity analysis is based on the Sobol’ indices that are defined from the variance of the surrogate model as:
      Si[%]=100Var[u(ξ)]i=0,iIisNcoeffαi2,
      (11)


      where Iis is the set of multi-indices mk that include i-th order only. In principle a large value of the index Si means that the associated parameter xi is important in the variation of a certain QoI.

      2.6 Setup of UQ analysis

      The analyses in this work are performed using the UQ Toolkit (UQTk) developed at the Sandia National Laboratories [].
      The development of an accurate high-fidelity model is the first step to exploit the full potential of non-intrusive UQ methods. For our application, an accurate and complete characterization of the beam properties along the transport lines in a cyclotron-based proton therapy facility requires a combination of multi-particle tracking codes with Monte Carlo simulations of particle-matter interaction. It has been shown in [
      • Rizzoglio V.
      • Adelmann A.
      • Baumgarten C.
      • Frey M.
      • Gerbershagen A.
      • Meer D.
      • Schippers J.M.
      Evolution of a beam dynamics model for the transport line in a proton therapy facility.
      ] that the OPAL (Object Oriented Parallel Accelerator Library) framework has such capability [

      Adelmann A, et al. The OPAL (Object Oriented Parallel Accelerator Library) Framework. Technical Report PSI-PR-08-02, Paul Scherrer Institut, 2008-2019.https://gitlab.psi.ch/OPAL/src/wikis/home.

      ]. For this reason, OPAL is used as black-box required by the non-intrusive UQ methods.
      Applying UQ techniques on a transport line means assuming the strength of magnetic fields, proton beam kinetic energy and transport element alignments as source of uncertainty. These factors constitute the input parameters (xj) of the UQ analysis, while the main QoIs are beam size, emittance, intensity and centroid position along the transport line and at the isocenter.
      After defining the input parameters and QoIs, the next step in UQ analysis is to choose the order p of the surrogate model. Depending on the number of input parameters, a low-order surrogate model (e.g. p = 2) has the advantage to require only a few training points to be developed. However, in case of complex large-scale models, this may return in an inaccurate surrogate model and wrong response of the QoI variation. An example of the inaccuracy of low-order surrogate model is shown in Fig. 2 in comparison with the expected results of the high-fidelity model evaluated at the training points.
      Figure thumbnail gr2
      Fig. 2Second-order surrogate model and confidence level of a generic QoI. The red dots represent the expected results of the high-fidelity model at the training points. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
      To properly interpret PCE plots, it has to be remarked that they show the different QoI outcomes for each training point. Since the training points have random inputs, there is no specific order of these points. Therefore, in the PCE plots these points have been ordered following, for example, an increasing value of QoI. To improve the precision of the PCE in Fig. 2, higher-order surrogate model, e.g. p = 4 is needed as displayed in Fig. 3.
      Figure thumbnail gr3
      Fig. 3Forth-order surrogate model and confidence level of a generic QoI. The red dots represent the expected results of the high-fidelity model at the training points. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
      Table 2 summarises the consequences of higher-order PCE in terms of increased number of training points and coefficients and therefore required computation resources.
      Table 2Number of coefficients and training points for low- and high-order surrogate model and different number of input parameters.
      PCE orderNtrainNcoeff
      d  = 3
      p = 22710
      p = 36420
      p = 412535
      d  = 5
      p = 224321
      p = 3102456
      p = 43125126

      3. UQ analysis for beam transport lines

      In this section non-intrusive UQ methods are applied on a beam transport line in the proton therapy facility at Paul Scherrer Institute (PSI) in Switzerland. Since 1984 the facility at PSI is treating patients with proton beams. Today the facility (PROSCAN) is equipped with a dedicated superconducting cyclotron (COMET) and consists of four clinical treatment rooms and one experimental area PIF (Proton Irradiation Facility). Three over the four treatment rooms are equipped with a rotating gantry, all supporting active spot-scanning technique. The forth treatment room is dedicated to eye-tumor treatment with an horizontal fixed beam line [,
      • Schippers J.M.
      • et al.
      The SC cyclotron and beam lines of PSI’s new protontherapy facility PROSCAN.
      ].
      In particular, the UQ analysis focuses on the beam transport lines (sketched in Fig. 4) that connects the cyclotron COMET to Gantry 3, that turned into clinical operation in 2018.
      Figure thumbnail gr4
      Fig. 4Beam transport line from COMET to Gantry 3. Named in blue the most important components of the beam line. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
      A 250 MeV proton beam is extracted from the superconducting cyclotron COMET and focused by a quadrupole triplet onto a multi-wedge graphite degrader. Based on the relative position of the wedges, the beam energies in the range 230–70 MeV are obtained. The growth of phase space and emittance is limited within the acceptance of the following transport line by a set of collimators [
      • Rizzoglio V.
      • Adelmann A.
      • Baumgarten C.
      • Meer D.
      • Snuverink J.
      • Talanov V.
      On the accuracy of Monte Carlo based beam dynamics models for the degrader in proton therapy facilities.
      ]. The first achromatic bending section that the beam encounters works as Energy Selection System (ESS). At the maximum dispersion peak, a horizontal slit is installed. The aperture of the slit determines the mean energy and spread of the beam that reaches the patient. Downstream of ESS, an energy-dependence beam intensity modulation is performed by means of the Intensity Compensation Collimator (ICC) [

      Koschik A, et al. Gantry 3: further development of the PSI PROSCAN proton therapy facility. In Proceedings of 6th International Particle Accelerator Conference, Richmond (USA); 2015. pp. 2275–2277.

      ]. The beam size after the second bending section is imaged at the coupling point to Gantry 3 through a magnetic telescope. For each energy in the range 230–70 MeV, the magnetic field strength of the transport line elements is optimized to achieve required beam properties (e.g. size, position and intensity) at the coupling point. At this location a symmetric double beam waist is needed to eliminate optics dependency on the gantry rotation. Additionally, a complete dispersion suppression is aimed to remove energy-dependent corrections on the beam position.
      Depending on the location, errors in the transport line setting lead to different undesired effects on the beam. Two examples are discussed in the next sections: variation in the quadrupole gradients in Section 3.1 and in the beam energy in Section 3.2.

      3.1 Variation in the quadrupole gradients

      A non-optimal quadrupole gradient leads to over-focusing/defocusing effects with a consequent variation of the beam size, symmetry and losses. To control such variation, the tolerance on the magnetic fields (ΔB/B) is of the order of 10-4 [

      Schippers JM. Beam transport systems for particle therapy. Notes for CERN specialised accelerator school on medical accelerators; 2016. pp. 321–333.

      ] and the optics layout has intermediate images between the quadrupole doublets.
      In this section we assume a variation on the gradients (K) of the last three quadrupoles (QMD10, QMD11 and QMD12) before the coupling point to Gantry 3. As starting condition, the high-fidelity OPAL model is developed for this last part of the transport line with the quadrupole setting optimized for 70 MeV proton beam. The 70 MeV beam envelope is quite broad before the coupling point (as shown in [

      Rizzoglio V. Precise Beam dynamics models for transport lines in a cyclotron-based proton therapy facility [Ph.D. thesis]. ETH Zürich, 2018. (25268).

      ]) and variation in the quadrupole gradients can lead to severe deterioration of the quality of the beam entering in the gantry.
      The UQ analysis focuses on the impact that the gradient variation has on the beam properties at the coupling point. Quantities of interest are the beam losses at the collimator installed at the coupling point and the symmetry factor in terms of ratio between the beam sizes in the two transverse planes. In the first test the variation domain of the three gradients is set by the tolerance level of 10-4. Table 3 summarizes the required parameters to build a third-order surrogate model in terms of training points (Eq. (2)) and coefficients (Eq. (5)).
      Table 3UQ parameters to build the surrogate models.
      ParameterValue
      Order of surrogate modelp = 3
      Number of inputsd = 3
      Input vectorξ=[KQMD10,KQMD11,KQMD12]
      Number of coefficientsα = 20
      Number of training points64
      The surrogate models for the two QoIs are shown in Fig. 5 with respect to a reference value from the high-fidelity model with the nominal 70 MeV quadrupole settings.
      Figure thumbnail gr5
      Fig. 5Surrogate models and confidence levels for two QoIs at the coupling point under quadrupole strength variation. The red dashed lines indicate the reference values for the nominal 70 MeV settings. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
      In Fig. 5a a deterioration of the beam symmetry is observed in the majority of the training points. The discrepancy with the reference value (1.017) is, for each training point, very small up to a maximum of 0.1%. For the same set of training points, the beam losses oscillate between 9.34% and 9.39% (Fig. 5b). A maximum increment of 0.14% is observed in the beam losses with respect to the reference value.
      The situation shown in Fig. 5 changes drastically if the quadrupole gradient variation goes outside the tolerance of 10-4. This scenario might happen in case of problems with power supplies or cooling system of the quadrupoles. To include also these aspects, a second surrogate model is developed assuming 1% variation domain of the three quadrupole gradients. The corresponding variation for beam symmetry and losses at the coupling point are displayed in Fig. 6.
      Figure thumbnail gr6
      Fig. 6Surrogate models and confidence levels for two QoIs at the coupling point under 1% quadrupole strength variation. The red dashed lines indicate the reference values for the nominal 70 MeV settings. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
      Here one can recognize that the beam losses and beam symmetry increase at the same training points, which indicates the correlation of beam losses with a beam size increase in a certain direction (i.e. beam symmetry >1).
      The analysis should be extended to all quadrupoles in the beam line for a more precise understanding of the effects of quadrupole gradient variation. However the results discussed above remain valid since the last three quadrupoles have direct influence on the beam quality at the coupling point. The surrogate models in Fig. 5 confirm how the strict tolerance of the magnetic field assures to maintain the required beam quality. The second surrogate models allow considering a wider spectrum of possible source of errors for quadrupole gradients, including also the variation within the required tolerance.

      3.2 Uncertainties in the beam kinetic energy

      A correct setup of the beam energy is of primary importance since it is related to the penetration depth into the patient as well as to the focussing of the quadrupoles, which are then set at the “wrong” strengths. At the PROSCAN facility the mean beam energy at the isocenter is set by the degrader and its final value and the energy spread are selected at the ESS. A small fraction (tenth of MeV) will be lost later, due to crossing of foils and air in the beam transport line [

      Rizzoglio V. Precise Beam dynamics models for transport lines in a cyclotron-based proton therapy facility [Ph.D. thesis]. ETH Zürich, 2018. (25268).

      ]. Uncertainties in the beam energy can arise from a wrong position of the degrader wedges, errors in the extracted beam due to the instability of the cyclotron or errors in the magnetic field of dipoles at ESS.
      The reference configuration for this analysis is the 230 MeV degrader setting where a variation in the beam energy is expected to have a bigger impact due to the small energy spread at ESS. As source of uncertainty we assume a variation in the thickness of the degrader that mimics errors in the wedge position and in the extracted beam properties from the cyclotron COMET. In fact the mean energy of the extracted beam can be influenced by small variation of the magnetic field and the energy spread by the multi-turn extraction. During the commissioning of the accelerator, both quantities have been measured and their stability verified within medical specification [
      • Schippers J.M.
      • et al.
      The SC cyclotron and beam lines of PSI’s new protontherapy facility PROSCAN.
      ,
      • Baumgarten C.
      • et al.
      Isochronism of the ACCEL 250 MeV medical proton cyclotron.
      ]. A summary of the input parameters with their variation domain used in this analysis is reported in Table 4.
      Table 4Input parameters for UQ analysis on energy error. The tolerances are from medical protocols and specifications.
      ParameterReferenceTolerance
      COMET energy (ECOMET)250.0 MeV±0.5 MeV
      COMET energy25 keV±250 keV
      spread (ΔECOMET)
      Half-wedge thick. (HWt)2.5 mm±0.2 mm
      For this analysis we setup a third-order surrogate model, whose number of training points and coefficients are the same as in Table 3. The aim is to understand the impact that the input parameters of Table 4 have on the transverse beam size at the coupling point to Gantry 3. The surrogate models for the two QoIs are shown in Fig. 7 where the beam size spread in the x-direction is clearly shown under energy errors. The x-direction coincides with the bending plane, so that dispersion (i.e. energy change) effects will contribute to the beam size in the x-plane. The spread in the y-beam size is mainly due to a wrong focussing of the quadrupoles, which are not matched to the energies taken.
      Figure thumbnail gr7
      Fig. 7Surrogate models and confidence levels of the transverse beam size at the coupling point under energy uncertainties. The black-dashed line indicates the reference value for the nominal 230 MeV setting.
      From the Sobol’ indices (Section 2.5), it is possible to identify which input parameter has the bigger influence on the variability of the QoIs. The sensitivity analysis associated with the surrogate models above is shown in Fig. 8, where the major contribution to the beam size variation arises from the degrader thickness and extracted beam energy from the cyclotron. The energy spread from COMET can be considered an insensitive parameter.
      Figure thumbnail gr8
      Fig. 8Sensitivity analysis: relative contribution of the input parameters to the two QoIs.

      4. UQ optimization in beam dynamics

      In this section we demonstrate the potential of UQ methods in developing optimized beam optics solutions for high-dimensional problems. The sensitivity analysis and surrogate models offer a fast way to exclude unimportant parameters from complex optimization problems such as the design of beam optics for a superconducting (SC) gantry project under development at PSI.

      4.1 Superconducting gantry project at PSI

      Research in proton therapy aims to reduce cost and size of accelerators and gantries and move towards single-room facilities. For this purpose the development of the superconducting technology in other research fields was of advantage. With SC magnets, the new generation gantries are expected to be lighter and more compact with an easier integration in local hospitals or existing cancer centers. Additionally, the possibility to increase the magnetic field strength and to use combined-function magnets opens the road for new treatment solutions [
      • Gerbershagen A.
      • Calzolaio C.
      • Meer D.
      • Sanfilippo S.
      • Schippers J.M.
      The advantages and challenges of superconducting magnets in particle therapy.
      ]. A major advantage of the new gantry generation is an increased momentum acceptance to avoid or, at least, reduce the time for dipole cycling. Normal conducting gantries have a momentum acceptance limited to ±0.5–1% Δp/p, defined by the aperture of ESS slit. Achieving larger momentum acceptance (10%) is one of the main potential offered by the SC gantries [
      • Nesteruk K.P.
      • Calzolaio C.
      • Meer D.
      • Rizzoglio V.
      • Seidel M.
      • Schippers J.M.
      Large energy acceptance gantry for proton therapy utilizing superconducting technology.
      ]. Since few years PSI researchers are developing the idea of a single-room proton therapy facility equipped with a 250 MeV cyclotron and a novel concept of isocentric gantry as sketched in Fig. 9.
      Figure thumbnail gr9
      Fig. 9Scheme of the single-room proton therapy facility proposed by PSI
      [
      • Gerbershagen A.
      • Meer D.
      • Schippers J.M.
      • Seidel M.
      A novel beam optics concept in a particle therapy gantry utilizing the advantages of superconducting magnets.
      ]
      .
      The PSI SC gantry (Fig. 9) consists of a degrader, series of collimators (to shape the phase space of the degraded beam) and a SC achromatic bending section. The transverse tumor scanning is performed by scanning magnets downstream of the SC bending section. The SC bending section is designed to deflect the proton beam by 135°. As shown in Fig. 10, it consists of three SC combined-function magnets: two 67.5° dipole magnets combined with quadrupole and sextupole components (SCD1 and SCD2) and a central quadrupole combined with sextupole component (SCQ). Outside of the SC bending section, two tunable normal conducting quadrupoles (Q1 and Q2) provide fine adjustments to the beam optics [

      Sanfilippo S, Calzolaio C, Anghel A, Gerbershagen A, Schippers JM. Conceptual design of superconducting combined-function magnets for the next generation of beam cancer therapy gantry. In Proceedings of 25th Russian Particle Accelerator Conference, St. Petersburg (Russia); 2016. pp. 138–140.

      ].
      Figure thumbnail gr10
      Fig. 10Projection of the magnetic field with reference particle trajectory along the mid-plane of the SC bending section. Magnet names are highlighted.
      In its preliminary phase, the project of the PSI SC gantry is mainly focused on the design of the bending section. The beam optics aims to an achromatic bending section with a complete suppression of the dispersion at the exit of the second dipole (SCD2) and a symmetric double waist of 5 mm FWHM at the isocenter. The gantry aims to satisfy a very large momentum acceptance from −11.5% to +13.5% Δp/p. The acceptance range is not large enough to cover all energies from 230 to 70 MeV and imposes to create two sub-energy ranges: 230–150 MeV and 150–70 MeV. Within each energy range the magnetic fields of the SC bending section are optimized for a reference energy only (185 MeV for the first setup and 115 MeV for the second). The other energies in each sub-range are transported through the gantry thanks to the large momentum acceptance [
      • Nesteruk K.P.
      • Calzolaio C.
      • Meer D.
      • Rizzoglio V.
      • Seidel M.
      • Schippers J.M.
      Large energy acceptance gantry for proton therapy utilizing superconducting technology.
      ,

      Sanfilippo S, Calzolaio C, Anghel A, Gerbershagen A, Schippers JM. Conceptual design of superconducting combined-function magnets for the next generation of beam cancer therapy gantry. In Proceedings of 25th Russian Particle Accelerator Conference, St. Petersburg (Russia); 2016. pp. 138–140.

      ]. The development of the beam dynamics model of the SC bending section requires a combination of three codes: COSY-Infinity for a preliminary lattice optimization [], OPERA-3D to design a realistic configuration of the SC magnetic field components [] and OPAL for a precise particle tracking through the 3D-magnetic field maps exported from OPERA-3D. The beam envelope at 185 MeV (reference energy) from the OPAL model, based on optimized magnetic field maps from OPERA-3D, is presented in Fig. 11 and benchmarked against COSY-Infinity.
      Figure thumbnail gr11
      Fig. 11Beam envelope of the SC bending section: comparison between OPAL (red and blue lines) and COSY-Infinity (green dots). On top the lattice is displayed. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
      The agreement between OPAL and COSY-Infinity envelope is acceptable up to the first SC dipole (SCD1) in both transverse planes. Downstream of SCD1 a deviation appears especially in the horizontal-bending plane, probably due the approximation of the magnetic field definition in COSY-Infinity. Different from OPAL, COSY-Infinity model does not use the same field maps exported from OPERA-3D, but only an approximation based on the field harmonics. Therefore, only a qualitative comparison between OPAL and COSY-Infinity is meaningful. However, to improve the optics model in OPAL and reproduce the results of COSY-Infinity, it is possible to adjust the relative position of five magnets (i.e. 3D field maps) and scale their magnetic fields. Each map is described by three input parameters: two coordinates (x, y) to place it in the lattice and a scaling factor of the magnetic field. Considering all five magnets, the total number of variables to be adjusted is 15. In this scenario, a fast evaluation of possible optics solutions can be obtained using UQ methods as described in Section 4.2.

      4.2 Beam envelope optimization with UQ

      With 15 input parameters, building a surrogate model requires millions of high-fidelity runs. However, looking at Fig. 11, the OPAL envelope deviates to COSY-Infinity model starting from the SC quadrupole (SCQ) till the isocenter. Therefore, only position shifts and field scaling of the last three magnets (SCQ, SCD2 and Q2) are considered in this analysis, reducing the number of input parameters to 9.
      At this point the sensitivity analysis can be used to identify which parameter and magnet field map contributes to improve the beam size at the isocenter. Around the reference configuration, small variation of the map position and field scaling is assumed as starting condition for the UQ analysis. The results of the sensitivity analysis are shown in Fig. 12 for the 3D map positions and field scaling.
      Figure thumbnail gr12
      Fig. 12Sensitivity analysis of the input parameters for the SC bending section.
      Fig. 12a suggests that the vertical position of SCQ (SCQy), horizontal and vertical positions of SCD2 (SCD2x and SCD2y) have a larger impact on the beam size at the isocenter. From Fig. 12b the field of SCD2 is the most influential, while surprisingly less important is the contribution of Q2. From the sensitivity analysis, 4 influential input parameters out of 9 are identified for the surrogate model: SCQy, SCD2x, SCD2y and SCD2B. Using Eqs. ((2), (5)), the parameters to build the surrogate model are summarised in Table 5.
      Table 5UQ parameters to build the surrogate models for the SC gantry optics.
      ParameterValue
      Order of surrogate modelp = 3
      Number of inputsd = 4
      Input vectorξ=[SCQy,SCD2x,SCD2y,SCD2B]
      α = 35
      Number of training points256
      Without previous knowledge, a larger variation of the input parameters is assumed to build the surrogate model. The resulting surrogate model is shown in Fig. 13 together with a band of accepted values of the beam size at the isocenter in both transverse planes.
      Figure thumbnail gr13
      Fig. 13Surrogate models for transverse beam size (σ) at the isocenter of the SC gantry. The violet band represents the accepted values. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
      The wide variation domain of the input parameters reduces the accuracy of the surrogate models with respect to the high-fidelity results. Fig. 14 shows a zoom of the complete surrogate models (Fig. 13) where the beam size values satisfy the requirements expressed by the violet band.
      Figure thumbnail gr14
      Fig. 14Zoom of around the surrogate model values within the accepted violet band. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
      Only two training points (IDs: 11 and 13) among 256 fit inside the accepted band. The beam sizes from the surrogate model and from the high-fidelity OPAL model are compared at the two training points as reported in Table 6.
      Table 6Beam size values (σ in mm) at the isocenter of the SC gantry for two selected training points.
      Train.Surrogate modelHigh-fidelity model
      point IDx-planey-planex-planey-plane
      112.382.411.724.88
      132.672.572.013.35
      The high-fidelity values for the training point number 13 seems promising to obtain an optimized optics solution. A small variation of the same input parameters (Table 5) around this configuration is performed and a second surrogate model is built, as shown in Fig. 15. In this second case a higher accuracy of the response of the surrogate model with respect to the high-fidelity runs is expected.
      Figure thumbnail gr15
      Fig. 15Surrogate models for a small variation around the training point n. 13 in . Variation of transverse beam size (σ) at the isocenter of the SC gantry.
      The SC gantry envelope, corresponding to the training point that satisfies the beam size requirements is shown in Fig. 16.
      Figure thumbnail gr16
      Fig. 16Optimized beam envelope from the surrogate models: comparison between OPAL and COSY-Infinity (green dots). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
      In this example, the sensitivity analysis results to be a fundamental instrument to remove unimportant parameters from the possible inputs. This application shows the potential of UQ analysis to optimize optics solutions, especially in high-dimensional problems.

      5. Conclusions

      In this work the UQ methods are applied, for the first time, on proton therapy beam lines and their flexibility is demonstrated in the two different applications. In Section 3, two common sources of uncertainty have been considered and the corresponding variations of the beam properties analysed. For proton therapy facilities already in operation UQ is a general and fast tool for error studies and the surrogate model is a powerful instrument to support machine learning. As further development, a more clear understanding of the transport system evolution could arise applying UQ techniques directly on measurements and not only on simulations. In the second application (Section 4) we have demonstrated the UQ potential to find optimized solution for the SC gantry optics. Dealing with 9 input variables, the sensitivity analysis helped removing unimportant parameters reducing the input variables to 4. The corresponding surrogate models provided an optimized optics solution of the beam parameters at the isocenter of the SC gantry. For the development of new proton therapy concepts, UQ reveals its potential for beam optics studies. In particular, unimportant parameters are excluded with the sensitivity analysis and improved solutions can be found running an optimizer on the surrogate model that is fast to evaluate.

      Acknowledgment

      The authors thank Dr. Ciro Calzolaio for providing the 3D magnetic field maps with OPERA used in the optimization of the SC gantry.

      References

        • Paganetti H.
        Proton therapy physics.
        CRC Press, Taylor & Francis Group2011
        • Paganetti H.
        Range uncertainties in proton therapy and the role of Monte Carlo simulations.
        Phys Med Biol. 2012; 57: R99-R117
      1. ICRU Bethesda. Clinical Proton Dosimetry Part I: beam production, beam delivery and measurement of absorbed dose. Technical Report ICRU-59, 1998.

        • Mairani A.
        • Böhlen T.T.
        • Schiavi A.
        • Tessonnier T.
        • Molinelli S.
        • Brons S.
        • Battistoni G.
        • Parodi K.
        • Patera V.
        A Monte Carlo-based treatment planning tool for proton therapy.
        Phys Med Biol. 2013; 58: 2471-2490
        • Hartman C.M.H.
        Sensitivity and uncertainty analysis in proton therapy treatment plans using Polynomial Chaos Methods.
        ([Master’s thesis]) Delft University of Technology, 2014
        • van der Voort S.R.
        Application of polynomial chaos in proton therapy. Dose distributions, treatment parameters, robustness recipes and treatment planning.
        ([Master’s thesis]) Delft University of Technology, 2015
        • Smith Ralph C.
        Uncertainty quantification. Theory, implementation, and applications.
        SIAM Computational Science Engineering, 2014
        • Adelmann A.
        On nonintrusive uncertainty quantification and surrogate model constructionin particle accelerator modeling.
        SIAM/ASA J Uncertainty Quantification. 2019; 7: 383-416
        • Schoutens W.
        The askey scheme of orthogonal polynomials.
        in: Stochastic processes and orthogonal polynomials. Lecture notes in statistics. vol. 146. Springer, New York, NY2000
        • Taylor J.R.
        An introduction to error analysis.
        University Science Books, California1996
      2. Uq toolkit.http://www.sandia.gov/uqtoolkit/.

        • Rizzoglio V.
        • Adelmann A.
        • Baumgarten C.
        • Frey M.
        • Gerbershagen A.
        • Meer D.
        • Schippers J.M.
        Evolution of a beam dynamics model for the transport line in a proton therapy facility.
        Phys Rev Accel Beams. 2017; 20124702
      3. Adelmann A, et al. The OPAL (Object Oriented Parallel Accelerator Library) Framework. Technical Report PSI-PR-08-02, Paul Scherrer Institut, 2008-2019.https://gitlab.psi.ch/OPAL/src/wikis/home.

      4. PROSCAN.https://www.psi.ch/protontherapy/center-for-proton-therapy-cpt.

        • Schippers J.M.
        • et al.
        The SC cyclotron and beam lines of PSI’s new protontherapy facility PROSCAN.
        Nucl Instr Meth B. 2007; 261: 773-776
        • Rizzoglio V.
        • Adelmann A.
        • Baumgarten C.
        • Meer D.
        • Snuverink J.
        • Talanov V.
        On the accuracy of Monte Carlo based beam dynamics models for the degrader in proton therapy facilities.
        Nucl Instr Meth A. 2018; 898: 1-10
      5. Koschik A, et al. Gantry 3: further development of the PSI PROSCAN proton therapy facility. In Proceedings of 6th International Particle Accelerator Conference, Richmond (USA); 2015. pp. 2275–2277.

      6. Schippers JM. Beam transport systems for particle therapy. Notes for CERN specialised accelerator school on medical accelerators; 2016. pp. 321–333.

      7. Rizzoglio V. Precise Beam dynamics models for transport lines in a cyclotron-based proton therapy facility [Ph.D. thesis]. ETH Zürich, 2018. (25268).

        • Baumgarten C.
        • et al.
        Isochronism of the ACCEL 250 MeV medical proton cyclotron.
        Nucl Instr Meth A. 2007; 570: 10-14
        • Gerbershagen A.
        • Calzolaio C.
        • Meer D.
        • Sanfilippo S.
        • Schippers J.M.
        The advantages and challenges of superconducting magnets in particle therapy.
        Supercond Sci Technol. 2016; 29083001
        • Nesteruk K.P.
        • Calzolaio C.
        • Meer D.
        • Rizzoglio V.
        • Seidel M.
        • Schippers J.M.
        Large energy acceptance gantry for proton therapy utilizing superconducting technology.
        Phys Med Biol. 2019; 64 (13): 175007
        • Gerbershagen A.
        • Meer D.
        • Schippers J.M.
        • Seidel M.
        A novel beam optics concept in a particle therapy gantry utilizing the advantages of superconducting magnets.
        Z Med Phys. 2016; 26: 224-237
      8. Sanfilippo S, Calzolaio C, Anghel A, Gerbershagen A, Schippers JM. Conceptual design of superconducting combined-function magnets for the next generation of beam cancer therapy gantry. In Proceedings of 25th Russian Particle Accelerator Conference, St. Petersburg (Russia); 2016. pp. 138–140.

      9. COSY-Infinity.http://www.bt.pa.msu.edu/index_cosy.htm.

      10. OPERA-3D.https://operafea.com/.