Biological systems are constantly under environmental and genetic changes, yet their phenotypes tend to remain invariant. It has been argued that robustness to changes arises as a buffering of the developmental process. Gene regulatory networks often exhibit robustness and are found at the core of the developmental process. Our main interest is to decipher the molecular mechanisms granting robustness to biological systems. Therefore, we adopt an in silico approach to modeling gene regulatory networks that integrates two levels: local sequence information and network architecture. Using a set of 10 TF-DNA complexes from the PDB as templates we modeled TF binding to all possible nucleotide sequences of 8bp in length by generating all 4^8 models of TF-DNA complexes. A distance dependent statistical pair-potential is used to determine the free energy of each modeled structure. For each TF we obtained free energy estimates for all possible interacting DNA binding sites. This allows an explicit sequence-level representation of upstream regulatory regions that determines the architecture of a gene regulatory network model. In the model, mutations at the DNA level cause changes on two levels: (a) at the local sequence level in individual motifs, e.g. by changing the binding affinity, and (b) at the network architecture level by creating and destroying motifs or binding sites, which results in dynamically rewiring network connections. The regulatory network model is used to determine gene expression dynamics while at a higher, population dynamics, level the networks undergo cycles of reproduction, mutation and selection. We use this multilevel hierarchical model to investigate the molecular mechanisms underlying the evolution of robustness in gene regulatory networks.