Analytical probabilistic modeling for radiation therapy treatment planning

This paper introduces the concept of analytical probabilistic modeling (APM) to quantify uncertainties in quality indicators of radiation therapy treatment plans. Assuming Gaussian probability densities over the input parameters of the treatment plan quality indicators, APM enables the calculation of the moments of the induced probability density over the treatment plan quality indicators by analytical integration. This paper focuses on analytical probabilistic dose calculation algorithms and the implications of APM regarding treatment planning. We derive closed-form expressions for the expectation value and the (co)variance of (1) intensity-modulated photon and proton dose distributions based on a pencil beam algorithm and (2) the standard quadratic objective function used in inverse planning. Complex correlation models of high dimensional uncertain input parameters and the different nature of random and systematic uncertainties in fractionated radiation therapy are explicitly incorporated into APM. APM variance calculations on phantom data sets show that the correlation assumptions and the difference of random and systematic uncertainties of the input parameters have a crucial impact on the uncertainty of the resulting dose. The derivations regarding the quadratic objective function show that APM has the potential to enable robust planning at almost the same computational cost like conventional inverse planning after a single probabilistic dose calculation. Beneficial applications of APM in the context of radiation therapy treatment planning are feasible.


Introduction
Quantitative science requires an adequate analysis of uncertainties. Both experimental and theoretical data are only meaningful in combination with the associated error bars.
In radiation therapy, uncertainties regarding treatment plan quality indicators may originate from setup uncertainties, range uncertainties, intrafractional motion, and interfractional motion, among others. Goitein (1985) was one of the first to explicitly discuss this issue in a 1985 paper, where he quantified the variability of a dose distribution by simulating three different dose calculation scenarios.
Despite the ever increasing complexity of dose distributions formed by multiple intensitymodulated beams, current clinical practice does not include a patient-specific uncertainty analysis of treatment plan quality indicators (Bortfeld 2006). Instead, uncertainties are accounted for with standardized margin concepts focusing on tumor coverage (van Herk et al 2004).
Approaches to quantify the uncertainties of treatment plan quality indicators for individual patients, such as worst case simulations and variance calculations based on dose blurring or sampling, remain computationally challenging or suffer from inherent limitations. Worst case simulations (Goitein 1985) provide only point estimates of a probability density and the worst case in the space of the input parameters does not necessarily correspond to the worst case in the space of the treatment plan quality indicators. Dose blurring (van Herk et al 2002, Baum et al 2004 may conceal the true variability of treatment plan quality indicators (Ploquin et al 2006). Reliable variance estimates usually require time consuming sampling processes, i.e., repeated simulations of different scenarios. The number of samples ranges from 10 4 , as applied by Maleike et al (2006) who consider uncorrelated motion of individual voxels of the patient anatomy, to 10 2 , as applied by Sobotta et al (2012) who consider rigid shifts of the patient anatomy and use Gaussian Process regression to decrease the number of samples.
Incorporating uncertainties in more than three translational parameters, accommodating arbitrary correlation models of the uncertainties in these parameters, and accounting for the different nature of random and systematic uncertainties in the context of fractionated radiation therapy additionally compromises the efficiency and validity of these approaches.
Furthermore, it is impossible to efficiently obtain information about the uncertainty induced by different intensity-modulations, i.e., different bixel weights, with dose blurring, sampling, and worst case simulations. Consequently, approaches which incorporate uncertainty information into the planning process in order to obtain robust treatment plans require repeated sampling processes or worst case simulations for the optimization (Unkelbach et al 2009, Unkelbach and Oelfke 2004, Pflugfelder et al 2008, Fredriksson et al 2011, Liu et al 2012. We introduce the concept of analytical probabilistic modeling (APM) to quantify uncertainties in radiation therapy treatment planning. Based on explicit assumptions about the probability density of the most important input parameters that are subject to uncertainty, we show that it is possible to derive closed-form expressions for the expectation value and the (co)variance of treatment plan quality indicators. APM is a bottom-up approach propagating the probability density over the underlying parameters into a probability density over the resulting quality indicators by analytical integration. It eliminates the need for sampling entirely, while being generally applicable to a wide range of problems in radiation therapy treatment planning. It accommodates complex correlation models of high dimensional uncertain input parameters, accounts for the different nature of random and systematic uncertainties in fractionated radiation therapy, and generalizes to arbitrary intensity-modulations.
In this paper, we focus on analytical probabilistic dose calculations incorporating range and setup uncertainties for photons and protons. We derive analytical expressions for the expectation value and the (co)variance of the resulting dose distributions based on a conventional pencil beam algorithm (Hong et al 1999) and show how APM may accelerate robust optimization in combination with the standard quadratic objective function used in inverse planning (Oelfke and Bortfeld 2001).
The overarching idea of APM is introduced in section 2, the analytical probabilistic dose calculation and its foundation are laid out in section 3, the derivations regarding the quadratic objective function are presented in section 4, and the implications of APM are discussed in section 5 where we also elaborate on the resulting novel possibilities in the context of probabilistic planning.

Analytical probabilistic modeling
We consider a treatment plan quality indicator I(X ) as a function of uncertain input parameters X, e.g. the dose distribution as a function of the patient position. If the uncertainty in the input parameters X is described with a probability density p(X ), the mth moment of the probability density over I(X ) is given by (1) In its most general form, this integral has to be solved with Monte Carlo integration. However, parameterizing the uncertainties in X by a multivariate normal distribution N (X;X, X ) with mean vectorX and covariance matrix X and replacing the treatment plan quality indicator I(X ) with a functional approximationÎ(X ) that can be integrated analytically against Gaussian densities renders the following integral analytically tractable: Hence, it is possible to compute the mean and (co)variance of I(X ) in closed form. The challenge is to find descriptions for N (X;X , X ) andÎ(X ) which allow for high approximation quality and low computational cost at the same time.
Modeling uncertainties with Gaussian densities is a widely accepted practice in radiation therapy even though Gaussian densities are not necessarily the best approximation of reality. We admit to use the Gaussian family mainly out of algebraic convenience as marginalization and conditioning, exercised repeatedly in section 3, correspond to simple linear algebraic computations. If required, more complex probability densities can be constructed as a weighted sum of Gaussians.
In general, APM may be applicable to various treatment plan quality indicators and sources of uncertainty. In this paper we focus on analytical probabilistic proton and photon dose calculations under the influence of range and setup uncertainties and the implications for robust treatment planning. Alternative applications of APM are outlined in section 5.

Nominal dose calculation
Let V be the number of voxels i and B be the number of bixels j. Then, the nominal dose d ∈ R V can be described as a linear function of the dose influence  Bortfeld (1997) (crosses) and the suggested analytical approximation (solid) composed by a weighted superposition of ten Gaussians (dashed). (b) Measured lateral profile for a 1 cm × 1 cm 6 MeV photon beam (crosses) at 50 mm depth and the suggested analytical approximation (solid) with two error functions (solid).
with the bixel weight vector w ∈ R B . For our derivations, we postulate that D i j factorizes into a depth part Z i j and a lateral part L i j , which in turn factorizes in x-and y-direction: 3.1.1. Protons. For a proton pencil beam, the lateral dose profile L can be parameterized with a Gaussian (Gottschalk et al 1993, Soukup et al 2005. Consequently, it can be integrated analytically against Gaussian probability densities. The Bragg curve describing the depth dose profile Z, in contrast, is a function of general nonlinear form (Bortfeld 1997) that cannot be integrated analytically against Gaussian probability densities. Hence, we suggest to model Z with a weighted superposition of multiple Gaussians, as (1) every Lipschitz-continuous function can be approximated arbitrarily well with a weighted sum of Gaussians (Schölkopf and Smola 2001) and (2) this particular choice for Z makes the subsequent integrals analytical. We found that ten Gaussians, each with a distinct weight ω jk , mean μ z jk , and width δ jk yield very good approximations for proton depth dose curves with ranges up to R max j = 350 mm, as exemplary shown in figure 1(a) for one selected range. The proton dose d p+ i can then be computed according to where μ x/y j denote the lateral positions of bixel j and x/y i j are the lateral coordinates of voxel i projected to the beam orientation of bixel j. The width of the lateral Gaussian proton beam profile λ i j is a function of the radiological depth z i j . For our simulations we calculate λ i j (z i j ) according to Gottschalk et al (1993) assuming initial beam widths of the proton beam ranging from 2.3 to 5.2 mm for beam energies from 230 to 70 MeV (Safai et al 2012). The described proton pencil beam algorithm corresponds to the implementation of Hong et al (1999) using a custom parameterization for Z.
3.1.2. Photons. For photons we use a sum of two error functions for the lateral parts L x/y and an exponential function for the depth part Z to describe the dose d γ i , as both error functions and exponential functions can be integrated analytically against Gaussians: x/y i j denote the lateral coordinates of voxel i projected to the beam orientation of bixel j. The integration limits a x/y i j and b x/y i j are given by the projection of the aperture outline onto the geometric depth z geo i j of voxel i from the virtual photon source of bixel j. λ i j is a function to the radiological depth z i j . We found that a linear model λ i j (z i j ) = τ · z i j + ζ is a good approximation to account for the blur-out of the photon penumbra. τ and ζ are found with a least squares fit of measured lateral dose profiles; A and c are found with a least squares fit of measured depth dose profiles. Comparisons of the analytical approximations to measured data for L x/y and Z are shown in figures 1(b) and 3(b), respectively.

Uncertainty model
The calculated nominal dose d does not correspond to the dose that is actually delivered during the treatment because of uncertainties regarding the input parameters of the dose calculation algorithm and uncertainties regarding the actual treatment scenario (Lomax 2008a(Lomax , 2008b. Among the most important sources of uncertainty are setup uncertainties, i.e., lateral offsets x , y of the patient relative to the beam, and-for particles-range uncertainties, i.e., offsets z of the calculated radiological depths z. It is possible to incorporate these offsets into the analytical expressions for the proton and photon dose: Note that we assume x ∈ R B , y ∈ R B , and z ∈ R B which corresponds to one offset in x, y, and z per bixel j. Alternative approaches are feasible within our framework, e.g. a voxel and bixel dependent range uncertainty z i j , but beyond the scope of this paper. The parameters x j , y j , and z j in equations (7) and (8) are subject to uncertainty. For our analytical probabilistic model, we assume that the underlying probability distribution p( x , y , z ) of the offsets x , y , and z corresponds to a product of multivariate normal distributions with the mean zero: (9) With the covariance matrices x ∈ R B×B , y ∈ R B×B , and z ∈ R B×B it is possible to explicitly specify the correlation between different uncertainties. E.g. setup uncertainties in x/y j of bixels belonging to the same beam orientation might be considered perfectly correlated, as they are applied almost simultaneously; range uncertainties z j of bixels belonging to the same beam and impinging at the same lateral position but with different energies may be considered perfectly correlated, as these bixels penetrate the same tissue on one ray. In these simple scenarios, the matrices x , y , and z exhibit block structure. With the standard deviation σ a/b of the uncertainty in bixel a/b we have ab = σ a σ b if the uncertainties in bixels a and b are perfectly correlated and ab = 0 if the uncertainties in bixels a and b are uncorrelated. More sophisticated correlation assumptions, such as distinct motion patterns of the patient, are feasible within our framework but not discussed in this paper.
In equation (9) we implicitly postulate that the uncertainties x , y , and z are independent from one another. This will allow for a clear structure of the following derivations but it is not a prerequisite for our calculations to remain analytically tractable.
The probability density p( x , y , z ) induces a probability density p(d) over the resulting dose distribution according to equations (7) and (8). p(d) will be non-Gaussian due to the sum over all bixels j and the nonlinear form of the lateral and depth dependent components of the dose. Our functional parameterizations of d p + and d γ , however, allow for the analytical calculation of all moments of the dose distribution, in particular of its expectation Furthermore, the components of the lateral parts L x/y and the depth part Z of the mean and (co)variance can be calculated separately as we postulated that the probability density of the considered uncertainties and the dose calculation factorize in x, y, and z and x , y , and z , respectively. The expectation value of the dose is given by For the covariance of the dose, which is defined as Hence, the calculation of the first moment of the probability density of the dose distribution is reduced to the calculation of L x/y i j and Z i j ; the calculation of the second moment is reduced to the calculation of the matrix elements ϒ x/y i jlm and i jlm . In the following sections 3.3 and 3.4 we derive analytical expressions for these quantities separately for protons and photons. When following our derivations, it is worthwhile to keep in mind that the indices i, l always refer to voxels, j, m always refer to bixels, and k, n always refer to the components of the ten Gaussians parameterizing the proton depth dose profile Z.

Analytical probabilistic proton dose calculation
3.3.1. Setup uncertainties. In order to calculate the expectation value of the proton dose incorporating setup uncertainties in x-direction 3 we use appendix A and appendix B to solve (13) x j j denotes the variance of the setup uncertainty for bixels j.μ andσ can be computed according to appendix B but their exact values are mute at this point as the integration over x j yields 1.
To compute the covariance we define and solve  The matrix elements ϒ i jlm correspond to a two-dimensional normal distribution whose covariance matrix is given by i jlm + jm . Figure 2(a) shows results of the analytical calculations in comparison to sampled data for a random lateral proton dose profile.

Range uncertainties.
In order to calculate the expectation value of the proton dose incorporating range uncertainties we have to solve To compute the covariance we define and perform analogous derivations as in equations (15) to obtain For proton range uncertainties, the matrix elements ϒ i jlm correspond to a quadratic form over two-dimensional normal distributions whose covariance matrices are given by jkmn + jm . Figure 2(b) shows results of the analytical calculations in comparison to sampled data for a random proton depth dose profile.
The histograms of the sampled dose values around x * show that the induced probability density over the dose p(d(x * )) is usually of general, non-Gaussian form, as explained in section 3.2. Due to the perfect correlation of the uncertainties of the underlying pencil beams, the dose maximum and minimum around x * are accumulation points in p(d(x * )).

Analytical probabilistic photon dose calculation
3.4.1. Setup uncertainties. In order to calculate the expectation value of the photon dose incorporating setup uncertainties we have to solve To compute the covariance we apply the vector notation introduced in equation (14), define and solve For photon setup uncertainties, the matrix element ϒ i jlm corresponds to an integral over a bivariate normal distribution which can be evaluated numerically (Genz 2004). Figure 3(a) shows results of the analytical calculations in comparison to sampled data for a random lateral photon profile. Again, we see that the induced probability density p(d(x * )) around x * is of general, non-Gaussian form.

Radiological depth uncertainties.
In order to calculate the expectation value of the photon dose incorporating uncertainties in the calculated radiological depth z i j we have to solve To compute the covariance we apply the vector notation introduced in equation (17), define the column vector c = (c, c), and solve In the last step we have used that (A −1 + I) −1 = A(A + I) −1 . For photons, the matrix element i jlm is given by the product of Z i j and Z lm multiplied with a correction factor e c 2 z jm that depends only on the covariance of the uncertainties of the bixels j and m. Figure 3(b) shows results of the analytical calculations in comparison to sampled data for a photon depth dose profile. For photons, uncertainties in the calculated radiological depth are negligible in the clinic yet still we wanted to include our derivations to demonstrate that our approach is flexible enough to accommodate any function that can be integrated against a multivariate normal distribution. In the remainder of this paper we do not account for uncertainties in the radiological depth for photons.

Preparation and execution errors
So far we have only considered one random error j which corresponds to a radiation treatment delivered in one fraction. In radiation therapy, however, we typically have multiple fractions with a systematic (preparation) error sys j which is the same for all fractions f = 1, . . . , F and a random (execution) error rand j f which is different for every fraction f (van Herk et al 2002). With APM it is possible to account for the different nature of preparation and execution errors. For demonstration purposes we consider a lateral proton dose profile delivered in F fractions. Incorporating fractionation effects, the factor L x/y i j from equation (4) becomes where we have dropped the superscript x. With the two probability distributions over sys j and rand j f we have to calculate L i j for the expectation value of equation (24) according to To calculate the second moment we have to solve With the standard variable definitions for x i jlm , μ jm , sys jm , rand j f mg , and i jlm the matrix element ϒ i jlm can be expressed as a sum of 2D Gaussians: jmS is given by For the definition jmR , however, we have to distinguish two different cases. While systematic errors are also correlated over all fractions, random errors are only correlated within the same fraction: (31) Consequently we have

Results
We implemented the analytical probabilistic dose calculation for two-dimensional test cases in MATLAB 4 . Figure 5 shows the nominal dose, the expectation value, and standard deviations assuming different fractionation schemes for a photon and proton treatment plan for a Cshaped target volume surrounding an organ at risk. Our calculations assume a systematic setup error of 1 mm and a random setup error of 2 mm for both protons and photons. A 3% range uncertainty (Yang et al 2012) is only considered for protons. Setup uncertainties of bixels impinging from the same beam orientation and range uncertainties of bixels impinging from same beam orientation at the same lateral position are considered to be perfectly correlated; other uncertainties are considered to be uncorrelated. The computation time is 75 s for the standard deviation of the proton plan (∼200 bixels, ∼2000 voxels) and 20 s for the standard deviation of the photon plan (∼100 bixels, ∼10 000 voxels).
The prescribed target dose is 60 Gy. For delivery in one fraction, we observe standard deviations of up to 8 Gy at the end of the range of the proton beam. For photons, we observe standard deviations of up to 9 Gy right after entering the patient for beams that penetrate a lot of healthy tissue before reaching the target volume. For delivery in 30 fractions, the standard deviation of the photon plan is significantly reduced because of the averaging out of the effects induced by random setup errors. For protons, however, the systematic range errors prevent an averaging out over multiple fractions.

Analytical probabilistic modeling for robust planning
We have shown that it is possible to calculate the first two moments, i.e., the expectation value and the covariance of the probability density p(d) of the dose distribution with an analytical d [mm] [  Figure 5. Nominal dose d (first row), expectation value E[d] (second row), and standard deviation σ for delivery in 1 (third row), 10 (forth row), and 30 fractions (fifth row) for 9 photon beams (left column) and 3 proton beams (right column). We assume a range uncertainty of 3%, a systematic setup error of 1 mm, and a random setup error of 2 mm. The color scale in the first row applies to the nominal dose and the expectation value and the color scale in the third row to the standard deviation.
probabilistic dose calculation. Consequently, p(d) can be approximated with a multivariate normal distribution with the mean E[d] and the covariance matrix d : Hence, the deviation of the dose from a prescribed dose d − d * also follows a multivariate normal distribution with the same covariance matrix but a shifted mean: With the penalty matrix P = diag(p 1 , . . . , p V ), the standard quadratic objective function F (Oelfke and Bortfeld 2001) If (1) d − d * follows a multivariate normal distribution with the mean vector E[d] − d * and a non-singular covariance matrix d and (2) P is a V × V non-negative definite matrix, then F follows a generalized χ 2 distribution with expectation value (Imhof 1961, Liu et al 2009) and variance (Imhof 1961, Liu et al 2009 The evaluation of equation (37) is computationally challenging because it requires the full V × V matrix d . However, equation (36), which corresponds to the objective function introduced by Unkelbach and Oelfke (2004) for robust planning, can be evaluated efficiently for different bixel weight vectors w. The second term corresponds to the standard objective function using the expectation value of the dose instead of the nominal dose, and the first term requires only the diagonal elements of d : Hence, equation (36) may be evaluated-and optimized-at almost the same cost as the standard quadratic objective function after the precalculation of the B × B matrix .

Discussion
This paper introduces APM to quantify uncertainties for radiation therapy treatment planning. We focus on analytical probabilistic dose calculations and derive closed-form expressions for the calculation of the expectation value and the covariance of dose distributions for intensitymodulated photon and proton radiation therapy. Thereby we explicitly account for the different nature of systematic and random errors in fractionated radiation therapy and accommodate complex correlation models of the underlying uncertainties. Currently these aspects are often neglected, especially for robust planning with worst case optimization (Unkelbach et al 2009, Pflugfelder et al 2008, Liu et al 2012. Figures 4 and 6, however, show that different correlation assumptions and an adequate handling of both random and systematic errors have a critical impact on the variance of dose distributions. APM is ideally suited for further investigations regarding the consequences of random and systematic errors as well as varying correlation models. [mm] [mm] Using more sophisticated correlation models it is also possible to simulate spatial motion of the patient for 4D treatment planning. If the uncertainty model is transferred from correlated uncertainties of bixel positions to correlated uncertainties of voxel positions even dose calculations on deforming geometries are analytically feasible.
As laid out in section 3, APM enables the efficient quantification of uncertainties in intensity-modulated dose distributions. In future, we want to analytically propagate these uncertainties into composite treatment plan quality indicators such as mean doses, maximum doses, dose volume histograms, equivalent uniform doses, normal tissue complication probabilities, or tumor control probabilities to improve the evidence for clinical decision making.
Naturally, the main application of APM will be in particle therapy but conventional photon radiation therapy may also benefit from fast computations of dose uncertainties. First, inverse planning with APM may ensure adequate coverage of the clinical target volume without definition of a planning target volume using direct constraints on the dose variance within the tumor. Just like conventional margin recipes (van Herk et al 2004), APM may consider random and systematic errors separately but it will also consider patient-specific factors such as the relative weight of certain bixels. Second, inverse planning with APM will also account for uncertainties in the dose applied to organs at risk. This may be important for photons where we observe pronounced uncertainties of the dose distributions right where the photon beams enter the patient, as shown in figure 5.
Beyond photons and protons, an analytical probabilistic dose calculation is also feasible for carbon ions where the fragmentation tail could be modeled with a linear function. For carbon ions, however, APM should include an analytical probabilistic survival calculation based on the linear quadratic model (Joiner and van der Kogel 2009, chapter 4). Here, it would be possible to account for uncertainties not only in the delivered dose D but also in the parameters α and β of the linear quadratic model. We see two main practical challenges for APM in radiation therapy treatment planning. First, the validity of the uncertainty calculations depends on the accuracy of the functional approximation of the treatment plan quality indicators. For the analytical probabilistic proton dose calculation as described in section 3.3, the shortcomings of conventional pencil beam algorithms compromise the accuracy in heterogeneous media. However, this limitation could be overcome with more elaborate dose calculation algorithms based on pencil beam decomposition (Soukup et al 2005) or an explicit measurement of the correlation of range and setup uncertainties with multiple ray tracings. For the analytical probabilistic photon dose calculation as described in section 3.4, the approximation of the lateral dose profile with error functions is good in the peak and penumbra region but poor in the low dose region, as shown in figure 1(b). This limitation could be overcome with an additional underlying broad Gaussian component in the lateral photon dose. As the approximation of the photon depth dose with an exponential does not include the build up effect and the dosimetric uncertainties induced by the uncertainties in the radiological depth were found to be negligible, we suggest to exclude the uncertainties in the radiological depth for photons. For both particles and photons a detailed comparison with more sophisticated dose calculation algorithms is part of ongoing research.
Second, APM may be computationally more demanding than nominal calculations. For the analytical probabilistic dose calculation, the most expensive computation is that of the covariance terms, which requires O(B 2 × V ) operations, while the computation of the expectation value only requires O(B×V ) operations-just like a conventional dose calculation. Currently, an analytical probabilistic dose calculation including variance calculation takes about 1 min for two-dimensional examples in MATLAB. The variance estimates do not come for free but APM may enable a substantial acceleration compared to sampling techniques, which require about the same time to simulate ∼10 scenarios in MATLAB due to multiple ray tracings.
The main advantage of APM regarding computational speed will manifest for probabilistic planning based on the quadratic objective function F, as explained in section 4. If we use APM to obtain estimates of the expectation value and covariance of the dose, we only need a single APM run before the optimization. Conventional approaches to robust planning, in contrast, use sampling to obtain estimates of the expectation value and covariance of the dose requiring 10 2 -10 3 samples (Sobotta et al 2012) in every iteration of the optimization process. While the estimates of the expectation value and the standard deviation stemming from sampling techniques are subject to statistical uncertainty, the expectation value and the standard deviation calculated with APM are exact; they are only compromised by the accuracy of the underlying model.
We acknowledge that the minimization of the expectation value of the objective function as outlined in section 5 does not necessarily guarantee the coverage of the tumor volume with a certain dose at a desired probability-which is the ratio behind margin concepts (van Herk et al 2004). However, using the functional form for the dose covariance as derived with APM, it may be feasible to incorporate such a constraint into probabilistic treatment planning approaches in the future. As a first step we could try to specify the target coverage probability with constraints on the expectation value of the dose minus a multiple of the standard deviation, or we could restrict the overall variability of the dose distribution within the target with constraints on the average standard deviation. In its current form the variance calculations can already serve as valuable guidance for clinicians to assess the robustness of the dose distribution.

Conclusion
Treatment plan quality indicators in radiation therapy are sensitive to uncertainties. We introduce the concept of analytical probabilistic modeling (APM) to quantify the influence of these uncertainties.
Based on a Gaussian uncertainty model and functional parameterizations that can be integrated analytically against Gaussian densities, it is possible to derive closed-form expressions for the expectation value and the covariance of intensity-modulated photon and proton dose distributions. Thereby we explicitly accommodate arbitrary correlation models of the uncertainties and account for the different nature of systematic and random uncertainties in fractionated radiation therapy.
Our calculations indicate that the correlation model and the interplay of random and systematic uncertainties have a critical impact on the overall uncertainty of the dose distribution.
Beneficial applications of APM in the context of treatment planning are feasible-in particular for robust planning with the standard quadratic objective function.

Appendix B. Products of multivariate normal distributions
The product of two multivariate normal distributions N 1 (x;x 1 , 1 ) and N 2 (x;x 2 , 2 ) is again of Gaussian shape (Rasmussen and Williams 2006, A.2), yet not normalized: