Abstract: Stochastic models of biochemical reaction networks are widely used to capture intrinsic noise in cellular systems. The typical formulation of these models are based on Markov processes for which there is extensive research on efficient simulation and inference. However, there are biological processes, such as gene transcription and translation, that introduce history dependent dynamics requiring non-Markovian processes to accurately capture the stochastic dynamics of the system. This greater realism comes with additional computational challenges for simulation and parameter inference. We develop efficient stochastic simulation algorithms for well-mixed non-Markovian stochastic biochemical reaction networks with delays that depend on system state and time. Our methods generalize the next reaction method and τ-leaping method to support arbitrary inter-event time distributions while preserving computational scalability. We also introduce a coupling scheme to generate exact non-Markovian sample paths that are positively correlated to an approximate non-Markovian τ-leaping sample path. This enables substantial computational gains for Bayesian inference of model parameters though multifidelity simulation-based inference schemes. We demonstrate the effectiveness of our approach on a gene regulation model with delayed auto-inhibition, showing substantial gains in both simulation accuracy and inference efficiency of two orders of magnitude. These results extend the practical applicability of non-Markovian models in systems biology and beyond.