Understanding the intricate dynamics of brain activities necessitates models that incorporate causality and nonlinearity. Dynamic causal modelling (DCM) presents a statistical framework that embraces causal relationships among brain regions and their responses to experimental manipulations, such as stimulation. In this study, we perform Bayesian inference on a neurobiologically plausible generative model that simulates event-related potentials observed in magneto/encephalography data. This translates into probabilistic inference of latent and observed states of a system driven by input stimuli, described by a set of nonlinear ordinary differential equations (ODEs) and potentially correlated parameters. We provide a guideline for reliable inference in the presence of multimodality, which arises from parameter degeneracy, ultimately enhancing the predictive accuracy of neural dynamics. Solutions include optimizing the hyperparameters, leveraging initialization with prior information and employing weighted stacking based on predictive accuracy. Moreover, we implement the inference and conduct comprehensive model comparison in several probabilistic programming languages to streamline the process and benchmark their efficiency. Our investigation shows that model inversion in DCM extends beyond variational approximation frameworks, demonstrating the effectiveness of gradient-based Markov chain Monte Carlo methods. We illustrate the accuracy and efficiency of posterior estimation using a self-tuning variant of Hamiltonian Monte Carlo and the automatic Laplace approximation, effectively addressing parameter degeneracy challenges. This technical endeavour holds the potential to advance the inversion of state-space ODE models, and contribute to neuroscience research and applications in neuroimaging through automatic DCM.