As whole-genome protein interaction network datasets become available for a wide range of species, evolutionary biologists have the opportunity to address some of the unanswered questions surrounding the evolution of these complex systems. Protein interaction networks from divergent organisms may be compared to investigate how gene duplication, deletion and ‘re-wiring’ processes may have shaped the evolution of their contemporary structures [1,2]. However, current approaches to aligning observed networks from multiple species are generally lacking the phylogenetic context necessary for meaningful conclusions to be drawn regarding network evolution. Here we show how probabilistic modeling can provide a platform for the quantitative analysis of multiple protein interaction networks. We apply this technique to the reconstruction of ancestral networks for the bZIP family of transcription factors [3] and find that excellent agreement is obtained with an alternative, sequence-based method for the prediction of leucine zipper interactions [4]. Further analysis shows our probabilistic method to be significantly more robust to the presence of noise in the observed network data than a simple parsimony-based approach [5]. In addition, the integration of evidence over multiple species means that the same method may be used to improve the quality of noisy interaction data for extant species. This is the first time that ancestral states of a protein interaction network have been reconstructed using an explicit probabilistic model of network evolution. We anticipate that it will form the basis of more general methods for probing the evolutionary history of biochemical networks.