Arbitrarily complex python code can be wrapped into pytensor Ops for inclusion in PyMC models, see here for an example focusing on likelihoods, or perhaps here for an example that wraps a scipy optimization routine.
Depending on the nature of the simulator, you might even be able to express it in pure pytensor, which would allow you to differentiate through the simulations (and compile it to numba/JAX for accelerated sampling). Some comments (and maybe a similar problem?) here.