Just to add, I think the real way to handle this issue is by having a multiscale algorithm that will resolve low population species with exact jumps and high population species with tau-leaping, but we don’t have that implemented yet.
Someone had a similar issue recently and posted a callback approach to having both methods, so maybe you could give that a try: