A dynamic system showing stable rhythmic activity can be represented by the dynamics of phase oscillators. This would provide a useful mathematical framework through which one can understand the system's dynamic properties. A recent study proposed a Bayesian approach capable of extracting the underlying phase dynamics directly from time-series data of a system showing rhythmic activity. Here we extended this method to spike data that otherwise provide only limited phase information. To determine how this method performs with spike data, we applied it to simulated spike data generated by a realistic neuronal network model. We then compared the estimated dynamics obtained based on the spike data with the dynamics theoretically derived from the model. The method successfully extracted the modeled phase dynamics, particularly the interaction function, when the amount of available data was sufficiently large. Furthermore, the method was able to infer synaptic connections based on the estimated interaction function. Thus, the method was found to be applicable to spike data and practical for understanding the dynamic properties of rhythmic neural systems.