Intercity traveling has been recognized as a leading cause for the continuation of the COVID-19 global pandemic. However, there lacks credible prediction of the spatiotemporal spread of COVID-19 with humans traveling between metropolitan areas. This study attempts to establish a novel framework to simulate human traveling and the spread of virus across an intercity population mobility network. A Markov process was introduced to capture the stochastic nature of travelers’ migration. A backward derivation algorithm was adopted and the Nelder-Mead simplex optimization method applied to overcome the limitation of existing deterministic epidemic models, including the difficulties in estimating the initial susceptible population and the optimal hyper-parameters required for simulation. We conducted two case studies with data from 24 cities in China and Italy. Our framework yielded state-of-the-art accuracy while being modular and scalable, indicating the addition of population mobility and stochasticity significantly improves prediction performance compared to using epidemic data alone. Moreover, our results revealed that transmission patterns of COVID-19 differ significantly with different population mobility, offering valuable information to the understanding of the correlation between traveling activities and COVID-19 transmission.