This paper presents a new method for regression based on the evolution of a type of feed-forward neural networks whose basis function units are products of the inputs raised to real number power. These nodes are usually called product units. The main advantage of product units is their capacity for implementing higher order functions. Nevertheless, the training of product unit based networks poses several problems, since local learning algorithms are not suitable for these networks due to the existence of many local minima on the error surface. Moreover, it is unclear how to establish the structure of the network since, hitherto, all learning methods described in the literature deal only with parameter adjustment. In this paper, we propose a model of evolution of product unit based networks to overcome these difficulties. The proposed model evolves both the weights and the structure of these networks by means of an evolutionary programming algorithm. The performance of the model is evaluated in five widely used benchmark functions and a hard real-world problem of microbial growth modeling. Our evolutionary model is compared to a multistart technique combined with a Levenberg-Marquardt algorithm and shows better overall performance in the benchmark functions as well as the real-world problem.