Human carboxylesterases (CES) are enzymes that play an important role in the metabolism and biotransformation of diverse substances. The two more relevant isoforms, CES1A1 and CES2A1, catalyze the hydrolysis of numerous approved drugs and prodrugs. The elucidation of CES isoform substrates specificity constitutes a very relevant medicinal chemistry issue. The general role pointed that the selectivity towards CES1A1 or CES2A1 depends on the size of the acyl and alkyl moieties present in the structure of the substrate, but several exceptions regarding substrate promiscuity towards both CES have been reported. In this work, a combination of classical molecular dynamics (MD) and hybrid quantum mechanics/molecular mechanics (QM/MM) simulations were applied with the purpose of studying the substrate selectivity of CES1A1 and CES2A1 on two sets of selected ligands: p-nitrophenyl ester derivatives (NPE) and pyrethroid stereoisomers (Pyr). The classical molecular modeling studies showed that the van der Waals (VDW) component of interaction, with the hydrophobic residues present on CES1A1 and CES2A1 subpocket 1 and subpocket 2, showed a significant contribution to the substrates-CES affinity properties. The hybrid QM/MM simulations exhibited that the rate-limiting step for the studied substrates reactions were related to the transition state (TS) with the higher steric hindrance molecular structure. In conclusion, it was possible to observe that the studied substrates generate the best possible interaction pattern with the residues from subpocket 1 and 2 in order to produce the corresponding affinity constant with the enzyme. Then, this interaction pattern drives the catalytic turn-over reaction through the presence or absence of a high steric hindrance center in the molecular structure of the rate-limiting reaction.