Molecular analyses of closely related taxa have increasingly revealed the importance of 10 higher-order genetic interactions in explaining the observed pattern of reproductive isolation 11 between populations. Indeed, both empirical and theoretical studies have linked the process of 12 speciation to complex genetic interactions. Gene Regulatory Networks (GRNs) capture the 13 inter-dependencies of gene expression and encode information about an individual's phenotype 14 and development at the molecular level. As a result, GRNs can-in principle-evolve via natural 15 selection and play a role in non-selective, evolutionary forces. Here, we develop a network-based 16 model, termed the pathway framework, that considers GRNs as a functional representation of 17 coding sequences. We then simulated the dynamics of GRNs using a simple model that included 18 natural selection, genetic drift, and sexual reproduction and found that reproductive barriers can 19 develop rapidly between allopatric populations experiencing identical selection pressure. Further, 20 we show that alleles involved in reproductive isolation can predate the allopatric separation of 21 populations and that the number of interacting loci involved in genetic incompatibilities, i.e., the 22 order, is often high simply as a by-product of the networked structure of GRNs. Finally, we discuss 23 how results from the pathway framework are consistent with observed empirical patterns for 24 genes putatively involved in post-zygotic isolation. Taken together, this study adds support for the 25 central role of gene networks in speciation and in evolution more broadly. 26 27 31Through this work, it is widely accepted that divergent selection on de novo mutations in geo-32 graphically isolated populations can facilitate speciation, as originally theorized by (Bateson, 1909; 33 Dobzhansky, 1936; Muller, 1942). Despite well-established examples from Drosophila (Brideau34 et al., 2006), Xiphophorus (Wittbrodt et al., 1989; Powell et al., 2020), Oryza (Yamamoto et al., 35 2010), Arabidopsis (Bikard et al., 2009), and Mus (Davies et al., 2016), the genetics and evolutionary 36 history of incompatibilities are typically far more complex and/or less well understood than what 37 is suggested by classical models (Noor and Feder, 2006; Lowry et al., 2008; Presgraves, 2010; Wolf 38 et al., 2010; Nosil and Schluter, 2011; Seehausen et al., 2014; Marques et al., 2019; Dagilis and 39Matute, 2020). 40 Post-zygotic, genetic isolation is thought to occur due to epistatic interaction between loci, where 41 alleles arise and fix in allopatry prior to secondary contact, e.g., the Bateson-Dobzhansky-Muller 42 1 of 26 Manuscript submitted to eLife (BDM) model (Bateson, 1909; Dobzhansky, 1936; Muller, 1942). However, many incompatibilities 43 uncovered using high-throughput molecular analyses (Castillo and Barbash, 2017; Kuzmin et al., 44 2018; Vaid and Laitinen, 2019) and quantitative trait locus (QTL) mapping (Moyle and Nakazato, 45 2008; Turner et al., 2014; Ch...