Natural deep eutectic solvents (NADES) are gaining significant attention in analytical chemistry due to attractive physico-chemical properties associated with sustainable aspects. They have been successfully evaluated in different fields, and applications in sample preparation have increased in the last years. However, there is a limited knowledge related to chemical interactions and mechanism of intermolecular action with specific analytes. In this regard, for the first time, this study exploited a computational investigation using molecular dynamics (MD) predictions combined with experimental data for the extraction/determination of steroidal hormones (estriol, β-estradiol, and estrone) in urine samples using NADES. The ultrasound-assisted liquid-liquid microextraction (UALLME) approach followed by high-performance liquid chromatography with diode array detection (HPLC-DAD) was employed using menthol:decanoic acid as extraction solvent. Experimental parameters were optimized through multivariate strategies, with the best conditions consisting of 3min of extraction, 150 μL of NADES, and 3mL of sample (tenfold diluted). According to molecular dynamics predictions confirmed by experimental data, a molar ratio that permitted the highest efficiency consisted of menthol:decanoic acid 2:1 v/v. Importantly, computational simulations revealed that van der Waals interactions were the most significant contributor to the interaction energy of analytes-NADES. Using the optimized conditions, limits of detection (LOD) ranged from 3 and 8μg L-1, and precision (n = 3) varied from 8 to 19%. Intraday precision was evaluated at 3 concentrations: low (LOQ according to each analyte), medium (100μg L-1), and high (750μg L-1). Accuracy was successfully assessed through recoveries that ranged from 82 to 98%. In this case, molecular dynamics simulations proved to be an important tool for in-depth investigations of interaction mechanisms of DES with different analytes.