A long-standing premise in theoretical neuroscience holds that the complex interactions between neurons or, at a different systems level, between networks in the nervous system can be expressed in terms of an underlying system of nonlinear dynamics. A classical approach has been dynamical systems theory, whereby biophysical processes in the brain can be naturally formalized in terms of differential equations (ODEs) that may be capable of offering an explicit and interpretable account of how different biochemical and physiological processes in the brain evolve to produce the complex repertoire of network dynamics. In parallel with these traditional approaches, popular machine learning architectures with vastly more parameters (recurrent neural networks – RNNs, variational autoencoders – VAEs) have been successfully implemented to study underlying nonlinear neural dynamics despite their lack of interpretability and direct analogy to biological systems. In this study, we combine the power of deep learning methods like VAEs with the interpretability and parsimony of ODEs to study the network dynamics given by Calcium activity in the nematode C. elegans during movement. VAEs extract important features of the data and generate time series of the coefficients for a system of ODEs (including the control parameter) describing the first two principal components of the neural data. To discover the ODE, we use the SINDy (sparse identification for nonlinear dynamics) framework, whereby the coefficients output from the VAE are combined with a pre-selected library of polynomial basis functions, along with sparsity regularization for the coefficients, to ensure parsimony of the system. The result is an ODE with control that switches between phases of the nematode behavior and otherwise has constant coefficients throughout. The resulting parsimonious system matches important features of the data: two stable states corresponding to the forward and backward crawling behavior, along with switching trajectories that describe turning. This work describes a new method to characterize complex nonlinear dynamics with control that combines the power of deep learning tools and the interpretability of differential equations, with the potential to generalize to other organisms and behaviors.