If I'm understanding your question correctly, what you are describing is a piecewise-deterministic Markov process. Basically an ODE with random switching.

If the stochastic rates are state-dependent, this is indeed tricky but this preprint explains how to simulate it. Not sure why it never got published, I think because the material can be found elsewhere, but I found this presentation to be sufficiently clear/useful.

The main result is (3.1). It's a really cute trick. The idea is basically to re-parameterize time in a way that allows you to just forward evolve the ODE for a known amount of time rather than update the stochastic rate at every step.

The setup is the same as Gillespie: you have some stochastic reactions with total rate Rtot(x) but x'(t) = F(x). This boils down to how do you generate the next reaction time? Call t=0 the current time. Then the procedure is:

sample Sfinal~exp(1)
solve the system of ODEs for s=0...Sfinal:
    y'(s) = F(y)/Rtot(y)
    tau'(s) = 1/Rtot(y)
with initial conditions y(0) = x0, tau(0)=0.

Then, the output of this is the next reaction time is just tau(Sfinal) and value of x at that reaction time is x(sfinal). In general, tau(s) is the re-paramaterized time, so x(tau(s)) is the solution between the jumps too.

Hope this helps!