REVIEW 3 major objections 3 minor 3 references
Probabilistic Proton Treatment Planning: a novel approach for optimizing underdosage and overdosage probabilities of target and organ structures
T0 review · 3 major / 3 minor · reviewed 2026-08-06 · deepseek-v4-flash
Pith's one-line read This paper claims that proton treatment plans can be optimized to meet explicit per-voxel underdosage and overdosage probability targets by iterating a simple $E[d] \pm \delta \cdot SD[d]$ percentile surrogate, and that on phantom…
desk verdict Solid proof-of-principle for a probabilistic proton planning method that adapts per-voxel δ-factors to non-Gaussian dose distributions; the outer-loop convergence is unproven, but the phantom evidence and PCE checks are strong enough to justify a serious referee. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The $\delta$-factor is the signed number of standard deviations from the expected voxel dose to a target percentile, defined per voxel as $\delta_i^{\alpha} = (E[d_i] - d_i^{\alpha})/SD[d_i]$. It converts a probabilistic objective into the deterministic dose level $E[d_i] \pm \delta_i \cdot SD[d_i]$ used in quadratic penalties with analytical gradient and Hessian. The outer loop updates $\delta_i$ after each inner optimization using percentiles sampled from the PCE meta-model, with damping $\kappa = 0.2$ and a moving-average convergence check, so that the surrogate tracks the true percentile across iterations.
What would settle it
Take a geometry with a strongly non-Gaussian voxel dose distribution, such as a CTV voxel at the field edge under range uncertainty alone, run the outer loop to convergence, then evaluate the final beam weights by direct Monte Carlo sampling of $10^4$ error scenarios and count the fraction of scenarios with $d_i \le \gamma_i$. If that fraction exceeds the requested $\alpha$ for a non-negligible set of voxels, or if the $\delta$-sequence cycles instead of settling, the probability-control claim fails.
Extended reading notes
Core claim
The central claim is that voxel-wise dose percentiles can be controlled during treatment-plan optimization by writing them as $E[d_i] \pm \delta_i \cdot SD[d_i]$ and updating the per-voxel factors $\delta_i$ in an outer loop, using Polynomial Chaos Expansion to sample the uncertainty distribution cheaply. With this machinery, the inner optimization replaces a probabilistic goal such as $P(d_i \le \gamma_i) \le \alpha$ with a deterministic dose level that can enter a quadratic objective with analytical gradient and Hessian, and the outer loop re-estimates each $\delta_i$ from the sampled percentile so the surrogate tracks the true distribution. In homogenous phantom geometries, the paper demonstrates that the resulting probabilistic plans meet their stated probability targets and outperform composite-wise mini-max robust plans on the chosen trade-off: for matched CTV coverage, $P(D_{2\%} > 30\,\mathrm{Gy})$ dropped by 10\textendash15\% in the spherical cases and spinal overdosage probability dropped by 24\textendash28\%; for matched OAR dose, $P(D_{98\%} > 57\,\mathrm{Gy})$ increased by 67.5\textendash71\% in the spherical cases and by 10\textendash15\% in the spinal plans.
Load-bearing premise
The load-bearing premise is that the outer-loop update $\delta_i = (E[d_i] - d_{\alpha,i})/SD[d_i]$, damped with $\kappa = 0.2$ and stopped by a moving-average rule, converges for every voxel to the true percentile; the paper shows this empirically for one spherical plan but gives no convergence proof.
Editorial extensions
If this is right
- Treatment planners could prescribe explicit probability budgets, such as at most 10\% of uncertainty scenarios underdosing the CTV, instead of selecting a margin or a robust scenario set, and the optimizer would work toward that budget voxel by voxel.
- Because probabilistic objectives weight scenarios by their probability, plans become more conformal than composite-wise mini-max plans: the margin shrinks in directions where large shifts are unlikely, sparing OARs without losing CTV coverage.
- On the phantom cases, the method produces a plan at least as good as a tuned robust plan on one side of the trade-off and better on the other: matched CTV coverage with lower OAR overdosage probability, or matched OAR dose with higher target coverage probability.
- The same $E \pm \delta \cdot SD$ machinery can be pointed at other statistics, including CVaR, dose-coverage objectives, and biological metrics, because $\delta$ is just a mapping from a sampled statistic to a deterministic dose level.
- Plan evaluation can use the same probabilistic language as the optimization, with dose population histograms and scenario fractions replacing worst-case DVH bands.
Reading between the lines
- If the $\delta$-loop converges reliably on clinical anatomies, a natural extension is per-patient, per-fraction re-estimation of the uncertainty distribution, turning the optimizer into a component of adaptive proton therapy rather than a one-time planning step.
- The memory bottleneck from the $E[D_{ij}D_{ij'}]$ terms could be reduced by restricting probabilistic objectives to CTV edge voxels and OAR near-edge voxels; the paper mentions voxel sampling strategies, and the $\delta$-formulation would survive such a restriction.
- Replacing the percentile-based $\delta$-factor with a CVaR-based one would make the inner objective convex and might remove the need for damping and moving-average convergence checks; the paper explicitly identifies CVaR as future work.
- The probability-control claim is modality-independent, so the same framework is testable in photon VMAT, where uncertainty is smaller but inter-patient variation in PTV coverage is documented.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The paper proposes a probabilistic proton treatment planning framework in which voxel-wise dose percentiles are approximated as E[d] ± δ·SD[d], with a per-voxel and per-probability-level δ that is iteratively re-estimated in an outer loop from PCE-sampled percentiles. The inner loop minimizes a weighted sum of such E ± δ·SD objectives together with low-weight expectation-based terms, using fmincon with analytical gradient and Hessian. Uncertainties are modeled as Gaussian setup (3 mm) and range (3%) errors. The method is validated on a spherical CTV with and without an OAR and on a horseshoe-shaped CTV surrounding a cylindrical spine, and compared with composite-wise mini-max robust plans. The authors report improved OAR sparing at matched CTV coverage and improved CTV coverage at matched OAR dose, and a consistency check with the Van Herk margin recipe for systematic setup errors.
Significance. If the outer-loop convergence gap is closed, this is a useful proof-of-principle contribution: it replaces the scenario-set selection of robust optimization with voxel-specific probability control and provides a tractable PCE-based inner/outer formulation. Strengths include the independent final-evaluation PCE, the Γ-analysis benchmark against the dose engine in Appendix A, explicit reporting of convergence behavior and computation times, and use of the open-source OpenGPC toolbox. The main limitation is that the central probability-control claim rests on an unverified fixed-point assumption for the δ iteration, and the reported probability metrics are not always the same quantities that are optimized voxel-wise.
major comments (3)
- [Section 2.2.2, Eqs. (25)–(27), Table 8] The central claim that the final plan realizes Eq. (4) requires the outer-loop δ-factors to match the percentiles of the delivered beam weights, but no convergence theorem or contraction argument is given for the δ update Eq. (12), and the only convergence display is for one spherical plan (Fig. 19). The stopping rule Eq. (27) compares moving averages of percentiles computed at the damped iterates x^k, not the current inner solution x^k_*, and the tolerances in Table 8 are loose (τ_OAR,ν = 0.1 for ν = 90% and 95%, and 0.05 for ν = 98%). A stale δ can therefore halt the outer loop while the inner objective Eq. (20) is still optimizing a surrogate that does not correspond to the target percentile. The paper should either provide a convergence argument or, at minimum, report for every geometry and probability level a final-plan check of Eq. (4) at x* for all voxels, evaluated with an independent PCE or dose-engine sampling, rather than only representative voxels and moving averages.
- [Appendix A.2 and Section 2.5.1] The final voxel-dose PCE used for evaluation has ΔD2% up to about 2 Gy over the test scenarios and mean dose differences up to about 0.9 Gy for 95% of the voxels. Given that the probabilistic thresholds are 57 Gy, 64.2 Gy, and 30 or 54 Gy and the target probabilities are 2–10%, a 1–2 Gy tail error can shift a reported probability by several percentage points. No analysis is given of how PCE approximation error propagates into the reported probabilities, and the PCE accuracy benchmark is only demonstrated for the spherical CTV-only setupXYrange plan, not for the spinal or OAR cases. Please add a sensitivity analysis, for instance by recomputing the final-plan probabilities with dose-engine sampling on a subset of scenarios, and report PCE error bounds for each treatment site.
- [Section 3.2, Tables 3 and Figures 9/10] For the spherical CTV+OAR cases, the reported OAR improvement is the DVH-metric probability P(D_2% > 30 Gy), whereas the optimized objective Eq. (14) controls the voxel-wise probability P(d_i > 30 Gy) ≤ 10% for each OAR voxel. These two quantities are not equivalent, and no voxel-wise probability map is given for the spherical OAR, in contrast to the spinal case in Fig. 13. The same remark applies to the CTV metric P(D_98% > 57 Gy) versus the optimized voxel-wise underdosage probability. To support the claim that the method controls under- and overdosage probability per voxel, please provide voxel-wise acceptance maps or summary statistics for the spherical plans as well.
minor comments (3)
- [Figure 3] In the submitted version, the flowchart contains placeholder text ('Lorem ipsum') and unreadable path tokens such as '/gid00035/...'; this figure must be replaced with a legible version.
- [Section 3.1.2] The 'verification against the Van Herk margin recipe' is a consistency check for one plan under ideal spherical and static-dose-cloud assumptions; the text should say 'consistent with' rather than 'probabilistically equivalent to,' since the derivation relies on a 2D Gaussian population formula.
- [Discussion] The statement that 'in some probabilistic optimizations, even one or two outer loop iterations were sufficient' is not quantified; please specify which cases and how this was determined, since it bears on the convergence and stopping-rule discussion.
Circularity Check
No significant circularity: delta-factors are re-estimated from PCE percentiles and the final evaluation uses an independent PCE, so the probability-control claim does not reduce to its fitted inputs; minor self-citations are not load-bearing.
full rationale
The central derivation is not circular. The inner objective (Eq. 10) uses d_alpha = E - delta*SD (Eq. 11), with delta defined as (E - d_alpha)/SD from PCE-sampled percentiles (Eq. 12). At each outer iteration this substitution is a definitional identity, not a fitted prediction: delta is the output of the percentile estimate, and the final probability claim is verified by scenario counting on an independently constructed voxel-dose PCE (Section 2.5.1), not by the E +/- delta*SD surrogate. The PCE itself is validated inside the paper against the dose engine via Gamma-analysis (Appendix A), so the cited OpenGPC toolbox and prior PCE work by the same group (Perko et al. 2014, 2016; Rojo-Santiago et al.) are supporting but not load-bearing. The robust-plan comparisons tune robust weights to match one metric (D10th_98% or D90th_2%) and then report the other, so the reported CTV/OAR probability differences are not forced by the matching procedure. The genuine weakness, acknowledged by the empirical rather than proven convergence of the outer loop (Eqs. 25-27; Fig. 19; Discussion notes that one or two iterations sometimes sufficed), is an unproven fixed-point assumption and hence a correctness/convergence risk, not a circularity.
Assumptions & free parameters
free parameters (4)
- Per-voxel δ-factor (δα, δβ, δν) =
Updated each outer iteration via Eq. 12; converged values not tabulated
- Probabilistic objective weights (πα_CTV, πβ_CTV, πν_OAR, πlow_CTV, πlow_OAR, πtissue) =
Table 8: e.g., CTV-only 15,15,1,1; spinal OAR 750,15
- Outer-loop hyperparameters (κ, ΔW, Δk, τ) =
κ=0.2; ΔW=15-20; Δk=5-10; τ=5e-4 to 0.1
- Robust comparison weights (ωCTV, ωOAR, ωOARmax, ωnomCTV, ωtissue) =
Table 9: e.g., {120,1,1,160} and {100,10,10,100} for XZ cases
assumptions (6)
- domain assumption Gaussian independent setup (σ=3 mm in x,y) and range (σ=3%) uncertainties, truncated at the 99% confidence ellipsoid.
- domain assumption PCE of the dose-influence matrix (GO7E8PO8) and of voxel dose (GO6E7PO7) accurately represents the dose engine, including the lower and upper tails used for 2-10% percentiles.
- ad hoc to paper The outer-loop δ iteration converges to the true percentile for every voxel.
- ad hoc to paper During an inner optimization with fixed δ, improving E[d] ± δ·SD[d] moves the true target percentile in the intended direction.
- standard math Dose is a linear function of pencil-beam weights with a precomputed dose-influence matrix.
- domain assumption The static dose cloud approximation holds for the spherical margin comparison.
Cite this review
Pith. "Pith review of Probabilistic Proton Treatment Planning: a novel approach for optimizing underdosage and overdosage probabilities of target and organ structures." pith.science (2026). https://pith.science/paper/DXSHR4D2
@misc{pith2026250701763,
author = {Pith},
title = {Pith review of: Probabilistic Proton Treatment Planning: a novel approach for optimizing underdosage and overdosage probabilities of target and organ structures},
year = {2026},
howpublished = {\url{https://pith.science/paper/DXSHR4D2}},
note = {Machine review of arXiv:2507.01763}
}
abstract
Treatment planning uncertainties are typically managed using margin-based or robust optimization. Margin-based methods expand the clinical target volume (CTV) to a planning target volume, generally unsuited for proton therapy. Robust optimization considers worst-case scenarios, but its quality depends on the uncertainty scenario set: excluding extremes reduces robustness, while too many make plans overly conservative. Probabilistic optimization overcomes these limits by modeling a continuous scenario distribution. We propose a novel probabilistic optimization approach that steers plans toward individualized probability levels to control CTV and organs-at-risk (OARs) under- and overdosage. Voxel-wise dose percentiles ($d$) are estimated by expected value ($E$) and standard deviation (SD) as $E[d] \pm \delta \cdot SD[d]$, where $\delta$ is iteratively tuned to match the target percentile given Gaussian-distributed setup (3 mm) and range (3%) uncertainties. The method involves an inner optimization of $E[d] \pm \delta \cdot SD[d]$ for fixed $\delta$, and an outer loop updating $\delta$. Polynomial Chaos Expansion (PCE) provides accurate and efficient dose estimates during optimization. We validated the method on a spherical CTV abutted by an OAR in different directions and a horseshoe-shaped CTV surrounding a cylindrical spine. For spherical cases with similar CTV coverage, $P(D_{2\%} > 30 Gy)$ dropped by 10-15%; for matched OAR dose, $P(D_{98\%} > 57 Gy)$ increased by 67.5-71%. In spinal plans, $P(D_{98\%} > 57 Gy)$ increased by 10-15% while $P(D_{2\%} > 30 Gy)$ dropped 24-28%. Probabilistic and robust optimization times were comparable for spherical (hours) but longer for spinal cases (7.5 - 11.5 h vs. 9 - 20 min). Compared to discrete scenario-based optimization, the probabilistic method offered better OAR sparing or target coverage depending on the set priorities.
Figures
Figures from the paper (25 more)
Reference graph
Works this paper leans on
-
[1]
An, Y., Liang, J., Schild, S. E., Bues, M. and Liu, W. (2017). Robust treatment planning with conditional value at risk chance constraints in intensity-modulated proton therapy, Medical Physics44(1): 28–36. Biggs, S., Jennings, M., Swerdloff, S., Chlap, P., Lane, D., Rembish, J., McAloney, J., King, P., Ayala, R., Guan, F., Lambri, N., Crewson, C. and Sob...
work page 2017
-
[897]
Wu, Q. and Mohan, R. (2000). Algorithms and functionality of an intensity modulated radiotherapy optimization system,Medical Physics27(4): 701–711. Xiu, D. and Karniadakis, G. E. (2002). The wiener–askey polynomial chaos for stochastic differential equations,SIAM Journal on Scientific Computing24(2): 619–644. 33 A Cost-accuracy analyses of the dose approx...
work page 2000
-
[4555]
Bortfeld, T. (1997). An analytical approximation of the bragg curve for therapeutic proton beams,Medical Physics24(12): 2024–2033. Chan, T. C., Mahmoudzadeh, H. and Purdie, T. G. (2014). A robust-CVaR optimization approach with application to breast cancer therapy,European Journal of Operational Research238(3): 876–885. Chen, W., Unkelbach, J., Trofimov, ...
1997
Reviewed August 6, 2026 · model on record in the stance chip above.
Discussion (0). Sign in to comment.