Evolvability of feed-forward loop architecture biases its abundance in transcription networks

Background Transcription networks define the core of the regulatory machinery of cellular life and are largely responsible for information processing and decision making. At the small scale, interaction motifs have been characterized based on their abundance and some seemingly general patterns have been described. In particular, the abundance of different feed-forward loop motifs in gene regulatory networks displays systematic biases towards some particular topologies, which are much more common than others. The causative process of this pattern is still matter of debate. Results We analyzed the entire motif-function landscape of the feed-forward loop using the formalism developed in a previous work. We evaluated the probabilities to implement possible functions for each motif and found that the kurtosis of these distributions correlate well with the natural abundance pattern. Kurtosis is a standard measure for the peakedness of probability distributions. Furthermore, we examined the functional robustness of the motifs facing mutational pressure in silico and observed that the abundance pattern is biased by the degree of their evolvability. Conclusions The natural abundance pattern of the feed-forward loop can be reconstructed concerning its intrinsic plasticity. Intrinsic plasticity is associated to each motif in terms of its capacity of implementing a repertoire of possible functions and it is directly linked to the motif's evolvability. Since evolvability is defined as the potential phenotypic variation of the motif upon mutation, the link plausibly explains the abundance pattern.


Background
Evolutionary adaptability in biological systems is often the result of trade-offs between flexibility and specialization [1]. In this context, buffering mutations and noise seem an important requirement for stability. This can be achieved by a robust response to parameter changes and correlates with the degree of specialization of the given structure. A given network insensitive to mutations will always perform the same function. On the other hand, adaptation and evolvability requires flexible structures that can be re-used to perform different (potential) functions and thus provide plasticity [2,3]. The problem here is often understanding why some particular structures are so common and what their (if any) functional meaning is. This is closely tied to the map-ping f between structure S and function F, namely the relationship: which is usually dubbed as the genotype-phenotype mapping problem [4]. Understanding the nature and origins of this mapping is at the core of many key questions concerning the evolution of complexity in nature.
Within the context of gene transcription networks, it has been suggested that the previous problem can be dissected by analyzing the frequency of some overabundant sub-networks of three or four elements, so called network motifs [5][6][7]. These sub-graphs only capture the topological pattern of connections and a dynamical description of their potential function requires a set of differential equations [8,9]. One particularly important example is provided by feed-forward loop (FFL) motifs [5]. Many genetic and biochemical systems, such as the Lac and Che systems in E. coli (responsible for lactose utilization and chemotaxis, respectively) involve FFL motifs [10][11][12]. Mounting evidence indicates that they have key roles in cell function [10] and morphogenesis [13,14].
However, the origin of a preferential bias towards given topologies remains under discussion.
The relative frequency of FFLs displays a well-defined pattern (figure 1c) dominated by two sub-graphs (C1 and I1). The uneven abundance of these graphs could be a fingerprint of their functional relevance [8,[15][16][17]. Such importance would be the blueprint of an evolutionary advantage, but it is not clear whether such functional connection really exists [15,[18][19][20][21] or if it resembles instead a byproduct of non-adaptive processes [22][23][24][25]. As shown below, motif structure does not directly relate to its frequency, but its plasticity in implementing different functions does.

Probability distribution of implementing different functions
Consider the FFL graphs Γ i from the set S = {C 1 , C 2 , C 3 , C 4 , I 1 , I 2 , I 3 , I 4 } shown in figure 1a. In previous studies [9,26] it has been shown that, the topology of a given FFL does not univocally define its function but it captures the probability distribution of implementing different functions. Our first goal is to identify an appropriate mapping f between FFL topology and each potential response j j (t), i.e. f ij : {Γ i , x(0), μ} j j (t) where we indicate as {Γ i , x(0),μ} the FFL graph together with its initial condition x(0) and the set of parameters used μ [26]. The six different responses (figure 1b) are triggered by an external input. These are either fast response (pulser) or delayed response (grader) considering the target concentration of the output. Here f ij indicates the likelihood that j j (t) is implemented by How likely is a motif to become part of a complex cellular network? Two extreme strategies can be envisioned. In the first, specific motifs play specific roles in a robust way and they are common because they are insensitive to mutational noise. In the second, the larger the variety of implementable functions, the more flexible the better. Such a scenario is feasible under the premise of ever-changing environments and comes with the cost of reduced robustness. In order to measure the plasticity of decision-making between these two strategies, let us first determine the (conditional) probability f ij = P(j j |Γ i ). These probabilities are normalized, i. e.
{φ j } P(φ j | i ) = 1 and can be systematically computed [26]. This set actually defines our structure-function map, namely and can be displayed (figure 2a) as a weighted, motiffunction bipartite graph (see Methods). The graph reveals that most motifs implement all functions, but the likelihood of each pair is case-dependent. Some motifs seem clearly more specialized (such as C4) whereas others are rather generalists (see for example I4). What influences the choice of a given topology over others? Since most of the functions can be implemented, it is not clear that a one-to-one, function-based argument will work. But we can go a step beyond and look at the structure of the probability distribution {P(j j |Γ i )} of each topology. This can be done by measuring the degree of homogeneity displayed, for each Γ i , by plotting P(j j |Γ i ). If the motif is highly specialized, some dominant peak(s) would be observed, whereas if it is very flexible no prominent peaks will appear. A simple way Figure 1 Structure and frequency of FFL motifs. In (a) we show the schematic representation of the FFL's genetic regulatory interactions ('+' represents activatory regulation and '-' represents inhibitory regulation). The external input I activates the signal protein X. Active X modulates expression of gene G Z directly and indirectly via regulation of Y expression, which in turn also regulates G Z . The dynamics of these regulatory interactions is described by a set of equations dy/dt = F(y, z), dz/dt = G(y, z) incorporating the nonlinearities associated to gene-gene interactions. In (b) we plot the general topology of FFL motifs and the six different functions j (t) represented by qualitative time-courses [Z(t)]. 'G' indicates grader dynamics, 'P' pulser dynamics. We specifically take into account the initial slope of the time-course ('+' or '-') and the concentration of the final target Z with respect to the non-induced protein concentration ('T+' and 'T-'). In (c) we display the relative abundance P obs (Γ i ) of these motifs in the transcription networks of yeast and E. coli (data from [27]). of measuring the homogeneity of the distribution is given by its kurtosis, defined as the fourth standardized moment by Kurtosis is the measure of the "peakedness" of a distribution. It quantifies the concentration of frequencies around the mean of the distribution. Higher kurtosis means that the variance is the result of infrequent, extreme deviations from the mean as opposed to frequent, modestly sized deviations resulting in low kurtosis. In order to define a measure characterizing the degree of plasticity of a given motif in terms of its specialization or its flexibility, we can consider two extreme cases, namely the most flexible graph Γ f equally likely to implement any function j j , and the most specialized graph Γ s implementing only one function j 1 . In the first case we would have P(j j |Γ f ) = 1/6 and the kurtosis associated is K(Γ f ) = -3.33, whereas in the second case we would have P(Γ s ) = 1 and 0 otherwise, with kurtosis K(Γ s ) = 6. Details on the calculation can be found in Methods. Any other FFL graph Γ i from the set S = {C 1 , C 2 , C 3 , C 4 , I 1 , I 2 , I 3 , I 4 } has kurtosis values locating within this interval, K(Γ i ) (K(Γ f ), K(Γ s )).
In order to measure the degree of plasticity in the decision-making between these extreme cases we introduce ψ(Γ i ). ψ(Γ i ) is the distance between the absolute value of the kurtosis K(Γ i ) and the origin K 0 [26] or in other words represents the intermediate level between specialization and flexibility. This transformation opens the way for a more intuitive biological interpretation: The values for ψ(Γ i ) range between high plasticity (low ψ(Γ i )) and high commitment towards one of the extreme strategies, i.e. maximal specialization or maximal flexibility (high ψ(Γ i )). The optimal solution here is likely to be strongly impacted by the predictability of the environment. As a first approximation we therefor place K 0 at the midpoint of the interval (K f ,K s ), i.e.
Finally we define the likelihood r(Γ i ) of a given motif Γ i to appear within a network as a function of its ψ(Γ i ). Assuming that a high degree of flexibility or specialization should be related with lower likelihood of appearance for a given motif Γ i , we write where a is a normalization coefficient defined as α = 8 j=1 ψ( j ) −1 . This function actually defines the expected probability of finding a given sub-graph and is thus mapping between the distributions associated to each motif and the expected abundance of motifs within networks. Figures 3a and 3b show the correlation between relative abundances of FFL motifs in E. coli and S. cerevisiae with respect to their expected probabilities r(Γ i ). The matching is striking (data concerning abundances obtained from [27]). The two most abundant graphs (C1 and I1) are consistent with our results and the actual distribution matches well the observed pattern.
Interestingly, the expected probabilities indicate a positive bias toward systems which show high plasticity as presented in figure 4. Intermediate values for kurtosis (figure 4a) or low ψ(Γ i ) (figure 4c) correlate with an increase in the likelihood of appearance. The values for kurtosis, ψ(Γ i ) and r(Γ i ) are collected in table 1. As an alternative measure for the homogeneity of a probability distribution the Shannon entropy is discussed in the Methods section, where the motif's entropy is used to characterize the degree of flexibility or specialization of of FFL motifs is displayed as a bipartite graph linking patterns (upper row) and processes (lower). The weight of the links indicates the relative probability P(j j | Γ i ) that a given motif Γ i implements a given function j j . In (b) the matrix of motif-function probabilities is displayed using a color scale. The plots highlight that some motifs look more specialized, whereas others display rather evenly distributed functional responses. a b Figure 3 Predicted probability and FFL abundance. In (a) we compare the natural abundance and its predicted counterpart r(Γ i ). S. cerevisiae (black box) is compared to E. coli (white box) and the predicted probabilities (black triangle). In (b) we present the correlation between r(Γ i ) and the natural abundances. The Pearson coefficient for the linear fit is r = 0.91 and r = 0.94 for E. coli and S. cerevisiae, respectively. a given motif Γ i [26]. Both, entropy (figure 4b) and kurtosis (figure 4a) yield similar results, ranking the most abundant motifs with intermediate values which translates to high functional plasticity.
For the less abundant motifs we see a more disordered trend in the two measures, as is the case for C3, C4, I3 (both measures) or C2 (entropy). The interpretation here is not straightforward. It is feasible that the disordered trends can be consequence of non-adaptive processes. An alternative hypothesis is related to the shape of the real distributions for the implementation of any function. We assume that for more and less frequent motifs the analytically deduced probability distributions does not fit equally well the real counterpart. The more abundant the motif, the better the underlying probability distribution is mirrored in its abundance, because the sampling space is covered more readily.
Our analysis of FFLs dynamics was performed considering single, isolated motifs. However, in real systems motifs are embedded in large networks allowing for the combination of motifs. The combination of more abundant motifs, such as C1 and I1, can cover the whole set of possible dynamics by that affecting the abundance of the rest of the motifs.

Evolvability
In order to have a relevant role in evolvability, the degree of plasticity of FFLs should correlate with the motif's capacity of generating phenotypic variation by exploring different functions under mutation. Two key aspects are of importance here, namely i) the reduction of mutational lethality and ii) the up-speeding of adaptational processes (reduction of the number of mutations needed to generate new phenotypes [1]). The evolvability of the circuits Γ k can be studied by calculating the transition probabilities ω m k (φ j |φ i ) of shifting from function j i to j j under m mutations. The matrix m (k) = (ω m k ) defines a flow graph (figure 5b) which allows us quantifying its evolvability ε( i ). We compared the robustness against single mutations versus sequential accumulation of multiple mutations. Mutations m are defined as single parameter changes. For better understanding of the procedure, we want to stress the conceptual difference of continuous real mutations and the here applied parameter changes m which result in a discretized observation pattern. In the presented   Sketch of the procedure. In (a) the update rule is shown. For a given string of conditional relations (BR) the associated dynamical pattern is calculated at time-step t -1. Next the entries of BR are mutated at time-step t and the (new) dynamical pattern is evaluated. Then the transition of pattern t-1 to pattern t is binned. This protocol is executed until no more changes in the bins occur as shown in plot (c). In (b) the graph Ω m (Γ i ) associated to the transitions between the possible types of dynamics is represented for C1. The thickness of the arrows correspond to transition probabilities obtained from procedure (a). From these graphs (see Methods) we can calculate the motif evolvability ε, which is found to positively correlate with r (Γ i ).
framework a mutation m can be both, numerically small or large without impact as long as it does not drive the system into another functional regime. Robustness is defined as the sum of the diagonal ele- give the probabilities of performing the same function after m mutations. A detailed description of the procedure is given in Methods. We have found that FFL motifs are very robust but some of them exhibit a high phenotypic variation under repeated mutations. The most abundant motifs, C1 and I1, show the highest phenotypic variation. In other words, C1 and I1 can widely change their function with greater ease than the rest of the circuits facilitated by their low ψ(Γ i ) (figure 6b). A network displaying little phenotypic diversity would give small values of ε whereas sub-graphs with high transition rates among states will have a high ε. As presented in figure 6a, ε correlates positively with the abundance of motifs, with C1 and I1 displaying the largest values. These results suggest that a proper degree of plasticity, in terms of a balance between flexibility and specialization, is the optimal strategy to increase evolvability providing the playground for adaptive responses without increasing mutational lethality.
Assuming that motif plasticity is a relevant trait, our analysis supports the idea that the observed FFL abundance pattern actually correlates with motif evolvability.
Our analysis suggests that neither a direct interpretation of motifs as functional modules [1,2,4] nor a purely non-adaptive view of their abundance [22][23][24] account for the uneven presence in transcription networks. Consistently with previous works [28,29] duplication-rewiring dynamics alone cannot explain the evolution of FFLs. The potential for evolvability associated to their topological structure might well be the missing ingredient connecting both views.

Conclusions
In this article we have interpreted a simplified, qualitative model of the FFL motif. The thorough analysis within the model framework allows to reconstruct its natural abundance pattern and provides insight in what might have shaped it. The argument leads to the very core of the genotype-phenotype mapping problem, since, due to its simplicity, a perfect mapping between the topology and all possible functions it can implement can be constituted. We claim, however, general applicability. FFL abundances are correlated with their plasticity and evolvability. Evolvability has been defined as a compromise between robustness against single mutations and the capability to modify the functional response upon increasing mutational pressure. The results indicate that a proper portion of intrinsic functional plasticity, which can be understood as a strategic trade-off between specialization and flexibility, is necessary to be abundant. Because only then one is suited to be readily evolvable in changing environments.
Future work should be devoted to analyzing how the coexistence of different motifs embedded in a large network affects their dynamics and abundance compared to the single motif analysis performed in this work.

Dynamical response of FFL motifs
Network motifs are recurrent interaction patterns, which are significantly more often encountered in biological interaction graphs than expected from random nets. It has been shown that feed forward loops (FFL) are capable of processing external signals by responding in a very specific, robust manner, either accelerating or delaying responses. They are composed of three genes. Firstly, gene G X that expresses the protein X. This protein X regulates the expression of the other two genes G Y and G Z encoding the proteins Y and Z, respectively. Additionally, protein Y regulates the expression of Z (see figure 1). Here, we assume that expression of X is unregulated and the protein is expressed in its inactive form, i.e. X does not regulate the expressions of Y and Z straightaway. Only upon the presence of an external signal (the input) X becomes active and regulation of Y and Z takes place, where Z resembles the output of the motif. The dynamics describing how the concentration of Z changes during time from the initial state (without input) to the final state (with input) are calculated from a set of two differential equations [26]: Abundance Evolvability a Psi Evolvability b Figure 6 FFL evolvability. In plot (a) we show the correlation between the motif's evolvability and its abundance. Evolvability ε, which is found to positively correlate with r(Γ i ), is highest for the most abundant motifs C1, I1. In black we show data point of S. cerevisiae, in white E. coli. In (b) we show the correlation between the FFL's evolvability and its ψ(Γ i ). We calculate a Pearson coefficient of r = -0.92 for the linear fit. The lower ψ(Γ i ), the higher the motif's plasticity and the higher its evolvability. The data points are developed from the motifs topology and thus are species independent (blue).
Here γ i describes the basal production of protein i, with i = {Y, Z}, subsuming the concentration of all biochemical elements which remain constant in time. The binding equilibrium of the regulators j with the gene G i are denoted by ω j i , with j = {X, Y, Z}. Parameters α x and b j define the type of regulatory interactions, i.e. activation or inhibition, for gene G Y and G Z , respectively, providing the regulatory rates with respect to the basal transcription. Values < 1 correspond to inhibitory regulation, whereas > 1 accounts for activation (denoted by '-' and '+' in figure 1a respectively). The parameter b xy accounts for the simultaneous regulation of G Z . The degradation rate of protein i is denoted as δ i . Finally, n and m are the degree of multimerization of the regulators.
If we consider the system in phase space, we find that in absence of input the system resides in a stable steady state determined by the crossing of the nullclines, Ẏ = 0 and Ż = 0, respectively. Upon external input, X is activated and hence the shapes of the nullclines change. They provide a new crossing and consequently a new steady state. Due to these changes in the nullclines' geometry the system must evolve from the initial state towards the new stable state. The evolution corresponds to a trajectory crossing phase space that depends on i) the location of the initial state, ii) the location of the final state, iii) the new shapes of the nullclines upon input. The specific dynamics implemented by a given motif is determined by this trajectory, which depends on the set of parameters. However, by analyzing the geometrical features of the nullclines it is possible to determine the so-called Backbone of Requirements for the FFL response (BR), i.e. a set of qualitative relationships between different geometrical features of the nullclines and the location of the initial and the final point that univocally determines the dynamics [26]. Therefore, for a given FFL motif all different sets of parameters satisfying the same BR implement the same function. Similarly, also different BRs may implement the same function.
Based on the analysis of the different BRs associated with a given motif and their impact on the motif's function, we are in the position to determine a distribution of probabilities for the implementation of any function (see [26] for details about quantification of the dynamical probabilities). Table 2 shows the conditional probability that each motif Γ i (rows) implements the function j j (columns).

Functional robustness and mutational perturbation
Parametric mutations have different impact on the motif's function, as in nature they can either be neutral or causing qualitative changes. For the system we present here, only those mutations cause functional change, which induce a qualitative alteration in the shape of the nullclines, represented in the Backbone of requirements for the FFL response (BR) [26]. However, the mutation will become visible only, if the resulting BR is actually associated with a different function.
To estimate and compare the degree of mutational robustness for the different FFL motifs, we carried out a numerical study calculating the frequency of functional shifts upon parametric perturbation of equation (6) as shown in figure 5. For a given motif, this can be done introducing random mutations in the parameters that define characteristically the dynamics of FFL motif (figure 5a). We restrict the analysis to that sort of mutations that does not change the topology of the FFL, i.e. mutations that do not change the qualitative type of regulations (activation or inhibition) described by a x , b x and b y .
The suffered mutations are reflected (or not) in a qualitative change of the BR. Here different scenarios are possible, i) mutations that do not affect qualitatively the nullclines' geometry, i.e. there are no changes in the BR (neutral), ii) mutations that are reflected in compatible changes in the BR, but the new BR is associated to the same dynamic than the previous one (neutral), and finally iii) mutations that are reflected in compatible changes in the BR, and the new BR is associated to a different dynamic (qualitatively changing mutations).
In our numerical study we have considered 1000 different sets of parameters for each FFL type (8000 circuits in total). For each FFL, the mutational process is repeated 10.000 times until the probabilities of functional shifts stabilize (figure 5c) and the effects of the accumulation of mutations can be analyzed. This evaluation of the transition probabilities does not depend on specific parameter values but on the conditional relations between them. Since the relations are of the sort a > 6, changes in the function may be achieved by very small or large parameter changes equally likely, Table 2 Conditional probabilities P(j j | Γ i )  Table 2 shows the conditional probability for each motif Γ i (rows) to implement the function j j (columns).These probability distributions are calculated from the association of a given motif, its possible BRs and the respective functions, as described in [26].
depending only on the conditional dependencies between the key-parameters and the corresponding values at time-step t -1. Tables 3 and 4 show the probabilities of transition between the different dynamics subsumed in a transition matrix. These matrices define a transition graph for each FFL motif shown in figure 5b.
The elements of the diagonal correspond to the functionally invariant mutations, i.e. changes in the BR without changes in the dynamics. Rows represent the initial, columns the final dynamics. In tables 3 and 4 we summarize the effects of single mutations, whereas in table 5 and 6 accumulated mutations are presented, i.e. multiple conditions of the BR changed in a successive manner.
The sum of the diagonal elements determines the fraction of the mutations without impact on the dynamics. It provides a measure of the robustness against perturbations. As our data show, all FFLs are highly robust against single mutations, they show similar values above 90%. In other words, upon single mutations the FFLs display low sensitivity against mutation, hence a low evolvability ε. However, there occurs a significant change in the degree of evolvability if the effect of accumulated mutations is studied. We quantify the evolvability ε( i ) of a given topology as: where i ω m k (φ i |φ i ) represents the robustness against accumulated and i 1 k (φ i |φ i ) represents the robustness against single mutations. Since

Entropy, kurtosis and the likelihood of appearance of FFL motifs
The degree of homogeneity of the distribution of probabilities for the implementation of any function can be Table 3 Transition probabilities for single mutations C1-C4  Table 3,4,5 and 6 list the transition probabilities between functions for every motif and its respective function probability distribution. The values are obtained from numerical simulation. Table 4 Transition probabilities for single mutations I1-I4 used to characterize the level of flexibility or specialization or, complementarily, the plasticity of a given motif Γ i . Different measures can be used to quantify the homogeneity of a distribution, namely kurtosis (see Results and Discussion) or entropy. Here we point out some more details on the calculation of the extreme cases and their kurtosis and how the Shannon entropy can be applied equally effective for this task. The expression for kurtosis can be written [30] as: where <P > is the mean value of the probability distribution and s is the standard deviation.
Considering the extreme case where the most specialized motif Γ s only implements a single function j 1 , i.e. P(j j | Γ s ) = 1 if j = 1 and 0 otherwise, expression (8) reduces to: On the other hand, the most flexible motif Γ f implements all possible functions with equal probability P(j j | Γ f ) = 1/n. Here expression (8) leads to a mathematical indetermination 0 0 , because P(j j | Γ f ) = <P >. This indetermination can be solved calculating the limit Applying L'Hopital's rule [31], expression (10) reduces to:   Table 6 Transition probabilities for accumulated mutations I1-I4 The calculated kurtosis values are K(Γ s ) = 6 and K(Γ f ) = -3.33, knowing that n = 6 (number of different possible functions).
Finally, we apply the Shannon entropy [32] to our data set to describe the qualitative differences in homogeneity of the motifs' probability distributions and develop a measure for functional specialization. It is defined as Again, we will first consider the extreme cases Γ s and Γ f to determine the range of all possible values. We find for the most flexible case P(j j | Γ f ) = 1/6 an entropy of H(Γ f ) = log 2 (6). For the most specialized case Γ s (P(j 1 | Γ s ) = 1 and 0 otherwise) the associated entropy is H(Γ s ) = 0. Any FFL motif will have entropy values residing within this range. When correlating the FFLs' entropy and their abundance in figure 4b we find that the most abundant motifs show intermediate values. This trend coincides with what has been found when applying kurtosis, where too, C1 and I1 show intermediate kurtosis values. They do neither exhibit high flexibility nor high specialization, indicating that a trade-off between both features is associated to the most abundant motifs. It can be understood in terms of adaptability or in other words the motifs' plasticity to explore the landscape of possible dynamics under mutational pressure.