@Luis can you point to any mathematical write ups of one of your models showing how they differ via the trait space from traditional Gillespie models? That might help in pointing you towards possible algorithmic improvements and/or seeing how to add your models into existing packages.
How many rates do you typically sum up each step of the algorithm to get the total rate? If it is many there are lots of possible optimized Gillespie methods one can use, typically designed around only updating individual rates that would have actually changed after the occurrence of a given jump.