Manufacturing, aerospace, energy and several other industries have witnessed a steep growth of increasingly complex, information rich, devices and systems of devices requiring simulation-based approaches. In fact, most modern systems have such complex behavior that their performance can only be evaluated through, usually computationally expensive, simulations. In such settings, it is of paramount importance to provide solutions with quality guarantees. In this manuscript, we focus on algorithms capable of identifying a level set of solutions in proximity of the global optimum, and specifically on the Probabilistic Branch and Bound (PBnB) method. We propose a new way to automate branching decisions by coupling this method with Gaussian process (GP) estimation. The result is PBnB-GP, where, at each iteration a collection of GPs is used to decide how to branch the input space. PBnB-GP not only returns an estimate of the regions with near-optimal reward (using the power of PBnB), but also a 'collection of Gaussian processes' that can produce point estimations for any location in the input space, thus harnessing the power of model-based black-box optimization. We present PBnB-GP for the first time together with preliminary numerical results.