A recent discovery of aromatase crystal structure triggered the efforts to design novel aromatase inhibitors for breast cancer therapy. While correlating docking scores with inhibitory potencies of known ligands, feeble robustness of scoring functions toward prediction was observed. This prompted us to develop new prediction models using stepwise regression analysis based on consensus of different docking and their scoring methods (GOLD, LIGANDFIT, and GLIDE). Quantitative structure-activity relationships were developed between the aromatase inhibitory activity (pIC(50) ) of flavonoid derivatives (n=39) and docking scores and docking descriptors. QSAR models have been validated internally [using leave-one-out cross-validated r(2)(cv) (LOO-Q2))] and externally to ensure the predictive capacity of the models. Model 2 [M2] developed using consensus of docking scores of scoring functions viz. ASP, potential of mean force and DOCK Score (r(2)(cv)=0.850, r(2) = 0.870, r(2)(pred) = 0.633, RMSE = 0.363 μm, r(2)(m(test)) =0.831, r(2)(m(overall)) =0.832) was found to be better in predicting aromatase inhibitory potency (pIC(50) ) compared to the Model 1 [M1] based on docking descriptors (r(2)(cv)= 0.848, r(2) = 0.825, r(2)(pred) =0.788, RMSE=0.421μm, r(2)(m(test)) =0.808, r(2)(m(overall)) =0.821). It has been observed that the natural flavonoids and their derivatives were less potent compared to these scaffolds with imidazolylmethyl substitution owing to the interaction of nitrogen atom of the imidazole ring toward the heme (Fe(3+) ) of the aromatase. Results confirm the potential of our methodology for the design of new potent non-steroidal aromatase inhibitors.