Change of gene regulation is an important factor of causing morphological difference. Recent studies show inter- and intra-specific phenotypic differentiation by mutations in cis-regulatory sequences such as enhancer and promoter. It is a long-standing problem in evolutionary biology to understand relative contributions of natural and adaptive changes to phenotypic evolution. Several methods are known to detect selection using cis-regulatory sequences, however they do not focus on changes of gene expression affecting directly the phenotypes. Therefore, we developed a method of detecting selection for gene expression level by implementing a thermodynamics-based model. This model is to predict gene expression driven by any DNA sequence, as a function of transcription factor concentrations and their DNA-binding specificities. Our method can detect selection as sequences accumulating impossible substitutions at random. We applied this method to segmentation network genes which expression level vary in blastoderm embryos among three fruit flies, Drosophila melanogaster, D. yakuba and D. pseudoobscura. These fruit flies also vary embryo morphology. We simulated whether these morphological variation is adaptive evolution of gene expression. Furthermore, we discuss the advantages and disadvantages of our method compared with other methods to detect selection.