Radiotherapy is among the most popular modalities used for cancer treatment. Proton therapy is a promising kind of radiotherapy, which uses protons characteristic maximum dose deposition to a specific tissue depth, for precise tumor irradiation. However, in comparison to conventional photon radiotherapy, proton therapy is sensitive to more treatment uncertainties, like errors in proton range and geometrical errors. To account for these uncertainties robust optimization and robustness evaluation have been developed, to obtain and guarantee the robustness of the treatment plans against potential error. Robust optimization uses a fixed number of scenarios, for which the plan is optimized using the worst-case scenario. However, a proper weighting of the sampled scenarios, with the corresponding probabilities is missing using this approach. Probabilistic treatment planning can resolve this limitation, but it requires significant time resources. To evaluate the quality of a treatment plan, one must quantify the effect of errors in dose distributions. However, dose distributions calculations are time-consuming, therefore in this work, polynomial chaos expansion (PCE) is used, as a dose meta-model, for fast and advanced dose analysis. This model uses a series of expansion in terms of polynomials, to evaluate the dose distribution when different errors occur. The first aim of this work was to improve the PCE construction speed and accuracy. To build the model a fixed number of dose scenarios from a dose engine are required. Currently, PCE dose meta-model is constructed using dose distributions with a 1% Monte Carlo (MC) noise level. The trade-off between the noise level (larger noise level allows a faster model) and the model's accuracy, was investigated. The accuracy for models built with larger noise levels (2 & 3%) was compared, to conclude that the default value is the most efficient choice. Additionally, PCE built using the dose differences (between a scenario with no errors and a shifted scenario), was used to improve models' accuracy for complex anatomies. However, probably due to the larger impact of MC noise when dose differences (smaller dose input value) are considered, the PCE accuracy was not high as expected. The second aim of this work is to evaluate robustness and trade-offs made in treatment planning in a clinically robust neuro-oncological patient, and a clinically complex patient case. As mentioned, robust optimization does not consider the occurrence probability of uncertainties. Therefore the treatment plans might be over-conservative. Using PCE robustness evaluations, we compare treatment plans with different robustness settings. From the results, we concluded that a further reduction of the settings was possible for the clinically robust patient. Third, for five robust skull base patients, we investigated the trade-off between homogeneity of target dose, robustness, and dose to healthy tissue. PCE robustness evaluations were used for comparison of inhomogeneous (120% maximum target dose) and homogeneous (107% maximum target dose) treatment plans, for different robustness settings. From these planning approaches, no intrinsic difference in the degree of robustness was observed. However, the inhomogenous treatment plan resulted in more healthy tissue sparing overall, at the expense of homogeneity of the target dose. The final aim of this project is to use PCE robustness evaluations towards a clinically feasible iterative probabilistic treatment planning. We attempted to find a linear relationship between the optimal robustness settings which meet a probabilistic goal, and a scaling factor, using iterations of PCE robustness evaluations. For our method, the overall scaling factor "α" was used, which is based on the assumptions that the actual margin recipe is approximately linear, and that the real underlying robustness recipe is a scaled version of the photon therapy margin recipe (van Herk's formula). The conclusion was that for one iteration an overall linear scale factor allows for a substantial gain. However, linear modeling and one iteration do not suffice to find the overall robustness settings optimum.