Let $X=(X_{t},\,t \in T)$ be a non-homogeneous, continuous time Markov process with a finite state space $S=\{1,...,K\}$.

Let $\alpha_{i,j}(t)$ be the hazard rates of some $\varGamma$-distributed random variables.

My question:How can I simulate random paths of the Markov process with a transition intensity matrix which is built with the above hazard rates $\alpha_{i,j}(t)$?

For example: $S=\{1,2,3,4\}$ with given $\alpha_{i,j}(t)$ with $i,j \in S$.

Any idea is greatly appreciated.