We discuss the identification of genetic networks based on a class of boolean gene activation rules known as hierarchically canalizing functions. We introduce a class of kinetic models for the concentration of the proteins in the network built on a family of canalizing functions that has been shown to capture the vast majority of the known interaction networks. The simultaneous identification of the structure and of the parameters of the model from experimental data is addressed based on a mixed integer parametrization of the model class. The resulting regression problem is solved numerically via standard branch-and-bound techniques. The performance of the method is tested on simulated data generated by a simple model of Escherichia coli nutrient stress response.