Biochemical reactions are often modelled as discrete-state continuoustime stochastic processes evolving as memoryless Markov processes. However, in some cases, biochemical systems exhibit non-Markovian dynamics. We propose a methodology for building stochastic simulation algorithms which model accurately non-Markovian processes in some specific situations. Our methodology is based on and implemented in Concurrent Constraint Programming (CCP). Our technique allows us to randomly sample waiting times from probability density functions (PDFs) not necessarily distributed according to a negative exponential function. In this context, we discuss an important case-study in which the PDF for waiting times is inferred from single-molecule experiments. We show that, by relying on our methodology, it is possible to obtain accurate models of enzymatic reactions that, in specific cases, fit experimental data better than the corresponding markovian models.