Is this adaption of the Gillespie algorithm using Michaelis constants justifiable?
The stochastic simulation algorithm by Gillespie makes the assumption that reaction propensity (probability) in an infinitesimally small interval of time is same as the macroscopic reaction rate. (For details, you can see my answer in Chemistry.SE).

You can model and simulate every step of enzyme catalysis or you can also use the Michaelis-Menten model. In any case, you can extend the same approximations to the stochastic model.

This is also applicable to other Michaelis-Menten like (or Hill) functions of the form $$\frac{x^n}{k^n+x^n} \qquad \text{or}\qquad \frac{k^n}{k^n+x^n}$$. 

There are many papers that have done this. Note that the non-linear nature of these functions results in deviation of the stochastic behaviour from the deterministic behaviour. Read about stochastic focussing. Modelling each step exactly (binding, unbinding and catalysis) may result in different behaviour depending on your parameters. Whether you should use Michaelis-Menten function to describe reaction rate totally depends on how reasonable your approximation is.