Today Time-of-Flight (TOF), in PET scanners, assumes a single, well-defined timing resolution for all events. However, recent BGO–Cherenkov detectors, combining prompt Cherenkov emission and the typical BGO scintillation, can sort events into multiple timing kernels, best described by the Gaussian mixture models. The number of Cherenkov photons detected per event impacts directly the detector time resolution and signal rise time, which can later be used to improve the coincidence timing resolution. This work presents a simulation toolkit which applies multiple timing spreads on the coincident events and an image reconstruction that incorporates this information. A full cylindrical BGO–Cherenkov PET model was compared, in terms of contrast recovery and contrast-to-noise ratio, against an LYSO model with a time resolution of 213 ps. Two reconstruction approaches for the mixture kernels were tested: 1) mixture Gaussian and 2) decomposed simple Gaussian kernels. The decomposed model used the exact mixture component applied during the simulation. Images reconstructed using mixture kernels provided similar mean value and less noise than the decomposed. However, typically, more iterations were needed. Similarly, the LYSO model, with a single TOF kernel, converged faster than the BGO–Cherenkov with multiple kernels. The results indicate that the model complexity slows down convergence. However, due to the higher sensitivity, the contrast-to-noise ratio was 26.4% better for the BGO model.