Monotonicity is a key feature of genotype-phenotype maps

It was recently shown that monotone gene action, i.e., order-preservation between allele content and corresponding genotypic values in the mapping from genotypes to phenotypes, is a prerequisite for achieving a predictable parent-offspring relationship across the whole allele frequency spectrum. Here we test the consequential prediction that the design principles underlying gene regulatory networks are likely to generate highly monotone genotype-phenotype maps. To this end we present two measures of the monotonicity of a genotype-phenotype map, one based on allele substitution effects, and the other based on isotonic regression. We apply these measures to genotype-phenotype maps emerging from simulations of 1881 different 3-gene regulatory networks. We confirm that in general, genotype-phenotype maps are indeed highly monotonic across network types. However, regulatory motifs involving incoherent feedforward or positive feedback, as well as pleiotropy in the mapping between genotypes and gene regulatory parameters, are clearly predisposed for generating non-monotonicity. We present analytical results confirming these deep connections between molecular regulatory architecture and monotonicity properties of the genotype-phenotype map. These connections seem to be beyond reach by the classical distinction between additive and non-additive gene action.


INTRODUCTION
Quantitative genetics is the major theoretical foundation for genetic studies in production biology, evolutionary biology, and biomedicine. A core concept in quantitative genetics is the genotypic value, the mean observed phenotype for a given genotype. It constitutes the basis for the genotype-to-phenotype (GP) map concept. The shape of a given GP map is typically described by the classical gene action terms: additivity, dominance, and epistasis. Together with genotype frequencies in a given population, the GP map is the basis for decomposing observed phenotypic variance into environmental variance and genetic variance components including additive variance, dominance variance and epistatic variance. This provides the basis for a very successful theory when it comes to predicting selection response and breeding values (Falconer and Mackay, 1996;Lynch and Walsh, 1998) and more recent statistical genetics methods for mapping Quantitative Trait Loci (QTL) (Neale et al., 2008). Quantitative genetics thus provides a mature machinery for predicting the population level consequences of a given GP map, but in order to understand several generic genetic phenomena there is a stated need for new tools for disclosing how the shape of the GP map is determined by underlying biology (Jaeger et al., 2012;Moore, 2012;Gjuvsland et al., 2013).
One such phenomenon is the resemblance between parents and offspring. An explanation in quantitative genetic terms is that the additive variance (V A ) makes up a substantial part of the phenotypic (V P ) and genetic variance (V G ). Hill et al. (2008) showed that in populations with extreme allele frequencies, high V A /V G ratios will arise regardless of the shape of the GP map. However, for populations with intermediate allele frequencies a much wider range of V A /V G ratios is observed (Wang et al., 2013). In such populations, high V A /V G ratios cannot be fully accounted for without considering properties of the GP map. Gjuvsland et al. (2011) showed that a key feature of GP maps that give high ratios of additive to genotypic variance (V A /V G ), is a monotone (or order-preserving) relation between gene content (the number of alleles of a given type) and phenotype. This led to the hypothesis that the regulatory circuitry of sexually reproducing organisms predominantly predisposes for highly monotone genotype-phenotype maps.
Here we address the above hypothesis by a two-step approach. First we provide methods and software tools for measuring monotonicity of generic GP maps (i.e., sets of genotypic values). Then we use these tools on the data generated by an extensive simulation study of a broad collection of gene regulatory network models. In these network models the steady state expression levels serve as phenotypes and genetic variation is introduced through parameters describing maximal production rates and the shape of the gene regulation function. Such causally cohesive genotype-phenotype (cGP) models [see Gjuvsland et al. (2013) and references therein] allow us to identify relationships between regulatory network architecture and properties of the resulting GP maps.
Our results confirm the prediction that the GP maps arising from a wide range of gene regulatory network motifs are in general highly monotone. In addition we show through numerical as well as mathematical analysis that regulatory motifs involving incoherent feed-forward or positive feedback stand out in their capacity to generate non-monotonicity. These relationships between molecular regulatory architecture and properties of the genotype-phenotype map-of substantial relevance to functional genomics in general-are beyond reach by the standard distinction between additive and non-additive gene action.
Our approach can be applied to cGP models of a wide range of biological systems at any level of model complexity. It opens for a systematic study of the monotonicity properties of molecular regulatory structures underlying the whole spectrum of physiological regulation. This suggests that the concept of monotonicity of GP maps can be used to build theory about heredity phrased in terms of molecular mechanism, something which standard genetic concepts and approaches appear to be incapable of.

BACKGROUND ON MONOTONICITY OF GP MAPS
To ease understanding we provide a brief recapitulation of the concept of monotonicity (or order-preservation) in GP maps introduced in (Gjuvsland et al., 2011). We consider a diploid genetic model with N biallelic loci (alleles indexed 1 and 2) underlying a quantitative phenotype. A genotype at a single locus k is denoted by g k ∈ {11, 12, 22}. In the case of two loci k and l there are 9 possible genotypes g kl = g k g l ∈ {1111, 1112, 1122, 1211, . . . , 2212, 2222}. The general N loci genotype space contains 3 N genotypes g 1 g 2 · · · g N (in condensed notation g 1:N ) constructed by concatenating single locus genotypes, = {g 1 g 2 · · · g N | g k ∈ {11, 12, 22}, k = 1, 2, . . . , N}.
We use the 2-allele content (i.e., the number of 2-alleles) of genotypes to define a partial order on the genotype space (see Figure 1, left panel for an illustration). For a particular locus k we order the three genotypes sharing the same background genotype g 1: k − 1 g k + 1: N as follows, We call this the partial genotype order relative to locus k, and it defines a strict partial order on .
A genotype-phenotype map is a mapping G that assigns to each genotype g ∈ a real-valued genotypic value G(g) (the mean trait value for a given genotype). We define monotonicity of G in terms of how it transforms the partial genotype order to the algebraic order of the genotypic values G(g). Without loss of generality we assume that the allele indexes at each locus have been chosen such that G(1111 · · · 11) is the smallest of all homozygote genotypic values. We call a genotype-phenotype map G monotone or order-preserving with respect to locus k if it preserves the partial genotype order relative to locus k, i.e., if, (2) FIGURE 1 | Examples of partial genotype order and genotype-phenotype maps. Left panel: The allele content defines a partial order on genotype space. A two-locus example is shown. The plot at the top displays the genotype at locus 1 (x-axis) and locus 2 (color) vs. the total number of 2-alleles (y-axis) in the two-locus genotype. The resulting partial ordering of genotypes is shown below. Right panel: Each lineplot shows the 9 genotypic values (y-axis) for a single GP map, coding of genotype are the same as in the left panel. GP maps that preserve the partial order of genotypes are called monotone. Examples shown are an intra-and interlocus additive map (A), a map showing partial dominance at both loci (PD), and duplicate dominant (DD) epistasis (see Table 1 in Phillips, 1998 for all genetic backgrounds of locus k. By allowing non-strict inequalities we include GP maps showing complete dominance and complete magnitude epistasis (Weinreich et al., 2005) in the class of order-preserving GP maps. Conversely we call a GP map non-monotone or order-breaking with respect to locus k if it does not preserve the partial genotype order relative to locus k for all backgrounds. Figure 1 (right panel) shows classical dominance and epistasis patterns, categorized into monotone and non-monotone GP maps.

STATISTICAL DECOMPOSITION OF GENOTYPE-PHENOTYPE MAPS
Given a genotype-phenotype map G as described above and a corresponding vector of genotype frequencies f in a population, quantitative genetic provides methods for orthogonal decomposition of genotypic values and resulting genetic variance in the population into additive and non-additive (dominance and epistasis) components (Lynch and Walsh, 1998). We performed such statistical decomposition with the function linearGPmapanalysis in the R package noia (http:// cran.r-project.org/package=noia; Le Rouzic and Alvarez-Castro, 2008) version 0.94.1. We assumed an idealized population where all genotype frequencies are equal (1/3 N ). In such a hypothetical population the NOIA (Alvarez-Castro and Carlborg, 2007) statistical and functional formulations and the unweighted regression model proposed by Cheverud and Routman (1995) are equivalent. Furthermore, the decomposition of genotypic values is equivalent to decomposing G into a sum of additive and nonadditive GP maps, and the genetic variance in this case is simply the variance of the 3 N genotypic values in G. We used the NOIA statistical formulation to decompose a GP map G into its additive and non-additive components, and computed the ratio of additive to total genetic variance V A /V G as a measure of how well the additive component describes the original GP map. In case of the illustrative GP maps depicted in Figure 1, this gives V A /V G = 1 for the fully additive GP map A, and V A /V G = 0 for the pure overdominance (OD) and the pure epistasis (Cheverud and Routman, 1996) maps A × A and D × D.

GENE REGULATORY NETWORK MODELS
Gene expression in eukaryotes is controlled through gene regulatory networks involving numerous regulatory mechanisms [see e.g., Latchman (2005), for details]. Modeling of such gene regulatory networks is well-established, and available modeling frameworks range from coarse-grained descriptions of the topology of genome-wide networks to very detailed mechanistic models describing the dynamics of small networks (De Jong, 2002;Schlitt and Brazma, 2007;Karlebach and Shamir, 2008). In line with a large number of authors we used ordinary differential equations (ODEs) to study a family of generic gene regulatory network models containing three diploid genes X 1 , X 2 , and X 3 , organized as a regulatory system where the rate of expression of each gene can be regulated by the expression level of one or both of the other genes. The wiring of the system is described by a 3 × 3 connectivity matrix A with elements A kl ∈ {−1, 0, 1}. The signs of the elements of A describe the mode of regulation, A kl = 0 indicates that X l is not a regulator of X k , if A kl = 1 then X l is an activator of X k , and if A kl = −1 then X l is a repressor of X k . Gene regulatory systems are often laid out visually as signed directed graphs. There is a one-to-one correspondence between a connectivity matrix and a signed directed graph, two examples are illustrated in Figure 4. We used the sigmoid formalism Plahte et al., 1998) in the diploid form (Omholt et al., 2000) where the expression the two alleles of gene k is described by the following ODEs, x k2 = α k2 R k2 (y 1 , y 2 , y 3 ) − γ k2 x k2 , Here α ki is the maximal production rate for allele i of gene X k , γ ki is the decay rate, while R ki is the gene regulation function (dose-response function). If X k has no regulators, we assume production is always switched on i.e., R ki = 1. If X k has a single regulator X l , the gene regulation function is given as R ki (y l ) = S(y l , θ lki , p lki ), where S(y, θ, p) = y p /(y p + θ p ) if X l is an activator and S(y, θ, p) = 1 − y p /(y p + θ p ) if it is a repressor. In both cases the parameter θ lki gives the amount of regulator needed to get 50% of maximal production rate, and p lki determines the steepness of the response. In the case of two regulators X l and X j we set R ki (y l , y j ) = S(y l , θ lki , p lki )S(y j , θ jki , p jki ), corresponding to the Boolean AND function. Modeling transcription regulation by means of Hill functions and Boolean composition has a long tradition in modeling of gene regulation and is widely used.
With three genes and up to two regulators per gene the number of possible connectivity matrices is 6859. We further required that the system is connected, and that X 3 is downstream to both X 1 and X 2 so either X 1 and X 2 both regulate X 3 directly (A 31 A 32 = 0), or one of them regulates X 3 directly and the other one indirectly (A 31 A 12 = 0 or A 32 A 21 = 0). This reduces the number of distinct connectivity matrices to 3724. Finally, we identified pairs of matrices that are symmetric with respect to interchanging X 1 and X 2 and picked just one matrix from each pair. The resulting 1881 connectivity matrices were used for our gene regulatory simulations.

IDENTIFYING FEEDBACK LOOPS AND FEEDFORWARD MOTIFS
Feedback and feedforward motifs appear recurrently as regulatory building blocks in transcription networks across all living organisms. These network motifs have several characteristic features (Alon, 2007), negative feedback can for example accommodate fast transcriptional responses and homeostasis, while positive feedbacks are utilized as biological switches. We went through all 1881 gene regulatory models and extracted information about their feedback and feedforward loop characteristics from their connectivity matrices. For each system we computed three autoregulatory feedback loop products FL 1 = A 11 , FL 2 = A 22 , FL 3 = A 33 , three two-gene feedback loop products: FL 12 = A 21 A 12 , FL 13 = A 31 A 13 , FL 23 = A 23 A 32 and two three-gene feedback loop products: FL 123 = A 32 A 21 A 13 , FL 213 = A 31 A 12 A 23 . Non-zero loop products indicate that the system contains the corresponding feedback loop, and the sign of the loop product gives the sign of the feedback loop. We also computed the products for two feedforward motifs: FFL 32 = A 32 (A 31 A 12 ), . Again non-zero products indicate that the system contains the corresponding feedforward motif, a positive value corresponds to a coherent feedforward while a negative value indicates incoherent feedforward. Figure 4 depicts the connectivity matrix and the signed digraphs of a system with a positive feedback loop as well as a system with incoherent feedforward. Spreadsheet S1 contains adjacency matrices and loop products for all 1881 motifs.

GENE REGULATORY NETWORK SIMULATIONS
The simulation were performed with the Python package cgptoolbox (http://github.com/jonovik/cgptoolbox), using the sigmoidmodel submodule, which contains an implementation of the gene regulatory network model (Equation 3) and the connectivity matrix A. A similar simulation setup is found in Gjuvsland et al. (2011) together with a discussion of gene regulation functions and the genotype-parameter map in molecular terms. We compared two different types of genotype-toparameter maps: • Genotype to parameter map without pleiotropy: biallelic genotypic variation for all three loci was introduced through the maximal production rates α ki . For each Monte Carlo simulation the allelic parameter values were sampled from U(100, 200). • Genotype to parameter map with pleiotropy: allelic parameter values were sampled for maximal production rates α ki (sampled from U(100, 200)), regulation thresholds θ lki (sampled from U(20, 40)), and regulation steepnesses p lki (sampled from U(1, 10)).
All decay rates γ ki were set equal to 10. We assembled parameter sets for all 27 diploid genotypes, and for each genotypic parameter set the system of Equation 3 was integrated numerically until convergence to a stable state. The equilibrium value of y 3 was recorded as phenotype. Datasets where the system failed to converge for one or more genotypes were discarded. For each of the 1881 motifs we performed 1000 Monte Carlo simulations. Some Monte Carlo simulations lead to very little phenotypic variation, in the sense that the span between the largest and smallest of the 27 genotypic values was small. In order to avoid artifacts arising from the numeric ODE solver tolerance, these essentially flat GP maps were discarded. Further analysis of monotonicity and variance components were only performed on GP maps where the absolute range (maximum genotypic value -minimum genotypic value) and relative range (absolute range/mean genotypic value) were both > 0.01.

MEASURING MONOTONICITY OF GP MAPS
In the following we present two numerical measures for quantifying monotonicity in a GP map G with N biallelic loci. The first quantifies the monotonicity for individual loci by comparing negative and positive allele substitution effects before weighting the individual loci into an overall measure. The second utilizes isotonic regression to quantify the distance between G and the closest fully monotone GP map.

Measure 1: quantifying non-monotonicity by substitution effects
We first develop a measure of monotonicity based on the effects of substituting a single allele at locus k, while keeping the background genotype g (k) = g 1: k + 1 g k + 1: N fixed. Monotonicity as defined by Equation 2 is equivalent to s i (g (k) ) ≥ 0 for i = 1, 2 across all genetic backgrounds of locus k. By taking into account also the magnitude of the substitution effects we can quantify the deviation from strict monotonicity. We start with the set S k = {s i (g (k) )} of single allele substitution effects for locus k for i = 1, 2 and across all genotypic backgrounds g (k) . The total number of elements in S k thus becomes 2 · 3 N−1 , and we split the set into two disjoint subsets reflecting their sign; S k We compute the sum of positive substitution effects and the sum of absolute values of negative substitution effects, and let T k = P k + N k denote the overall sum of absolute substitution effects. We then define the degree to which the GP map G is monotone with respect to locus k by, |s 1 (g (k) )| + |s 2 (g (k) )| .
The absolute value in the numerator ensures that the measure m k is invariant with respect to the choice of indexes for the two alleles of locus k. Interchanging the numbering of the alleles leads to the mappings s 1 (g (k) ) → −s 2 (g (k) ), s 2 (g (k) ) → −s 1 (g (k) ), which leaves the value of m k unchanged. By the triangle inequality m k ≤ 1. If m k = 1, then G is monotonic with respect to locus k, whereas m k < 1 implies that G is order-breaking w.r.t. locus k. If m k = 0, then the positive substitution effects equal the negative substitution effects in magnitude and we say that G is completely order-breaking w.r.t. locus k. This measure distinguishes well between the monotone and non-monotone maps in Figure 1.
Clearly m 1 = m 2 = 1 for the additive map (A) and GP maps showing partial dominance and duplicate dominance epistasis. In contrast, m 1 = m 2 = 0 for the maps showing pure OD and pure epistasis (A × A and D × D). In order to quantify the overall monotonicity of the GP map G we introduce the degree of monotonicity (m) which is a weighted mean of all m k , where the weights reflect the relative effect size of the loci in terms of T k , As shown in Figure 3A, the degree of monotonicity is accordingly 1 for the monotone maps in Figure 1 while it is 0 for the pure OD and pure epistasis maps. This definition of degree of monotonicity allows us to establish a vocabulary that is analogous to the classification of single locus dominance; i.e., a GP map is called For example, the degree of monotonicity of the GP map published by Cheverud and Routman (1995), with two loci underlying 10-week body-weight (in grams) at 10 weeks in a mouse F 2 cross, may be computed as follows. After renaming the two loci (B →1, A →2) and indexing alleles to conform to our notation, the nine genotypic values ( Table 1 in (Cheverud and Routman, 1995)) are G (1111)  For random GP maps (randomly sampled genotypic values as in (Gjuvsland et al., 2011)) there is a strong positive correlation between the degree of monotonicity and the size of the additive component (V A /V G ) ( Figure 3A). A similar relationship was observed for three-locus random GP maps ( Figure A1A). All GP maps in Figure 3A with m < 0.1 have V A /V G < 0.1. At the other end of the spectrum there is much more variation, for instance the most extreme completely monotone map (the duplicate dominant factors DD) has V A /V G as low as 0.375.

Measure 2: quantifying monotonicity by isotonic regression
This measure quantifies the monotonicity of a particular GP map G in terms of the least-squares distance to the closest monotone map. We build on the mathematical notation introduced in section "Background on monotonicity of GP maps" where is the genotype space for N biallelic loci and a GP map is a function that assigns a real-valued genotypic value G(g) to each genotype  Table 1 in Phillips (1998); duplicate recessive genes (DR) and recessive epistasis (RE). Red dots show 1000 random two-locus GP maps, while blue dots show the same 1000 GP maps after rearranging genotypic values to introduce order-preservation for 1 locus [see Model and Methods in Gjuvsland et al. (2011)].

FIGURE 2 | Decomposition of genotype-phenotye map into monotone and non-monotone components. Left panel:
Genotype-phenotype map G for two loci underlying 10-week body-weight at 10 weeks in a mouse F 2 cross. The GP map shown here is equivalent to the one in the original publication [see Figure   3A in Cheverud and Routman (1995)], but we have changed indexing of loci and alleles for consistency with the notation used here. The GP map G is decomposed with isotonic regression into a (middle panel) monotone component G M and a (right panel) non-monotone component G N .

www.frontiersin.org
November 2013 | Volume 4 | Article 216 | 5 g in . For any particular GP map G, we identify the monotone component of G as the map G M which minimizes the residual variance var(G − G M ), i.e., G M is the monotone GP map which is closest to G in the least-squares sense. For a given G the monotone component G M is unique (Barlow and Brunk, 1972) and can be computed numerically by isotonic regression (Leeuw et al., 2009) of G subject to the partial ordering of genotypes defined in Equation 1. Furthermore, the residual G N is orthogonal to G M in the sense that g ∈ G M (g)G N (g) = 0. This allows the orthogonal decomposition, of a genotype-phenotype map into a monotone component G M and a non-monotone component G N such that var(G) = var(G M ) + var(G N ). The orthogonality property allows us to measure monotonicity of G in terms of the coefficient of determination R 2 mono of the isotonic regression given by the ratio R 2 mono = var(G M )/var(G). In the case that G itself is monotone for all loci we have R 2 mono = 1, while order-breaking for one or more loci will result in R 2 mono < 1. The isotonic regression approach can be illustrated in a straightforward way on the two-locus GP map provided by Cheverud and Routman (1995) (see text above and left panel of Figure 2). The partial ordering of genotypes defined by Equation 1 is illustrated in Figure 1 (left panel). By isotone regression (Leeuw et al., 2009) on this partial genotype ordering, the original GP map is decomposed into a monotone and a non-monotone component (Figure 2, middle and right panels), and the coefficient of determination (R 2 mono ) is 0.97. Our simulation results for random GP maps show that R 2 mono is positively correlated to the size of the additive component ( Figure 3B for two-locus GPs maps and Figure A1B for threelocus GP maps) and that for a given V A /V G the lower bound for R 2 mono is close to a straight line from (0, 0.2) to (1, 1). However, due to the search for the closest monotone GP map, R 2 mono will not become zero even for purely overdominant or purely epistatic maps. As shown in Figure A2, the two monotonicity measures are highly correlated.

An R package for studying monotonicity in GP maps
We developed an R package gpmap for studying functional properties of GP maps. The package takes GP maps in the form of vectors of genotypic values as input, and provides functions for (i) determining whether the map is order-breaking or orderpreserving w.r.t. any given locus, (ii) the degree of monotonicity m, (iii) R 2 mono using isotonic regression from the isotone package (Leeuw et al., 2009), and (iv) plots of the original and decomposed GP maps. Code example 1 (Box 1) below illustrates the usage and functionality of the gpmap package. The package is available from CRAN http://cran.r-project.org/package=gpmap under GPLv3.

MONOTONICITY IN GP MAPS ARISING FROM GENE REGULATORY NETWORKS
To search for generic relationships between monotonicity and regulatory network structure, we used the above measures of monotonicity to characterize GP maps emerging from the gene regulatory network models (see Models and Methods). Based on earlier results (Gjuvsland et al., 2007(Gjuvsland et al., , 2011Wang et al., 2013) we hypothesized that incoherent feed forward (Figure 4, right panel) or positive feedback (Figure 4, left panel) would be necessary in order to obtain highly order-breaking GP maps, and we characterized all 1881 networks in terms of these two properties. Table 1 shows the number of motifs falling into the resulting four categories. We summarized the number of Monte Carlo simulations where all genotypic parameter sets gave convergence to a stable steady state, and where the resulting GP maps were not essentially flat (see Models and Methods for details). Motifs with less than 100 usable GP maps were discarded from further analysis. For the genotype-to-parameter maps without pleiotropy (in the sense Box 1 | Code example 1.

FIGURE 4 | Connectivity matrices and signed directed graphs.
Connectivity matrix A and the corresponding signed directed graph for two of the 1881 systems in the simulation study. The left panel depicts the connectivity matrix and the signed digraph of a system with a positive feedback loop between X 1 and X 2 while the right panel shows a system with incoherent feedforward from X 1 to X 3 .

Frontiers in Genetics | Genetic Architecture
November 2013 | Volume 4 | Article 216 | 6 that genetic variation at one locus influences only a single parameter, see Model and Methods) 868 motifs were discarded, while for the genotype-to-parameter map with pleiotropy (genetic variation at one locus influences three parameters) 791 motifs were discarded. All (but one) discarded motifs contained at least one positive feedback loop ( Table 1). A plausible explanation for this is that many motifs with positive feedback loops have a stable steady state at, or very close to 0 for one or more state variables regardless of parameter values, and this leads to essentially flat GP maps.
The introduction of pleiotropy in the genotype to parameter map has a marked effect on the monotonicity characteristics of the associated GP map (Figure 5). When genetic variation at a locus X i affects only its maximal production rate the GP maps come out as highly monotone (Figure 5A), with a large majority being fully monotone or order-breaking for just a single locus. When genetic variation at locus X i affects the threshold and steepness of the dose-response curve in addition to the maximal production rate (pleiotropy in the genotype-to-parameter map), the majority of GP maps still show order-breaking either for no loci or just one locus (Figure 5B). But a considerable number of GP maps are in this case order-breaking for two or three loci. Furthermore, by dividing the motifs into the four groups given in Table 1 it is evident that the regulatory anatomy of a network determines its predisposition for non-monotonicity in its associated GP map. Presence of incoherent feedforward or positive feedback loops appears to be prerequisites for the majority of the observed non-monotonic GP maps.
The class of motifs lacking both incoherent feedforward and positive feedback contains very few order-breaking GP maps, and with no pleiotropy in the genotype-to-parameter map we observe only fully order-preserving GP maps for this class (cyan in Figure 5A). In the Appendix we generalize this to an arbitrary number of nodes and formally prove that without pleiotropy in the genotype-to-parameter map, the presence of incoherent feedforward or positive feedback is indeed a necessary condition for non-monotone GP maps to arise from networks with monotone gene regulation functions.
The introduction of pleiotropy in the genotype-to-parameter map increases the frequency of order-breaking GP maps substantially ( Figure 5B). Motifs lacking both incoherent feedforward  Table 1 for the number of motifs in each class. A single boxplot summarizes, for all motifs in the given class, the proportion of the GP maps (y-axis) that are order-breaking with respect to a given number of loci (x-axis). For example, consider the red box at x = 0 in panel (A). This boxplot contains results for motifs with both incoherent feedforward and positive feedback and from Table 1 we find that the red boxplot summarizes results for 135 motifs. From the y-axis we find that at least half (box median at y = 1) of these 135 motifs result in only monotone GP maps, while for the most extreme (end of whisker) of the 135 motifs only 25% of the GP maps are monotone. Similarly, the cyan box is compressed into a line at x = 0, y = 1 indicating that all 251 motifs that lack both incoherent feedforward and positive feedback result in only monotone GP maps. and positive feedback may in this case lead to GP maps that are order-breaking for one or two loci, but never for all three loci. Using isotonic regression to quantify the overall monotonicity of the GP maps reinforces the finding that incoherent feedforward and positive feedback predispose for non-monotonicity (Figure 6). Figure 6 also shows that for all classes of motifs the majority of GP maps are fully monotone, while the most non-monotone GP maps (lowest R 2 mono values) are observed for motifs with positive feedback. The differences between classes www.frontiersin.org November 2013 | Volume 4 | Article 216 | 7

A B
FIGURE 6 | Empirical distribution functions for R 2 mono . Summary of R 2 mono values from isotone regression for all motifs for which at least 100 (out of 1000) Monte Carlo simulations lead to GP maps with non-negligible phenotypic variation (see Models and Methods section "Gene regulatory network simulations," for detailed criteria). Results are shown for 1013 motifs with a genotype-to-parameter map without pleiotropy (A) and 1090 motifs with a genotype-to-parameter map with pleiotropy (B). Each panel is divided into 4 subplots containing classes of motifs based on the presence/absence of incoherent feedforward and positive feedback loops, see Table 1 for the number of motifs in each class. Each curve shows, for a single motif, the empirical distribution function value (y-axis) of R 2 mono for all GP maps (x-axis).
of motifs are also evident when inspecting the additivity of GP maps (Figure A3), but since monotone GP maps can still be non-additive, the patterns are much more blurred than for monotonicity. Fisher's (1918) regression on gene content and the concepts derived from this, such as additive effects and dominance deviation, provide the theoretical basis for most of quantitative genetics (Falconer and Mackay, 1996;Lynch and Walsh, 1998). By regressing on gene content, including the extensions by Cockerham (1954), the genotype-phenotype map is decomposed into additive, dominant, and epistatic components. The use of gene content or the number (0, 1, or 2) of alleles with a particular index in a genotype implies the same partial ordering of genotype space as defined in Equation 1. Thus, our proposed definition of monotonicity of GP maps, and in particular the use of isotonic regression to quantify monotonicity, may be viewed as a relaxation of the linearity assumption underlying current quantitative genetics theory. In this perspective the positive correlation between monotonicity and additivity (Figure 3) is expected. We have addressed GP maps with 2 and 3 loci as we considered an in-depth study of the properties of GP maps with higher number of loci to be outside the scope of this study. Some general observations can be made, though. Since m is a weighted average, the m k of major loci (i.e., for which T k is large relative to T k ) will tend to dominate. For instance, in a case with a single major locus showing monotone gene action and several minor loci showing order-breaking, the GP map will overall be close to monotone (m close to 1). Conversely, order-preservation in a number of minor loci would have little influence on m if major loci have strongly non-monotone gene action. Isotonic regression gives an overall measure of monotonicity of a GP map, but provides no locus-specific measures corresponding to m k . Similar to the case for m, the gene action of major loci will have high influence on the value of R 2 mono . The observation that monotonicity is an important property of GP maps is in principle not new. For a single locus, non-monotone gene action appears in the form of over-or under-dominance, while complete and partial dominance as well as additivity exemplify monotone gene action. Weinreich et al. (2005) distinguished between sign epistasis and magnitude epistasis and showed that sign epistasis limits the number of mutational trajectories to higher fitness. As sign epistasis reflects a nonmonotone GP relationship and magnitude epistasis reflects a monotone one, this insight concords with our results. A similar distinction has been proposed (Wang et al., 2010) for statistical interactions where removable interactions are those that can be removed by a monotone transformation of the phenotype scale, while non-monotonicity in the GP map leads to essential interactions. Wu et al. (2009) developed a method to screen for and test the significance of essential interaction in genome-wide association studies. Isotonic regression has also recently been applied to link genotype and phenotype data (Beerenwinkel et al., 2011;Luss et al., 2012). Our treatment of monotonicity is more general than these earlier works in three major ways. First, we deal with monotonicity of the GP map as a whole rather than either intra-locus (dominance vs. overdominance) or inter-locus (magnitude vs. sign epistasis and removable vs. essential interactions). Second, where the earlier treatments have focused on classifying the type of gene action, we make use of quantitative measures of monotonicity. Third, our approach combining the concept of monotonicity with cGP models opens a direct link between genetics and the theory of dynamical systems in the wide sense.

DISCUSSION
Monotonicity is a property of the GP map separate from the allele frequencies, making it a physiological (Cheverud and Routman, 1995) or functional (Hansen and Wagner, 2001) descriptor rather than a statistical one. The distinction between physiological and statistical epistasis has lead to much debate (Phillips, 2008). Zeng et al. (2005) argued the distinction was unnecessary and potentially misleading. Although their arguments around orthogonality and variance components are valid, our results demonstrate very clearly that describing the properties of the GP map without reference to any particular study population is essential if we want to connect quantitative genetics with regulatory biology.
It is clear from our results that positive feedback and incoherent feedforward promote non-monotonicity. The clear-cut differences in monotonicity between different classes of regulatory networks, combined with the strong correlation between monotonicity and additivity of GP maps, appear therefore to explain the findings that regulatory systems with positive feedback give considerably more statistical epistasis than those without (Gjuvsland et al., 2007;Wang et al., 2013). Even though both incoherent feedforward and positive feedback predispose for non-monotone GP maps, the underlying mechanisms are different for the two regulatory motifs. In the case of incoherent feedforward the sum of direct and indirect effects may result in a non-monotone doseresponse relationship (Kaplan et al., 2008). That positive feedback loops can give non-monotonicity is intuitively less clear, but in the Appendix we show both results analytically. Positive feedback predisposes for multiple steady states, and order-breaking might also emerge from different genotypes corresponding to different states. It should be noted, however that positive feedback is only a necessary condition for multistationarity , and a positive loop in the connectivity matrix A of a system is not necessarily active at any point during the time course of the system.
Without any restrictions on the connectivity of a threegene system there are 3 9 = 19, 683 possible distinct networks. The main restriction we imposed (see Models and Methods for details) was a maximum of two regulators per gene, which allowed us to use Boolean gene regulation functions already established in the sigmoid formalism (Plahte et al., 1998). Other model formalisms allowing an arbitrary number of regulators are also available (Wagner, 1994(Wagner, , 1996Siegal and Bergman, 2002) and could be extended to diploid forms and used in later studies.
Although this study has focused on gene regulatory networks, the concept of monotone gene action applies to the propagation of genetic variation across the whole physiological hierarchy. One may therefore systematically use the concepts and methods presented here to study the orderpreserving and order-breaking properties of genotype-phenotype mappings that are associated with any regulatory structure amenable for mathematical modeling. Through this it will be possible to make a wide-ranging survey of which regulatory anatomies promote monotonicity and which promote nonmonotonicity. We foresee that this classification may become instrumental for predicting how phenotypic effects of genetic variation propagate across generations in sexually reproducing populations.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://www.frontiersin.org/journal/10.3389/fgene.2013.00216/ abstract Spreadsheet S1 | Excel spreadsheet with connectivity matrices and loop products for all 1881 gene regulatory networks.

APPENDIX
In this appendix we complement the simulation studies in the main text with some analytic results for GP maps emerging from ODE models of gene regulatory networks. We study a generalization of the gene network model in Equation (3) with an arbitrary number of loci and monotone gene regulation functions, but restrict the analysis to genotype-parameter maps without pleiotropy. In particular, we show that (i) if there are no positive feedback loops and no incoherent feedforward loops in the network, the resulting GP maps are always monotone, (ii) a positive feedback loop or an incoherent feedforward loop may lead to non-monotone GP maps. The results hold for phenotypes given as the stable concentration of the product of one of the genes, and under certain restrictions also for phenotypes given as a function of one or several stable gene product concentrations that is monotonic with respect to each of its arguments.

GENE NETWORK MODEL
We consider a dynamic system consisting of n mutually interacting diploid loci X j , j ∈ N = {1, . . . , n}, regulating each other's expression. The time dependent output of X j is denoted z j , and we define z = [z 1 , z 2 , . . . , z n ]. It goes without saying that z j in general depends on the genotypes of all the genes even though we will not always state this explicitly.
For a given genotype g = g j g (j) = a j b j g (j) , where g j ∈ {11, 12, 22} denotes the genotype and a j , b j ∈ 1, 2 denote the indexes of the two alleles of locus X j , the equations of motion for X j areż where z 1 j and z 2 j are the time-dependent outputs of the two homologous copies of X j . The two allele rate functions r 1 j (z) and r 2 j (z) have range [0, 1] so that α 1 j and α 2 j represent the maximum production rates of the two alleles. We assume that all dose-response functions in Equation (A1) are differentiable and monotonic with respect to each of its arguments, and that for each j, k, the signs of ∂r 1 j /∂x k and ∂r 2 j /∂x k in the stable point x are equal. This model generalizes Eq. (3) to an arbitrary number of loci and a broader class of gene regulation functions.
In the following we are only concerned with the steady states of Equation (A1), and assume for simplicity that they have just a single stable equilibrium point. Solving the equilibrium conditions of Equation (A1) with respect to z 1 j and z 2 j and adding gives j . Since our definition of monotonicity of GP maps does not depend on the numbering of alleles, we will without loss of generality assume μ 1 j ≤ μ 2 j for all j.
The network architecture can be read out from the structure of the system's Jacobian matrix in the stable state x. We define the elements of the Jacobian J for the set of functions f j defined in Equation (A2) by To the Jacobian J it is customary to assign a signed directed graph G in which each locus X k is represented by a node X k , and in which there is an arc from X j to X k if and only if J kj = 0, its sign given by the sign of J kj . A chain from X j to X k is a set of arcs in G leading from X j to X k in which all intermediate nodes are visited only once. The sign of a chain is equal to the product of the signs of the J ij corresponding to the arcs in the chain. If there is a chain from X i to X j and also a chain from X j to X i through a disjoint set of nodes, the two chains constitute a proper feedback loop (FBL).
To each FBL is associated a loop product L which is the product of the Jacobian elements corresponding to all the arcs in the loop. The sign of the loop is given by the sign of L. Two chains from X j to X i , i = j, with only the endpoint nodes in common, constitute a feedforward loop (FFL). If the two chains have opposite signs, the FFL is incoherent (IFFL), otherwise it is coherent (CFFL).
The system's phenotype could be any scalar quantity defined by its equilibrium value x. In the following we assume the genotype-phenotype map G(g) = x q (g), q ∈ N, for a given and fixed q, and investigate the monotonicity properties of G(g k g (k) ) with respect to genetic variation in any locus X k for different backgrounds g (k) . In the following sections we analyse the causes of order-breaking in G in the restricted case in which there is only genetic variation in μ 1 k and μ 2 k , not in the shape of the doseresponse functions r 1 k and r 2 k , implying r 1 . This is what we mean by a genotype-to-parameter map without pleoitropy.
In the next sections we prove the following result: At the end of this note we show that under some reasonable conditions this result is also valid for more general phenotypes depending on more than one x q .

NETWORKS WITHOUT LOOPS
We first consider networks containing no feedforward loop and no feedback loop. In these networks there is at most one chain from one node to another, and of course no autoregulatory loops. If there is a chain from X j to X k , there is no chain from X k to X j . Any node is either unregulated (constitutively expressed) or regulated by one or several other nodes.
of J composed of the rows and columns V. Because there is no feedback loop among the nodes represented in D (11) and D VV , only the diagonal degradation terms contribute to these two determinants. Hence D (11) = (−1) n − 1 i = 1 γ i . Similarly, D VV = (−1) n − m i∈V γ i , giving q ml = γ 1 C U / U , where U = i∈U γ i . Finally, we note that P = (∂R 1 /∂x m )C U is the loop product of the loop.
Solving Equation (A9) with respect to dx 1 /dx k and using all these expressions lead to The sign of ∂x 1 /∂x l is independent of the genotype of X k and the sign of q 1k is fixed. Genotypic variation in X k may change the magnitude of P, but its sign is fixed because all Jacobi elements have fixed sign independent of the system parameters. Thus, genotypic variation in X k does not alter the sign of dx 1 /dx k if the loop is negative (P < 0), while for a positive loop the sign of U − P may switch. In the latter case, an increase in x k due to genetic variation in X k may increase x 1 in some cases and decrease it in others, leading to order breaking. As there is only a single chain from X 1 to X q , no order breaking in X 1 implies no order breaking in X q , while order breaking in X 1 may propagate to X q . The same result follows if X q is downstream a node in the loop because order breaking in this node may propagate to X q .

FEEDFORWARD LOOPS (FFLS)
A feedforward loop (FFL) is a motif in the network in which there are two different chains C 1 and C 2 from one particular node to another particular node. To each chain C i is associated a chain product P i defined as the product of the Jacobian elements corresponding to the arcs in C i . If P 1 and P 2 have equal signs, the FFL is coherent, otherwise it is incoherent. In a network with a single feedforward loop and no feedback loops we now investigate the effect on G(g) = x q (x k (g)) of genetic variation in X k for varying background g (k) . Our starting point is again Equation (A7). We first let X k and X q be the initial and terminal nodes in the FFL. The two chains C 1 and C 2 leading from X k to X q comprise ρ 1 and ρ 2 nodes including X k and X q , respectively. Let the set of nodes in C 1 and C 2 be X U 1 and X U 2 , respectively, where U 1 and U 2 are the corresponding subsets of N, and let V 1 and V 2 be their complements.
Roughly speaking, the derivative of the propagation function p qk (x k ) can be expressed as a sum of terms, each term corresponding to one of the chains leading from X k to X q (Plahte et al., 2013).
To the chain C i is assigned the chain weight w i given by where D V i V i is the Jacobian subdeterminant for the nodes not included in C i , and D (kk) is the Jacobian subdeterminant for all nodes except X k . Because there are two chains from X k to X q , the derivative of p qk is a sum of two terms: dp qk dx k = w 1 P 1 + w 2 P 2 , where P 1 and P 2 are the two chain products, and w 1 and w 2 their weights (Plahte et al., 2013). When there is no feedback loop in the system, only the diagonal elements in J stemming from the term −γ i x i in Equation (A7) contribute to the determinants D V i V i and D (kk) : Altogether this gives where 1 and 2 are the products of the γ j in the two chains, respectively. The chain products P 1 and P 2 depend on the genotype g k of X k as well as on the genotypic background g (k) , but their signs S 1 and S 2 are invariant under genotypic variation. It is easy to see that a negative autoregulatory loop, which is a common feature in gene regulatory networks, would not invalidate the conclusion, but a positive autoregulatory loop might.
If the FFL is incoherent, P 1 and P 2 have opposite signs, implying that the sign of dx q /dx k may vary. If the FFL is coherent, however, no order-breaking can occur.
If X k is upstream relative to the initial node X init of the FFL, it follows from the above section on networks without loops that there will be no order-breaking in X init , and the above argument is still valid.

MORE GENERAL PHENOTYPES
In real life, relevant phenotypes are not direct gene products, but rather functions of the concentrations of one or several gene products. Let the phenotype G(g) be a function of x U (g), G = h(x U (g)), where U is a subset of N, and assume that for any u ∈ U, ∂h/∂x u has fixed sign for all genotypes. To analyse this case we extend the original system Equation (A2) to and apply our above results to this system, in which G(g) = x n + 1 (g), i.e. q = n + 1. If there are two nodes among X U which have a common predecessor X k , then there will exist two chains from X k to X n + 1 . These two chains constitute a feedforward loop with X n + 1 as final node. If this FFL is incoherent, order breaking due to genetic variation in X k may occur even if there is no order breaking in the original system comprising the nodes X 1 , . . . , X n . If the FFL is coherent, order breaking only occurs if it occurs in the original system. Red dots show 1000 random three-locus GP maps, blue dots show the same 1000 GP maps after sorting to introduce order-preservation for 1 locus while green dots show the same 1000 GP maps after sorting to introduce order-preservation for 2 loci [see Model and Methods in Gjuvsland et al. (2011)].

FIGURE A2 | Comparing measures of monotonicity GP maps.
Scatterplots showing degree of monotonicity (m) vs. R 2 mono . Black dots correspond to the maps shown in Figure 1. Red dots show 1000 random two-locus GP maps, while blue dots show the same 1000 GP maps after sorting to introduce order-preservation for 1 locus [see Model and Methods in Gjuvsland et al. (2011)]. and1090 motifs with a genotype-to-parameter map with pleiotropy (B). Each panel is divided into 4 subplots containing classes of motifs based on the presence/absence of incoherent feedforward and positive feedback loops, see Table 1 for the number of motifs in each class. Each curve shows, for a single motif, the empirical distribution function value (y-axis) of V A /V G from unweighted regression for all GP maps (x-axis).