We present an evolutionary modeling method for modeling genetic regulatory networks. The method features a hybrid algorithm of genetic programming with statistical analysis to derive systems of differential equations. Genetic programming and the least mean squares method were combined to identify a concise form of regulation between the variables from a given set of time series. Results of multiple runs were statistically analyzed to indicate the term with robust and significant influence. Our approach was evaluated in artificial data and real world data.