🔢 mathematics

Geometry, Energy and Sensitivity in Stochastic Proton Dynamics

This paper presents a geometrically consistent numerical framework for stochastic proton transport that combines a logarithmic Milstein scheme for energy loss and a Lie-group integrator for angular diffusion to achieve stable, first-order convergence and reliable pathwise sensitivity estimates for dose deposition.

Original authors: Veronika Chronholm, Tristan Pryer

Published 2026-04-22

📖 6 min read🧠 Deep dive

CC BY 4.0

Original authors: Veronika Chronholm, Tristan Pryer

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to hit a tiny target with a stream of marbles, but these aren't ordinary marbles. They are protons (tiny particles used in cancer therapy), and they are traveling through a complex, bumpy fog (human tissue).

Your goal is to stop these protons exactly inside a tumor, depositing all their energy there to kill cancer cells, while sparing the healthy tissue before and after. This is called Proton Therapy.

The problem? Protons are chaotic. They don't travel in straight lines like laser beams. They bump into atoms, lose energy, and get knocked sideways. Sometimes they stop a little too early; sometimes a little too late. This randomness is called stochasticity.

This paper is about building a better mathematical simulation to predict exactly where these protons will go and how much energy they will drop along the way. The authors, Veronika and Tristan, didn't just build a better simulator; they built a "smart" one that respects the laws of physics and can tell doctors how sensitive the treatment is to tiny changes in the setup.

Here is the breakdown of their work using simple analogies:

1. The Problem: The "Naive" Simulator

Imagine you are trying to simulate a drunk person walking on a tightrope (the proton moving through tissue).

The Old Way (Standard Math): If you use standard math to guess their next step, you might accidentally calculate that they stepped off the tightrope and fell into the void (negative energy) or that they are now floating in mid-air (direction vector no longer points the right way). In physics, energy can't be negative, and direction must always be a valid angle. Standard math often breaks these rules, leading to garbage results.
The Sensitivity Problem: Doctors also need to know: "If I change the energy of the beam by 1%, how much does the damage to the tumor change?" The old way of answering this is like guessing by trial and error: "Let's run the simulation, then run it again with a slightly different number, and see the difference." This is slow, noisy, and often gives the wrong answer because the randomness of the protons makes the comparison messy.

2. The Solution: "Geometry-Aware" Math

The authors built a new set of rules for the simulation that forces the math to obey the laws of physics at every single step.

A. The Energy Trick: The "Logarithmic Elevator"

Protons lose energy as they travel. Eventually, they stop.

The Issue: If you simulate energy directly, a random "jitter" in the math might make the energy dip below zero. That's impossible in real life.
The Fix: Instead of tracking the energy number directly, they track the logarithm of the energy (think of it like tracking the height of a building instead of the number of floors).
- Analogy: Imagine you are tracking a balloon's height. If you track the height directly, you might accidentally calculate it going underground. But if you track the logarithm of the height, the math naturally keeps it above zero. When you convert it back at the end, the balloon is guaranteed to be floating, never underground. This ensures the proton always has positive energy.

B. The Direction Trick: The "Globe Walker"

Protons move in 3D space, and their direction is always a vector pointing somewhere.

The Issue: Standard math updates might make the direction vector grow longer or shorter, effectively stretching or shrinking the arrow. But a direction arrow should always be exactly the same length (like a compass needle).
The Fix: They used a Lie-Group Integrator.
- Analogy: Imagine walking on the surface of the Earth. If you take a step using flat-map math, you might accidentally walk off the edge of the world or into the sky. Instead, the authors use a "curved step" (an exponential map) that keeps you glued to the surface of the sphere. No matter how many steps you take, you stay on the globe. This keeps the proton's direction perfectly accurate.

3. The "Super-Sensitivity" Feature

This is the paper's biggest innovation.

The Old Way: To find out how sensitive the treatment is, you run the simulation twice. Once with the original settings, once with a tiny tweak. You subtract the results. Because protons are random, the "noise" in the two runs often drowns out the tiny difference you are trying to measure. It's like trying to hear a whisper in a hurricane by shouting twice.
The New Way (Pathwise Sensitivity): The authors developed a method where the simulation tracks the "sensitivity" alongside the proton's path, using the exact same random jitters.
- Analogy: Imagine you are walking a dog on a leash. You want to know how much the dog's path changes if you pull the leash slightly harder. Instead of walking the dog twice (once with a loose leash, once with a tight one), you walk the dog once, but you carry a second, invisible leash attached to a sensor that measures exactly how the dog's path would have changed if you had pulled harder.
- Because the "sensitivity" and the "proton" experience the exact same random bumps, the noise cancels out. The result is a crystal-clear, high-precision answer about how the treatment will react to changes.

4. Why This Matters

Safety: By ensuring the math never breaks the laws of physics (like negative energy), the simulations are more reliable.
Precision: The new method gives doctors a much clearer picture of where the "Bragg Peak" (the spot where the proton dumps its energy) will land.
Optimization: Because the sensitivity calculation is so fast and accurate, doctors can use it to automatically tweak treatment plans to maximize tumor killing and minimize damage to healthy organs, something that was previously too slow or noisy to do effectively.

Summary

The authors took a chaotic, messy problem (protons bouncing through tissue) and built a physics-respecting, high-precision simulator. They used clever math tricks (logarithms and curved surfaces) to keep the simulation honest, and they invented a "shadow tracking" method to calculate sensitivity with incredible clarity. This helps doctors aim proton beams with the precision of a sniper, rather than a shotgun.

1. Problem Statement

The paper addresses the challenges in numerically simulating stochastic proton transport, a critical component in high-precision applications like proton therapy. Current deterministic models (e.g., Continuous Slowing Down Approximation) fail to capture statistical fluctuations such as range straggling (variations in stopping distance) and angular scattering (lateral beam broadening).

While Monte Carlo and Stochastic Differential Equation (SDE) based methods can model these fluctuations, they face three specific numerical difficulties:

Positivity Violation: Standard discretization schemes (like Euler-Maruyama) can produce negative energy values due to stochastic noise, which is physically impossible.
Geometric Drift: Angular dynamics occur on a manifold (the unit sphere $S^2$ ). Standard updates often cause direction vectors to drift off the sphere, violating the unit-norm constraint.
Inefficient Sensitivity Analysis: Calculating dose sensitivities with respect to physical parameters (e.g., stopping power, straggling coefficients) is essential for treatment planning and uncertainty quantification. Standard finite-difference methods applied to stochastic systems suffer from a severe bias-variance tradeoff, resulting in high variance and computational inefficiency.

2. Methodology

The authors propose a unified framework combining structure-preserving geometric integrators with pathwise sensitivity analysis.

A. Mathematical Formulation

The proton trajectory is modeled by the state $(X_t, E_t, \Omega_t)$ , representing position, energy, and direction. The dynamics are governed by a coupled SDE system:

Position: $dX_t = \Omega_t dt$
Energy: Includes deterministic stopping power $S(E)$ and a stochastic term for range straggling $\sqrt{T(E)} dW^E_t$ .
Direction: Modeled as Brownian motion on the sphere with diffusion coefficient $\epsilon(E)$ .

To handle the singularity at $E=0$ and ensure positivity, the authors apply a logarithmic transformation $Y_t = \log(E_t)$ . This converts the energy SDE into one defined on $\mathbb{R}$ , ensuring that the reconstructed energy $E_t = \exp(Y_t)$ remains strictly positive.

B. Numerical Schemes

The paper introduces two primary discretization schemes:

Geometric Euler Scheme (Strong Order 1/2):
- Energy: Uses the Euler-Maruyama method on the log-transformed variable $Y_t$ .
- Angle: Uses the Riemannian exponential map ( $\text{Exp}_{\Omega_n}$ ) to update the direction vector. This ensures the updated vector remains exactly on the unit sphere $S^2$ , preserving the geometric constraint without post-hoc projection.
High-Order Structure-Preserving Scheme (Strong Order 1):
- Energy: Applies a Milstein scheme to the log-energy SDE. This includes a correction term involving the derivative of the diffusion coefficient, achieving strong order 1 convergence while maintaining positivity.
- Angle: Employs a Lie-group integrator based on the Runge-Kutta-Munthe-Kaas (RKMK) method. This lifts the dynamics to the Lie algebra (tangent space), performs a higher-order update including Lévy area corrections (essential for strong order 1 in non-commutative vector fields), and maps back to the sphere via the exponential map.

C. Pathwise Sensitivity Analysis

Instead of finite differences, the authors derive coupled forward sensitivity equations.

They differentiate the SDE system with respect to parameters $\theta \in \{\alpha, p, \kappa\}$ (stopping power and straggling parameters).
This yields linear SDEs for the sensitivities $J^\theta_t = \partial_\theta E_t$ , driven by the same Brownian motion as the state variables.
Regularization: Since the stopping time (when energy hits a threshold) is not differentiable with respect to parameters in the stochastic case, the authors introduce a regularized dose functional using a smooth transition function (hyperbolic tangent) instead of a hard cutoff. This allows for valid pathwise differentiation.

3. Key Contributions

Structure-Preserving Discretization: The development of a numerical framework that strictly enforces physical constraints (energy positivity and spherical geometry) at the discrete level, avoiding the artifacts common in standard schemes.
High-Order Convergence: Proving that the proposed Milstein-based log-energy scheme and RKMK-based angular integrator achieve strong order 1 convergence, a significant improvement over the standard order 1/2.
Efficient Pathwise Sensitivities: Deriving and implementing coupled sensitivity equations that provide unbiased gradient estimators with significantly lower variance than finite-difference methods.
Regularized Observable Framework: Formulating dose deposition as a regularized path-dependent functional, enabling the application of pathwise sensitivity methods to stochastic stopping-time problems.

4. Results

Numerical experiments validate the theoretical findings:

Convergence Rates: The schemes achieve the expected strong convergence rates (0.5 for Geometric Euler, 1.0 for Geometric RKMK/Milstein).
Geometric Fidelity: Unlike naive Euler schemes which exhibit norm drift and bias in angular distributions, the geometric integrators preserve the unit-norm constraint and reproduce correct statistical distributions even over long time intervals.
Dose Distribution: Simulations successfully reproduce physical phenomena such as the Bragg peak, range straggling (broadening of the peak with increasing $\kappa$ ), and lateral beam spreading. The inclusion of energy-dependent angular diffusion (Molière theory) was shown to significantly affect the distal edge of the dose distribution.
Sensitivity Performance:
- Pathwise estimators captured sharper sensitivity profiles compared to finite-difference benchmarks.
- Finite-difference methods were shown to underestimate gradient magnitudes and introduce artificial smoothing (numerical diffusion).
- The pathwise method avoided the bias-variance tradeoff inherent in finite differences, providing stable estimates with lower empirical variance.

5. Significance

This work bridges the gap between rigorous stochastic analysis and practical medical physics applications.

Clinical Relevance: By providing accurate, low-variance sensitivity estimates, the framework facilitates uncertainty quantification and machine-learning-guided treatment optimization in proton therapy, where small errors in dose prediction can have critical clinical consequences.
Theoretical Advancement: It demonstrates the power of geometric integration (Lie groups and manifold constraints) in solving complex, coupled SDE systems arising in particle transport, offering a template for other physical systems with strict geometric or positivity constraints.
Computational Efficiency: The pathwise sensitivity approach offers a computationally superior alternative to finite differences for stochastic inverse problems and optimization, potentially reducing the computational cost of calibration and robustness analysis in radiotherapy planning.