Phase sensitivity analysis is a powerful method for studying (asymptotically periodic) bursting neuron models. One popular way of capturing phase sensitivity is through the computation of isochrons-subsets of the state space that each converge to the same trajectory on the limit cycle. However, the computation of isochrons is notoriously difficult, especially for bursting neuron models. In [W. E. Sherwood and J. Guckenheimer, SIAM J. Appl. Dyn. Syst., 9 (2010), pp. 659-703], the phase sensitivity of the bursting Hindmarsh-Rose model is studied through the use of singular perturbation theory: cross sections of the isochrons of the full system are approximated by those of fast subsystems. In this paper, we complement the previous study, providing a detailed phase sensitivity analysis of the full (three-dimensional) system, including computations of the full (twodimensional) isochrons. To our knowledge, this is the first such computation for a bursting neuron model. This was made possible thanks to the numerical method recently proposed in [A. Mauroy and I. Mezíc, Chaos, 22 (2012), 033112]-relying on the spectral properties of the so-called Koopman operator-which is complemented with the use of adaptive quadtree and octree grids. The main result of the paper is to highlight the existence of a region of high phase sensitivity called the almost phaseless set and to completely characterize its geometry. In particular, our study reveals the existence of a subset of the almost phaseless set that is not predicted by singular perturbation theory (i.e., by the isochrons of fast subsystems). We also discuss how the almost phaseless set is related to empirically observed phenomena such as addition/deletion of spikes and to extrema of the phase response of the system. Finally, through the same numerical method, we show that an elliptic bursting model is characterized by a very high phase sensitivity and other remarkable properties.