Stachyose is a typical prebiotic that can be utilized by the probiotic strain Bacillus licheniformis. Pioneering X-ray crystallography has determined the structure of stachyose in complex with the solute-binding protein MsmE in B. licheniformis (BIMsmE). The present work describes a combined strategy for the identification of putative BIMsmE-specific ligands, which can be used for the development of prebiotics. After a ligand-based virtual similarity screening of a large ZINC database containing similar to 22 M compounds, we identified 3575 ligands. A total of 600 structures for which the Tanimoto coefficient's value was larger than a cutoff of 0.23 were selected for molecular docking. Based on the docking scores, we identified 100 top-scoring ligands, followed by molecular dynamics (MD) simulations. During simulations, 35 candidates were abandoned because of serious steric clashes in the complexes. Finally, the top 10 ligands with free energies below an energy threshold of -50.84 kcal/mol were selected. The top two ligands were stachyose and raffinose, which have proved their health benefits as prebiotics and their safety. The remaining eight ligands were further analyzed by the in silico ADME tool; only galactinol did not violate any of the criteria required for a lead compound. These three ligands were further analyzed for understanding their binding to BIMsmE. Isothermal titration calorimetry analysis suggested that stachyose, raffinose, and galactinol bound strongly to BIMsmE with K-d values of 299, 170, and 134 nM, respectively. Microsecond MD simulations suggested significant conformational changes of BIMsmE upon ligand binding. Our results provide new insight into the thermodynamics of sugars and MsmE, which would promote the development of novel prebiotics.