Objective. 4D dose reconstruction in proton therapy with pencil beam scanning (PBS) typically relies on a single pre-treatment 4DCT (p4DCT). However, breathing motion during the fractionated treatment can vary considerably in both amplitude and frequency. We present a novel 4D dose reconstruction method combining delivery log files with patient-specific motion models, to account for the dosimetric effect of intra- and inter-fractional breathing variability. Approach. Correlation between an external breathing surrogate and anatomical deformations of the p4DCT is established using principal component analysis. Using motion trajectories of a surface marker acquired during the dose delivery by an optical tracking system, deformable motion fields are retrospectively reconstructed and used to generate time-resolved synthetic 4DCTs (‘5DCTs’) by warping a reference CT. For three abdominal/thoracic patients, treated with respiratory gating and rescanning, example fraction doses were reconstructed using the resulting 5DCTs and delivery log files. The motion model was validated beforehand using leave-one-out cross-validation (LOOCV) with subsequent 4D dose evaluations. Moreover, besides fractional motion, fractional anatomical changes were incorporated as proof of concept. Main results. For motion model validation, the comparison of 4D dose distributions for the original 4DCT and predicted LOOCV resulted in 3%/3 mm gamma pass rates above 96.2%. Prospective gating simulations on the p4DCT can overestimate the target dose coverage V_95% by up to 2.1% compared to 4D dose reconstruction based on observed surrogate trajectories. Nevertheless, for the studied clinical cases treated with respiratory-gating and rescanning, an acceptable target coverage was maintained with V_95% remaining above 98.8% for all studied fractions. For these gated treatments, larger dosimetric differences occurred due to CT changes than due to breathing variations. Significance. To gain a better estimate of the delivered dose, a retrospective 4D dose reconstruction workflow based on motion data acquired during PBS proton treatments was implemented and validated, thus considering both intra- and inter-fractional motion and anatomy changes.