Heterogeneous Markovian Susceptible-Infected-Susceptible (SIS) epidemics with a general infection rate matrix $\widetilde{A}$ are considered. Using a non-negative matrix factorization to approximate $\widetilde{A}$, we are able to identify when a metastable state can be expected, and that the metastable distribution, under certain conditions, will feature a normal distribution with known expectation and covariance. Furthermore, we model a heterogeneous Markovian SIS epidemic, that starts with a fraction of initially infected nodes different from that in the metastable state, by approximating its behaviour by a standard linear stochastic differential equation (SDE) in sufficiently high dimensions. By exploiting the knowledge of the covariance matrix from the SDE, we demonstrate significant accuracy improvements over the first-order mean-field approximation NIMFA.