Prediction of binding and unbinding kinetics from all-atom molecular simulations could one day be an important component in the toolkit for structure-based drug design. Unfortunately, in contrast to binding affinity prediction, algorithms for estimating on- and off-rates in relatively complicated systems remain relatively undeveloped. Here we demonstrate one approach that has been applied to successfully compute the entry and exit rates of small gas molecules into globular proteins. The method is an application of transition-path theory and involves performing milestoning MD simulations in cells forming a Voronoi tesselation of some chosen state-discriminating collective-variable space. Particular attention is paid to handshaking the simulations with a continuous approximation of the solute flux to a protein target surface as a function of ligand diffusivity and bulk concentration. The method is applicable beyond small gas molecules and extensions involving more complicated molecules will be discussed.