From gene regulatory networks to population dynamics: robustness, diversity and their role in progression to cancer
11 From gene regulatory networks to population dynamics:robustness, diversity and their role in progression to cancer
Tom´as Alarc´on , ∗ and Henrik Jeldtoft Jensen ∗ E-mail: Corresponding author [email protected]
Abstract
The aim of this paper is to discuss the role of robustness and diversity in population dynamics in particularto some properties of the multi-step from healthy tissue to fully malignant tumours. Recent evidenceshows that diversity within the cell population of a neoplasm, a pre-tumoural lession that can developinto a fully malignant tumour, is the best predictor for its evolving into a tumour. By studying thedynamics of a population described by a multi-type, population-size limited branching process in termsof the evolutionary formalism, we show some general principles regarding the probability of a residentpopulation to being invaded by a mutant population in terms of the number of types present in thepopulation and their resilience. We show that, although diversity in the mutant population poses a barrierfor the emergence of the initial (benign) lession, under appropiate conditions, namely, the phenotypesin the mutant population being more resilient than those of the resident population, a more variablenoeplastic population is more likely to be invaded by a more malignant one. Analysis of a model ofgene regulatory networks suggest possible mechanisms giving rise to mutants with increased phenotypicdiversity and robustness. We then go on to show how these results may help us to interpret some recentdata regarding the evolution of Barrett’s oesophagus into throat cancer.
Author Summary
Recent results by Maley and co-workers [1] regarding progression of the Barrett’s esophagus, a benign,neoplastic throat condition, into a malignant tumour have revealed that clonal diversity within diseasedtissue is the best risk predictor outperforming usual genetic markers. So far, there is no detailed expla-nation as to why this is so. What is the role of diversity in the progression to cancer and what are themechanisms involved? In this paper, we address these issues and put forward some generic mechanismsthat may help to gain a better understanding of the empirical results obtained by Maley et al. Our ap-proach to this problem is in two steps. First, we propose a simple model of gene regulatory networks andanalyse their behaviour upon gene silencing. We observe that, following gene silencing, both phenotypicdiversity and robustness increase. We then proceed to analyse how these two effects affect the dynamicsof the cell populations and find that more diverse pre-malignant lessions are likely to evolve into a moremalignant form only if phenotypic robustness increases, hence suggesting that robustness is essential tothe results put forward by Maley and coworkers.
Introduction
Complex diseases such as cancer pose an enormous scientific challenge. Cancer involves an extra level ofcomplexity as its development involves evolutionary processes driven by Darwinian natural selection [2].The traditional view of the emergence of cancer consists of a series of mutations with each of them a r X i v : . [ q - b i o . P E ] M a r increasing the fitness of the mutated cells with respect to that of their normal counterparts, which areeventually taken over by succesive rounds of clonal expansion [3, 4]. Recent experimental analysis on thegenetic landscape of tumours [5], however, show an even more complex situation whereby tumours havebeen found to be much more genetically diverse than expected where low-frequency mutations are likelyto play a leading part in the evolution of tumours. Furthermore, recent analysis of data from patients ofBarrett’s oesophagus, a neoplastic condition that may evolve into throat cancer, has found that the bestpredictor for its evolution into a fully malignant lession is diversity in the cell population, outscoring allof the usual genetic markers [1]. In particular, Maley et al. have shown that quantities normally usedin ecology to quantify diversity, such as the Shannon entropy, H S , and the number of clones detectedin tissue samples extracted from patients of Barrett’s oesophagus, are excellent predictors. Accordingto [4], these findings bring into question the vision of tumour evolution as a series of succesive roundsof mutation/clonal evolution to favour a model in which premalignant lessions sustain a number of co-existing cellular types.In order to explore some of the issues put forward by results of Maley et al. [1], in particular theirresults concerning the higher probability of pre-malignant lessions with larger diversity to evolve intotumours, we first analyse the robustness process of a simple model of gene regulatory networks. Wefocus on the study of the effects of gene silencing on both the number of stable phenotypes and theirrobustness [6–11]. We will show using a simple model, similar to the ones used in previous relatedworks [7–9, 12], that such modification on the gene regulatory network leads to an increase on both thenumber of stable (robust) phenotypes and their ability to sustain further mutations (robustness).We then analyse how these changes induced by gene silencing have on the dynamics of the corre-sponding population, in particular on the ability of one mutant carrying such modifications to invadethe resident cell population. To so doing, we consider a model of population dynamics consisting of amulti-type branching process with population-size limited proliferation probability [13]. Our aim is toput forward a number of general properties of multi-type populations which could shed some light onthe mechanisms involved in the transition from neoplasm to malignant tumours. Our analysis is car-ried out in terms of the so-called evolutionary formalim developed by Demetrius and coworkers [14–16].This framework extends the thermodynamic formalism developed within ergodic theory and statisticalphysics [17] and it allows us to study the evolutionary behaviour of structured populations in terms ofthe evolutionary entropy , which is a measure of the dynamical diversity of the population.The evolutionary formalism, whose main results concerning the present problem are summarisedin Materials & Methods, has been recently applied to analyse the competition between resident andmutant populations modelled in terms of age-structured population models [16] and multi-type branchingprocesses [13]. Demetrius et al. [15,16] have shown that the probability of a mutant structured populationto invade the resident one is determined by the rate at which a population returns to its correspondingsteady state after a random perturbation of the parameters that characterise the population dynamics.In turn, it has been proofed that such rate is given by the evolutionary entropy [15]. They have alsoshown that the evolutionary stability of a population, that is, its resilience to invasion by a mutant, isdetermined by the evolutionary entropy [16].Our aim in this paper is to understand the dynamical mechanisms involved in the findings of Maley etal. [1]. In other words, we intend to find how the Shannon entropy, which is an index of diversity withinthe population, relates to the evolutionary entropy, which determines the evolutionary stability (or lackthereby) of the corresponding population. Results
Robustness of gene regulatory networks
In order to model the dynamics of gene regulatory models, we consider a network with N G nodes, eachrepresenting a gene. To each of the nodes we associate a state g i = ± , i = 1 , . . . , N G , with state g = 1( −
1) corresponding to the gene being activated (inactivated). We generate a directed network,defined in terms of the corresponding adjacency matrix, according to a prescribed degree distribution(in the present case the in- and out-degrees are exponentially distributed). The links of this networkcorrespond to the interactions between genes, which can be activatory or inhibitory. If a gene i activatesgene j , the corresponding link, represented in a matrix W = ( w ij ) whose entries are zero if there is nolink between nodes i and j and non-zero otherwise, carries a weight w ij = 1. If, on the contrary, gene i inhibits gene j , the corresponding weight w ij = −
1. In our model, positive and negative weights arerandomly assigned with probability p + and 1 − p + , respectively.The dynamics of the gene regulatory netowrk is defined in terms of the above elements as follows. Aninitial condition g i ( t = 0) , i = 1 , . . . , N G and a matrix W are given. For each node and at each time stepthe following quantity is computed: I i ( t ) = (cid:80) j ∈(cid:104) i (cid:105) in w ji g j ( t ). Then at each time step: g i ( t + 1) = I i ( t ) > − I i ( t ) < g i ( t ) if I i ( t ) = 0 (1)which is ran until the system reaches a steady state. The phenotype, φ , is defined as the vector φ =( g , . . . , g N G ) at steady state. This steady state may be a periodic solution. How we deal with this typeof solution will be explained in more detail later.In our model, we consider two type of mutations. On the one hand, following [7, 9, 12], we introducerandom changes on how genes regulate each other (i.e. w ij → − w ij with a certain probability). Thesecond type of mutation we consider is a form of gene silencing where one particular gene is chosen (inprinciple, at random) and its state is fixed at g i = − φ , which we defineas the probability that a mutation of the first type according to the description in the paragraph aboveinduces a phenotype switch (defined as a change in sign of at least one of the components of the vector φ ),we carry out a numerical experiment using the algorithm defined in the Materials & Methods section. Notethat the population dynamics defined by this algorithm is completely neutral, since all the individualsyield offspring with the same probability and we have not imposed any viability conditions as assume e.g.in [7]. The corresponding results are shown in Fig. 1. We see that even in neutral conditions the diversityof phenotypes within the population decreases, while the extant types exhibit a greater resilience againstmutations. In other words, the number of different phenotypes decreases with time but the survivingones are such that their robustness increases, i.e. we observe the emergence of canalisation [18,19] withinour simple model.These results can be understood in terms of the concept of the neutral network as introduced in [20,21].The neutral network is a meta-network where each node corresponds to a different choice of W and twonetworks are connected if they differ only in one of the values of the corresponding entries w ij . It hasbeen shown [20] that phenotypic robustness can be understood in terms of the structure of such network:robust phenotypes are those whose subnetwork is highly interconnected, i.e. mutation events are likelyto produce a network with the same phenotype. This means that robust phenotypes act like traps withinthis network making it harder for mutations to drive the system out of the corresponding phenotype. Thesame topological mechanism causes non-robust phenotypes to disappear at early stages in the evolution ofthe population. This is because mutations will drive them towards network configurations representingthe more stable phenotypes. Thus, we can qualitatively explain both why the number of phenotypesdecreases in time and also how the remaining phenotypes become more resilient.We now repeat the same simulatons with the difference that at a given point in time we randomlyselect one of the nodes, say node i , and introduce a gene silencing event, i.e. g i = − Results for multi-type branching process with resource limitation.
In order to advance further our discussion, we focus on a multi-type branching process with resourcelimtation. In the context of this paper, the different types in which the population divides correspondto the robust phenotypes. We further assume that all the (pheno)types yield offspring with the sameprobability, namely that of an individual in our original population model. For the sake of example, wewill consider two different populations. The incumbent population is defined by an strategy whereby thepopulation splits in d I types. The average number of offspring to be produced by a given indvidual isgiven by: A = 2 e − µN − ν νd I − · · · νd I − νd I − − ν · · · νd I − νd I − · · · · · · νd I − νd I − · · · · · · − ν (2)where ν ij is the mutation probability per generation from type i to type j . For simplicity we will assume ν ij = ν ∀ i, j . Here N = N I + N M , i.e. the total population (both resident and mutant). In all thesimulations shown in this paper the initial mutant population is N M = 1.The mutant population adopts a strategy with a larger variety of available types ( d M ≥ d I differenttypes) but which may have an increased robustness (smaller type mutation probability). Otherwise, theproliferation probability is the same as for the incumbent population. Thus: B = 2 e − µN e ϕ (1 − ν + ρ ) e ϕ ν − ρd M − · · · e ϕ ν − ρd M − ν − ρd M − − ν + ρ · · · ν − ρd M − ν − ρd M − · · · · · · ν − ρd M − ν − ρd M − · · · · · · − ν + ρ (3)In Eq. (3), ρ is a measure of the increase in resilience of the corresponding phenotype. The parameter φ corresponds to an increment in the growth rate of one of the phenotypes, which can be either positiveor negative. The remaining phenotypes produce offspring at the same rate as the phenotypes of theincumbent population.With relation to the results of Maley et al. on the Barrett’s esophagus [1], the competition between thetwo populations in our model could represent the initial emergence of the neoplasm (mutant population)which overtakes the normal tissue (incumbent population).There are a number of results that we can obtain from numerical simulation of the correspondingbranching process (see Materials & Methods for details). First, increased diversity ( d M > d I ) withoutincreased resilience results in loss of the ability of the mutant to invade the resident population, this isbecause P F decreases when d M increases (see Fig. 3). In other words, increasing the phenotypic diversityinduces an error catastrophe-like behaviour. Furthermore, increasing resilience (i.e. increasing ρ ) rescuesthe more diverse mutant populations from the error catastrophe, so that they can now invade populationswith less diversity, as shown in Fig. 3. Both these results can be explained in terms of the evolutionaryformalism [13]. We can also observe (Fig. 4) that H S and P F are negatively correlated: the largerthe Shannon entropy of the mutant population the less likely fixation is. In other words, high-diversitymutant populations are less likely to invade and reach fixation.We consider now the stability and resilience to invasion of the mutants that have achieved fixation.Within the context of the problem we are dealing with, namely, the transition from neoplastic lesions tofully malignant tumours, this transition requires further mutations that perturbs the neoplastic cells anddrives them into the malignant state. We model this type of mutations as perturbations in the parametersthat determine the population dynamics. In particular, we consider that a second phenotype acquires aslight advantage in its proliferation probability quatified by a factor e ϕ . We assume that perturbationsin the resilience of the phenotypes is represented by an additive term ∆ ρ . The mean-field dynamics (seeAppendix) of this second mutant population is thus described by the matrix:¯ B = 2 e − µN e ϕ (1 − ν + ρ + ∆ ρ ) e ϕ ν − ρ − ∆ ρd M − · · · e ϕ ν − ρ − ∆ ρd M − e ϕ ν − ρ − ∆ ρd M − e ϕ (1 − ν + ρ + ∆ ρ ) · · · e ϕ ν − ρ − ∆ ρd M − ν − ρ − ∆ ρd M − · · · · · · ν − ρ +∆ ρd M − ν − ρ − ∆ ρd M − · · · · · · − ν + ρ − ∆ ρ (4)We consider now essentially the same problem as above, namely, the likelyhood of a population whose(deterministic) dynamics is described by ¯ B to take over a resident population whose dynamics is given by B (Eq. (3)). Using the methods described at length in our previous work [13], the relative fitness s (seeMaterials & Methods for the precise definition of this quantity) of ¯ B with respect to B , which characterisesthe fixation probability of the former when competing with the latter, can be computed. The results areshown in Fig. 5. We observe that in the case in which the original mutation (that is, the one givingrise to the neoplastic lesion) increases the diversity of the population, i.e. d M , without a concomitantincrease in the phenotypic resilience, i.e. ρ = 0 (see Fig. 5, red dashed line), the ability of the neoplasmfor further evolution decreases with increasing d M . In other words, the robustness of the neoplastic lesionincreases with d M , as the probability of it being invaded is a monotonically decreasing function of d M . If,on the contrary, the original mutation, in addition to increasing the number of cellular types within thepopulation, also increases their resilience, i.e. ρ >
0, the robustness of the resulting population decreaseswith increasing d M , as the corresponding relative fitness, and therefore the fixation probability, is amonotinically increasing function of d M (see Fig. 5, black solid line). Numerical results correspondingto this situation are presented in Fig. 6 and show agreement with the theoretical predictions of theevolutionary formalism. Discussion
In this paper we have analysed the effects of gene silencing by means of gene silencing on a model ofgene regulatory networks and we have studied how it affects phenotypic diversity and robustness. Wehave then investigated how these changes can modify the ability of a resident population to withstandinvasion by mutants, and, last, we have discussed how our general results may shed some light on recentresults regarding the evolution from neoplastic lesions to fully-malignant tumours.We have shown that, within our model, gene silencing leads to both an increase in phenotypic diversityand also enhances the robustness, i.e. in the populations ability to sustain gene mutations. Furthermore,we have done this in a completely neutral model with no viability conditions [7]. Our neutral modelallows us to compute a lower bound on the robustness.We have then proceeded to analyse the consequences of the above observations for the competitionbetween resident and mutant populations and on the ability of the mutant population to further evolveinto more malignant variants. We have done this within the framework of multitype branching processewith resource limitation and by applying the evolutionary formalism. We have shown that increased di-versity, as measured by the corresponding Shannon entropy, correlates with decreased fixation likelyhoodfor the initial invasion (i.e. the emergence of the neoplasm, e.g. Barrett’s oesophagus). If, in addtionto an increased number of phenotypes, there is a concomitant increase in the phenotypic resilience, i.e. ρ (cid:54) = 0, the likelyhood of invasion increases. However, even the latter case where invasion is more likely,the same general pattern of anti-correlation between Shannon entropy and fixation probability remains(see Figs. 3 and 4). Furthermore, we have shown that both increased phenotypic resilience and increaseddiversity correlate with lower robustness of the population and therefore they are more likely to be furtherinvaded.These general results on population dynamics can help us to make sense of some recent data regardingthe progression of the Barrett’s oesophagus into a fully malignant throat cancer [1]. First, less than 5%of the cases of Barrett’s oesophagus actually evolved into cancer. This fact can be understood in termsof our results: lesions with increased diversity, which are the more likely ones to evolve into a tumour,are going to be less frequent, as their fixation probability is smaller. Furthermore, we predict that thefact that more diverse lesions are more likely to evolve into cancer implies that the event leading to thepre-malignant lession must also increase the phenotypic resilience, as that induces loss of robustness inmore diverse populations. Materials and Methods
Dynamics of the population of networks
In order to analyse the robustness properties of thecorresponding phenotypes, φ , which we define as the probability that a mutation of the first type accordingto the description in the paragraph above induces a phenotype switch (defined as a change in sign ofat least one of the components of the vector φ ), we carry out a numerical experiment defined by thefollowing algorithm:1. An initial population of N W -matrices, i.e. directed networks, are randomly generated as describedabove. Each of these matrices is assumed to define a particular individual within the population.2. The gene regulatory network dynamics described above is iterated for a number of time steps (fixedto a much larger value than typically needed to reach a steady state) for each of the N individuals.The vector ( g , . . . , g N G ) at the end of iteration round is taken to be the phenotype correspondingto each individual within the population.3. The population is next binned according to the corresponding phenotype, with N φ j being thenumber of individuals with a particular phenotype φ j = { g i , i = 1 , . . . , N G }
4. Each individual produces two offspring with probability p o = e − µN ( t ) , where µ is the carryingcapacity5. Upon division, each regulatory interaction is mutated ( w ij → − w ij ) with probability p m = r m /E ,where r m is the mutation probability per generation per edge and E is the number of non-zeroentries in the matrix w ij
6. If no mutation occurs, two new individuals are incorporated in the next generation with the same w ij and φ i
7. If mutations happen, steps 1 to 3 are repeated and new individuals are created accordingly
Summary of the evolutionary formalism.
Here we present a brief summary of the main results andformulae of the evolutionary formalism. For a full account of the details we refer the reader to [14–16].For specific results concerning the particular model we are analysing here, see [13].Consider a population with d I different types whose dynamics is given by iteration according to u ( t + 1) = A ( t ) u ( t ), where u ( t ) = ( n ( t ) , . . . , n d I ( t )) T is the vector consisting of the population sizes ofeach of the types and A = ( a ij ) is the average number of offspring of type i produced per generation byindividuals of type j . Note that, in general, the entries of this matrix depend on the total population N I = (cid:80) d I i =1 n i . We will consider the long-time behaviour of the population, i.e. once it has settled in asteady state. We will further assume that A is an irreducible matrix.The long time behaviour of the population is determined by the dominant eigenvalue of the matrix A , which will denoted by λ . Let u and v be the corresponding left and right eigenvectors, respectively: Au = λ u and vA = λ v . We can now define an associated Markov matrix, P , such that: p ij = a ji v j λ v i (5)One of the main results consists in a variational principle, which states that the following relationbetween the growth rate r = log λ , the evolutionary entropy, H E , and the so-called proliferative potential, F , [14]: r = H E + FH E = d I (cid:88) i,j =1 π i p ij log p ij F = d I (cid:88) i,j =1 π i p ij log a ij (6)where π i is the stationary distribution corresponding to the Markov process defined by P . The requirementthat A should be irreducible ensures the uniqueness of this distribution.A number of interesting results follow from Eq. (6) and the corresponding variational principle, forexample, a proof concerning large deviations states that the rate of relaxation to the steady state aftera perturbation of the entries of the matrix A , specifically, a perturbation of the type A ( δ ) = ( a δij ), ispositively correlated with the corresponding change in the evolutionary entropy H E ( δ ) − H E ( δ = 0), thusproving that this quantity is a measure of robustness of the population [15].However, we are more interested in a (related) result concerning the evolutionary stability of anincumbent population against an invading mutant population when both compete for a common limitedresource. By means of a diffusion approximation, Demetrius et al. [16] have showed that the fixationprobability of the mutant, P F is given by: P F ( y ) = 1 − (cid:16) − ∆ σ σ M y (cid:17) (cid:104) N (cid:105) s ∆ σ +1 − (cid:16) − ∆ σ σ M (cid:17) (cid:104) N (cid:105) s ∆ σ +1 (7)where the total population is assumed to be constant, y is the initial proportion of mutants, and s isdefined by: s = ∆ r − ∆ σ (cid:104) N (cid:105) (8)where (cid:104) N (cid:105) is the average stationary population. In the case of the branching process we consider in thispaper, (cid:104) N (cid:105) is given by 2 e − µ (cid:104) N (cid:105) = 1, i.e. µ (cid:104) N (cid:105) = log 2. The quantity ∆ r is defined as ∆ r ≡ r M − r I , wherethe subindices I and M denote the corresponding quantities for the incumbent and mutant populations,respectively. The quantity ∆ σ is defined as ∆ σ = σ M − σ I where [13]: σ = − dH E ( δ ) dδ (cid:12)(cid:12)(cid:12)(cid:12) δ =0 (9)The behaviour of P F can be analysed in terms of the quantity s . If s > P F ( y ) is convex and P F ( y ) > y . If, on the contrary, s < P F ( y ) is concave and P F ( y ) < y . This means that if s < s > s increases. Moreover, the larger the value of s the more likelyfixation is. Therefore s can be taken as a measure of fitness. Definition of the muitype branching processes
To precisely define the branching processes in-volved in our models of population dynamics, we use the corresponding generating function formalism [22]where a generating function for the probability distribution of the per individual offspring number is pre-scribed. In this article we consider three different populations (corresponding to the matrices A , B , and¯ B ). The set of generating functions corresponding to the population associated to the matrix A is: g i ( (cid:126)z ) = (1 − e − µN ) + (1 − ν ) e − µN z i + (cid:88) j (cid:54) = i νd I − e − µN z j (10)where i = 1 , . . . , d I and (cid:126)z is a d I -dimensional vector with 0 ≤ z ≤ B and ¯ B are, respectively: g ( (cid:126)z ) = (1 − e − µN + ϕ ) + (1 − ( ν − ρ )) e − µN + ϕ z + (cid:88) j (cid:54) =1 ν − ρd M − e − µN + ϕ z j g i ( (cid:126)z ) = (1 − e − µN ) + (1 − ( ν − ρ )) e − µN z i + (cid:88) j (cid:54) = i ν − ρd M − e − µN z j for i ≥ g i ( (cid:126)z ) = (1 − e − µN + ϕ ) + (1 − ( ν − ρ − ∆ ρ )) e − µN + ϕ z i + (cid:88) j (cid:54) = i ν − ρ − ∆ ρd M − e − µN + ϕ z j for i ≤ g i ( (cid:126)z ) = (1 − e − µN ) + (1 − ( ν − ρ − ∆ ρ )) e − µN z i + (cid:88) j (cid:54) = i ν − ρ − ∆ ρd M − e − µN z j for i > i = 1 , . . . , d M and (cid:126)z is now a d M -dimensional vector with 0 ≤ z ≤ i , with probability given by the coefficient of the y i - or z i -terms or of a differenttype, say type j , with probability equal to the coefficent of the y i - or z i -terms.Note that the entries matrices A , B , and ¯ B correspond to the branching ratios m ij = ∂ j g i ( ) as definedin Eqs. (10), (11) and (12), respectively, and therefore they determine the mean-field dynamics. Acknowledgments
TA and HJJ gratefully acknowledge the EPSRC for funding under grant EP/D051223.
References
1. Maley CC, Galipaeu PC, Finley JC, Wongsurawat VJ, Li X, et al. (2006) Genetic clonal diversitypredicts progression to esophageal adenocarcinoma. Nature Gen 38: 468-473.2. Nowak MA (2006) Evolutionary dynamics. Harvard University Press, Cambridge, Mass., USA.3. Attolini CSO, Michor F (2009) Evolutionary theory of cancer. Ann N Y Acad Sci 1168: 23-51.4. Shibata D (2006) Clonal diversity in tumour progression. Nature Gen 38: 402-403.5. Wood LD, Parsons DW, Jones S, Lin J, Sjoblom T, et al. (2007) The genomic landscapes of humanbreast and colorectal cancers. Science 318: 1108-1113.6. Rutherford SL, Lindquist S (1998) Hsp90 as a capacitor for morphological evolution. Nature 396:336-342.7. Siegal ML, Bergman A (2002) Waddington’s canalization revisited: Developmental stability andevolution. Proc Natl Acad Sci 99: 10528-10532.8. Bergman A, Siegal ML (2003) Evolutionary capacitance is a general feature of complex gene net-works. Nature 424: 549-552.9. Ciliberti S, Martin OC, Wagner A (2007) Robustness can evolve gradually in complex regulatorygene networks with varying topology. PLoS Comput Biol 3: e15.10. Wagner A (2008) Neutralism and selectionism: A network-based reconciliation. Nature ReviewsGenetics 9: 965-974.011. Levy SF, Siegal ML (2008) Network hubs buffer environmental variation in saccharomyces cere-visiae. PLoS Biology 6: e264.12. Wagner A (1996) Does evolutionary plasticity evolve? Evolution 50: 1008-1023.13. Alarc´on T, Jensen HJ (2009) Invasion in multi-type populations: The role of robustness andfluctuations. Submitted to Math Med Biol .14. Arnold L, Gundlach VM, Demetrius L (1994) Evolutionary formalism for products of positiverandom matrices. Ann Appl Prob 4: 859-901.15. Demetrius L, Gundlach VM, Ochs G (2004) Complexity and demographic stability in populationmodels. Theor Pop Biol 65: 211-255.16. Demetrius L, Gundlach VM, Ochs G (2009) Invasion exponents in biological networks. Physica A388: 651-672.17. Ruelle D (1978) Thermodynamic formalism. Encyclopedia of Mathematics and its applications,vol. 5. Addison-Wesley, New York, USA.18. Jablonka E, Lamb MJ (2005) Evolution in four dimensions. Massachusets Institute of TechnologyPress.19. Waddington CH (1942) Canalisation of development and the inheritance of acquired characters.Nature 150: 563-565.20. Ciliberti S, Martin OC, Wagner A (2007) Innovation and robustness in complex regulatory genenetworks. Proc Natl Acad Sci 104: 13591-13596.21. Wagner A (2007) Robustness and evolvability in living systems. Princeton University Press, Prince-ton, NJ, U.S.A.22. Kimmel M, Axelrod DE (2002) Branching processes in Biology. Springer-Verlag, New York, U.S.A.1 P h e no t yp i c S e n s iti v it y K = 50000 200 400 600 800Time (Generations)00.050.10.150.2 P h e no t yp i c A bund a n ce Figure 1.
Simulations of the evolution of robustness and phenotypic diversity in a model of generegulatory networks (see text for details). The upper plot corresponds to the phenotypic sensitivity tomutations, i.e. the probability that a mutation induces a change in phenotype, which is an inversemeasure of robustness: the more sensitive, the less robust the phenotype. The lower plot corresponds tothe phenotypic abundance defined as the number of phenotypes that are actually present within thepopulation at generation t divided by the total population. We observe that both quantities decreasewith time. Parameter values: N G = 10, µ = K − = 0 . r m = 12 P h e no t yp i c S e n s iti v it y K = 50000 200 400 600 800Time (Generations)00.050.10.150.2 P h e no t yp i c A bund a n ce Figure 2.
Simulations of the evolution of robustness and phenotypic diversity in a model of generegulatory networks (see text for details). The upper plot corresponds to the phenotypic sensitivity andthe lower plot to the phenotypic abundance (see caption of Fig. 1). In these simulations we haveintroduced a gene silencing event (see main text for details) at generation t = 250, which is released at t = 500. Parameter values: N G = 10, µ = K − = 0 . r m = 13 d M P F ρ P F Figure 3.
This figure shows the fixation probability, P F , as a function of d M corresponding to φ = 0 . µ = 0 . ρ = 0 (squares) and ρ = 0 .
49 (circles) and µ = 0 .
001 (squares) with δ = 0. The insetshows the fixation probability for ϕ = 0 .
01 and µ = 0 . ρ .4 P F S h a nnon e n t r opy , H S Figure 4.
This plot shows the negative correlation between the average Shannon entropy of themutant population and the corresponding fixation probability. Circles correspond to simulation resultsfor ρ = 0 and squares to ρ = 0 . d M s Figure 5.
This plot shows the invasion fitness, s , as a function of d M . The red dashed line and theblack solid line correspond to ρ = 0 and ρ = 0 .
49, respectively. ∆ ρ = 0 . d M P F Figure 6.
This plot shows the fixation probability of a population described by ¯ B over a B -population, s , as a function of d M . Squares correspond to ρ = 0 and circles to ρ = 0 .
49, respectively. ∆ ρ = 0 ..