The traditional approach of molecular biology research was on examining and collecting data on a single gene or a single reaction. However, recently, there has been much interest on the dynamics of gene regulatory networks (Klipp et al., 2005). We applied mathematical approach for modeling of gene network. The models depict the reaction kinetics of the constituent parts and the functions are ultimately made from basic principle of simple expressions derived from Michaelis-Menten enzymatic kinetics, and the functional forms are usually chosen as Hill functions that serve as an approximation for the real molecular dynamics (Klipp et al., 2005). These dynamics depend on many parameters and the parameters strongly influence the behavior of the resulting gene network. Thus, we used simulated annealing algorithm to calculate a high fitness and optimal parameters of the gene network. The simulated annealing algorithm is suitable for calculating many degree of freedom (Tomshine and Kaznessis, 2006), and is the closest analogy with the shifting balance theory of populations (Kirkpatrick et al., 1983). We developed 3 different models that have two genes and experience two different environments, and simulated to describe the behavior of evolutionary gene networks. From simulation, we could obtain a high fitness of each gene network model, and we could indicate how gene network is evolved in evolutionary time from tracks of parameters and a fitness. Also, we analyzed the relations of a high fitness and parameters. We think we can apply to design and optimize other gene network, and these findings are useful to analysis of the evolutionary gene network.