We derive local sensitivities of statistical quantities of interest with respect to model parameters in dynamical systems.Our main contribution is the extension of adjoint-based a posteriori analysis for differential operators of generic dynamical systems acting on states to the Liouville operator acting on probability densities of the states. This results in theoretically rigorous estimates of sensitivity and error for a broad class of computed quantities of interest while propagating uncertainty through dynamical systems.We also derive Monte-Carlo type estimators to make these estimates computationally tractable using spatio-temporal normalizing flows and exploiting the hyperbolic nature of the Liouville equation. Three examples demonstrate our method.First, for verification of the theoretical results, we use a 2D linear dynamical system with an initial multivariate Gaussian density. Then, we apply our method to the challenging task of propagating uncertainty in a double attractor system to illustrate sensitivities in bimodal distributions. Finally, we show that our method can provide sensitivities with respect to the parameters of Neural Ordinary Differential Equations (here, in the context of classification).