Nonlinear model predictive control (NMPC) is an efficient approach for the control of nonlinear multivariable dynamic systems with constraints, which however requires an accurate plant model. Plant models can often be determined from first principles, parts of the model are however difficult to derive using physical laws alone. In this paper a hybrid Gaussian process (GP) first principles modeling scheme is proposed to overcome this issue, which exploits GPs to model the parts of the dynamic system that are difficult to describe using first principles. GPs not only give accurate predictions, but also quantify the residual uncertainty of this model. It is vital to account for this uncertainty in the control algorithm, to prevent constraint violations and performance deterioration. Monte Carlo samples of the GPs are generated offline to tighten constraints of the NMPC to ensure joint probabilistic constraint satisfaction online. Advantages of our method include fast online evaluation times, possibility to account for online learning alleviating conservativeness, and exploiting the flexibility of GPs and the data efficiency of first principle models. The algorithm is verified on a case study involving a challenging semi-batch bioreactor.