Introgressive hybridization is an important evolutionary process and new analytical methods provide substantial power to detect and quantify it. In this study we use variation in the frequency of 657 AFLP fragments and DNA sequence variation from 15 genes to measure the extent of admixture and the direction of interspecific gene flow among three Heliconius butterfly species that diverged recently as a result of natural selection for Müllerian mimicry, and which continue to hybridize. Bayesian clustering based on AFLP genotypes correctly delineated the three species and identified four H. cydno, three H. pachinus, and three H. melpomene individuals that were of mixed ancestry. Gene genealogies revealed substantial shared DNA sequence variation among all three species and coalescent simulations based on the Isolation with Migration (IM) model pointed to interspecific gene flow as its cause. The IM simulations further indicated that interspecific gene flow was significantly asymmetrical, with greater gene flow from H. pachinus into H. cydno (2Nm = 4.326) than the reverse (2Nm = 0.502), and unidirectional gene flow from H. cydno and H. pachinus into H. melpomene (2Nm = 0.294 and 0.252, respectively). These asymmetries are in the directions expected based on the genetics of wing patterning and the probability that hybrids of various phenotypes will survive and reproduce in different mimetic environments. This empirical demonstration of extensive interspecific gene flow is in contrast to a previous study which found little evidence of gene flow between another pair of hybridizing Heliconius species, H. himera and H. erato, and it highlights the critical role of natural selection in maintaining species diversity. Furthermore, these results lend support to the hypotheses that phenotypic diversification in the genus Heliconius has been fueled by introgressive hybridization and that reinforcement has driven the evolution of assortative mate preferences.