Neuroscientific data exhibits complex, nonlinear dynamics, varying in time and over trials. A growing body of work has focused on capturing nonlinear dynamics, doing so in increasingly complex ways for modeling and inference. Meanwhile, capturing the dependency on experimental conditions surrounding the neural recordings remains understudied, and few methods leverage the correlation structure across experimental conditions. This poses difficulties for model estimation as well as any analyses comparing and interpolating between conditions. In this work, we develop the Gaussian Process Linear Dynamical System model class, which consists of a Linear Dynamical System (LDS) whose parameters vary smoothly as a function of conditions using Gaussian Processes. The model thus captures nonlinear dynamics while benefiting from the interpretability of time-varying LDS models, and we develop a closed-form parameter fitting methodology for fast Bayesian inference and learning. After validating the model and inference on a simulated dataset, we turned to the analysis of neural datasets. First, we analyzed recordings from the neural head direction (HD) system in mice foraging in open environments and investigated how the fixed points of the dynamics in the anterodorsal thalamic nucleus change with respect to head direction. Second, we examined dorsal premotor cortex recordings in macaques performing a center-out reaching task and highlight how the consistencies across reach angles can be leveraged to interpolate dynamics in unseen angles. These experiments showcase the model’s ease of deployment and ability to capture complex tasks with multiple covariates.