The quasipotential function allows for comprehension and prediction of the escape mechanisms from metastable states in nonlinear dynamical systems. This function acts as a natural extension of the potential function for non-gradient systems and it unveils important properties such as the maximum likelihood transition paths, transition rates and expected exit times of the system. Here, we leverage on machine learning via the combination of two data-driven techniques, namely a neural network and a sparse regression algorithm, to obtain symbolic expressions of quasipotential functions. The key idea is first to determine an orthogonal decomposition of the vector field that governs the underlying dynamics using neural networks, then to interpret symbolically the downhill and circulatory components of the decomposition. These functions are regressed simultaneously with the addition of mathematical constraints. We show that our approach discovers a parsimonious quasipotential equation for an archetypal model with a known exact quasipotential and for the dynamics of a nanomechanical resonator. The analytical forms deliver direct access to the stability of the metastable states and predict rare events with significant computational advantages. Our data-driven approach is of interest for a wide range of applications in which to assess the fluctuating dynamics.