Complex heritable phenotypic variation fuelling adaptive evolution citep{huxley1942evolution,

Complex living systems are inherently robust to internal and external changes. In particular, mutational phenotypic robustness is widespread across organisms at various level of organization, from molecules to individuals citep{wagner2011origins}. Phenotypic robustness is a key feature deriving from the genotype-phenotype map structure, namely from the fact that many genotypes map on the same phenotype, leading to a reduction of the probability that a genetic mutation produces a novel phenotype. Adaptive evolution requires heritable phenotypic variation for selection to act upon, and the standing paradigm that emerged from Modern Synthesis argued that random genetic mutations of fixed phenotypic effects are one of the most important source of heritable phenotypic variation fuelling adaptive evolution citep{huxley1942evolution, kutschera2004modern, pigliucci2010evolution}. Thus, phenotypic robustness may seem to act as a buffer on mutations and consequently on phenotypic variation and adaptive evolution citep{draghi2010mutational}. So far, phenotypic robustness has received little attention in standard population genetic models, and only recently long-term effects of robustness on evolvability and innovability have started to be studied more in depth. However, except for some few works, the mechanisms by which robustness might be established during evolution are far from clear and overall little explored citep{wagner2011origins}.In a previous study, we showed that including phenotypic robustness in standard population genetics and evolutionary models (the Quasi-species model and the Price equation, cite{price1970selection, eigen1989molecular}) can explain the widespread phenotype resistance to genetic and environmental perturbations. Specifically, and counterintuitively, that a certain (sizable) level of phenotypic robustness is not only a favourable but also a necessary condition (although not sufficient) for adaptation to occur. This appears as a threshold effect, i.e. as a critical level of phenotypic robustness under which evolutionary adaptation cannot occur, even in the case of significant positive selection coefficients and in absence of any drift effect. These theoretical results were obtained under the assumptions of fixed evolutionary parameters (e.g., mutation rates) and large, virtually infinite, population sizes. In other words, the models we used to obtain the analytical results are deterministic, and this can potentially affect the generality of previous results. However, evolution is certainly not a completely deterministic process and stochasticity has been shown to play a crucial role in evolutionary dynamics, leading to the formulation of many important stochastic models like the genetic drift model citep{kimura1968evolutionary} or the probabilistic version of the Price’s equation citep{rice2008stochastic}.Here we tested for the effects of stochasticity on the predictions of our deterministic derivations, to assess their generality and soundness with respect to assumption violation. We built and put to work a simple heuristic individual-based computational model, comparing adaptive evolution in populations of different sizes and fixed or changing parameters (variables) such as the mutation rate, phenotypic robustness and the selection coefficient. We examined how phenotypic robustness affects adaptive evolution under different population sizes, mutation rates, and different selection regimes, and then we examine how phenotypic robustness evolve under common scenarios in constant and variable environments. section{The Model}This model description follows the Overview, Design concepts and Details protocol for describing individual-and agent-based models citep{grimm2005pattern, grimm2006standard, gilbert2008agent}. The model is implemented in NETLOGO v. 5.0.3, citep{tisue2004netlogo}, and the code is available in the electronic supplementary material (Model.txt). The model is designed to simulate the effects of different levels of phenotypic robustness $
ho$, ($0< ho<1$) on adaptive dynamics. As proposed by Rigato and Fusco (in prep.) we adopted a narrow, quantitative definition of phenotypic robustness, that is the probability that, across one replication, mutation of a given genotype g takes to a genotype $g'$ that exhibits the same phenotype of $g$, in other words, the probability that a mutation has no phenotypic effect.The model consists of a population of entities (or individuals), each assigned with initial fitness and robustness values. Each entity produces offspring according to its fitness value. Each single offspring can exhibit a mutated genotype (not modeled) with probability , and, conditional on the mutated genotype, a mutated phenotype with probability $(1- ho)$. Magnitude of the phenotypic mutation can be fixed or modeled with a random variable, depending on the simulation. There are no fitness costs directly associated to phenotypic robustness.subsection{Purpose}The main purpose of the model is to explore the consequences of phenotypic robustness in adaptive evolution, testing the effects of stochasticity with respect to the deterministic predictions of Rigato and Fusco (in prep.). This is done by simulating finite population persistence and phenotypic evolution under constant or changing environmental conditions (according to the simulation). subsection{Entities, state variables and timing.}Entities of the model are asexual individuals. Each individual has a genotype (not modeled) and a phenotype that includes its absolute fitness $W_i$, and phenotypic robustness $ ho_i$. Different genotypes can map on the same phenotype. Depending on the simulation, state variables are: population average relative fitness ($w_{mean}$), relative fitness variance ($var(w_i)$) and average robustness ($ ho_{mean}$). Relative fitness is computed with respect to a hypothetical maximum fitness value attainable in a given environment ($W_{opt}$), thus, $w_i= W_i / W_{opt}$. The model is time-discrete, one time step corresponding to one generation, and generations do not overlap. subsection{Process overview.}See a schematic diagram in Figure 1. Each cycle (generation) start with a parent population of $N$ individuals. Parents reproduce according to their fitness and die. Offspring initially inherit their parent's genotype and phenotype, but immediately the genotype can mutate with probability $mu$, and this can have a phenotypic effect with probability $(1- ho)$. If the resulting offspring population is larger than $N$, this is reduced to size $N$ through random elimination of the exceeding entities. These are the parents of the succeeding generation. In simulation with changing environment, fitness values are updated before parents reproduce.egin{figure}Hincludegraphicswidth=1linewidth{schema.pdf}centeringcaption{Schematic diagram of the process overview (see text).}label{placeholder}end{figure}subsection{Design concepts.}Evolution (changes in population mean and variance values of phenotypes) and other population dynamics (e.g., stability, bottlenecks and extinction) emerge from the combined effects of heredity (mutation and phenotypic robustness), natural selection (differential fecundity fitness of individuals) and demographic (population size) processes. Stochasticity has effects on i) genotype and phenotype mutation (with their magnitude), ii) survival and reproduction and iii) environmental changes.subsection{Initialization}Initialized parameters and their initial values depend on the simulation (see below).subsection{Input}The model does not have any external input; parameters are updated according to internal model rules.subsection{Sub-model "environmental-change".}Environmental change is simulated as a periodic indiscriminate negative effect on entity's fitness. The magnitude of the negative effect is described as a fixed negative factor ($-X$). To simulate the negative environmental effect, the negative factor ($-X$), is added to the individual absolute fitness values ($W_i$) before reproduction. Changes happens with a period $T$, in generation time unit ($T_{min}=1 generation$).subsection{Sub-model "reproduction".}In each generation, each individual i produces $W_i$ offspring and die.subsection{Sub-model "mutation".}The genotype can mutate with probability $mu$ , and if a genetic mutation occurs, the phenotype can mutate with probability $(1- ho)$. Depending on the simulation, this sub-model can be activated for fitness only or phenotypic robustness as well. In the latter case they operate independently. If the phenotype mutates, the new phenotype value is updated according to the following rules:\For fitness values: \1) in simulations with fixed selection coefficients, $s$, the new mutated fitness is set to $W_{updated} = W_{mean} (1+s)$. s can be positive with probability $p = 0.17$ citep{allen2000adaptation}, and negative with probability $(1- ho)$.\2) in simulations with variable selection coefficients, $s$, the new mutated fitness value is set to a value ranging from 0 to the maximum attainable absolute fitness $W_{max}$, with equal probability. Accordingly, fitness mutant gains decrease as the population approach to the optimum, as implied by the Fisher's geometric model.\For robustness values:\1) the new mutated robustness value is set to a value ranging from 0 to 1 with equal probability.subsection{Sub-model "maximum population size".}In each generation, if offspring population is larger than $N$, the population is reduced to size $N$ through random elimination of the exceeding entities. The surviving entities are going to be the parents of the succeeding generation.section{Simulations}The deterministic model proposed by Rigato and Fusco (in prep.) predicted a minimum level of phenotypic robustness for adaptation to occur, i.e. for the mean population fitness to increase. This minimum level depends on the magnitude of the selection coefficient and on the genotypic mutation probability according to the following equation:\$ ho_c=((1+s)mu-s)/((1+s)mu)$\where $ ho_c$ is the critical amount of robustness required, $mu$ the genotype mutation probability, and $s$ the selection coefficient of a given mutant phenotype. To test the consistency of this prediction we run four simulations.Simulations went on for 100-500 generations. Each simulation included several runs, each characterized by different initialization parameters (equal for all individuals or not), and a number of replicas for each run.subsection{Simulation 1.} It tested the existence of a $ ho_c$, by fixing for all runs the initialization parameters of the selection coefficient and the genotype mutation probability. Different runs were initialized with different fixed parameters such as the population size $N$ and phenotypic robustness $ ho$. Different combinations of $N$ and $ ho$ constitute the separate runs, each replicated three times. Individual fitness is the only changing internal variable, and the average population relative fitness is the only state variable. This simulation allows to verify the existence of a maximum robustness level under which the mean population fitness (state variable) does not increase significantly during time (generations).subsection{Simulation 2.} It tested how stochasticity can affect the structure of the relation between $ ho_c$, genotypic mutation probability  and the selection coefficient s, predicted by Rigato and Fusco (in prep.), in populations of nearly constant finite size $N$. We initialized all runs with the same fixed population size of 500, but different fixed values of $s$,  and $ ho$, each run in three replicas. As a result of the simulation, for each combination of s and , an observed $ ho_c$ was selected as the minimum level of $ ho$ under which adaptation did not occur with that particular combination of $s$ and .subsection{Simulation 3.} This simulation allowed to test the existence of a critical robustness level under which the mean population relative fitness does not increase during time in a fixed environment, with constant population size, mutation rate but variable selection coefficients among individuals, in a more realistic scenario that considers a population adapting through phenotypic mutations of declining magnitude effect in approaching the optimum and adapting with different (not fixed) selection coefficients. Simulation 3 tested the existence of a $ ho_c$, by fixing for all runs the initialization parameter of genotype mutation probability. Among runs, different fixed initialization parameters are population size $N$ and phenotypic robustness $ ho$. Different combinations of $N$ and $ ho$ constitute the separate runs, each replicated ten times. Individual fitness is the only changing internal variable while the distribution of the selection coefficients is an emergent property of the system, namely a new state variable like the average population relative fitness. The latter derives from the random assignment of fitness to newly mutated phenotypes (see sub-model "mutation"). subsection{Simulation 4.} To assess the evolvability of phenotypic robustness under different conditions, we simulated the fixed initialization conditions of Simulation 3, but allowing phenotypic robustness to evolve during adaptation. This means that, like fitness, phenotypic robustness is a changing internal variable. Population average phenotypic robustness is the new state variable as the population average relative fitness. We simulated a constant and a variable environment (see "environment" sub-model). Evo