SDE-based Monte Carlo dose calculation for proton therapy validated against Geant4

Purpose: To assess the clinical impact of Monte Carlo (MC) dose calculations in proton therapy by analyzing uncertainties of analytical algorithms when predicting dose to target and critical structures and beam ranges in patient geometries. Methods: We used TOPAS to simulate double scattering proton treatments, which were compared to analytical treatment planning dose calculations (TPC). We investigated: 1)beam ranges, dose distributions, does‐rate profiles and output factors using various dedicated experiments.range differences caused by the patient heterogeneity considering 48 fields and 4 treatment sites. We calculated distal 2)dose surfaces composed of beam ranges and assessed average range difference (dR) and root‐mean‐square deviation (RMSD).3)dose delivery accuracy using a DVH based analysis considering patient treatment plans from 3 sites. 4)the accuracy of dose delivery for small fields (diameter below 7cm) based on patient heterogeneity for 38 head fields. We assess differences in D50. Results: 1)MC and measurements agree within statistical uncertainties.2)For 3 of 4 patient sites the RMSD between TPC and TOPAS was <2mm. Head treatments indicate significant RMSD of up to 6mm. 3)TPC DVHs showed underdosage for head cases by 2% compared to TOPAS, discrepancies in critical volumes can reach up to 50%, prostate patients show significant penumbra differences, liver patients show good agreement between TOPAS and TPC. The site‐specific findings can be explained by beam range and geometrical complexities.4)For small proton fields, discrepancies of up to 5.4% were found in D50 for single fields. A clear correlation between the accuracy and patient heterogeneity was established. Conclusion: This study highlights treatment sites and scenarios where MC simulations are recommended for proton therapy, i.e. small fields (due to scattering disequilibrium) and head treatments with bone/air/tissue interfaces affecting the ability of analytical algorithms to predict the correct beam range. Absolute dose differences where typically within 2% for large fields.

In proton therapy, pencil-beam algorithms (PBAs) are the most widely used dose calculation methods. However, the PB calculations that employ one-dimensional density scaling neglect the effects of lateral density heterogeneity on the dose distributions, whereas some particles included in such pencil beams could overextend beyond the interface of the density heterogeneity. We have simplified a pencil-beam redefinition algorithm (PBRA), which was proposed for electron therapy, by a spatial resampling technique toward an application for proton therapy. The purpose of this study is to evaluate the calculation results of the spatial resampling technique in terms of lateral density heterogeneity by comparison with the dose distributions that were measured in heterogeneous slab phantoms. The pencil beams are characterized for multiple residual-range (i.e., proton energy) bins. To simplify the PBRA, the given pencil beams are resampled on one or two transport planes, in which smaller sub-beams that are parallel to each other are generated. We addressed the problem of lateral density heterogeneity comparing the calculation results to the dose distributions measured at different depths in heterogeneous slab phantoms using a two-dimensional detector. Two heterogeneity slab phantoms, namely, phantoms A and B, were designed for the measurements and calculations. In phantom A, the heterogeneity slab was placed close to the surface. On the other hand, in phantom B, it was placed close to the Bragg peak in the mono-energetic proton beam. In measurements, lateral dose profiles showed a dose reduction and increment in the vicinity of x = 0 mm in both phantoms at depths z = 142 and 161 mm due to lateral particle disequilibrium. In phantom B, these dose reduction∕increment effects were higher∕lower, respectively, than those in phantom A. This is because a longer distance from the surface to the heterogeneous slab increases the strength of proton scattering. Sub-beams, which were generated from the resampling plane, formed a detouring∕overextending path that was different from that of elemental pencil beams. Therefore, when the spatial resampling was implemented at the surface and immediately upstream of the lateral heterogeneity, the calculation could predict these dose reduction∕increment effects. Without the resampling procedure, these dose reduction∕increment effects could not be predicted in both phantoms owing to the blurring of the pencil beam. We found that the PBA with the spatial resampling technique predicted the dose reduction∕increment at the dose profiles in both phantoms when the sampling plane was defined immediately upstream of the heterogeneous slab. We have demonstrated the implementation of a spatial resampling technique for pencil-beam calculation to address the problem of lateral density heterogeneity. While further validation is required for clinical use, this study suggests that the spatial resampling technique can make a significant contribution to proton therapy.

Advancing (Proton) Radiation Therapy

Background: Growth trajectories are highly variable between children, making epidemiological analyses challenging both to the identification of malnutrition interventions at the population level and also risk assessment at individual level. We introduce stochastic differential equation (SDE) models into child growth research. SDEs describe flexible dynamic processes comprising: drift - gradual smooth changes – such as physiology or gut microbiome, and diffusion - sudden perturbations, such as illness or infection. Methods: We present a case study applying SDE models to child growth trajectory data from the Haydom, Tanzania and Venda, South Africa sites within the MAL-ED cohort. These data comprise n=460 children aged 0-24 months. A comparison with classical curve fitting (linear mixed models) is also presented. Results: The SDE models offered a wide range of new flexible shapes and parameterizations compared to classical additive models, with performance as good or better than standard approaches. The predictions from the SDE models suggest distinct longitudinal clusters that form distinct ‘streams’ hidden by the large between-child variability. Conclusions: Using SDE models to predict future growth trajectories revealed new insights in the observed data, where trajectories appear to cluster together in bands, which may have a future risk assessment application. SDEs offer an attractive approach for child growth modelling and potentially offer new insights.

Stochastic differential equation (SDE) models are the foundation for pricing and hedging financial derivatives. The drift and volatility functions in SDE models are typically chosen to be algebraic functions with a small number ( < 5 ) of parameters which can be calibrated to market data. A more flexible approach is to use neural networks to model the drift and volatility functions, which provides more degrees-of-freedom to match observed market data. Training of models requires optimizing over an SDE, which is computationally challenging. For European options, we develop a fast stochastic gradient descent (SGD) algorithm for training the neural network-SDE model. Our SGD algorithm uses two independent SDE paths to obtain an unbiased estimate of the direction of steepest descent. For American options, we optimize over the corresponding Kolmogorov partial differential equation (PDE). The neural network appears as coefficient functions in the PDE. Models are trained on large datasets (many contracts), requiring either large simulations (many Monte Carlo samples for the stock price paths) or large numbers of PDEs (a PDE must be solved for each contract). Numerical results are presented for real market data including S&P 500 index options, S&P 100 index options, and single-stock American options. The neural-network-based SDE models are compared against the Black-Scholes model, the Dupire's local volatility model, and the Heston model. Models are evaluated in terms of how accurate they are at pricing out-of-sample financial derivatives, which is a core task in derivative pricing at financial institutions. Specifically, we calibrate a neural network-SDE model to market data for a financial derivative on an asset with price S t with a payoff function g ( s ) , and we then evaluate its generalization accuracy for a financial derivative on the same asset S t but with a different payoff function f ( s ) . In addition to comparing out-of-sample pricing accuracy, we evaluate the hedging performance of the neural network-SDE model.

We develop a new continuous‐time stochastic gradient descent method for optimizing over the stationary distribution of stochastic differential equation (SDE) models. The algorithm continuously updates the SDE model's parameters using an estimate for the gradient of the stationary distribution. The gradient estimate is simultaneously updated using forward propagation of the SDE state derivatives, asymptotically converging to the direction of steepest descent. We rigorously prove convergence of the online forward propagation algorithm for linear SDE models (i.e., the multidimensional Ornstein–Uhlenbeck process) and present its numerical results for nonlinear examples. The proof requires analysis of the fluctuations of the parameter evolution around the direction of steepest descent. Bounds on the fluctuations are challenging to obtain due to the online nature of the algorithm (e.g., the stationary distribution will continuously change as the parameters change). We prove bounds for the solutions of a new class of Poisson partial differential equations (PDEs), which are then used to analyze the parameter fluctuations in the algorithm. Our algorithm is applicable to a range of mathematical finance applications involving statistical calibration of SDE models and stochastic optimal control for long time horizons where ergodicity of the data and stochastic process is a suitable modeling framework. Numerical examples explore these potential applications, including learning a neural network control for high‐dimensional optimal control of SDEs and training stochastic point process models of limit order book events.

Stochastic Differential Equations (SDE) are often used to model the stochastic dynamics of biological systems. Unfortunately, rare but biologically interesting behaviors (e.g., oncogenesis) can be difficult to observe in stochastic models. Consequently, the analysis of behaviors of SDE models using numerical simulations can be challenging. We introduce a method for solving the following problem: given a SDE model and a high-level behavioral specification about the dynamics of the model, algorithmically decide whether the model satisfies the specification. While there are a number of techniques for addressing this problem for discrete-state stochastic models, the analysis of SDE and other continuous-state models has received less attention. Our proposed solution uses a combination of Bayesian sequential hypothesis testing, non-identically distributed samples, and Girsanov's theorem for change of measures to examine rare behaviors. We use our algorithm to analyze two SDE models of tumor dynamics. Our use of non-identically distributed samples sampling contributes to the state of the art in statistical verification and model checking of stochastic models by providing an effective means for exposing rare events in SDEs, while retaining the ability to compute bounds on the probability that those events occur.

The chance that a freeway will breakdown, transition from a free-flow to a congested state, is normally assumed to increase with an increase in traffic volume V (vehicles per unit time). In this paper, this assumption is challenged. Traffic density K (vehicles per unit length) proves to be a better predictor. Diffusion or stochastic differential equation (SDE) modeling is used to substantiate the claim. SDE modeling is especially useful in explaining the role that traffic noise (volatility) plays in breakdown. The SDE models take advantage of the unique properties of the geometric Brownian motion (gBM) and Ornstein-Uhlenbeck (OU) model structures. The breakdown probability model of π (K) and delay models provide accurate forecasts.

Background: Growth trajectories are highly variable between children, making epidemiological analyses challenging both to the identification of malnutrition interventions at the population level and also risk assessment at individual level. We introduce stochastic differential equation (SDE) models into child growth research. SDEs describe flexible dynamic processes comprising: drift - gradual smooth changes – such as physiology or gut microbiome, and diffusion - sudden perturbations, such as illness or infection. Methods: We present a case study applying SDE models to child growth trajectory data from the Haydom, Tanzania and Venda, South Africa sites within the MAL-ED cohort. These data comprise n=460 children aged 0-24 months. A comparison with classical curve fitting (linear mixed models) is also presented. Results: The SDE models offered a wide range of new flexible shapes and parameterizations compared to classical additive models, with performance as good or better than standard approaches. The predictions from the SDE models suggest distinct longitudinal clusters that form distinct ‘streams’ hidden by the large between-child variability. Conclusions: Using SDE models to predict future growth trajectories revealed new insights in the observed data, where trajectories appear to cluster together in bands, which may have a future risk assessment application. SDEs offer an attractive approach for child growth modelling and potentially offer new insights.

Background: Growth trajectories are highly variable between children, making epidemiological analyses challenging both to the identification of malnutrition interventions at the population level and also risk assessment at individual level. We introduce stochastic differential equation (SDE) models into child growth research. SDEs describe flexible dynamic processes comprising: drift - gradual smooth changes - such as physiology or gut microbiome, and diffusion - sudden perturbations, such as illness or infection. Methods: We present a case study applying SDE models to child growth trajectory data from the Haydom, Tanzania and Venda, South Africa sites within the MAL-ED cohort. These data comprise n=460 children aged 0-24 months. A comparison with classical curve fitting (linear mixed models) is also presented. Results: The SDE models offered a wide range of new flexible shapes and parameterizations compared to classical additive models, with performance as good or better than standard approaches. The predictions from the SDE models suggest distinct longitudinal clusters that form distinct 'streams' hidden by the large between-child variability. Conclusions: Using SDE models to predict future growth trajectories revealed new insights in the observed data, where trajectories appear to cluster together in bands, which may have a future risk assessment application. SDEs offer an attractive approach for child growth modelling and potentially offer new insights.

Background and ObjectivesLevodopa concentration in patients with Parkinson’s disease is frequently modelled with ordinary differential equations (ODEs). Here, we investigate a pharmacokinetic model of plasma levodopa concentration in patients with Parkinson’s disease by introducing stochasticity to separate the intra-individual variability into measurement and system noise, and to account for auto-correlated errors. We also investigate whether the induced stochasticity provides a better fit than the ODE approach.MethodsIn this study, a system noise variable is added to the pharmacokinetic model for duodenal levodopa/carbidopa gel (LCIG) infusion described by three ODEs through a standard Wiener process, leading to a stochastic differential equations (SDE) model. The R package population stochastic modelling (PSM) was used for model fitting with data from previous studies for modelling plasma levodopa concentration and parameter estimation. First, the diffusion scale parameter (σw), measurement noise variance, and bioavailability are estimated with the SDE model. Second, σw is fixed to certain values from 0 to 1 and bioavailability is estimated. Cross-validation was performed to compare the average root mean square errors (RMSE) of predicted plasma levodopa concentration.ResultsBoth the ODE and the SDE models estimated bioavailability to be approximately 75%. The SDE model converged at different values of σw that were significantly different from zero. The average RMSE for the ODE model was 0.313, and the lowest average RMSE for the SDE model was 0.297 when σw was fixed to 0.9, and these two values are significantly different.ConclusionsThe SDE model provided a better fit for LCIG plasma levodopa concentration by approximately 5.5% in terms of mean percentage change of RMSE.

To investigate the feasibility of proton therapy planning using stopping power ratio (SPR) maps directly generated from spectral CT raw data, and to perform a comparative evaluation of dose calculation uncertainties between SPR maps derived from conventional CT Hounsfield Unit (HU) conversion and direct spectral CT SPR generation. A retrospective analysis was conducted on 30 patients with mid-thoracic esophageal squamous cell carcinoma (ESCC) who underwent pre-treatment spectral CT imaging. Target volumes and organs at risk (OARs) were delineated on contrast-enhanced CT images and subsequently registered to both non-contrast CT and SPR maps. Three treatment plans were generated: Intensity-modulated radiotherapy (IMRT) plan based on conventional CT, Intensity-modulated proton therapy (IMPT) plan using HU-SPR conversion, IMPT plan utilizing direct SPR maps (IMPT-SPR) from spectral CT. Dose-volume parameters for target volumes and OARs (lungs, heart, spinal cord) were systematically analyzed. Comparative dosimetric analyses were performed among the three plans and between paired groups. All plans met clinical radiotherapy requirements. For OARs (lungs, heart), IMPT plans demonstrated significantly lower dose-volume parameters compared to IMRT, except for maximum dose (Dmax). Between the two IMPT approaches, no statistically significant differences were observed in dose-volume parameters (p>0.05), except for the gradient index which was significantly higher in the HU-converted IMPT plan (p<0.05). No significant differences were detected in heart, lung and spinal cord dosimetric parameters between IMPT approaches. IMPT demonstrated superior OAR sparing compared to IMRT. For mid thoracic ESCC patients under proton therapy, dose calculations based on CT-HU conversion was showed comparable dosimetric impact to DECT-derived SPR in terms of target coverage and OAR protection. These findings support the clinical feasibility of conventional CT-based proton therapy planning and dose calculation.

While cone beam computed tomography (CBCT) is able to provide patient anatomical information, its image quality is severely degraded due to scatter contamination, which degrades the accuracy of CBCT-based dose distribution estimation in proton therapy. In this work, we combined two existing scatter kernel correction methods: the point-spread function (PSF)-based scatter kernel derivation method and the fast adaptive scatter kernel superposition (fASKS) model, and evaluated the impact of the modified fASKS (mfASKS) correction on the accuracy of proton dose distribution estimation. To evaluate feasibility of the mfASKS approach using accurate scatter distributions, both Monte Carlo simulations and experiments were performed for an on-board CBCT machine integrated with a proton therapy machine. We developed a strategy to modify central intensity, constant intensity, and amplitude of the scatter kernels derived from PSFs for the fASKS model. A parameter required for the fASKS model was derived by optimizing uniformity in the mfASKS-corrected reconstructed images. Subsequently, the mfASKS model was used to remove scatter in CBCT imaging. We quantitatively compared the Hounsfield Unit (HU) and proton stopping power ratio (SPR) images for five different phantoms. To assess improvement of dose calculation accuracy, a series of proton treatment plans were produced using the CBCT images with and without the mfASKS correction. The accuracies of both HU and SPR intensity quantifications are improved as a result of the mfASKS correction. Mean absolute water-equivalent path length difference to the true value decreases from 10.3 to 0.934mm for the Gammex phantom (simulation). At the same time, mfASKS is able to offer more accurate dose distributions, especially at the distal fall-off region where noticeable dose overestimation is observed in the uncorrected scenario. Mean absolute relative error of proton range in the pelvic phantom improves from 5.03% to 2.57% (experiment). mfASKS enables more accurate CBCT-based proton dose calculation. This technique has significant implications in image-guided radiotherapy and dose verifications in adaptive proton therapy.

Accurate dose calculation is a critical step in proton therapy. A novel machine learning-based approach was proposed to achieve comparable accuracy to that of Monte Carlo simulation while reducing the computational time. Computed tomography-based patient phantoms were used and three treatment sites were selected (thorax, head, and abdomen), comprising different beam pathways and beam energies. The training data were generated using Monte Carlo simulations. A discovery cross-domain generative adversarial network (DiscoGAN) was developed to perform the mapping between two domains: stopping power and dose, with HU values from CT images incorporated as auxiliary features. The accuracy of dose calculation was quantitatively evaluated in terms of mean relative error (MRE) and mean absolute error (MAE). The relationship between the DiscoGAN performance and other factors such as absolute dose, beam energy and location within the beam cross-section (center and off-center lines) was examined. The DiscoGAN model is found to be effective in dose calculation. For the abdominal case, the MRE is found to 1.47% (mean), 3.30% (maximum) and 0.67% (minimum). For the thoracic case, the MRE is found to ~2.43% (mean), 4.80% (maximum) and 0.71% (minimum). For the head case, the MRE is found to ~2.83% (mean), 4.84% (maximum) and 1.01% (minimum). Comparable accuracy is found in the independent validation dataset (different CT images), achieving a mean MRE of ~1.65% (thorax), 4.02% (head) and 1.64% (abdomen). For the energy span between 80 and 130MeV, no strong dependency of accuracy on beam energy is found. The results imply that no systematic deviation, either over-dose or under-dose, occurs between the predicted dose and raw dose. The DiscoGAN framework demonstrates great potential as a tool for dose calculation in proton therapy, achieving comparable accuracy yet being more efficient relative to Monte Carlo simulation. Its comparison with the pencil beam algorithm (PBA) will be the next step of our research. If successful, our proposed approach is expected to find its use in more advanced applications such as inverse planning and adaptive proton therapy.

This article examines the critical role of fast Monte Carlo (MC) dose calculations in advancing proton therapy techniques, particularly in the context of increasing treatment customization and precision. As adaptive radiotherapy and other patient-specific approaches evolve, the need for accurate and precise dose calculations, essential for techniques like proton-based stereotactic radiosurgery, becomes more prominent. These calculations, however, are time-intensive, with the treatment planning/optimization process constrained by the achievable speed of dose computations. Thus, enhancing the speed of MC methods is vital, as it not only facilitates the implementation of novel treatment modalities but also leads to more optimal treatment plans. Today, the state-of-the-art in MC dose calculation speeds is 106–107 protons per second. This review highlights the latest advancements in fast MC dose calculations that have led to such speeds, including emerging artificial intelligence-based techniques, and discusses their application in both current and emerging proton therapy strategies.