Abstract Cellular phenotypes underpinned by regulatory networks need to respond to evolutionary pressures to allow adaptation, but at the same time be robust to perturbations. This creates a conflict in which mutations affecting regulatory networks must both generate variance but also be tolerated at the phenotype level. Here, we perform mathematical analyses and simulations of regulatory networks to better understand the potential trade-off between robustness and evolvability. Examining the phenotypic effects of mutations, we find an inverse correlation between robustness and evolvability that breaks only with nonlinearity in the network dynamics, through the creation of regions presenting sudden changes in phenotype with small changes in genotype. For genotypes embedding low levels of nonlinearity, robustness and evolvability correlate negatively and almost perfectly. By contrast, genotypes embedding nonlinear dynamics allow expression levels to be robust to small perturbations, while generating high diversity (evolvability) under larger perturbations. Thus, nonlinearity breaks the robustness-evolvability trade-off in gene expression levels by allowing disparate responses to different mutations. Using analytical derivations of robustness and system sensitivity, we show that these findings extend to a large class of gene regulatory network architectures and also hold for experimentally observed parameter regimes. Further, the effect of nonlinearity on the robustness-evolvability trade-off is ensured as long as key parameters of the system display specific relations irrespective of their absolute values. We find that within this parameter regime genotypes display low and noisy expression levels. Examining the phenotypic effects of mutations, we find an inverse correlation between robustness and evolvability that breaks only with nonlinearity in the network dynamics. Our results provide a possible solution to the robustness-evolvability trade-off, suggest an explanation for the ubiquity of nonlinear dynamics in gene expression networks, and generate useful guidelines for the design of synthetic gene circuits.

Citation: Steinacher A, Bates DG, Akman OE, Soyer OS (2016) Nonlinear Dynamics in Gene Regulation Promote Robustness and Evolvability of Gene Expression Levels. PLoS ONE 11(4): e0153295. https://doi.org/10.1371/journal.pone.0153295 Editor: Stephen R. Proulx, UC Santa Barbara, UNITED STATES Received: April 13, 2015; Accepted: March 28, 2016; Published: April 15, 2016 Copyright: © 2016 Steinacher et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability: All relevant data are within the paper and its Supporting Information files. Funding: This work was funded by the UK Engineering and Physical Sciences Research Council, grant number EP/I017445/1. Competing interests: The authors have declared that no competing interests exist.

Introduction Biological systems are subject to random mutations as well as noise arising from internal and external stochastic fluctuations. It is important that the potential perturbing effects of noise and mutation are buffered at the phenotype level, i.e. biological systems are expected to display a phenotype that is robust to these perturbations [1–3]. However, changes in fitness pressures over evolutionary time make it similarly important that biological systems are able to produce heritable phenotypic variants that are adaptive. This ability has often been termed evolvability [4–6]. How can biological systems generate phenotypes that are robust to mutations and noise, yet also able to evolve through the effects of these same perturbations [7]? A suggested solution to this trade-off is that robustness of a phenotype to mutations could allow the accumulation of genotypic diversity, which could then translate into phenotypic diversity under subsequent mutations or changes in selective pressures [8]. In support of this idea, analysis of computational models of diverse biological systems has shown that there exist mutationally linked genotypes that display the same phenotype, forming a so-called neutral network [9–11], and that these genotypes can still have access to high phenotypic diversity [8, 12–14]. It has also been found that the size of accessible neutral neighborhood for genotypes determines speed of adaptation [15]. Confirming these theoretical findings, experiments on several biological systems have found these systems to display phenotypes that are robust to most mutations [1, 16, 17], but at the same time able to display high levels of phenotypic diversity under certain mutations [6, 16, 18]. Furthermore, it is shown experimentally that a period of neutral evolution of RNA enzymes under one selective pressure increased the speed of evolution under a different selection pressure [19]. While these findings suggest that robustness can increase evolvability in a population context, they do not provide any mechanistic understanding of how the potential evolvability-robustness trade-off can be broken at the level of genotypes. A clearer understanding of the evolvability-robustness relationship at the genotype level requires the formulation of definitions of robustness and evolvability that permit their quantification across a large class of genotypes and mutations. So far, most computational studies have used a discretized set of mutations to measure the robustness of genotypes and the phenotypic diversity available to them [9, 10, 14, 20–23]. When the phenotypic diversity is used as a proximate measure for evolvability, these studies have found an inverse correlation between robustness and evolvability [21–24]. It is not clear, however, how this result depends on the discretisation method used. In particular, a number of studies indicate that genotypes’ robustness to noise (i.e. intrinsic variation of gene expresssion levels related to infinitesimally small perturbations of system parameters) and mutation (i.e. large perturbations of system parameters, affecting gene expression variability) are interlinked with their evolvability [25–27]. Other studies have used in silico evolution, and measured evolvability as the change in evolutionary fitness with respect to mutation size [24, 28–30] or as the speed or frequency of the emergence of specific phenotypes [9, 31]. Some of these studies have indicated a positive relation between robustness and evolvability [31], but such conclusions are likely to be dependent on the choice of fitness function, and other details of the in silico simulations. In order to overcome these limitations and provide a comprehensive analysis of the relationship between robustness and evolvability at the genotype level, here we develop several mathematically rigorous measures for robustness and evolvability, and apply these to genotypes defined for two common gene regulation network architectures involving a single gene under self-regulation or under the regulation of an upstream transcription factor. Evaluating several million genotypes for each of these systems and using the developed measures, we show that increased nonlinearity in the system dynamics breaks the inverse correlation between robustness and evolvability. We analyze the genotype-phenotype mapping in these circuits by varying the degree of nonlinearity in the equations governing gene expression (see Methods). Thus, we distinguish between genotypes encoding for “linear” and “nonlinear” dynamics and their resulting phenotypes as steady state levels of gene expression. Our finding holds irrespective of the mutational distributions considered for measuring these quantities. We find that robust and evolvable genotypes display low expression levels and occupy a special region, presenting sudden changes in phenotype with small changes in genotype. These findings suggest that nonlinear system dynamics in gene regulation are crucial for maintaining robustness and evolvability of expression levels. Furthermore, they predict that the empirically found correlation between gene expression noise and plasticity [25, 32] results from nonlinearity in gene regulatory systems.

Material and Methods Gene circuit models To study robustness and evolvability in the context of gene regulatory networks, we consider here two network architectures that are commonly observed in nature. This analysis considers the system dynamics of these networks in isolation from other cellular components. The two circuits we consider are described in detail in the following sections. Circuit I—Auto-activation model. In this network architecture, it is assumed that translated protein positively regulates the transcription of its own gene by binding to its cis-regulatory module. Such regulation is common in biology; moreover, synthetic implementations of auto-activation feedback motifs have demonstrated experimentally that they can give rise to bistable system dynamics in which two distinct steady state expression levels are possible [33, 34]. The particular model of auto-activation feedback that we considered comprises four reaction processes: transcription, translation, mRNA degradation and protein degradation. The equations governing the time evolution of the concentrations of mRNA (denoted M) and protein (denoted P) are shown below: (1) Here, transcription is controlled by the maximum transcription rate a, the basal transcription rate k 1 , and the rate of feedback-mediated transcription k 2 . The parameter k 3 denotes the translation rate, while k 4 and k 5 are the mRNA and protein degradation rates respectively. Finally, the Hill coefficient N indicates the degree of cooperativity in the feedback loop and k D quantifies the protein concentration at which activation is half maximal. It is straightforward to show that nondimensionalising the steady state protein level via (2) leads to the following simple equation for steady state expression: (3) In the above, the function is defined by (4) where (5) It follows that steady state expression is determined by the composite parameters α and K, together with the level of nonlinearity N. A detailed stability analysis of the model can be found in the supplementary information (S1 File), where it is shown that the circuit can exhibit both monostability and bistability, depending on the values of α, K and N. Circuit II—Simple-activation model. In this network, gene expression is driven by an external transcription factor (TF) which is assumed to be at steady state (a schematic diagram of this circuit is shown in Fig 1b). This TF could represent, for example, the final component of a signalling pathway, such as the MAPK cascade ([35]). The model equations are (6) where mRNA and protein are represented by M and P respectively, T denotes the concentration of the TF, and all other parameters represent the same processes that they did previously for the auto-activation circuit in Eq (1). Gene activation is again modeled using a Hill-type function in which the Hill coefficient N determines the nonlinearity of the signal response. PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 1. Genotype-phenotype mapping of gene regulation networks. (a): The auto-activation circuit (circuit I) consists of a single gene, whose protein product (P) binds the promoter region to activate the expression of its own mRNA (M). The sink sign indicates mRNA and protein degradation/dilution. (b): A simple-activation circuit (circuit II) in which a transcription factor (purple) binds to the promoter region of a gene to activate the expression of mRNA (M). The protein product (P) is produced by translation of M. (c): A schematic representation of the G-P mapping. Genotypes are represented as discrete entities in the multi-dimensional parameter space of kinetic rates (indicated by the axes k 1 → k n ) and are mapped to their corresponding phenotypes, which are steady state solutions in the two-dimensional space spanned by M and P expression levels. https://doi.org/10.1371/journal.pone.0153295.g001 At steady state, the protein level is determined by (7) while the corresponding steady state mRNA level is (8) There is therefore always only one steady state of the system , and it can be shown that this is always stable (see S1 File). For simplicity, we introduce the following composite parameters (cf. Eq (5)): (9) Note that γ represents the nondimensionalized TF concentration. Also, as in the analysis of the auto-activation circuit, K is the ratio of the TF-mediated and basal transcription rates. The steady state protein level in this case can thus be written in the form (10) where f(γ) is the function defined previously in Eq (4), and is thus a function of β, K, N and γ. ‘Degrees’ of nonlinearity. Both gene regulation networks analyzed here are characterized by nonlinear response dynamics. Our analysis provides arguments for the extent to which nonlinearity shapes the robustness-evolvability relationship. Since such statements could easily be interpreted in different ways, we feel the need to explain in more detail what we mean by ‘nonlinearity’ in the context of this study. In both gene regulation networks at hand, expression of mRNA (M) is dependent on levels of transcription factor (P or T, respectively, for circuits I and II), to the power of N. The parameter N is the Hill coefficient in the kinetic equation, directly relating to the level of cooperativity in transcriptional regulation. From a mathematical viewpoint, expression of M is always nonlinear unless N = 0. It is well known that nonlinear response dynamics can result in ultrasensitive relationships, characterized by a sigmoidal function between signal and response [36, 37]. A system exhibiting nonlinearity leading to ultrasensitivity is characterized by insensitivity of the system response to a stimulus of a certain range of concentration, whereas outside this range the system response might be dramatic. Ultrasensitivity is one specific manifestation of nonlinear response, which can in more general terms be described as a deviation of system response towards stimuli from a perfect straight line [38]. It is this notion of nonlinearity that we refer to when we characterize our system parameter N in its ability to increase the nonlinearity in its regulatory response—it increases the deviation of M expression as a response to increase in transcription factor levels from a straight line, as quantified by the nonlinearity measure L in [38]. Since the regulatory networks featured in this study are specifically showing ultrasensitive response dynamics at high levels of N, the quantification of response coefficient R in [37] might be an even more relevant equivalent to our notion of increasing or decreasing levels of nonlinearity, captured by levels of N. Genotype-phenotype mappings For both circuits, the genotype-phenotype mapping was constructed by discretizing the parameter space, and then calculating the steady state expression level of the system (phenotype) for each possible parameter combination (genotype) (Fig 1c). For circuit I, over 10 million genotypes were simulated. A large subset of these, represented by 7.77 × 106 genotypes were tested based upon parameters in [29]. Of these genotypes, 6.28 × 106 mapped to monostable phenotypes and 1.48 × 106 mapped to bistable phenotypes. Another subset of 3.38 × 106 genotypes were simulated, using parameter ranges based on experimental measurements of the lac operon in Escherichia coli. These yielded 2.85 × 106 monostable phenotypes and 5.2 × 105 bistable phenotypes. For circuit II, 2 million genotypes were considered. Parameter values were obtained by randomly sampling within the specified bounds. For parameters k 1 , k 2 , T and k D , samples were taken from a lognormal distribution; the remaining parameters were sampled from a uniform distribution. The parameter ranges used to construct the G-P mappings are presented in Tables A–C in S1 File. Quantifying robustness and evolvability We use several sets of measures to quantify robustness and evolvability, with the overall aim to capture the effects of mutations of different sizes on expression levels. Before we characterize each measure in detail, we would like to point out that one reason to use several measures for both robustness and evolvability is to show that our outcomes are consistent, ie. not just artifacts of chosen measures. As some measures might only be valid for a certain size of mutations, the finding that genotypes which are robust to mutations of a certain size could be evolvable to mutations of a different size would not provide a conclusive answer to the question whether these genotypes are effectively both robust and evolvable. Therefore, we here present two measures each for robustness and evolvability: one which is valid for infinitesimally small changes in genotype, hence for mutations with a small effect, and another one which is valid for arbitrarily large changes in genotypes. We then provide results based upon combinations between these measures to establish which genotypes are both robust and evolvable irrespective of the size of a mutational effect. Quantifying evolvability. Evolvability can be thought of as the ability to produce variation by mutation. For a gene regulation network, looking at steady state expression levels as phenotype, we can relate this to the ability of a genotype to give rise to different phenotypes, given mutations upon the genotype. Since we model gene expression network by ordinary differential equations based upon mass-action kinetics, genotypes are directly represented as reaction rates and hence parameter sets in differential equations. Hence, mutations can in this framework be represented as parameter changes in the set of equations characterizing a gene regulation network. Within this framework, we can formulate measures for evolvability that capture the variation of steady state expression levels (phenotypes) with regard to changes in the underlying system parameters (mutational effects). Quantifying evolvability for small mutations. Gene expression noise is one way of accomplishing phenotypic variation, and it is due to very small fluctuations in genotypes (ie. Reaction rates). Hence, we use gene expression noise as a measure of evolvability, arising from small fluctuations. We measured levels of intrinsic noise by applying the Linear Noise Approximation (LNA) to genotypes that met the necessary condition for the LNA to hold, namely their Lyapunov stability; this condition is met by genotypes that are monostable in circuit I, and by all genotypes in circuit II. LNA was performed using the approach of [39]. For our analyses, the noise measure was taken to be the covariance of the protein levels P, scaled (normalized) by the respective steady state expression level . Quantifying evolvability for arbitrarily large mutations. While infinitesimally small fluctuations on the genotype level, resulting in expression noise, are able to produce variation on the phenotype level, the LNA framework for quantification of this variational effect is not suitable for larger mutations, since LNA is defined only for infinitesimally small changes in parameters. For larger mutational effects, one way of quantifying evolvability would be to characterize how large the ‘spread’ of expression levels from mutated genotypes is, with regard to the overall size of mutation causing phenotypic changes. The coefficient of variation provides a straightforward quantification of the variation in expression levels, regardless of the expression level of the focal genotype, i.e. the genotype which is mutated. Since this measure is not bounded by size of mutational effect, we employ it for quantifying the effect of arbitrarily large mutations. Following these considerations, we consider, for each focal genotype, a set of mutated genotypes g n such that every member of this set differs by one parameter change from the focal one. Writing P n for the corresponding set of phenotypes (i.e. steady state protein expression levels), evolvability E is defined as the coefficient of variation of this set, scaled by a normalizing factor m: (11) whereby σ and μ represent the standard deviation and the mean of the distributions of protein expression levels, respectively. We used two approaches to apply mutational changes: in the first approach, the set of mutated genotypes g n consists of the 1-mutant neighbors of each focal genotype in the discretized parameter space of the genotype-phenotype mapping. In the second approach, different mutation sizes were considered, such that the set of mutated genotypes consists of fixed-percentage perturbations of single parameters. e.g. perturbations of ±5%, ±10% or ±20%. For all evolvability measurements presented in this study, the factor m in Eq (11) was defined as the number of mutated genotypes that comprise g n . This number can vary between 8 and 16 in circuit I, depending on whether mutated parameters are at the boundaries of the tested parameter vectors or not. We also considered alternative definitions of m, such as the norm of all relative parameter changes (12) where p refers to the parameter values of the focal genotype, p n is the parameter value of its neighbors, and ϵ 1 , ϵ 2 are small numbers (ϵ 1 = ϵ 2 = 10−10) which prevent potential divisions by zero. The results obtained using this measure were qualitatively equivalent. Quantifying robustness. Robustness can be understood as the ability to withstand mutational effects, i.e. to sustain a phenotype amid changes on a genotype level. Considering our framework of kinetic equations describing gene regulation networks, robustness as described above would then relate to the insensitivity of the output of a dynamical system to changes in underlying parameters, reflecting mutational changes. As in the case of evolvability, fluctuations on the genotype level, represented as parameter changes, can be of diverse magnitude and the quantification of robustness will depend on the scale of such changes. Considering further our introductory remarks in this section, to compare robustness and evolvability for a given genotype, we need to ensure that such a comparison is not thwarted by potential artifacts emanating from the usage of different scales for evolvability and robustness measures. Thus, following the above definitions of evolvability, we here present quantifications of robustness on two scales: one on the scale of infinitesimally small genotypic changes, and one to arbitrarily large mutations. Quantifying robustness for small mutations. For mutations of small effect, we quantify the robustness as the inverse of the global sensitivity of the linearized dynamical system that it encodes. Hence, robustness to small mutations is represented as insensitivity to parameter changes in the underlying dynamical system. By focussing on insensitivity, we are able to capture the ability of the system to retain the phenotype of a corresponding focal genotype subject to mutational effects. For infinitesimally small changes at the genotypic level, a commonly used measure is the global sensitivity of the linearized system, and given our interest in insensitivity of this system, we use the inverse of global sensitivity as a quantification of robustness. Following [3], we defined the robustness of a given genotype to mutations of small effect as the reciprocal 1/S of the sensitivity S of the steady state protein level to parameter changes: (13) The global sensitivity S, given in this expression, is the extent of change in steady state protein expression level, given a change in model parameters. Here, ‖⋅‖ 2 represents the standard Euclidian 2-norm, and k = (k i ) is the vector of model parameters. Given an ensemble of zero-mean, independent, identically distributed scaled parameter perturbations, the variance of the corresponding scaled protein levels is approximately given by . The sensitivity S thus quantifies the extent to which protein levels can be adjusted by small bounded fluctuations that affect biochemical reaction rates [40, 41]. The smaller the value of S, the smaller the relative change in under parameter variations, and hence the greater the robustness of the circuit. Throughout all figures of our manuscript, the robustness measure was normalized to have a maximum value of 1. Detailed derivations of S for the circuits considered in this study can be found in the SI. In the case of circuit I, S can be expressed in the form shown below: (14) It follows from Eqs (3) and (4) that the robustness of circuit I only depends on α, K and N. The sensitivity of circuit II is given by: (15) Eqs (4) and (10) therefore imply that the robustness of circuit II only depends on β, K, N and γ. Quantifying robustness for arbitrarily large mutations. To ensure robustness is not just capturing small perturbations around the steady state, but also larger perturbations, representing mutations of larger effect (for instance dramatic reduction of complex stability during protein degradation, captured by the protein degradation rate), we developed an additional measure that does not depend on linearization of the system, but takes a more general approach of how close expression levels of a mutated genotype lies with respect to expression levels of a focal genotype. Congruent with our formulations of evolvability measures, as described above, the definition of global sensitivity in a dynamical system is bounded by the magnitude of parameter changes. And, similarly to our rationale for developing an evolvability measure for arbitrarily large mutations, we developed a similar statistical expression for robustness towards large mutations. In short, this measure captures how similar the expression level of mutated genotypes remain to a focal phenotype. This is quantified by imposing a weighting function around the steady state value of the focal genotype, with perturbations around the genotype assigned probability values specified by this function. Robustness is defined as the mean value of the set of parameter perturbations around the focal genotype mapped onto a normal distribution , with μ r set equal to the steady state value, and a σ r value that reflects the standard deviation of the parameter fluctuations used to calculate evolvability (see Figure A in S1 File). Thus, this alternative robustness measure scores phenotype neighbors by their distance to a focal phenotype under fluctuations in the genotype, such that neighbors with the same or nearly unchanged steady state expression level get a high score, and those further from a focal phenotype get a lower score, based on the mapped Gaussian. The robustness of a focal genotype with respect to a certain perturbation size is then the mean of all phenotypic neighbor scores. One assumption that needs to be made for this alternative robustness measure to be consistent concerns the width of the score function, i.e. the value of σ r . If σ r is small, then the robustness measure is strict and only genotypic fluctuations that lead to phenotypes that are extremely close to the focal phenotype result in high robustness. On the other hand, if σ r is large, nearly all genotypes are relatively robust, since even large deviations from a focal phenotype would earn a nonzero score. Assuming that the magnitude of fluctuations m f applied to a genotype correlates with the magnitude of deviation from a corresponding phenotype, we set σ r = σ(m f ). Thus, if the measure is based upon 10% perturbations to a genotype, we would use σ r = 0.1 in the scoring normal distribution; with 20% perturbations we would use σ r = 0.2. Evolutionary simulations Evolutionary simulations of the auto-activation model followed the implementation described in [29]. A population of 1000 cells was considered. At the start of each simulation, the population was taken to be homogeneous with initial parameters set to the following values: a = 1; k 1 = 0.02; k 2 = 0.2; k 3 = 0.1; k 4 = 0.1; k 5 = 0.002; k D = 50 and N = 0. The population was modeled in a fluctuating environment that switched between selection for high and low protein levels at a constant rate v = 0.05 (corresponding to an environmental switch every 20 generations) for a total of 5000 generations. Fitness under the two environments (w high and w low ) was defined as (16) where P end is the amount of P in a cell at the end of a generation, evaluated using Gillespie’s stochastic simulation algorithm run for 2000 time units. Initial conditions for the simulations were M = P = 0. New populations were produced by randomly drawing a cell from the population, and then cloning it into the new population if it had a fitness value above a random number drawn from . This procedure was repeated until the new population consisted of 1000 cells, thereby maintaining a constant population size. Mutations were assumed to occur at a rate μ, which for all simulations was set to the value 0.01. Mutations were performed by adding a normally distributed random variable to the current parameter value, with the exception of mutations to the Hill coefficient N which were implemented by adding or subtracting 0.5 with equal probability. All parameter values were restricted to be positive. Robustness and evolvability for genotypes arising in the evolutionary simulations were computed using the same genotypic fluctuation sizes as in the G-P mapping they were compared to: robustness against small fluctuations versus evolvability based on 20% perturbations around focal genotypes. Figure G in S1 File shows the relationship between evolvability and robustness for a set of typical evolutionary simulations generated using this method. Simulations and software Genotype-phenotype mappings were constructed using custom software developed in Scientific Python. Steady state values were computed using iterative approximation algorithms, as implemented in the Scientific Python module scipy.optimize. Evolutionary simulations were coded in the C language, and simulation code was taken from [29]. All computations were carried out on desktop computers and a Sun Grid Engine high-performance computing cluster.

Discussion We have performed an extensive analysis of genotype robustness and evolvability using mathematical models of two common gene regulatory networks involving a single gene under auto- (or transcription factor) mediated regulation. We have used the expression levels of this gene as a phenotype and the system parameters controlling expression level as the genotype. Defining several complementary measures for genotype robustness and evolvability under mutations of different size, we have evaluated these properties for several million genotypes for each network architecture. This analysis revealed that for most genotypes, robustness and evolvability display a negative correlation, but there exist a significant number of genotypes for which this trade-off can be broken. This observation holds for all the combinations of the different measures utilized. Furthermore, the identified robust and evolvable genotypes using these fitness-independent measures are also found to emerge under in silico evolution when selection schemes that are shown to facilitate adaptation time are used. This suggests that our fitness-independent measures applied to a genotype-phenotype map are then able to identify genotypes that are evolvable in a population dynamics context and using fitness functions based on that same phenotype (such as adaptation time, or performing well in fluctuating environments, section 4 of the Results). Thus we conclude that the fitness-dependent and the fitness-independent view on evolvability need not be mutually exclusive. We show that the robust and evolvable genotypes found in this analysis are characterized by parameter combinations that confer nonlinearity and ultrasensitivity in system dynamics. Among the system parameters that can determine whether a given genotype confers these characteristics, we find the strongest effect to come from the parameter controlling nonlinearity. This effect extends to the point that the breaking of the robustness and evolvability trade-off is only observed when a certain threshold level of nonlinearity is exceeded. We find that this strong effect comes from the fact that nonlinear system dynamics allow for the emergence of phenotypic cliffs in the genotype-phenotype map, and thereby enable the presence of genotypes that can be robust and evolvable. This result is corroborated by experimental studies on transcriptional circuits controlled by LexA in E. coli, which show that nonlinearity stemming from negative feedback brings about the ability to withstand mutational effects (i.e. robustness), while at the same time enabling capacitance of the system (i.e. evolvability)[43]. This finding is dependent on nonlinear dynamics: upon removal of negative feedback, both mutational robustness and evolvability are reduced. Here, we show the mechanistic basis of the relationship between these two features by directly tuning the degree of nonlinearity in a given architecture. It is evident that the results presented here need to be considered in the context of the robustness and evolvability measures we applied. These measures were chosen to capture the properties of individual genotypes (rather than populations) under different mutational effects, and without invoking reference to organismal fitness. The relationships we found between the different measures of robustness and evolvability are in line with experimental findings. In particular, the correlation we observe between the evolvability measure defined for mutations of large effect and intrinsic noise fits with the empirical finding of a high correlation between environmental plasticity and gene expression noise in yeast [25, 32]. Gene regulation circuits can and frequently do exhibit more complex architectures than the two architectures covered in our study. The latter have, however, gained much attention previously, for instance as model systems for implementation in synthetic biology [44, 45]. Restricting our analysis to these well-researched architectures has the valuable benefit of tractability concerning their mechanistic underpinnings. By focusing on such minimal architectures, we are able to restrict the search space for potential mechanisms which can affect the relationship between robustness and evolvability. Thus, we find that even the simplest cases of nonlinear gene expression circuits are capable of solving the robustness-evolvability trade-off. Our findings provide a genotype-based resolution of the robustness-evolvability trade-off, and do not contradict previous suggestions based on population dynamics [6, 8, 12–16, 18]. In particular, our results relate to a recent population genetics study suggesting that there can be an optimal level of robustness that promotes evolvability, depending on the properties of the fitness landscape [15]. Using fitness landscapes with specified statistical properties, that study found that while robustness is negatively correlated with evolvability when mutations allow access to any of the possible phenotypes, there can be a positive correlation between robustness and evolvability when mutations can allow access to only a small fraction of all phenotypes. In light of those findings, it is interesting to see that our biologically well-defined genotype-phenotype mapping contains robust and evolvable genotypes only at increasing levels of nonlinearity in the equations governing gene expression levels. It is possible that increasing nonlinearity re-shapes this genotype-phenotype mapping in a way that is in line with the statical assumptions made in [15]. Furthermore, we find that robust and evolvable genotypes identified in our analysis occupy specific regions of the genotype-phenotype map, characterized by sudden changes in phenotype with small changes in genotype. These genotypes also present a specific level of nonlinearity and other kinetic parameters. This could relate to them being “tuned” to have a specific level of robustness as suggested in [15]. This proposition is also reminiscent of the idea of biological systems being in a critical state that enhances their evolvability [21–23, 46], but does not rely on the presence of chaos in system dynamics. It is also possible to draw an analogy between the findings presented here and changes observed when considering a system passing through distinct dynamical regimes due to parameter alterations [47]. It must be noted, however, that in our study the analysis is restricted to a single dynamical regime. In particular, the finding that the ability of nonlinearity to break the robustness-evolvability trade-off extends to systems that do not allow for bistability shows that proximity to a bifurcation surface is not a necessary condition for genotypes to be robust and evolvable at the same time. The mechanistic insights gained from this study could provide an improved explanation for the emergence of evolvability in laboratory evolution experiments [48]. In particular, our findings would suggest that genetic mutations identified in these studies could relate to their effects on noise and nonlinearity in gene regulation. Future experimental studies towards confirming this suggestion will significantly improve our understanding of how natural gene regulatory systems achieve robustness and evolvability, and allow us to better design robust synthetic gene regulatory circuits [49].

Conclusion The question of how biological systems are able to withstand mutational changes (robustness), yet still remain able to produce phenotypic variation (evolvability) has gained considerable interest in recent years. Despite important insights towards the understanding of scenarios that allow biological networks to be both robust and evolvable, the mechanistic underpinnings of such phenomena are still elusive. By performing mathematical analyses of common gene regulatory network motifs, we provide a mechanistic explanation of the robustness-evolvability trade-off. Using several measures across scales of perturbation, we find that nonlinearity consistently breaks the predominant negative correlation found between robustness and phenotypic variability. This holds for both small and large perturbations in genotypes with low expression levels. Our results provide a potential link between the abundance of nonlinearity in biological regulatory networks and their apparent ability to be both robust and evolvable at the same time.

Supporting Information S1 File. Supplementary information containing detailed mathematical analyses of gene regulation networks, as well as supplemental figures and tables. https://doi.org/10.1371/journal.pone.0153295.s001 (PDF)

Acknowledgments We acknowledge insightful discussions with Joanna Masel, Frank Bruggeman and members of the Soyer lab. This work has made use of the resources provided by the University of Exeter Science Strategy and resulting Systems Biology initiative. Primarily these include HPC facilities managed by Konrad Paszkiewicz of the College of Environmental and Life Sciences and Pete Leggett of the University of Exeter Academics services unit. We thank the two anonymous reviewers for their valuable comments that helped to improve the manuscript.

Author Contributions Conceived and designed the experiments: AS DGB OEA OSS. Performed the experiments: AS OEA. Analyzed the data: AS OEA. Contributed reagents/materials/analysis tools: AS OEA OSS. Wrote the paper: AS DGB OEA OSS.