  Efficient stochastic simulation algorithms are of paramount importance to the
study of spreading phenomena on complex networks. Using insights and analytical
results from network science, we discuss how the structure of contacts affects
the efficiency of current algorithms. We show that algorithms believed to
require $\mathcal{O}(\log N)$ or even $\mathcal{O}(1)$ operations per
update---where $N$ is the number of nodes---display instead a polynomial
scaling for networks that are either dense or sparse and heterogeneous. This
significantly affects the required computation time for simulations on large
networks. To circumvent the issue, we propose a node-based method combined with
a composition and rejection algorithm, a sampling scheme that has an
average-case complexity of $\mathcal{O} [\log(\log N)]$ per update for general
networks. This systematic approach is first set-up for Markovian dynamics, but
can also be adapted to a number of non-Markovian processes and can enhance
considerably the study of a wide range of dynamics on networks.
