Background: Tracking the changes of neural dynamics based on neuronal spiking activities is a critical step to understand the neurobiological basis of learning from behaving animals. These dynamical neurobiological processes associated with learning are also time-varying, which makes the modeling problem challenging. New method: We developed a novel multiwavelet-based time-varying generalized Laguerre-Volterra (TVGLV) modeling framework to study the time-varying neural dynamical systems using natural spike train data. By projecting the time-varying parameters in the TVGLV model onto a finite sequence of multiwavelet basis functions, the time-varying identification problem is converted into a time invariant linear-in-the-parameters one. An effective forward orthogonal regression (FOR) algorithm aided by mutual information (MI) criterion is then applied for the selection of significant model regressors or terms and the refinement of model structure. A generalized linear model fit approach is finally employed for parameter estimation from spike train data. Results: The proposed multiwavelet-based TVGLV approach is used to identify both synthetic input-output spike trains and spontaneous retinal spike train recordings. The proposed method gives excellent the performance of tracking either sharply or slowly changing parameters with high sensitivity and accuracy regardless of the a priori knowledge of spike trains, which these results indicate that the proposed method is shown to deal well with spike train data. Comparison with existing methods: The proposed multiwavelet-based TVGLV approach was compared with several state-of-art parametric estimation methods like the steepest descent point process filter (SDPPF) or Chebyshev polynomial expansion method. The conventional SDPPF algorithm, or SDPPF with B-splines wavelet expansion method was shown to have the poor performance of tracking the time varying system changes with the synthetic spike train data due to the slow convergence of the adaptive filter methods. Although the Chebyshev polynomial basis function method gave the good parametric estimation results, it requires prior parameter estimation. It was shown that the proposed multiwavelet-based TVGLV method can track the time-varying parameter changes rapidly and accurately. Conclusions: The multiwavelet-based TVGLV modeling framework developed in this paper can not only provide a computational modeling scheme for investigating such nonstationary properties, track more general forms of changes in time-varying neural dynamics, and but also may potentially be applied to investigate the spatial-temporal information underlying biomedical spiking signals. (C) 2016 Elsevier B.V. All rights reserved.