As the repetition rates of ultra-high intensity lasers increase, simulations used for the prediction of experimental results may need to be augmented with machine learning to keep up. In this paper, the usage of Gaussian process regression in producing surrogate models of laser-plasma interactions from particle-in-cell (PIC) simulations is investigated. Such a model retains the characteristic behaviour of the simulations but allows for faster on-demand results and estimation of statistical noise. A demonstrative model of Bremsstrahlung emission by hot electrons from a femtosecond timescale laser pulse in the 10 20 − 10 23 Wcm − 2 intensity range is produced using 800 simulations of such a laser-solid interaction from 1D hybrid-PIC. While the simulations required 84 000 CPU-hours to generate, subsequent training occurs on the order of a minute on a single core and prediction takes only a fraction of a second. The model trained on this data is then compared against analytical expectations. The efficiency of training the model and its subsequent ability to distinguish types of noise within the data are analysed, and as a result error bounds on the model are defined.