Objective: We present a novel modeling framework for identifying time-varying (TV) couplings between time-series of biomedical relevance. Methods: The proposed methodology is based on multivariate autoregressive (MVAR) models, which have been extensively used to study couplings between biosignals. Contrary to the standard estimation methods that assume time-invariant relationships, we propose a modified recursive Kalman filter (KF) to track changes in the model parameters. We perform model order selection and hyperparameter optimization simultaneously using Genetic Algorithms, greatly improving accuracy and computation time. In addition, we address the effect of residual heteroscedasticity, possibly associated with event-related changes or phase transitions during a given experimental protocol, on the TV-MVAR coupling measures by using Generalized Autoregressive Conditional Heteroskedasticity (GARCH) models to fit the TV-MVAR residuals. Results: Using simulated data, we show that the proposed framework yields more accurate parameter estimates compared to the conventional KF, particularly when the true system parameters exhibit different rate of variations over time. Furthermore, by accounting for heteroskedasticity, we obtain more accurate estimates of the strength and directionality of the underlying couplings. We also use our approach to investigate TV hemodynamic interactions during exercise in young and old healthy adults, as well as individuals with chronic stroke. We extract TV coupling patterns that reflect well known exercise-induced effects on the underlying regulatory mechanisms with excellent time resolution. Conclusion and Significance: The proposed modeling framework can be used to efficiently quantify TV couplings between biosignals. It is fully automated and does not require prior knowledge of the system TV characteristics.