A Family of Fitness Landscapes Modeled through Gene Regulatory Networks

Fitness landscapes are a powerful metaphor for understanding the evolution of biological systems. These landscapes describe how genotypes are connected to each other through mutation and related through fitness. Empirical studies of fitness landscapes have increasingly revealed conserved topographical features across diverse taxa, e.g., the accessibility of genotypes and “ruggedness”. As a result, theoretical studies are needed to investigate how evolution proceeds on fitness landscapes with such conserved features. Here, we develop and study a model of evolution on fitness landscapes using the lens of Gene Regulatory Networks (GRNs), where the regulatory products are computed from multiple genes and collectively treated as phenotypes. With the assumption that regulation is a binary process, we prove the existence of empirically observed, topographical features such as accessibility and connectivity. We further show that these results hold across arbitrary fitness functions and that a trade-off between accessibility and ruggedness need not exist. Then, using graph theory and a coarse-graining approach, we deduce a mesoscopic structure underlying GRN fitness landscapes where the information necessary to predict a population’s evolutionary trajectory is retained with minimal complexity. Using this coarse-graining, we develop a bottom-up algorithm to construct such mesoscopic backbones, which does not require computing the genotype network and is therefore far more efficient than brute-force approaches. Altogether, this work provides mathematical results of high-dimensional fitness landscapes and a path toward connecting theory to empirical studies.


Introduction
Since its introduction by Wright [1], the concept of fitness landscapes has grown and matured into a cornerstone of biology [2][3][4]. A fitness landscape consists of a space of genotypes that are mutually accessible through mutations and a fitness value associated with the phenotype each genotype encodes. In this context, fitness describes the evolutionary potential of each genotype, and the set of navigable genotypes on these landscapes is termed the genotype network [5]. Continuing with this metaphor, the evolution of a population can be depicted as a trajectory wandering on the fitness landscape. As a consequence, the topography of a fitness landscape sheds light on various evolutionary processes, including constraints on adaptation [6][7][8][9], speciation via genetic incompatibilities [10,11], (dis)advantages of sexual reproduction and recombination [12][13][14], the repeatability/reversibility (or not) of evolutionary trajectories [15][16][17][18], and the role of neutral networks-components of the genotype network with the same fitness-in epochal evolution [19][20][21][22][23][24].
Despite being introduced by Wright [1], Fisher's 1930 geometric model of adaptation is the first mathematical model of evolution on what we now call fitness landscapes [15,25,26]. Later work by Kingman [27] and Kauffman and Levin [28] constructed what they termed a "house of cards" (HoC) model where fitness values for each genotype are drawn independently from a specified probability distribution. Building on the HoC model, Kauffman and Weinberger [29] introduced the NK model, which forces each locus to interact with a fixed number of other loci and where a genotype's fitness becomes the sum of the fitness contributions of every interaction group. More recently, the "rough Mount Fuji model" [30,31] combines the HoC landscape with an additional field penalizing a genotype's Hamming distance away from a referenced genotype with the optimal fitness. The dependence of a genotype's fitness on that of neighboring genotypes is thought to be a key feature of empirical fitness landscapes.
Over the past three decades, the fitness landscapes for various organisms, including bacteria [32][33][34], fungi [35,36], and fruit flies [37], have been empirically reconstructed. While the number of genotypes included in these early landscapes was limited, modern sequencing techniques and high-throughput analyses have enabled the construction of many large landscapes. Notable studies have been conducted in HIV [38,39], yeast [40], E. coli [41], jellyfish [42], human cancers [43], human stem cells [44], and DNA/RNA networks [45][46][47][48][49]. Comprehensive landscapes for multiple eukaryotic species have also been analyzed based on the binding affinity of transcription factors [50] and after accounting for the ecological context the species experiences [51]. What emerged from these studies is a set of prominent topographical features conserved across diverged taxa [52].
Together, empirical and modeled fitness landscapes exhibit three key topographical features. First, fitness landscapes are more often "rugged" than smooth [52]. The degree of ruggedness can be assessed via a variety of measures, such as the roughness of the slope ratio [53,54] and the number of local fitness maxima [55], which are often strongly correlated with each other [4]. Empirical studies typically show moderate ruggedness in the observed fitness landscapes [32][33][34]36,50,56]. The degree of ruggedness in these empirical landscapes is less than the HoC model assumes and comparable to a fine-tuned NK model or rough Mount Fuji model [4]. Second, fitness landscapes reveal mutational trajectories from one genotype to another where the fitness is non-decreasing, which implies accessibility (typically to a fitness optimum) across the landscape [32,[57][58][59][60]. Lastly, whereas the inaccessible region in the HoC model expands when distant from the fitness optimum [61], other models find accessible trajectories despite high genotypic dimensionality [62,63].
Due to the often pervasive interaction between loci, determining phenotype from genotype can have a high degree of computational complexity [64,65]. Many existing fitness landscape models have dealt with this complexity by strongly constraining the state-space of possible genotypic interactions and/or reducing the complexity of how information is processed when mapping genotype to phenotype. For example, studies have focused on the folded structure of short RNA sequences, where the resulting stability or affinity is a fitness proxy [66][67][68][69][70], networks of molecular/genetic pathways whose expression pattern or homeostasis determines fitness [71][72][73][74][75], and modular mutational effects at different loci in Fisher's geometric model [76].
Here, we model the genotype-phenotype map using the pathway framework of gene regulatory networks (GRNs), where mechanistic knowledge of how phenotypes are computed from genotypes is encoded in the GRN (see [77,78] for a more formal introduction). To study the fitness landscapes induced by GRN evolution, we integrate the pathway framework into a family of fitness landscape models where the fitness value is uniquely determined by the phenotype corresponding to the regulatory outcome of a genotype. For a fitness landscape of GRNs, we first prove the existence of two key topographical features: (a) GRNs with the same phenotype are themselves connected in the underlying genotypic network, and (b) there exists accessible trajectories between all pairs of GRNs with similar phenotypes. Second, utilizing the idea of symmetries and automorphisms in the genotype network, we coarse-grain GRNs into groups with equivalent roles in the fitness landscapes and deduce an underlying mesoscopic structure with which we can predict the trajectory of evolution with minimal complexity. Lastly, using this coarse-graining, we develop a bottom-up algorithm for constructing the underlying fitness landscape of GRNs, which does not require computing the genotype network and is thus more efficient than the conventional brute-force approach.

Methods
Here, we introduce a family of fitness landscape models where the genotype-phenotype mapping is constructed from regulatory interactions. We first summarize a modeling framework of GRNs proposed in our previous work [77,78], termed the pathway framework, and then fitness landscape models of GRNs built upon the pathway framework.

Pathway Framework of GRNs
Genotypes in the pathway framework of GRNs contain all necessary information to construct a regulatory network [77,78]. More specifically, alleles at each locus include both a transcription activator and a protein product, which means the regulatory interactions among the loci can be deduced by connecting genes whose expression product corresponds to the activator of another. Compared to existing work on regulatory circuits-where mutations are modeled as rewiring a single interaction between genes [79,80]-the pathway framework considers a mutation as changing the activator/product of a gene. Lastly, the phenotype is determined by the set of loci reached in a regulatory cascade induced by external stimuli. These stimuli could be completely external to the individual or simply come from another regulatory network in the organism. For additional details on the pathway framework, see [77,78] and Figure 1 for illustration.  [77,78]. Under our four simplified assumptions, a GRN (genotype) consists of a fixed number of proteins as nodes and a constant number of directed edges depicting the activator/product pairs of genes. The phenotype is modeled as the Boolean states of proteins (colored), which are determined by their reachability from the external stimulus (lightning icon).
In this work, when building the family of fitness landscape models, we restrict the pathway framework with four assumptions. First, we consider a fixed set of genes underlying the genotypes, i.e., gene duplication and deletion events are excluded. Second, we assume a fixed underlying collection of proteins that can possibly exist in the organism. Third, we consider the case where a gene's expression is activated by a specific protein, and it generates only one protein product. Fourth, we assume that the associated chemical state of each protein is modeled as a Boolean/binary variable (present or absent), and external environmental signals stimulate the existence of specific proteins in the organism. As a consequence, the Boolean state of a phenotype-related protein is determined by whether it is reached by a regulatory cascade starting from an initial stimulus.
While the above assumptions seem naive, as we will show in Section 3, this simplified model still predicts the topographical features observed in empirical landscapes (see Section 1). As a result of these assumptions, we are able to derive rigorous theoretical insights into GRN evolution and obtain fitness landscapes consistent with far more complicated models. We believe these assumptions are conservative with respect to the biology and a justified starting point for modeling fitness landscapes, and we discuss the implications of these assumptions and possible extensions to the model in Section 4.

Fitness Landscape of GRNs under the Pathway Framework
Let Γ and Ω be the fixed, underlying collection of loci and proteins, respectively. A genotype is represented by its GRN g such that every locus γ ∈ Γ is associated with a protein activator/product pair e g (γ) = (u, v), u, v ∈ Ω. Equivalently, any GRN is a directed graph with |Ω| nodes labeled by the proteins Ω and |Γ| edges labeled by the loci Γ. In the rest of this paper, we will use the terminology "source/target" node of edge γ interchangeably to refer to the protein activator/product of locus γ. We also write G to be the set of all GRNs with the underlying loci Γ and proteins Ω.
The backbone of a fitness landscape of GRNs, i.e., the genotype network, is an undirected network of networks encoding the mutational relationship between the GRNs. Let G be the genotype network, and we denote its mega-nodes by V(G) = G and its edges by E(G). There is an edge (g 1 , g 2 ) ∈ E(G) between any GRNs g 1 , g 2 ∈ G when they only differ by the allele of a single locus γ, e g 1 (γ) = e g 2 (γ). In other words, g 1 and g 2 are connected in G when they can be transformed into each other through one edge rewiring.
Furthermore, we write x ω to be the binary state of protein ω ∈ Ω, where x ω = 1 indicates the presence of ω, and x ω = 0 designates its absence. We also partition Ω into three disjoint groups: (a) proteins Ω 0 whose presence is externally stimulated by the given environment, (b) proteinsΩ whose states influence the fitness value, and (c) the remaining ones, which we call the dummy proteins Ω since their specific identities are irrelevant to the external environment and the resultant phenotypes/fitness. (In this paper, we assume that the stimuli Ω 0 must be proteins that cannot be produced by expression, and we leave no constraint to the fitness-relevant and dummy proteinsΩ and Ω ).
A phenotype is then treated as a vector of zeros and ones, where each entry corresponds to the binary state of a protein inΩ. The resultant phenotype xΩ(g) of a GRN g is determined by the reachability in g: For any ω ∈Ω, x ω = 1 if and only if there is a stimulus ω 0 ∈ Ω 0 and a path from ω 0 to ω in g, which represents a chain of sequentially expressed genes that generates protein ω. Finally, the fitness f is simply a function of the phenotype xΩ(g).
Combined, a fitness landscape of GRNs is characterized by three key elements: the genotype network G, the external stimuli Ω 0 , and the fitness function of phenotype f (which implicitly identifies the fitness-relevant proteinsΩ). The genotype network G serves as the skeleton of the fitness landscape, whereas the environment-dependent stimuli Ω 0 and fitness function f determine the phenotypes of GRNs and their selective advantages.

Results
In this work, we derive three theoretical insights into fitness landscape models using GRNs as the embedded genotype-phenotype mapping. First, we show that the resulting family of fitness landscapes must always contain two topographical properties: connectivity, i.e., GRNs with the same phenotype can be mutually reached via mutations, and accessibility, i.e., that any GRN can be reached from an arbitrary less-fit GRN (once certain similarity criterion is met). Second, we propose a mesoscopic coarse-graining for fitness landscapes, which is a more compact alternative to analyzing evolutionary processes than the original landscape. This mesoscopic backbone recognizes "symmetries" in the genotype network, and it aggregates GRNs with the same role in the fitness landscape into a single representative genotype. Third, we provide a bottom-up approach to algorithmically construct this mesoscopic backbone and demonstrate its efficiency over coarse-graining the genotype network using brute force.

Connectivity and Accessibility in a Fitness Landscape of GRNs
A fitness landscape model of GRNs features a handful of properties that have either been discovered in empirical fitness landscapes or investigated mathematically. First, its underlying space, i.e., the genotype network G, presents immense dimensionality. Second, the fitness function f is flexible and can effectively tune the ruggedness of the fitness landscape. For example, a highly rugged "holely" landscape can be modeled by a binary f such that any GRN g ∈ G has high fitness once some single protein ω ∈Ω is present, x ω = 1, and otherwise, g has low/zero fitness. Because one can always find several mutational neighbors of g whose phenotype shows an opposite state x ω , the resultant fitness landscape is inevitably rugged. In what follows, we further show that fitness landscape models of GRNs must hold the characteristics of connectivity and accessibility.
Let y be a phenotype and denote by G y the set of all GRNs with phenotype y, i.e., xΩ(g) = y for g ∈ G y , under the given external stimuli Ω 0 . We also write Ω + y to be the required-present proteins in the phenotype y, so x ω = 1 for ω ∈ Ω + y and x ω = 0 for any other ω ∈Ω \ Ω + y . Note that the number of required-present proteins |Ω + y | is bounded from above by the number of loci |Γ| since any present protein that is not a stimulus must be triggered by the expression of some locus.
We observe that some GRNsG y ⊂ G y play a "central" role among GRNs with the same phenotype y. Specifically, for anyg ∈G y , all the edges ing point from the stimuli Ω 0 to the required-present proteins Ω + y , and each ω ∈ Ω + y is targeted by at least one edge ing. We demonstrate an example of suchg in Figure 2a. TheseG y are deemed central because they can be reached by any GRN g ∈ G y through mutations among G y themselves: First, for every edge in g that points to an ω ∈ Ω + y , we rewire the edge such that it still points to ω but now from an ω 0 ∈ Ω 0 . Arbitrarily rewiring the remaining edges between Ω 0 and Ω + y then leads to some central GRN inG y (see Figure 2a). In addition, if the phenotype y has strictly less required proteins than the number of loci, the central GRNs G y are mutually reachable by edge rewiring amongG y . There is always a redundant edge whose rewiring makes no change to the phenotype, and it helps us rewire each edge to any desired source/target pair between Ω 0 and Ω + y (see Figure 2b), which subsequently creates a chain of mutations between anyg,g ∈G y . These results implicate that, for any phenotype y with |Ω + y | < |Γ| and any g 1 , g 2 ∈ G y , there is always a mutational trajectory between g 1 and g 2 that only traverses over GRNs in G y , especially through the central onesG y (see Figure 2c). In the extreme case where |Ω + y | = |Γ|, however, G y fragments into multiple connected components (detailed in Appendix A).
Next, we turn to accessibility between GRNs of different phenotypes y and y , where without loss of generality f (y ) ≥ f (y). We observe that, if |Ω + y ∪ Ω + y | ≤ |Γ| + 1, there are always two "peripheral" GRNsĝ ∈ G y andĝ ∈ G y , which only differ by one edge rewiring. To be more specific, there are two independent chains inĝ, one of which begins with a stimulus ω 0 ∈ Ω 0 and sequentially connects the proteins required to be present in y but not in y , i.e., Ω + y \ Ω + y , while the other consecutively joins Ω + y \ Ω + y . The rest of the edges inĝ merely point from Ω 0 to Ω + y ∩ Ω + y , and each ω ∈ Ω + y ∩ Ω + y is targeted by at least one edge (see example in Figure 3, left). The other GRNĝ only differ from g by the first edge in the chain of Ω + y \ Ω + y , which is rewired such that it points from the stimulus ω 0 to the first node in the chain of Ω + y \ Ω + y (Figure 3, right). Example of peripheral GRNs connecting two different phenotypes-a peripheral GRN g of phenotype y in this example. There is a chain that triggers the presence state of proteins Ω + y \ Ω + y = {3, 4}. However, the other peripheral GRNĝ of phenotype y contains a chain of proteins Ω + y \ Ω + y = {5, 6}.ĝ andĝ are mutational neighbors since they only differ by rewiring the dark green edge, i.e., the first edge in either chain.
Our observation suggests that there is a sequence of mutations with non-decreasing fitness from any GRN g ∈ G y to any GRN g ∈ G y , as long as |Ω + y ∪ Ω + y | ≤ |Γ| + 1. In particular, when |Ω + y | < |Γ|, the mutational trajectory starting at g first traverses within G y to a peripheral GRN and then transitions into G y to reach g . An analogous trajectory exists even under the extreme scenario |Ω + y | = |Γ| (see Appendix A). We also note that if the number of fitness-relevant proteins is |Ω| ≤ |Γ| + 1, then the condition |Ω + y ∪ Ω + y | ≤ |Γ| + 1 is assuredly satisfied for any two phenotypes y and y . As a corollary, if |Ω| ≤ |Γ| + 1, the fitness optimum will always be accessible.

Mesoscopic Skeleton Derived from "Symmetries" in the Genotype Network of GRNs
Because the number of possible GRNs grows super-exponentially as the underlying loci and proteins expand, constructing the genotype network becomes extremely challenging beyond a small Γ and Ω. Here, we present a more compact skeleton of the fitness landscape of GRNs based on "symmetries" in the genotype network.
As the underlying space of a fitness landscape of GRNs, the genotype network G appears to contain redundant information. On the one hand, GRNs leading to the identical phenotype are deemed to have equal fitness. On the other hand, given any GRN, for example, the mega-node rounded by orange in Figure 4a, one can always find some other GRN such that their neighborhoods in G are locally similar, e.g., the mega-node rounded by blue. This simple demonstration suggests that the structure of the genotype network G is not arbitrary; instead, some structural symmetries exist. Figure 4. The genotype network has symmetry such that multiple GRNs have similar local neighborhoods, as we demonstrate in (a) since the corresponding GRNs only differ by exchanging the role of loci A and B. More formally, these GRNs constitute an equivalence class under phenotype-preserving automorphisms, which can be found by graphical operations of (b) permuting loci, (c) permuting dummy proteins (circles), (d) exchanging edges pointing from two different stimuli (squares), and (e) exchanging self-loops at two different nodes.
In graph theory, symmetries in a network are formally described through the network's automorphisms. An automorphism of a graph is a way to shuffle the labels of its nodes such that the graph remains identical before and after shuffling. For instance, in Figure 5b, exchanging nodes 2 and 3 generates the same network and is thus an automorphism, whereas exchanging nodes 2 and 4 is not because there is an edge from 2 to 3 after shuffling. Formally, an automorphism of the genotype network G is a permutation σ of all plausible GRNs G = V(G) such that, for any g 1 , g 2 ∈ G, (σ(g 1 ), σ(g 2 )) ∈ E(G) if and only if we also have (g 1 , g 2 ) ∈ E(G). (A permutation of G is a mapping σ : G → G where no two GRNs are mapped to the same GRN, i.e., σ(g 1 ) = σ(g 2 ) if g 1 = g 2 for any g 1 , g 2 ∈ G.) Once two GRNs g and g are related through an automorphism σ of G, e.g., g = σ(g), they share the same mega-node properties that are fully determined by the connections in the genotype network (see Proposition A1).
Furthermore, automorphisms partition the GRNs by their roles in the genotype network through the mathematical concept of equivalence classes. For a high-level and general description, imagine a set of elements and a group of operations acting on them. Each operation turns one element into another, and these two elements are related by the operation, which describes the similarity between them. An equivalence class consists of elements that are mutually related by any operation, and the set of elements is said to be partitioned into equivalence classes under the action of the operations (see Figure 5a for an illustrative example). For automorphisms Σ(G ) of a graph G , the equivalence classes of nodes V(G ) under the action of Σ(G ) then gather nodes with a similar "structural position" in G (Figure 5c). However, to reveal GRNs with identical roles in a fitness landscape, these automorphisms also need to preserve the phenotype. Denote by Σ x (G) the set of such automorphisms of G, i.e., for any σ ∈ Σ x (G), and GRN g ∈ G, σ(g) and g have the same phenotype. The equivalence classes of mega-nodes V(G) under the action of phenotype-preserving automorphisms Σ x (G) then unite GRNs that (a) show similar mutational relationships with others and (b) lead to the same fitness due to their identical phenotype. We will mildly abuse the terminology to call them the equivalence classes of GRNs, which we denote by Θ, and each θ ∈ Θ is a set of GRNs related through Σ x (G). Crucially, since the mutational relationship and the resultant phenotype are the two components that characterize a GRN in the fitness landscape, GRNs in a θ ∈ Θ are deemed equivalent semantically, and they can be reduced to an arbitrary representative among them. Therefore, the equivalence classes of GRNs provide an efficient way to depict the underlying space of the fitness landscape.
However, what exactly composes the phenotype-preserving automorphisms Σ x (G) of the genotype network? From a sufficiency direction, we show that there exist a few graphical operations on the GRNs that produce phenotype-preserving automorphisms. These graphical operations involve permuting/shuffling different sorts of elements in a GRN: The identities of loci Γ, e.g., exchanging edge labels of loci A and B in Figure 4b; (ii) The identities of dummy proteins Ω , e.g., exchanging node labels of proteins 3 and 4 in Figure 4c.
Then, potentially rewiring a given edge (see details in Definitions A1 and A2): (iii) Change the source node of an edge from one stimulus to another stimulus and vice versa, e.g., in Figure 4d, moving an edge pointing from node 1 to node 3 to pointing from node 2. (Note that this operation is not necessarily equivalent to permuting the identities of stimuli since at most only the single focal edge will be affected.) (iv) Move a self-loop at one node to another node and vice versa, for example, re-allocating a self-loop at node 3 to node 4 in Figure 4e.
For the formal proofs, we point the reader to Theorem A1. Additionally, from a necessity direction, one can computationally obtain a partitionΘ of the GRNs G that is coarser than the equivalence classes Θ. (A partition P is coarser than another partition P if any group in P is included in some group in P.) Specifically, start with a partition Θ 0 where GRNs with the same resultant phenotype are grouped together. We create a sequence of partitions of G through the following iterative procedure: Given the partition Θ i , the next partitionΘ i+1 is obtained by further dividing groups intoΘ i (if needed) such that for each group θ ∈Θ i and θ ∈Θ i+1 , any two GRNs in θ have the same number of neighbors among θ. This iterative procedure is terminated when no further division is required, i.e.,Θ k+1 =Θ k for some integer k (see Figure 6a for an illustrative cartoon of the iterative procedure). We then haveΘ =Θ k to be our desired partition of GRNs. To see why the proposed iterative procedure generates a coarser partitionΘ than the equivalence classes Θ of GRNs, we stress that the equivalence classes under automorphisms always form an equitable partition. A partition P = {P i } m i=1 of nodes of a graph is equitable [81] if for every P i , P j ∈ P, any two nodes u, v in group P i have the same number of neighbors in P j (Figure 6b). Since GRNs in an equivalence class θ ∈ Θ must have the same amount of neighbors for each different phenotype, we inductively show that any two GRNs g 1 , g 2 ∈ θ are never separated during the iterative procedure that generateŝ Θ (see Theorem A2). Therefore, any equivalence class θ ∈ Θ must be included in a computationally acquired groupθ ∈Θ. Figure 7 demonstrates the coarser partitionΘ generated by the iterative procedure for an arbitrary toy example. The obtainedΘ contains 154 groups of GRNs, and the size of groups ranges from 2 to 96. We also count the number of different kinds of GRNs that can not be transformed through graphical operations (i) and (ii), and this number varies from 1 to 4 in our exampleΘ. Moreover, for every group inΘ, we observe that those different kinds of GRNs can be related by changing the stimulus that an edge is pointing from and re-allocating self-loops (e.g., see Figure 7b).Θ is thus not simply a coarser partition than the equivalence classes; according to (i)-(iv), we know that groups inΘ are exactly the equivalence classes Θ. This arguably general toy example implicates that there is no need for other graphical operations to determine the equivalence classes of GRNs.
As a result, we conjecture that all the phenotype-preserving automorphisms Σ x (G) of the genotype network can be generated by combining graphical operations (i) to (iv) on the GRNs. In other words, two GRNs g 1 and g 2 belong to the same equivalence class if and only if, after removing all the self-loops and merging stimuli Ω 0 into a single node, there exist permutations of loci Γ and dummy proteins Ω that jointly transform g 1 into g 2 . This condition reconciles with the concept of isomorphisms between graphs. Whereas an automorphism is a mapping of nodes such that a graph preserves itself, an isomorphism is a mapping of nodes that transform one graph into another. We will borrow the terminology and call the two permutations of Γ and Ω together a phenotype-preserving isomorphism from g 1 to g 2 .

Algorithmic Construction of the Mesoscopic Backbone of GRN Fitness Landscape
Next, we investigate algorithmic approaches to construct the mesoscopic backbone of a fitness landscape based on equivalence classes, where a representative GRN replaces all other GRNs in an equivalence class due to their identical role. In particular, the desired algorithm must (a) acquire the equivalence classes Θ from scratch and (b), for a representative GRN in any equivalence class, count the number of its mutational neighbors in other equivalence classes and also within the class it belongs to.
To avoid any confusion, we emphasize that, although drawing mutational connections between equivalence classes Θ can be achieved by grouping mega-nodes in the genotype network G, this naive exercise is unsuitable. First and foremost, grouping mega-nodes demands prior knowledge of the genotype network itself, but its construction is computationally heavy. Second, in contrast to coarse-graining nodes in a graph where the groups of nodes are pre-specified, listing all GRNs in an equivalence class requires examining pairs of GRNs and assuring a phenotype-preserving isomorphism between them after removing self-loops and merging stimuli. Determining the equivalence classes Θ from all the GRNs G = V(G) can thus be costly as well. These reasons again show the value of the equivalence classes Θ, which consolidate GRNs into their equivalent representatives.
Here, we present a bottom-up approach that enumerates each equivalence class of GRNs and simultaneously computes the number of mutational connections among them. To begin, recall from Section 2.2 that a mutation from a GRN g 1 ∈ G to another g 2 ∈ G corresponds to rewiring a single edge in g 1 , where g 1 may rewire a self-loop/non-self-loop edge to a self-loop/non-self-loop edge in g 2 . We observe that the number of non-self-loop edges in mutational neighbors g 1 and g 2 differ at most by one. We denote by Γ (g) the loci representing the non-self-loop edges in the GRN g, and |Γ (g)| the number of those non-self-loop edges. In other words, given equivalence classes θ, θ ∈ Θ and representative GRNs g ∈ θ and g ∈ θ , g has no mutational neighbors in θ if ||Γ (g)| − |Γ (g )|| > 1.
We can therefore build the mesoscopic backbone by incrementally examining each equivalence class with an increasing number of non-self-loop edges in the representative GRN. This strategy is envisioned in Figure 8, where the backbone can be viewed as "layers" of equivalence classes of GRNs. Let Θ k be the set of equivalence classes where for every θ ∈ Θ k , the representative GRN g ∈ θ has exactly k non-self-loop edges, |Γ (g)| = k. We start with layer Θ 0 , which consists of the only equivalence class with no non-self-loop edges. Then, with layers Θ 0 , Θ 1 , . . . , Θ k and all the mutational connections among them, we will find the equivalence classes in the next layer Θ k+1 and their mutational connections with layer Θ k and within themselves up until k = |Γ|, where all the edges are non-self-loops. To be more precise, we introduce the concept of M + neighborhood: For any GRN g ∈ G, denote by M + (g) the mutational neighbors of g that have one more non-self-loop edge than g. M + neighborhoods are sufficient to capture the relationship between two mutational neighbors g and g : • If g has one more non-self-loop edge than g, then g ∈ M + (g); • If g has one less non-self-loop edge than g, then we have g ∈ M + (g ); • If g has the same number of non-self-loop edges as g, and then they share a common mutational neighbor g , where the only different edge between g and g is rewired to a self-loop and thus g, g ∈ M + (g ). The mutational connections between equivalence classes can hence be uncovered by examining the M + neighborhood of the representative GRNs. Moreover, the M + neighborhood of representative GRNs in layer Θ k reveals the equivalence classes in layer Θ k+1 because any GRN must have a mutational neighbor with one less non-self-loop edge. All that remains is to join different M + neighbors into equivalence classes. In particular: (A) For an equivalence class θ ∈ Θ k and its representative GRN g ∈ θ, under what condition will g 1 , g 2 ∈ M + (g) belong to the same equivalence class in layer Θ k+1 ? (B) For two distinct equivalence classes θ 1 , θ 2 ∈ Θ k and their representative GRNs g 1 ∈ θ 1 and g 2 ∈ θ 2 , under what condition will g 1 ∈ M + (g 1 ) and g 2 ∈ M + (g 2 ) belong to the same equivalence class in layer Θ k+1 ?
For our ease of illustration, we hereafter choose the GRNs g, g 1 , g 2 , g 1 and g 2 such that only one stimulus node is incident to out-going edges.
To address (A), let g 1 , g 2 ∈ M + (g) belong to the same equivalence class, so there is a phenotype-preserving isomorphism π from g 1 to g 2 after self-loop removal. Recalling from Section 2.2, e g (γ) = (u, v) denotes that "the source-target pair of edge γ is (u, v) in GRN g." Furthermore, we write e g 1 (γ 1 ) = (u 1 , v 1 ) and e g 2 (γ 2 ) = (u 2 , v 2 ), where γ 1 and γ 2 are the non-self-loop edges "added" to g that forms g 1 and g 2 , respectively. A few observations follow:
We detail the reasoning behind these observations in Lemma A1-A3. Critically, our fifth observation implies that, after self-loop removal, the isomorphism π between g 1 and g 2 is in fact a phenotype-preserving automorphism of a subgraphḡ of the GRN g. In addition, observations 3. and 4. show that those edges in g-but not inḡ-are sequentially mapped from one to another via this automorphism π, i.e., Γ (g) \ Γ (ḡ) = π k (γ 1 ) q−1 k=1 , and they bridge the newly added edges γ 1 and γ 2 = π q (γ 1 ). We show that the converse is also true (see Theorem A3): After self-loop removal, if we find a phenotype-preserving automorphism π of a subgraphḡ of g where γ 1 is consecutively mapped to γ 2 through the edge differences Γ (g) \ Γ (ḡ), π is guaranteed a phenotype-preserving isomorphism from g 1 to g 2 .
The sufficient and necessary condition for two M + neighbors of g to be in the same equivalence class, intriguingly, lies in the phenotype-preserving automorphisms of subgraphs of the representative GRN g. Here, we demonstrate a few simple examples in Figure 9a. In the top row, an automorphism of g directly maps between the two additional edges (u 1 , v 1 ) = (3, 5) and (u 2 , v 2 ) = (1, 4). In the middle row, the two edges (u 1 , v 1 ) = (1, 2) are consecutively mapped to (u 2 , v 2 ) = (3, 4) through edge (2, 3), and (u 2 , v 2 ) is consecutively mapped back to (u 1 , v 1 ) through the non-edge (4, 1), so we have q = 2 and p = 4. As a mixture of both, in the bottom row, (u 1 , v 1 ) = (2, 5) is consecutively mapped to (u 2 , v 2 ) = (3, 6) through edge (1,4), and this isomorphism is exactly an automorphism of a subgraphḡ of g where edge (1, 4) is removed.
Switching gears to the remaining question (B), suppose that g 1 and g 2 are the representative GRN in two different equivalence classes where |Γ (g 1 )| = |Γ (g 2 )| and that g 1 ∈ M + (g 1 ) and g 2 ∈ M + (g 2 ) belong to the same equivalence class. Let γ 1 and γ 2 be the newly added edges to g 1 and g 2 that generate g 1 and g 2 , respectively, where e g 1 (γ 1 ) = (u 1 , v 1 ) and e g 2 (γ 2 ) = (u 2 , v 2 ), and let π be a phenotype-preserving isomorphism from g 1 to g 2 after self-loop removal. We observe that applying the permutation π on g 1 transforms it into another GRNg 1 in the same equivalence class. Since g 1 simply has one less edge γ 1 than g 1 , andg 1 and g 2 only differ by a missing edge π(γ 1 ). Namely, we have g 2 ∈ M + (g 1 ) with the additional edge e g 2 (π(γ 1 )) = (π(u 1 ), π(v 1 )). Moreover, since g 2 also belongs to the M + neighborhood of g 2 with the additional edge e g 2 (γ 2 ) = (u 2 , v 2 ), by removing both the extra edges from g 2 , we find a GRN g such thatg 1 , g 2 ∈ M + (g ).
Our observation resolves the necessary condition of (B): For the representative GRNs of two different equivalence classes g 1 ∈ θ 1 and g 2 ∈ θ 2 , if their M + neighbor g 1 ∈ M + (g 1 ) and g 2 ∈ M + (g 2 ) belong to the same equivalence class, then we can always find two GRNs g 1 and g such that (a)g 1 falls into the equivalence class of g 1 , and (b)g 1 and g 2 are M + neighbors of g . Moreover, the converse is true as well (Theorem A4). Therefore, whether the M + neighborhood of g 1 and g 2 reveal a common equivalence class depends on the existence of a GRN g that both the equivalence classes θ 1 and θ 2 are rooted from. Our strategy to build the mesoscopic backbone is now complete, and here, we detail our algorithm that incrementally generates the equivalence classes Θ of GRNs and establishes the mutational connections among them. Suppose that we have already built layers of equivalence classes Θ 0 , Θ 1 , . . . , Θ k and determined the mutational connections among them. For each representative GRN g in layer Θ k and every g ∈ M + (g), we will view g as the combination of g and an additional, non-self-loop edge e g (γ) = (u, v), for which we write g = g ⊕ (γ, u, v). All such combinations form a collection of M + neighbors of the representative GRNs in layer Θ k , for which we abuse the notation M + (Θ k ).
The joining operations Φ further provide useful information to count the number of mutation neighbors that a representative GRN g ∈ θ in layer Θ k has among any equivalence class θ , which we will denote by A g (θ ). Let us first consider θ ∈ Θ k+1 . For anỹ g ∈ M + Φ (θ ),g is a mutational neighbor of g if it can be viewed as a combination of g and an arbitrary extra non-self-loop edge, and hence Note that, in this case, A g (θ ) is easily acquired when building up the layer Θ k+1 through Φ. Second, for θ ∈ Θ k−1 , A g (θ ) can be computed given A g (θ), where g is the representative GRN of θ . Since the equivalence classes Θ generate an equitable partition of the genotype network G (see Section 3.2), we have A g (θ ) × |θ| = A g (θ) × |θ | equal to the total number of mutational connections between θ and θ . Moreover, the size of the equivalence class θ is (see Appendix D) where (a) we denote by Π the set of all permutations of dummy proteins Ω and denote by Σ (g) the set of automorphisms of the representative GRN g after self-loop removal that only permutes Ω ; (b) n l (|Γ| − k) is the number of ways to allocate |Γ| − k labeled self-loops among the proteins Ω; (c) m s (g) is the number of ways to re-distribute the edges pointing from stimuli Ω 0 in g; and (d) r(g) is the number of ways to divide loci Γ into self-loops, non-self-loop edges pointing from stimuli, and others. As a result, Third, we turn to the case where θ ∈ Θ k but θ = θ. Recall that, if anyg ∈ θ is a mutational neighbor of g, then there is a GRNg in layer Θ k−1 , where g,g ∈ M + (g ), and suchg is unique up to arbitrary self-loop re-allocation. Additionally, the extra edge in g andg must correspond to the same locus, so in which we use g to be the representative GRN of equivalence class θ . Lastly, if θ = θ, we also need to include the scenario that the mutational neighborg of g is generated by rewiring a self-loop to another self-loop. Therefore, In Algorithm 1, we summarize our proposed approach that constructs the mesoscopic backbone. It is apparent that the core of our algorithm is determining the joining operations Φ for a given layer Θ k . This task can be achieved by pre-computing the phenotypepreserving automorphisms of every representative GRN once it is chosen. In addition, since these joining operations reflect the mutational neighbors and the phenotype-preserving isomorphisms in previous layers, the type-(III) Φ for layer Θ k is generated as a composition of the already uncovered operations. Furthermore, the remaining Φ of type (II) consists of combinations of the uncovered joining operations and the newly computed automorphisms of representative GRNs in layer Θ k . As a result, the only prerequisite in our proposed algorithm is producing the phenotype-preserving automorphisms of a GRN.

Conclusions
In this work, we integrate mechanistic knowledge of how phenotypes are computed from genotypes via regulatory interactions into fitness landscape models. The resulting family of fitness landscape models features flexibility for tunable ruggedness and accessibility among phenotypes. Furthermore, we introduce the concept of equivalence classes of GRNs, where GRNs of the same phenotype and with similar structural positions in the genotype network are coarse-grained into a group. These equivalence classes of GRNs lead to a compact and informative description of the fundamental space of a fitness landscape. Using this coarse-graining, we develop a bottom-up, efficient algorithm for constructing the underlying space of a fitness landscape based on the equivalence classes. Critically, this algorithm does not require pre-computing the genotype network and therefore permits the exploration of substantially larger GRNs.
Naively, ruggedness and accessibility would seem to be contradictory characteristics of a fitness landscape. Indeed, reciprocal sign epistasis has been shown to yield a strong influence on a landscape's ruggedness and was regarded as an impediment to evolutionary accessibility when first introduced [2,32,55]. Nevertheless, recent studies suggest that fitness landscape models most closely aligned with empirical observations show that sign epistasis (and thus ruggedness) can co-exist with accessibility [63,82]. In addition to demonstrating that ruggedness and accessibility are not mutually exclusive, our model is compatible with three additional empirical observations. First, GRNs result in high dimensional genotype-phenotype maps [63]. Second, selection acts on the superposition of mutations and the background GRN rather than a few pairs of mutations [60]. Third, and perhaps most importantly, a GRN may experience a series of neutral mutations and then evolve into a nearby phenotype [3,8,83,84]. The accessibility induced in fitness landscapes of GRNs via neutral evolution agrees with the phenomenon of punctuated equilibrium/ epochal evolution [23,85,86].
Our derived equivalence classes for GRNs provide a novel, mesoscopic, and optimally descriptive skeleton of a fitness landscape. Neither the genotypic space nor the phenotypic space alone fully characterize a fitness landscape; however, models with even a relatively simple genotype-phenotype map are computationally intensive because they must retain all plausible genotypes [70,[73][74][75]. Intuitively, the complexity of a genotype-phenotype map can be reduced by combining similar phenotypes into high-level descriptors [87]. The equivalence classes of GRNs, on the other hand, serve as an intermediate level between the genotypic and phenotypic space, which provides an optimal coarse-graining that encodes all necessary information to predict the evolutionary trajectory on the fitness landscape.
We argue that our proposed algorithm for coarse-graining GRN fitness landscapes is more efficient than brute-force approaches. First, because we consolidate an equivalence class into a single representative GRN, our method is less costly in memory and requires fewer computations when finding mutational neighbors. Second, suppose all plausible GRNs were organized into layers by the number of non-self-loop edges (see Section 3.3), every layer would still super-exponentially contain many GRNs. Our algorithm instead finds the equivalence classes in each layer iteratively. To construct the (k + 1)-th layer, we only have to exhaust the representative GRNs in the k-th layer-along with any plausible additional non-self-loop edge(s)-this amount will be significantly fewer than the number of GRNs in the (k + 1)-th layer. Lastly, existing heuristics for graph automorphisms [88,89] can be used to produce the phenotype-preserving automorphisms of the representative GRNs, which is the only prerequisite when joining together different GRN-edge pairs. Because the set of automorphisms becomes more limited as the complexity of GRNs increases, we expect only a minor overhead in the joining procedure as compared to the exhaustive, brute-force approach.
Despite our model being constrained to the pathway framework of GRNs [77,78] and a few naive assumptions described in Section 2, we believe our methodology to be flexible and, in what follows, we outline some potential directions to extend the framework. First, when GRNs are modeled through more complex computation, e.g., with different logic gates connecting multiple expression activators/suppressors/products, those GRNs that only consist of naive interactions are never excluded. Thus, the current model represents a subset of the complete landscape built by more complex gene regulation. The derived connectivity and accessibility among the naive GRNs still hold, and we expect these topographical features to manifest for complex GRNs if mutations between the simple and complex expressions are permitted. Second, hypergraphs [90] could be used to describe the expression behavior of genes where multiple activators/products appear. Third, stable motif identification [91] and target control [92] for Boolean network models could be used to explore the phenotypes of mutational neighbors of a focal complex GRN. Lastly, our methodologies are likely applicable to other classes of genotype-phenotype maps [93,94]. In particular, once the mapping and the genotype network are determined, one can simply follow the proposed iterative procedure (Figure 6) to obtain a genotype partition coarser than the equivalence classes.
More broadly, this work showcases the potential of combining biological computation across different scales along the hierarchy of living systems. Computing biological functionality on the organism level with genotype-phenotype mapping provides a blueprint of the overall fitness landscape, where evolutionary processes occur/compute on the population level. Furthermore, several intriguing perspectives arise from the proposed mesoscopic backbone if we consider evolution to be a random walk on the fitness landscape. The process of evolution not only manifests genotypes with higher fitness values but also reveals genotypes whose mutational neighbors are more fit [19,23,78]; in other words, the prevalence of different genotypes would reflect the connection counts between equivalence classes of GRNs. In addition, these "connection counts" could become associated with an analogous theory of computation in evolution that addresses questions such as how likely a genotype in an equivalence class is to evolve into a specified phenotype, as well as how likely it is to "reset" to another genotype in the same equivalence class and recover its position in the fitness landscape.

Data Availability Statement:
The authors affirm that all data necessary for confirming the conclusions of the article are present within the article, gures, and tables.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Central and Peripheral GRNs Where No Regulation Presents
Here, we show the connectivity of "central" and "peripheral" GRNs mentioned in Section 3.1 under the extreme scenario where no gene regulation appears, in particular, when |Ω + y | = |Γ|. We observe that, since |Ω + y | = |Γ|, each node ω ∈ Ω + y is incident to one and only one incoming edge, and thus, each central GRNg ∈G y corresponds to a bijective mapping from Ω + y to edge labels Γ. First, we show that if two central GRNsg 1 ,g 2 ∈G y correspond to different mappings between Ω + y and Γ, there is no mutational trajectory connecting them among G y . Start with assuming that such a mutational trajectory does exist. Due to the different associated mappings ofg 1 andg 2 , there is an ω ∈Ω with distinct labels of the incident edge γ 1 and γ 2 ing 1 andg 2 , respectively. Because mutatingg 1 intog 2 requires rewiring both γ 1 and γ 2 , there must exist a GRN g where either none or both the edges γ 1 and γ 2 point to ω. Nevertheless, g contradicts our observation following the constraint |Ω + y | = |Γ| so g / ∈ G y . As a result, under this extreme scenario, G y fragment into multiple connected components when only mutations among themselves are considered.
Next, for any phenotype y for which |Ω + y ∪ Ω + y | ≤ |Γ| + 1, we show that there is a mutational trajectory among G y connecting an arbitrary central GRNg ∈G y and a peripheral GRNĝ ∈ G y at the boundary of G y and G y . Specifically, take ω ∈ Ω + y \ Ω + y and its incident edge pointing from ω 0 ∈ Ω 0 . For each ω ∈ Ω + y \ Ω + y where ω = ω, one can sequentially rewire the incident edge of ω to form a chain of Ω + y \ Ω + y initiated by ω 0 , which leads to a resultant GRNĝ. Moreover, since |Ω + y | = |Γ|, we have |Ω + y \ Ω + y | = 1. Recall from Section 3.1, thisĝ is indeed a peripheral GRN between G y and G y .

Appendix B. Phenotype-Preserving Automorphisms of the Genotype Network of GRNs
In this section, we demonstrate (a) why GRNs mapped by automorphisms of the genotype network G are equivalent, (b) four graphical operations that generate phenotypepreserving automorphisms of G, and (c) the correctness of our iterative procedure to obtain a coarser partition than the equivalence classes of GRNs.
Proposition A1. Given an automorphism σ of the genotype network G and a mega-node function f G that depends on the adjacency matrix A of G, for any GRNs g 1 , g 2 ∈ G where g 2 = σ(g 1 ), f G (g 2 ) = f G (g 1 ).
Let π and π be a permutation of the loci Γ and the dummy proteins Ω , respectively. It is not hard to see that π and π also generate a permutation of the GRNs G. For g ∈ G, we abuse the notation π(g) to be the GRN mapped through the locus permutation π, where an edge with e g (γ) = (u, v) is transformed into e π(g) (π(γ)) = (u, v). Similarly, in the GRN π (g) mapped through the dummy protein permutation π , an edge with e g (γ) = (u, v) is transformed into e g (γ) = (π(u), π(v)).
Furthermore, we have two more types of graphical operations on GRNs: Definition A1. For a locus γ and two stimuli ω, ω ∈ Ω 0 , ρ γ,ω,ω : G → G transforms a GRN g into g such that edge γ becomes Definition A2. For a locus γ and two nodes ω and ω , γ,ω,ω : G → G transforms a GRN g into g such that the edge γ becomes Note both ρ γ,ω,ω and γ,ω,ω are permutations of G as well, where pairs of GRNs are mutually mapped from one to the other.
These four graphical operations introduced above are, more importantly, automorphisms of the genotype network G that also preserve the phenotype of GRNs: Theorem A1. The transformations π, π , ρ γ,ω,ω , and γ,ω,ω are phenotype-preserving automorphisms of G.
Furthermore, since π, π and γ,ω,ω simply change the labels of edges, labels of the intermediate nodes, and the location of a self-loop, they maintain any path from a stimulus ω 0 ∈ Ω 0 to a fitness-relevant proteinω ∈Ω. ρ γ,ω,ω may alter the path betweenω 0 and ω, but the reachability from Ω 0 toω remains. Therefore, the four transformations also preserve the phenotype of GRNs.
Finally, we turn to the computationally acquired partition that can be shown to be coarser than the equivalence classes of GRNs. Recall from Section 3.2 that our iterative procedure starts from a partitionΘ 0 where GRNs with the same phenotype are grouped together. Given the partitionΘ i , the next partitionΘ i+1 is obtained by further dividing groups intoΘ i (if needed) such that for each group θ ∈Θ i and θ ∈Θ i+1 , any two GRNs in θ have the same number of neighbors among θ in the genotype network G. The procedure terminates when a stationary partitionΘ is reached. We then have: Theorem A2. Every equivalence class θ ∈ Θ is included in a groupθ ∈Θ, θ ⊂θ.
Proof. Recall that GRNs in an equivalence class have the same phenotype, so for each θ ∈ Θ, there is some θ 0 ∈ Θ 0 where θ ⊂ θ 0 . Suppose that θ ⊂ θ i for each θ ∈ Θ and some θ i ∈ Θ i . Since Θ forms an equitable partition, every g ∈ θ has the same number of neighbors in each θ ∈ Θ and thus also in each θ i ∈ Θ i . Consequently, no two GRNs in θ will be separated into two different groups in Θ i+1 , and the theorem is proved by induction.

Appendix C. Combining Mutational Neighbors into Equivalence Classes
In this section, we tackle the two questions raised in Section 3.3: (A) For an equivalence class θ ∈ Θ k and its representative GRN g ∈ θ, under what condition will g 1 , g 2 ∈ M + (g) belong to the same equivalence class in layer Θ k+1 ? (B) For two distinct equivalence classes θ 1 , θ 2 ∈ Θ k and their representative GRNs g 1 ∈ θ 1 and g 2 ∈ θ 2 , under what condition will g 1 ∈ M + (g 1 ) and g 2 ∈ M + (g 2 ) belong to the same equivalence class in layer Θ k+1 ?
Furthermore, recall that for ease of demonstration, we constrain the GRNs g, g 1 , g 2 , g 1 and g 2 where only one stimulus node is incident to out-going edges.
The following theorem resolves the necessary and sufficient condition in question (A): Theorem A3. Let g 1 , g 2 ∈ M + (g) with the additional non-self-loop edge e g 1 (γ 1 ) = (u 1 , v 1 ) and e g 2 (γ 2 ) = (u 2 , v 2 ), respectively. g 1 and g 2 belong to the same equivalence class if and only if there exist two integers q < p and a phenotype-preserving automorphism σ of a subgraphḡ of g such that i.
Next, the necessary and sufficient condition of question (B) is described in the theorem below: Theorem A4. Let g 1 and g 2 be of two different equivalence classes, and let g 1 ∈ M + (g 1 ) and g 2 ∈ M + (g 2 ) with the additional non-self-loop edge e g 1 (γ 1 ) = (u 1 , v 1 ) and e g 2 (γ 2 ) = (u 2 , v 2 ), respectively. g 1 and g 2 belong to the same equivalence class if and only if there exist GRNsg 1 and g such that i.g 1 and g 1 belong to the same equivalence class; ii.g 1 , g 2 ∈ M + (g ); iii. g 2 ∈ M + (g 1 ).

Appendix D. Size of an Equivalence Class of GRNs
Here, given the representative GRN g of an equivalence class θ, we calculate the number of GRNs in θ in Equation (2). Observe that |θ| is proportional to (a) the number of GRNs to which there is a phenotype-preserving isomorphism from g, (b) the number of ways to arbitrarily allocate self-loops on g, and (c) the number of ways to arbitrarily rewire source nodes of edges among the stimuli Ω 0 .
First, for part (a), if we temporally ignore labels on the edges, the set of permutations over dummy proteins Ω is partitioned into groups of isomorphisms from g to different GRNs. In addition, the size of each group is exactly the number of automorphisms of g since the composition of an automorphism of g and an isomorphism from g to g also generates an isomorphism from g to g . Thus, there are in total |Π | |Σ (g)| different GRNs isomorphic to g, where Π is the set of permutations over Ω , and Σ (g) is the set of automorphisms of g that only permutes Ω .
Third, we write k s (g) as the number of incident edges to stimuli Ω 0 in g. Any possibility in part (c) chooses a source node among Ω 0 for each of the k s (g) incident edges. Providing that the labels of incident edges of the stimuli are already known, we have m s (g) = |Ω 0 | k s (g) .