Evolution of phenotypic fluctuation under host-parasite interactions.

Robustness and plasticity are essential features that allow biological systems to cope with complex and variable environments. In a constant environment, robustness, i.e., insensitivity of phenotypes, is expected to increase, whereas plasticity, i.e., the changeability of phenotypes, tends to diminish. Under a variable environment, existence of plasticity will be relevant. The robustness and plasticity, on the other hand, are related to phenotypic variances. As phenotypic variances decrease with the increase in robustness to perturbations, they are expected to decrease through the evolution. However, in nature, phenotypic fluctuation is preserved to a certain degree. One possible cause for this is environmental variation, where one of the most important "environmental" factors will be inter-species interactions. As a first step toward investigating phenotypic fluctuation in response to an inter-species interaction, we present the study of a simple two-species system that comprises hosts and parasites. Hosts are expected to evolve to achieve a phenotype that optimizes fitness. Then, the robustness of the corresponding phenotype will be increased by reducing phenotypic fluctuations. Conversely, plasticity tends to evolve to avoid certain phenotypes that are attacked by parasites. By using a dynamic model of gene expression for the host, we investigate the evolution of the genotype-phenotype map and of phenotypic variances. If the host-parasite interaction is weak, the fittest phenotype of the host evolves to reduce phenotypic variances. In contrast, if there exists a sufficient degree of interaction, the phenotypic variances of hosts increase to escape parasite attacks. For the latter case, we found two strategies: if the noise in the stochastic gene expression is below a certain threshold, the phenotypic variance increases via genetic diversification, whereas above this threshold, it is increased mediated by noise-induced phenotypic fluctuation. We examine how the increase in the phenotypic variances caused by parasite interactions influences the growth rate of a single host, and observed a trade-off between the two. Our results help elucidate the roles played by noise and genetic mutations in the evolution of phenotypic fluctuation and robustness in response to host-parasite interactions.

Introduction

Robustness and plasticity are two important properties of biological systems. To maintain function and high fitness against internal noise, environmental variation, and genetic changes, the fitted state must be robust to such disturbance, whereas the phenotype needs to possess plasticity to adapt to environmental variation. Indeed, the evolution of robustness and plasticity has been investigated extensively both theoretically and experimentally [1-5].

In this context, phenotypes are shaped by dynamic processes involving several genetically determined variables, as well as external and internal noise. Hence phenotypes can vary not only as a consequence of genetic variation, but even between isogenic individuals [6-10]. On the other hand, the dynamical process varies by the genetic changes, so that there also exists phenotypic fluctuations by the genetic variation [11-13]. The sum of both sources of phenotypic variation defines the genotype–phenotype mapping [14-16]. The evolution of genotype–phenotype mapping, thus, is essential to understand the evolution of robustness and plasticity, and has been explored extensively both in numerical [4, 17] and laboratory experiments [18, 19].

The magnitude of these phenotypic fluctuations generally decreases as a result of evolution under a given environment, because, with this decrease, the system deviates less from the fitted state [20-23]. In fact, experiments and simulations of adaptive evolution under fixed conditions have shown that fluctuations decrease over the course of evolution [24], and evolvability, i.e., the rate of increase in the fitness or the phenotypic change per generation, declines. Accordingly, as evolution progresses, the robustness of the phenotype increases, while the phenotypic fluctuations decrease. The decreased phenotypic fluctuation will then be accompanied by the decreased adaptability to environmental variation, as discussed by generalizing the fluctuation-response relationship in statistical physics [17, 25, 26].

In nature, however, phenotypic fluctuations and evolutionary potentiality persist. Even after evolution in a constant environment, the phenotype in question is not necessarily concentrated on its optimal value, instead its variance often keeps a rather large value. For instance, even for isogenic cells (clones), fluctuations in the concentration of each protein remain sufficiently high. This leads to the following questions: Why are phenotypic fluctuations not reduced further, which, in principle, would be possible by evolving appropriately negative feedback processes for stabilization? How are phenotypic fluctuations or plasticity sustained?

One possible cause for the preservation of plasticity or fluctuation is environmental fluctuation [27-29]. In natural populations, environment is not fixed; instead, it is variable. If phenotypic fluctuation is reduced excessively, the plasticity required to adapt to new environmental conditions would be lost. The plasticity of a biological system, would then decrease its capacity to cope with environmental changes.

For every species, interactions with other species are one of the most important “environmental” factors that affect its existence. Even if a certain “external” environmental condition itself is fixed, the types and populations of other species may change because of species–species interactions [30]. For example, in the interaction between hosts (prey) and parasites (predators), the former may change their phenotype to escape the attack by latter’s, which, in turn, will change the distribution of the latter phenotypes to continue attacking the former [31, 32]. Thus, each species may retain evolutionary plasticity to cope with dynamic changes in other species. Furthermore, phenotypic variation in one species may influence the other species, which can result in a dynamic change in inter-species interactions. If sufficient species–species interaction exists, the dynamic variation of phenotypes may be mutually sustained across multiple species.

As a first step to investigate the phenotypic fluctuations originating from such species–species interactions, we present a study of a simple two-species system that comprises hosts and parasites. Within this system, distinct parasite types that can attack hosts by possessing a specific phenotype exist; the host phenotypes are determined by genotype–phenotype mapping that incorporates possible fluctuations. Hosts are expected to evolve to achieve a phenotype that optimizes fitness and increases the robustness of the successful phenotype by reducing phenotypic fluctuations. in addition, they are expected to evolve plasticity, i.e., phenotypic adaptability, to cope with parasite attacks. We use a host gene-expression dynamics model to study the evolution of genotype–phenotype mapping and the resultant phenotypic fluctuations, and explore the evolution of phenotypic robustness and plasticity in response to interactions with parasites.

Here, we address the following questions: first, how are robustness and plasticity of phenotypes, which tend to exhibit opposite trends, maintained through evolution caused by interactions between the host and the parasite species? Second, are phenotypic fluctuations sustained to cope with parasite attacks? As mentioned before, there are two sources of phenotypic fluctuations: genetic variation and stochasticity in genotype–phenotype mapping (i.e., noise in the dynamics that shape the phenotype). Then, we address a third question: which of these sources of phenotypic fluctuation is dominant in sustaining the level of phenotypic fluctuation? Finally, we discuss the dependence of the evolution of phenotypic fluctuation and plasticity on the strength of the host–parasite interaction and the sensitivity of phenotype dynamics to noise.

Models

Gene regulatory networks

We employed a simple model for gene expression dynamics based on a sigmoidal function to study the evolution of the genotype–phenotype map. There are dynamics of protein expression level for each protein (gene) i. This model comprises a network consisting of genes that mutually activate or inhibit each other’s expression. In this model, there are M genes whose gene expression level, (i = 1, 2, …, M), corresponding to the host genotype l at time t is described by where J denotes the matrix representing the influence of gene j on the expression of gene i. When j activates (represses) the expression i, it takes the form J = 1(−1), while J = 0 if j does not influence i. The sigmoidal function f(z) is given by with β = 25; this implies that f(z) closely resembles a step function. Furthermore, θ represents the threshold of the input required for the expression of gene i. Therefore, whether the gene(protein) is expressed (f(z) ≈ 1) or not (f(z) ≈ 0) depends on whether the sum of inputs from the expression of other genes(proteins) is larger than θ. In Eq (1), the value θ is selected randomly from a uniform distribution between 0.05 and 0.3. The term representing the interaction with other genes is scaled by , where α denotes the path density, that is, the fraction of nonzero values of J. We adopted this scaling so that x is comparable with the order of θ regardless of total gene number M(= 64). (Because the connection is sparse, many of J = 0, and the average number of inputs for gene is αM. We normalized the interaction with other genes by the square root of the path number αM so that the variance of J is independent of αM.) Moreover, there are lout(= 8) “output genes” i = 1, 2, …, lout, which determine the fitness as defined below. To eliminate the potential for the output genes to influence others, the summation is taken only for j > l. The term η((t) denotes a Gaussian white noise originating from molecular fluctuations in chemical reactions.

Here, the noise is applied to the expression of each gene, and it is not correlated across genes. The value of σ represents the strength of this noise (In the case where the gene expression of x is negative, we set x to zero.) For constant f(z), Eq (1) defines an Ornstein-Uhlenbeck process, and to give the equilibrium variance of σ2/(2γ). Initially, all gene expression levels are set smaller than the threshold θ. Furthermore, ϵ denotes the spontaneous expression level, which is smaller than θ. I denotes the external input for linp(= 8) “input genes,” with I = 1 for j = M − linp + 1, M − linp + 2, …, M, and I = 0 for j ≤ M − linp. The input genes are always/constitutively expressed. In order for x to be expressed, x must receive external or internal inputs from other x values to advance beyond θ.

Initially, we select J at random to obtain J = 1(−1) with a probability of 0.15 so that α = 0.3. We set the initial number of host cells as N and total number is constant at N = 300 throughout all simulations. Here, we set M = 64, lout = 8, and linp = 8. The model parameters used in this study are set as in Table 1.

Table 1

Model parameters and variable definitions.

Symbol	Description	value
M	Total number of host genes	64
N	Total number of host cells	300
l inp	Number of input genes	8
l out	Number of output genes	8
I j	Input strength	1
γ	Time constant	0.1
α	Average path density	0.3
ϵ	Spontaneous expression level	0.01
β	Sensitivity in the expression	25
μ 0	Basal growth rate	0.0002
m	Mutation rate	0.03
l p	Total number of parasite genes	3
N p	Number of parasite types	8
T p	Parasite generation time	100
D p	Mutation rate of parasite	0.01

Raw fitness and reproduction

The raw fitness of the host alone, i.e., in the absence of parasites, is determined by the fraction of output genes that are expressed as where x denote the expression of the output genes and μ0 denotes the basal fitness. Therefore, the state with all output genes expressed (x = 1, for j = 1, 2, …, lout) is an optimal phenotype for achieving maximal growth.

We introduce parasite attacks into the model. There are parasite genotypes that are coded as with i = 0 or 1. For example, for l = 3, there are 23 = 8 genotypes coded by 000, 001, …, 111. Each parasite attacks hosts that possess the corresponding target gene expression, , that is, each parasite attacks the host whose gene expression pattern, x (1 ≤ j ≤ l), matches (). Here, we term the l genes as “target genes.” The target condition for the attack of each parasite is given by s(x) = 1 for x > 0.5 and s(x) = 0 otherwise (as the threshold function is close to a step function, x takes either ∼0 or ∼1). For instance, a host whose expression pattern (s(x1), s(x2), s(x3)) ≈ (0, 1, 0) is attacked by the parasite (i1, i2, i3) = (0, 1, 0). Let m () be the parasite type whose matches with host for i = 1, …, l. The growth rate of each host genotype l decreases owing to the effect of the parasite m as here P denotes the population density of the m-th parasite that matches with the host gene expression for (i = 1, …, l) as mentioned and c denotes the strength coefficient of the attack by the parasite. The evolution of hosts is modelled by an individual-based approach with a fixed population size akin of the Moran process. Each host has its own genotype (gene regulatory networks) and the growth rate determined by the gene expression levels calculated by numerical simulations in Eq (1). The volume of the host cells grows according to

When the volume exceeds a threshold value of 2 (v(0) set to 1), the host cell is divided into two parts. The total number of hosts N is set to be constant over time, and thus, whenever a cell divides, another cell is randomly selected and removed from the original population.

Here, after division, the initial expression of x is reset to take a random value smaller than the threshold θ. The division time of the host is determined by the individual growth rate that ranges between 500 and 1,000 time units. At this time scale, the host gene expression reaches a steady state. The gene expression of hosts reaches a steady state after around 200 time units. Mutation is introduced in the division process and added to the network J with a low mutation rate m = 0.03. An (i, j) path in the network matrix J is selected at random and its value is changed to one of the others to preserve the total path number. For example, if J = 1(−1) and J = 0, then they are changed to J = 1(−1) and J = 0. Here, we randomly select a gene k ≠ i, i ≠ j, k ≠ j. Then, the path j → i is replaced by j → k to represent the mutation of a path in gene regulatory networks. We prohibited direct connections between the “input” and “output” genes. Initially, J is selected randomly with ±1 with the rate α = 0.3, and α remains constant throughout the simulation. (However, the result is independent of this choice.)

It should be noted that once cells settle down to different fixed-point attractors, no more switching occurs in time. The on/off difference in the gene by cells occurs due to the developmental noise from the initial condition, which subsequently results in phenotypic variance in the population.

Parasite population dynamics

Unlike host phenotype (gene expression levels), which is partially stochastic, the parasites’ phenotype is determined uniquely by its genotype. Hence parasites of a given genotype share a common phenotype and parasite population dynamics can be modelled by deterministic equations (instead of individual-based simulations).

The population of each parasite P increases in proportion to the host population with the corresponding phenotype, where the host density H to be attacked by the parasite strain i is determined by the host population that satisfies s(x) = i. Further, there is a mutation from other parasite types . Because parasite types are represented by a binary string of length l, they are represented in a l-dimensional hypercube space. Each distinct genotype is represented by a vertex on the cube. Mutations in the parasite genotypes occur by flipping each gene i = 0 ↔ i = 1 at the mutation rate D. This is represented by diffusion in the hypercube (e.g., 010 ↔ 000) determined by a diffusion constant D. Therefore, the population dynamics is expressed as where κ takes unity only if the binary string of j can change to that of i via a single mutation, otherwise it equals zero. We assume competition within all the parasite genotypes, such that the total number of populations is assumed to be fixed, to avoid the divergence of parasite populations. Accordingly, we normalize the density P such that . Thus, the reduction term is introduced ensuring that and remain satisfied. We computed the population dynamics of parasites with a given generation time T (s). We set the initial density P to take the same value for each type, i.e., (The result below, however, is independent of this initial distribution). The generation time of parasites is shorter than that of hosts. Therefore, the fittest type would quickly become dominant.

Phenotypic variances attributed to noise and mutation

The model includes a noise component that allows gene expression level to vary even among individuals sharing the same gene regulatory network. We define the phenotypic variance as the variance of the phenotype over the entire host population with different genotypes. Here, the variation of to give arises both from noise and from network mutations. In addition, we compute the isogenic phenotypic variance of the gene expression level within each host genotype l. Because of the noise term ση in Eq (1), the expression level x for the same host genotype (network) varies between individuals and this variance is positive. We define two phenotypic variances to distinguish between these sources: V(i) and V(i)[23, 25]. The variance is defined as the variance of gene expression levels in an isogenic population, i.e., the phenotypic variance of x caused by noise in the gene expression dynamics within the clonal population of hosts, i.e., sharing the common J. In contrast, is defined as the phenotypic variance due to genetic mutation, as computed by the variance of within a population of initial genotype i in which each individual receives one random mutation, where is the average over the clonal distribution of host l. Details describing the calculation of variances V and V are provided in the Methods section. Here, as we are interested in the variance of the expression of output genes, we term the average variances of the output genes (i = 1, 2, …, 8) as V, V, and V, thereby omitting the l notation. In addition, we define the average variances of the target genes corresponding to parasite attacks as and .

Note that developmental noise (the ση term in Eq (1)) has no significant effect on host fitness once gene expression dynamics have reached a stable state with all x close to 0 or 1. This is because host growth depends only on mean gene expression levels and parasite attack depends on rounded values (whether x is larger or smaller than 0.5). However, noise is important, because when amplified in the GRN, it may cause the initial dynamics to converge to a different stable state (i.e., noise-induced switches in gene expression). Regarding phenotypic variances, we discuss two types of robustness: mutational robustness and developmental robustness. Mutational robustness describes the “ability” of a GNR to prevent changes in gene-expression patterns due to mutations in the network structure, and developmental robustness describes the ability to prevent changes in expression patterns due to noise-induced fluctuations in gene expression level during the initial network dynamics.

Results

Parasite interaction accelerates host diversity

We started simulations with N host cells of an initial gene-regulatory network created by creating random paths with probability α = 0.3. We took 50 replicated populations to obtain the statistics. The initial condition for each cell is given by a state in which only the input genes are expressed: all other genes are initialized with expression levels drawn uniformly between 0 and 0.05. This same initial condition was also applied to the daughter cells after each cell division. In the following results, gene regulatory networks (GRNs) that went extinct as a result of parasite interactions are excluded immediately from the sample at the initial time t = 0.

Fig 1 shows the population dynamics and fitness values of the host in three parasite environments (c = 0, 2, and 4) in the absence of noise (σ = 0). In the absence of parasite interaction, i.e., at c = 0, the host population is dominated by GRNs that generate a phenotype expressing all target genes, such as type “111” (black line), which dominate the host population (see Fig 1A). Then, if sufficient interactions with the parasites occur, the host population evolves into multiple groups with different phenotypes (Fig 1B and 1C). This is explained by the selection pressure to avoid parasite attacks. When the host population is concentrated on the fittest type “11111111,” which includes the target gene expressions “111,” the parasite population is also concentrated on the corresponding type matching the target “111,” which results in suppression of the population of the fittest host types by the interaction with the parasite. Consequently, other host types with less preferred/less fit target patterns have more chances to survive. In short, parasites induce negative frequency-dependent selection pressure, and the effective fitness (Eq (4)) varies over time because of the change in the distribution of parasite genotypes. Thus, multiple phenotypes and genotypes can coexist within a host population. Note that in the absence of noise, the expression dynamics of isogenic individuals always converge to the same pattern.

Fig 1

Dynamics of host phenotypes under interaction with parasites.

The populations of hosts H corresponding to the parasite type () are plotted against the generation time. The number of parasite types is N = 23 = 8, which corresponds to the on/off expression of the three target genes of the host. The population of each host type i for H represents {s(x1), s(x2), s(x3)} = “000,” “001,”…, and “111” is indicated by different colors. The interaction strength is c = 0 (A), 2 (B), and 4 (C). The noise level is σ = 0.0. In the absence of parasites (c = 0), the population is dominated by the “111” type in which all output gene expression levels are turned on, whereas other phenotypic types with low growth rates (“110,” “101,” and “011”) also coexist at c = 2 and 4. In addition, N = 300, M = 64, linp = 8, lout = 8, and l = 3 (see Table 1).

Next, we computed the phenotypic variance V that is the variance of output gene expression levels over the host population with polymorphic genotypes (i.e., network J) and plotted it as a function of the interaction strength c to determine the extent to which host–parasite interactions enhance phenotypic diversity (Fig 2). Independent of the level of developmental noise, Fig 2 shows a marked increase in the phenotypic variance at a threshold interaction strength c = c ≈ 1. This diversification indicates that multiple phenotypes coexist within the population, as shown in Fig 1B and 1C. In contrast, below the transition point c, the variance decreases and the host population is concentrated on the fittest phenotype, as shown in Fig 1A.

Fig 2

Total phenotypic variance V of the output genes against the parasite interaction strength.

The variance V of the output genes plotted against the parasite interaction strength c (bars represent standard deviation). V is computed from the expression levels of the output genes over the host population (N = 300 individuals) and plotted as the average V over 2500–3000 generations. Each color represents a different noise strength: σ = 0.03 (blue), 0.02 (orange), 0.01 (black), and 0.001 (red).

This critical interaction strength c ≈ 1 follows from the fitness function in Eq (4). The host growth rate is maximized when the expression level of all output genes is close to x = 1. If the parasite population is concentrated on the type that attacks this fittest phenotype, the population growth is reduced to cP ∼ c, with P as the population fraction of the parasite type m(= 1, …, l). To allow the phenotype to escape from the parasite attack, one target gene should be switched off, which reduces the raw growth rate by one unit. Then, the raw fitness loss needs to be smaller than the gain by the escape c. Since, P ≤ 1, the total fitness gain is positive if c ≳ 1. Indeed, Fig 2 shows that, for c > c, the variance increases substantially. When the interaction is strong (e.g., c = 4), the number of coexisting host genotypes approaches the maximal number 2 = 8, which corresponds to the total number of parasite types. (see Fig 1C).

In our model, phenotypic variance can have two different origins: genetic variation and noise-induced phenotypic variation. These correspond to two strategies to sustain phenotypic diversity. We investigated how the choice of each of the two adaptive strategies depends on c and σ, by computing the noise-induced isogenic phenotypic variation, V, and the variance induced by mutations, V.

Evolution of robustness and phenotypic fluctuation

Fig 3 shows evolutionary time courses for the variance components V and V (for both output and target genes) for various combinations of c and σ. The temporal averages of V and V after evolution of the system reaches a stable state are shown in Figs 4 and 5, respectively. These results are independent of the choice of the initial network. Initial values of V are high because random networks are not canalized and are highly sensitive to mutation. In addition to the critical interaction strength c, we can see (e.g., from Figs 4 and 5 against c and σ) the existence of a critical noise level σ ≈ 0.02, above which V increases markedly. Based on these two thresholds, behaviors of host gene expression states can be classified into four regimes:

Fig 3

Evolutionary time course of the average phenotypic variations V and V.

The time course of the average phenotypic variations V and V (i.e., and ). The plot is from the first generation rather than 0th generation. Note that host types that went extinct during the simulation have been excluded. This explains the different initial variance levels between parameter combinations. (A) Average V and V over all output genes are plotted against the generation for different values of the noise level σ and interaction strength c, as indicated by different colors. Red: Example of evolution with no host–parasite interactions. Hosts are occupied by individuals with the fittest phenotype, thereby losing phenotypic diversity (σ = 0.02). Both the isogenic phenotype variation V and the genetic variance V decrease. Orange: Case with c = 2 and σ = 0.01, which shows that V decreases, whereas V maintains high values, implying the evolution of genotypic diversity. Green: c = 2 and σ = 0.03. Both V and V increase. Blue: c = 4 and σ = 0.05. V decreases slightly, whereas V maintains low values. (B) and for the target genes. The line colors represent the equivalent conditions as in (A).

Fig 4

Dependencies of V and V on the noise level σ and interaction strength c.

V and V for all output genes upon noise level σ (horizontal axis) and interaction strength c (vertical axis). The variance values for each parameter are displayed using color maps. The variance is computed by taking the average of 2500–3000 generations. The color maps indicate that the host–parasite interaction enhances V and V. Parameter values of σ and c for each square are shown at their left and bottom edge (e.g., the bottom-leftmost square is the results for σ = 0.001 and c = 0) Each parameter regime is categorized into one of four regimes based on the values of V and V (see the text for details).

Fig 5

Transition of V and V with noise strength.

Dependence of the average variances V (red) and V (blue) for output genes (solid lines) and (red) and (blue) for target genes (dotted lines) upon the noise level σ for a strong host–parasite interaction (c = 3). Bars represent the σ standard deviation. exceeds at σ ≈ 0.02. Based on the data in Fig 4, the phenotypic diversification occurs for all parameter values here, but the main origin of the phenotypic variances changes from the mutation to the noise at σ ≈ σ.

When the noise is below the threshold σ = 0.02 and the interaction is weak (c < c = 1), V is low, whereas V is maintained at a high level (Fig 4). The phenotype varies considerably by mutation, which indicates that mutational robustness does not evolve. (A large V implies that the phenotypes are not particularly robust against mutations.) Host phenotypes (i.e., gene expression) fluctuate around the fittest type, as shown in Fig 1A, where most of the target genes are expressed, whereas the variances are sustained at moderate values because of genetic variation in the network structure.

The evolution of developmental robustness against the noise in gene expression dynamics occurs when the noise is above a threshold of σ = 0.02 and the interaction is weak (c < c = 1). In this region, the strength of the interaction is weaker than at the transition point c = 1, which means that the attack by the parasite is not sufficient to cause phenotypic diversification (Fig 2). Fig 4 shows that the values of both V and V are low relative to the case with σ < σ and c < c. The phenotype is not significantly changed by genetic variation. V and V both decrease while satisfying the inequality V > V (see Fig A in S1 Appendix), and this leads to the evolution of robustness to both noise and mutation. This supports previous studies [17, 23, 25]. The host population loses both genetic and phenotypic inhomogeneity, which indicates that the output genes are all switched on and exhibit minimal variation with respect to the fittest type.

The evolution of genetic diversification occurs when the noise level is below the threshold σ < σ and the interaction strength exceeds the transition point c. V remains at high levels, whereas V remains small, as shown in Fig 5. Phenotypic diversity is generated by genetic diversity. For each given genotype, the generated phenotype is unique and exhibits almost no variation. Multiple groups with different genotypes coexist because of interactions with diversified parasites. The stronger the interaction, the more phenotypes arise, which leads to a larger V and higher genetic diversity in the population. For moderate c > c, the genotypes with the fittest type “111,” coexists with genotypes having one target gene switched off (“110,” “101,” “011,” etc.), whereas for c = 4, even more host types appear, leading to coexistence, including expression patterns where only one target gene switched on (“100,” “010,” “001,” etc.)

As exemplified by the case of c = 3 and σ = 0.02, the phenotypic variation of the output gene expression (V) increases in the initial stage of evolution (up to 500 generations; see Fig 3). Both V and V maintain high values, while satisfying V > V. Host phenotypes are diversified without resorting to genetic diversification. Here, isogenic hosts can have more than one phenotype. Gene expression levels are diversified by noise, which implies the phenotypic variation of the isogenic population. Variances in expression levels are highest for target genes, whereas those of the other output genes are maintained at high levels (see Fig B in S1 Appendix). Accordingly, the fitness is reduced.

Transition between Regime III and IV against noise strength

We now explore the transition between the last two regimes. Fig 5 shows the dependence of the variances V and V of the output genes on the noise level σ for c = 3. Below the threshold noise level σ, V exceeds V, indicating that there is little robustness against mutation; the diversity of the target gene expression levels is based primarily on genetic variation. In contrast, when the noise level exceeds σ, V > V and the diversification of gene expression patterns is achieved over isogenic phenotypic variation, i.e., due to noise.

Note that, in Fig 5 V decreases slightly with σ once σ exceeds the critical value σ. This reflects previous results for c = 0 (i.e., in the absence of parasites), which showed that mutational robustness decreases as a correlated response to the evolution of developmental robustness against increasing noise levels (Fig A in S1 Appendix). This is because, for a larger noise strength, the acquisition of developmental robustness to noise is needed which also leads to an increase in robustness to genetic mutations. These factors have been investigated in detail in previous studies [23]. When c is nonzero, V is larger than the case for c = 0, but its decrease against σ for the case with c = 0 still remains, so that V shows slight decrease against σ.

Evolution of the noise-induced on/off switches

As mentioned above, once expression levels of host target genes have settled down close to one of the fixed points x = 0 or 1, continued variation due to the noise term in Eq (1) has no effect on the interaction with the parasite. The noise is important, however, for the initial network dynamics that determine which of these fixed points is approached in the first place. In order to separate the contributions of noise-induced on/off-switches from those of the noise itself, we define a variance which only reflects the former. More precisely, is the variance of the on/off variable z = tanhβ(x − 0.5) with β = 100 and defined as the variance of z, computed in the same manner as V(i).

Fig 6 shows the dependence of the variances (output genes) and (target genes) on the noise level σ for c = 3. The key difference to Fig 5 is that for σ > 0.02, remains essentially constant for both target and output genes, suggesting that the variance due to noise-induced switches cannot be increased arbitrarily, and that the continued increase of V seen in Fig 5 is merely a reflection of the increase in white noise. The reason that does not increase much beyond σ > 0.02 will be because the level of phenotypic variation reached by the noise, in conjunction with genetic variation, is sufficient to produce the maximal phenotypic variation. Indeed, V with c = 3 in Fig 2 is the same regardless of the noise strength. Evolutionary time courses for both V and are shown in Fig E in S1 Appendix. For target genes, both measures increase from their initial value and then stay at a high level, with no evidence for evolution of developmental robustness. For non-target output genes, in contrast, V and decrease over time relative to the values in the initial random networks. This shows that the decrease in the variance has not been observed in the target gene implying that the developmental robustness has not evolved. Also, as discussed in the previous section, V decreases with σ, but not more than the case with c = 0, suggesting that there exists some coupling between the and . On the other hand, the variance of the non-target output gene () decreases compared to the initial network, and the developmental robustness evolves in correlation with mutational robustness (V). At around the noise level (σ = σ ∼ 0.02), the variances and show remarkable increase. Albeit less significant as compared to , also increases. Although only the on/off of the target gene is relevant to protect against the parasite attack, the noise-induced switch is increased for other genes, due to the coupling in GRN. In the figure, is of the same level beyond σ = 0.02. The required switching frequency is of the same level, once at σ ∼ σ it occurs. (It increases with the interaction strength c, though).

Fig 6

Transition of with noise strength.

Dependence of the average variances for output genes (solid lines) and for target genes (dotted lines) due to noise-induced switches plotted against the noise level σ for a strong host–parasite interaction (c = 3). We used the function z = tanhβ(x − 0.5) with β = 100 and defined as the variance of z. Bars represent standard deviation.

Dependencies of the variances on genetic change and noise

In the absence of parasites, the variances V and V decrease as the host phenotype adapts to the environment. This indicates the evolution of robustness to mutation and noise (regime II). Under the pressure of the host–parasite interactions, the gene expression dynamics evolved to diversify the phenotypes in the manner described above. Here, we discuss the increase in V and V through evolution.

First, we plot the evolutionary time course of V and V starting from the initial random networks J setting with parameters c = 3, σ = 0.03 in Fig 7A. The variances and of target genes both increase rapidly because of interactions with parasites. (GRNs that went extinct due to parasite infections of the 0th generation are excluded. Therefore, the variance V tends to be larger at c = 3 initially.) As shown in Fig 7A(b), the two variance terms exhibit considerably closer correlation between V and V for the expression of the target gene. Moreover, both variances maintain large values, and the host can escape parasite attacks by producing a different phenotype via mutation. In contrast, Fig 7A(a) shows that the variances of other output genes decreased and developmental robustness increased for those phenotypes that were not attacked by parasites. This decrease, however, is much smaller than that observed in the absence of the host–parasite interactions (see Fig F in S1 Appendix). The decrease implies that the expression of target and non-target output genes is partially decoupled, but the reduction of the decrease suggests some coupling. The evolution of robustness and phenotypic fluctuation occurs in the same GRN’s.

Fig 7

Time course of the variances (V, V) after switching interaction strength.

(A) The time course of the variances (V, V) for output genes (a) and the variance (, ) for target genes (b) plotted over generations for c = 3 when starting with initial random networks (circle; dotted lines). The variances are computed from the isogenic variance over 100 iterations. The time course over generations is plotted for the variances of gene expression. The plots cover 3500 generations. We set the noise level to σ = 0.03 > σ. The two variances of the non-target output genes decrease, while V < V is maintained throughout the evolutionary course (a). Conversely, under host–parasite interactions, and show a correlated increase (b). (a) shows that the variances of other output genes not attacked by parasites decreased. However, this decrease is much smaller than that corresponding to the absence of host–parasite interactions. (B) The time course of the variances (, ) for target genes (a) and variance (V, V) for output genes (b) over generations for σ = 0.03 when switching interactions c from 0 to 3 (asterisks; bold lines). For the first 2500 generations, the evolution of the host was modeled without parasites (c = 0). Up to this generation, the evolution of developmental and mutational robustness was complete V and V decreases. Then, we introduced the interaction with the parasites by changing c from 0 to 3. The variances were computed from the isogenic variance over 20 iterations. The time course of (V, V) and () over generations after this switch is plotted for the output and target genes. The color bar indicates the number of generations since the switch. The variances and V both increase because of interactions with the parasites, which results in the evolution of phenotypic variation. Although this increase is prominent for target genes (b), those for other output genes also show a slight increase up to 2000 generations (a). Variances and finally reach the same values as those in (A). Other model parameters: N = 300, M = 64, linp = 8, lout = 8, and l = 3.

Coping with parasites by phenotypic variation

Next, we investigated whether high variability in phenotype (i.e., phenotypic fluctuation) can evolve in response to the introduction of parasites, even after developmental robustness has evolved. To investigate this issue, we first set c = 0 and let the system evolve, and once the host had adapted to the environment and acquired robustness (i.e., after V and V had decreased from their high initial values), we switched on the interaction by increasing c from 0 to 3. The time course of the variance is shown in Fig 7B. After switching the interaction strength, the host growth rate decreases temporarily before increase in phenotypic variation evolves (see Fig G in S1 Appendix). Then, both the variances and increase. Thus, notwithstanding the prior evolution of developmental robustness, phenotypic fluctuations can evolve in response to parasite infection. The phenotypic fluctuations (V) in the host population show a notable increase within 2000 generations. In addition to the target genes, the fluctuations in the expression level of the output genes also increase. Hence, the output and target genes are not completely decoupled.

A comparison of Fig 7A and 7B shows that the final variances are nearly identical. This means that the variances in both target and output genes after evolution are independent of the initial conditions. To increase the phenotypic variance of the target genes, variances of other gene expressions need to be maintained to a certain degree because of gene interactions through the GRN. This trend of a weak increase in the phenotypic variance for all genes holds for c > c (see Figs C and D in S1 Appendix).

Trade-off between growth and tolerance

We examined how the increase in the phenotypic variances caused by parasite interactions influences the growth rate of a single host. Parasites concentrate their attacks on the hosts with the most frequent phenotype in the population. Fig 8 plots the relationship between V and the raw fitness, and it shows the trade-off between the two. Individuals exhibiting high plasticity tend to have lower growth rates because of the uncertainty in the expression of output genes; however, they are attacked less by the parasites. As the variance V is larger, the parasite attacks are reduced. Consequently, the host genotypes with larger V have a larger chance of survival, even if their growth rate as a single cell is lower. In contrast, in region III with large V, host populations are adapted to parasite infection through genetic variation. In this case, low phenotypic fluctuations and high raw fitness are compatible. Here, target genes and non-target output genes are well decoupled and raw fitness tends not to decrease even if the expression of target genes is variable (Fig 8 σ = 0.01 and σ = 0.001 case). Therefore, the fitness costs are higher for the variance by noise (V) than that by the genetic variation (V). The lack of fluctuations, however, makes the host population more susceptible to parasite infections when genetic variation is suddenly reduced because of extinction.

Fig 8

Trade-off between phenotypic variation (V) and fitness.

Trade-off between isogenic phenotypic variation (V) of the output genes and raw fitness . Each value is computed by the average from 2500–3000 generations. Each color represents a different interaction strength: c = 0 (black), c = 1 (green), c = 2 (red), c = 3 (blue), and c = 4 (yellow). Each marker shape represents a different noise strength: σ = 0.05(⋄), σ = 0.04(⋆), σ = 0.03(⊳), σ = 0.02(∘), σ = 0.01(△), and σ = 0.001(⊲).

Discussion

In this study, we investigated the evolution of phenotypic variances using host gene expression dynamics within a regulatory network in the presence of a host–parasite interaction. If the interaction is weak, the host with the fittest phenotype evolves to reduce phenotypic variances. In contrast, if the interaction is sufficiently strong, the phenotypic variance of the host increases as strains evolve adaptation to decrease attack by specific parasites. We identified two strategies—either to increase the noise-induced phenotypic variation (V) or to increase the genetic variance (V)—depending on the strength of the noise in the stochastic gene expression. If the noise strength is below the threshold, the diversification is primarily genetic in origin, whereas above the threshold noise-induced phenotypic variation dominates, leading to the increase in variances V and V. In the latter case, both variances increase in correlation, and thus, which enhances phenotypic variation and helps avoid parasite attacks.

In a fixed environment without inter-species interactions, both variances V (due to noise) and V (due to mutation) tend to decrease, which causes loss of phenotypic variation and evolvability. Under host–parasite interactions, the GRN evolves to increase these variances, even after robustness has evolved and the variances have diminished. We classified the diversification of phenotypes into four regimes based on the interaction strength c and noise level σ. The regimes were classified according to the respective degrees of phenotypic variance caused by phenotypic noise (V) and genetic variation (V).

It should be of interest to investigate how the structure of evolved GRN depends on each regime. For example, we examined the distribution of 3-genes motifs [33] for each regime. We detected the increase of motifs that activate target genes in all regimes; however, so far we could not detect a clear difference in the motif distribution among the four regimes yet. The characterization of network structure inherent to each regime remains as a future problem.

Whether biological populations deal with environmental changes by genetic variation or phenotypic variation has been discussed both theoretically and experimentally [3, 34–36]. Indeed, the relevance of species–species interactions to phenotypic variation and genetic diversification has been discussed extensively [37-39]. Therefore, it is interesting to examine how the two strategies for phenotypic diversification studied here are adopted therein.

Here, we observed a trade-off between phenotypic variation and raw fitness. Hosts that evolved the increase in phenotypic variation to deal with parasite interactions tended to exhibit a decrease in raw fitness (i.e., the growth rate). Interestingly, such a trade-off has been observed for predator-induced phenotypic plasticity [39-42] and in bacteria-phage experiments [43], where the bacteria gain resistance to the phage by reducing the competition for resources. Furthermore, it has been suggested that excessively strong phenotypic plasticity may reduce adaptability to climate change [44].

In contrast, host populations that exhibit phenotypic diversification by increasing genetic variation (V) to reduce the rate of infection do not show a significant decrease in fitness. However, in this case, their phenotypes are less robust against noise or mutation.

The generation of diverse phenotypes against uncertain environmental changes is known as a “bet-hedging strategy” [45-51]. The question of whether to cope with environmental changes by genetic evolution or by phenotypic plasticity has been discussed in relation to the speed of environmental change. The diversification strategies originating from either V or V, as discussed here will provide the basis for the increase in phenotypic variance, which is required for “bet-hedging.”

Moreover, the time scale of environmental change is expected to be an important factor in determining which adaptation strategy is selected. Nevertheless, in the present model, the hosts do not receive parasite information explicitly. The effect of parasite infection on the host is considered only as a negative effect on population growth. The influence of the parasite does not have any effect on the host gene expression dynamics. The characteristics of gene expression dynamics are determined by the noise level σ, which indicates that the transition point σ is independent of the time scale of the parasite population dynamics (see Fig H in S1 Appendix). Furthermore, because the population density is fixed, the interaction strength is constant in time, which means that the hosts do not become extinct even if they cannot cope with rapid changes in parasite dynamics. It would be interesting to study the evolution of a model in which the host receives the influence of the parasite directly as the input to phenotypic dynamics. For example, in the perceptron model that produces a favorable phenotype for an input, different strategies emerge depending on the environment [50]. By introducing such parasite inputs, host adaptation strategies are enriched.

In summary, we demonstrated the evolution of phenotypic variation and robustness under host–parasite interactions based on the change in phenotypic variances. Depending on the strength of the interaction and the noise level in the gene expression dynamics, either phenotypic variation or genetic diversification evolves. A stronger manifestation of genetic diversification results in speciation. In our forthcoming paper, we aim to show how the isogenic phenotypic variation induced by host–parasite interactions will lead to genetic speciation.

Methods

We define the two phenotypic variances V and V as follows.

Gene expression can take different values even among individuals of the same genotype because of noise during the development process. The variance, denoted as , is defined by the variance of for gene i of the l-th host genotype in the isogenic population, i.e., among those with the same genotype (network ).

is defined as the phenotypic variance caused by the genetic variation resulting from mutation for hosts of the genotype l. We first computed , i.e., the average of over noise, for each genotype l to distinguish the variance by mutation and that by noise. Then, we computed the variance of over a heterogenic population created by adding random one-step mutations to the interaction matrix of this genotype l. In the simulation, there are up to N individuals with different genotypes , each of which potentially takes a different gene expression pattern. To calculate the and , L individuals with different genotypes are selected at random from the population. For each genotype individual (with genotype l), we created N clones. Then, we carried out the simulation for each clone under noise. We computed the variance in gene expression caused by noise for each gene at t = 1000, after gene expression reaches a steady state. This gives rise to the value , that is, the noise-induced variance of gene i for genotypes l. Besides the variance, the average is obtained for each clone l. Subsequently, we created genetic variation in by adding a single random mutation to the network. From each of mutated genotype l′, we computed the average over N clones. From L different mutated genotypes, was computed as the variance of the over L different mutants. We used L = 30 and N = 50 for all simulations.

Supplementary figures.

Fig A: Evolutionary change of the variances for output genes for c = 0. Fig B: Average variances of V(i) and V(i) (non-target output genes i = 4, 5, …, 8). Fig C: Average variances of and (target genes i = 1, 2, 3). Fig D: Average variances of V(i) and V(i) (other genes i = 9, 10, …, 64). Fig E: Evolutionary time course of the average phenotypic variations V and . Fig F: Evolutionary change of the variances for target genes and output genes for c = 0 and c = 3. Fig G: Time course of the variances and the growth rate μ and . Fig H: Timescale of the parasite changes. Fig I: Phenotypic variance V and V of the output genes against the parasite interaction strength.

(PDF)

Click here for additional data file.

22 Mar 2021

Dear Prof. Kaneko,

Thank you very much for submitting your manuscript "Evolution of Phenotypic Plasticity under Host-Parasite Interactions" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Jacopo Grilli

Associate Editor

PLOS Computational Biology

Natalia Komarova

Deputy Editor

PLOS Computational Biology

***********************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: See the attachment

Reviewer #2: see attachment

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: None

Reviewer #2: No: Computer code and simulation results should be provided.

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Anton S. Zadorin

Reviewer #2: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

Submitted filename: Review-PCOMPBIOL-D-21-00026_Kaneko.pdf

Click here for additional data file.

8 Jun 2021

Submitted filename: Response to reviewers.docx

Click here for additional data file.

6 Jul 2021

Dear Prof. Kaneko,

Thank you very much for submitting your manuscript "Evolution of Phenotypic Fluctuation under Host-Parasite Interactions" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Jacopo Grilli

Associate Editor

PLOS Computational Biology

Natalia Komarova

Deputy Editor

PLOS Computational Biology

***********************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The authors have replied all my questions and have resolved all my concerns. I recommend the article for the publication. The reasons for my recommendation are explained in my first review.

Reviewer #2: General comments:

-----------------

I am still not too happy with the revised version, and I think another

major revision is necessary. My main points of critique are the

following:

1) The writing still needs lots of improvement. It might be good if

the authors could have the manuscript checked by a native English

speaker. I make some suggestions below, but they are far from

exhaustive.

2) I think I understand the model assumptions much better now, but

some points need to be made clearer for the reader. I also think that

there is a problem with the measure V_ip, as part of this variance

does not affect fitness (see below). I therefore suggest to look at a

version of V_ip that only looks at on/off switches of genes, not

small variation in expression level.

3) I am still having a hard time understanding some of the results.

While this is in part due to the "black-box nature" of the GNRs, I

think the manuscript might be improved by (i) adding a measure of V_ip

that focuses on on/off switches (see above) and (ii) adding

information about the relation between V_ip and V_g in random

networks.

I detail these comments in the following section.

Major comments:

--------------

- For a single gene, if the function $f$ is assumed to be constant,

the Langevin equation (1) describes an Ornstein-Uhlenbeck process,

in which gene expression levels x_i fluctuate around a mean of f +

epsilon. More precisely, the stationary distribution is Gaussian

with mean f + epsilon \\approx f and variance sigma^2/(2 gamma); this

should be mentioned in the manuscript; for the largest noise level

sigma = 0.05 (and gamma = 0.1), the variance of the OU process equals

0.0125, and the corresponding standard deviation is 0.112, meaning

that, for f close to either 0 or 1, genes can still be essentially

"switched on and off", with expression levels in the vicinity of 1

or 0 (but note that in the latter case, the actual x values can

become negative, unless this is somehow excluded in the simulations;

please clarify). Importantly, the variation around these values is

insignificant for fitness (the interaction with parasites depends

only on "rounded" or discretized values, and growth rate mu depends

on the mean expression level of output genes), even though it is

essential for the dynamics of the gene-regulatory network (i.e.,

for triggering noise-induced variation).

- I think it is useful to introduce the terms mutational robustness

(which prevents changes in gene-expression patterns due to mutations

in the network structure) and developmental robustness

(which prevents on/off switches in the face of the white-noise

fluctuations in gene expression level).

- The variance sigma^2/(2 gamma) of the Ornstein-Uhlenbeck

process (i.e., the variance of small-scale fluctuations in

expression levels) directly contributes to the measure V_ip (which

should be the sum of this variance and the variance due to on/off

switches). This blurs the interpretation of some results, and in

particular, makes it hard to see evidence for evolution of

developmental robustness. For example, in Fig. 4A, V_ip increases

with increasing sigma, but we cannot say how much of this effect is

due to the direct effect of expression-level fluctuations (increase

in sigma^2/(2 gamma)), and how much due to an increase in

noise-induced on/off switches. Similarly, in Fig. 5, the blue lines

(V_ip) beyond sigma = 0.02 appear to be roughly parallel to

sigma^2/(2gamma), which would suggest that noise-induced variation is

triggered around $sigma = sigma_c$, but that nothing much happens

beyond this value in terms of genes being switched on or off.

- Since for the interaction with the parasite only the on/off switches

matter, I suggest to look at a measure V_ip' that reflects only this

variation. Such a measure might be V_ip - sigma/(2gamma), the

variance of x_i after rounding to 0 or 1, or the probability

that a given gene is switched on (note that when x_i values are

discretized to 0 and 1, the resulting variance is the variance of

Bernoulli distribution, which is p(1-p), where p = P(x = 1)). If

this measure goes down (relative to the value in an unevolved

network), it is evidence for evolution of developmental robustness.

- I repeat that V_g is not standing genetic variation (this would be the

variance of $\\bar x$), and in contrast to what the authors claim in

the reply letter, V_g + V_ip does not equal V_p. One way to

see this is that V_ip is expressed within a generation, but V_g only

between generations. So be careful to not call V_g the variance "due

to genetic differentiation". A population might be genetically

homogeneous and still have positive V_g.

- I am still struggling to understand the initial variance values in

Fig. 3 and 6 (e.g., why is V_g so low initially in Fig. 3[i]?). You

write that genotypes that die out due to parasitism are excluded

from the sample, but do these extinctions happen immediately at t =

0?

- In addition, you write that in a previous paper you have shown that

V_ip and V_g are correlated (I suppose this means positively

correlated?). For readers of the present paper, I think it would be

useful to repeat these results, that is, to show the relation of

V_ip (or V_ip') and V_g in random networks for various values of

sigma.

- With the additional information mentioned above, I hope to better

understand some of the results. Currently, I still have the

following questions:

- Fig. 1(i): I don't understand why for c = sigma = 0, there are such

large fluctuations in host phenotype frequencies, with the optimal

phenotype sporadically dropping to about 50% of the total

population. As you notice this is evidence that mutational

robustness does not evolve in regime I, but why? And where do

these fluctuations come from exactly?

- As mentioned above, in Fig. 4A, it would be nice to see the

behavior of V_ip' (variance due to on/off switches only), and in

particular, whether it increases as a function of sigma. If

V_ip' does not increase (or increases less than in random

networks) this would be evidence for developmental robustness.

In addition, we see that V_g decreases as a function of sigma

(mutational robustness). This might be either a correlated

response to increased developmental robustness, or (in case

developmental robustness is absent) and kind of "compensatory"

evolution (but this would mean that evolution of developmental

and mutational robustness can be decoupled, or at least that one

evolves easier than the other).

- In Fig. 5, as also mentioned above, it seems to me that V_ip'

does not increase beyond sigma = 0.02. In this case, why does

V_g decrease? Again, is this a correlated response to the

evolution of developmental robustness (which keeps V_ip

constant), or is there a different cause? For example, I have

been wondering whether V_g automatically decreases with sigma

(maybe mutations have less effect if expression levels

fluctuated all the time anyway?, see my question about random

networks above).

Some more question:

- In the presence of noise-induced variations, do expression patterns

remain for individual cells once they have settled down to a stable

state, or do they occasionally change? In other words, does

noise-induced variation manifest only in the differences between

isogenic individuals, or also in single individuals over time?

- In the absence of noise (sigma = 0), do the expression dynamics of

isogenic individuals always converge to the same pattern, or are

can there be differences cause by differences in the initial random

expression levels?

Specific comments:

-----------------

- Equation (1): As mentioned above, for fixed $f$, this describes an

Ornstein-Uhlenbeck process (a random walk with a tendency to return

to the mean), in which $x_i$ fluctuates around $f + \\epsilon$ with

variance $\\sigma/(2\\gamma)$.

- 93: Make it easier for the reader to understand that (i) there

are are M (= 64) genes in total, (ii) the first $l_{out} = 8$ genes

are output genes affecting growth, (iii) of these the first $l_p =

3$ are target genes, and (iv) the last $l_{inp} = 8$ are input

genes.

- 96: Maybe the formula in this line is needlessly complicated

here, and it is enough to just say that $\\eta$ describes white noise

which uncorrelated across genes. This is especially true if you

state that equation (1) describes an Ornstein-Uhlenbeck process (and

maybe mention that it is a Langevin equation, even though most

biologists will not know this term).

- 94: k should be l_inp.

- 98: Does epsilon play any role at all for the results?

- 102: I = 1 means that input genes are always expressed.

(Indeed, the idea of external input seems unnecessary here; you

might just say that these genes are always/constitutively

expressed).

- 104: output genes have already been mentioned in line 94.

- 106: If I understand correctly, N = 300 is constant throughout,

as in the Moran model.

- Equation (3) and (4): I suppose that the growth rate is evaluated

instantaneously at each moment and can change over time (so it should be

$\\hat \\mu^{(l)}(t)$ in equation (5).

- 130-1: the part about parasite dynamics should be moved to the next

paragraph.

- 136: When there are 300 cells and one divides, the upper limit N

will always be exceeded, so it is more precise to say that whenever

a cell divides, another randomly selected dies.

- 140: What is s? Seconds? Maybe better to just say between 500 and

1000 time units.

- 141: When the gene-regulatory network reaches a steady state, but

there is phenotypic variation due to noise, do individual networks

settle into a stable state (which may vary between individuals), or

are genes occasionally switched on or off in response to

fluctuations in expression levels?

- 149: alpha = 0.3 is not just an initial choice, but a choice that

remains valid throughout the simulation.

- Equation (6): c is the attack parameter from equation (4). Here it

should probably be replaced by 1/T_p.

- 178: As mentioned above, V_ip has two components: One to the

fluctuations in eta, and one due to noise-induced on/off

switches.

- 185: V_g is not the variance due to genetic differences, but the

variance induced by mutations. In particular, unlike V_{ip}, V_g

is not expressed within a generation, but only between generations.

It is really important to make this clear, especially since readers with a

background in quantitative genetics will tend to assume that V_g is

the genetic variance (i.e., the variance in $\\bar x$), which is not

the case.

- Table 1: Is the mutation rate m = 0.0002 (as in the table) or 0.03

(as in line 143)?

- 200: From which sample? The one for t = 0?

- 241: Again, V_g is not the genotypic variance, but the variance

induced by new mutations.

- 247: Initial values of V_g are high because random networks are not

canalized and highly sensitive to mutation.

- 266-7: V_ip is higher, not lower, under large sigma for small c in Fig. 4

(regime 2).

- 268: Where can I see that V_ip and V_g both decrease (the purple

lines in Fig. 2?)

- 286: c = 4?

- 152: Unlike host dynamics, which are stochastic and

individual-based, parasite dynamics are modelled via a deterministic

model of genotype frequencies.

- 177-180: These two sentences seem redundant, and they disrupt the

connection between "both from noise and from network mutations"

(177) and "We define two phenotypic variances ..." (180).

- 235: the number of coexisting host genotypes approaches the maximal

number 2^{lp} = 8, corresponding to the total number of parasite

genotypes.

- 259: are not

- 312: GRNs that became extinct due to parasite infections are excluded

- 317: and robustness increased

- 320: Fig. S6

- 324: whether high variability in phenotype can evolve in response

to...

- 327: once the host had adapted

- 331: before increase phenotypic variation evolves

- 332: Thus, notwithstanding ..., phenotypic fluctuations can evolve

in response to parasite infection

- 338: A comparison of Fig. 6[i] and [ii] shows the final variances

are nearly identical. This means that the evolution variances in

both target and output genes is independent of the initial

conditions.

- 358: the expression of target genes is variable

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: No: I answer "No" to the question of the computational code availability, as it has not been yet published at the date of the review request, but the authors promised to do so on Github in the future.

Reviewer #2: No: Authors say they will upload the code to Github upon publication.

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Anton S. Zadorin

Reviewer #2: Yes: Michael Kopp

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

27 Aug 2021

Submitted filename: Response to reviewers.docx

Click here for additional data file.

7 Oct 2021

Dear Prof. Kaneko,

Thank you very much for submitting your manuscript "Evolution of Phenotypic Fluctuation under Host-Parasite Interactions" for consideration at PLOS Computational Biology. As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. The reviewers appreciated the attention to an important topic. Based on the reviews, we are likely to accept this manuscript for publication, providing that you modify the manuscript according to the review recommendations.

Please prepare and submit your revised manuscript within 30 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to all review comments, and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Thank you again for your submission to our journal. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Jacopo Grilli

Associate Editor

PLOS Computational Biology

Natalia Komarova

Deputy Editor

PLOS Computational Biology

***********************

A link appears below if there are any accompanying review attachments. If you believe any reviews to be missing, please contact ploscompbiol@plos.org immediately:

[LINK]

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: I have not issues with the current version of the manuscript and I, as before, recommend it for the publication.

Reviewer #2: In this revision, the authors have answered most of my questions. In

particular, the new Fig. is a real improvement, which shows that the

variance due to noise-induced on/off switches does not increase for

noise levels beyond sigma = 0.02. In would be good to discuss why this

is so. One potential explanation is that the attained level of

noise-induced phenotypic variation, in conjunction with genetic

variation, is enough to produce the maximal phenotypic variation with

8 coexisting target-gene expression patterns.

Editorial comments:

------------------

- Abstract, first "blue" sentences: two typos in plasticity

- Abstract, line 9: where one of the most important...

- Abstract, second "blue" part: By using a dynamic model of gene

expression for the host ... evolution of the genotype-phenotype map

and of phenotypic variances

- 8: Maybe: In this context, phenotypes are shaped by dynamic

processes involving several genetically determined variables, as

well as external and internal noise. Hence phenotypes can

vary not only as a consequence of genetic variation, but even

between isogenic individuals [6-10]. The sum of both sources of

phenotypic variation defines the genotype–phenotype mapping [14-16].

- 31: remain

- 46: latter's

- 81: I still think it would be useful for readers to say that, for

constant f, equation 1 defines an Ornstein-Uhlenbeck process, and to

give the equilibrium variance of sigma^2/(2gamma).

- 92: total gene number

- 110: constant at N = 300 throughout all simulations

- 125: Let m be the parasite type

- 142: after around 200 time units

- 145: value

- 145: one of the others

- 151: alpha = 0.3, and alpha remains valid...

- 151: of this choice

- 153: Unlike host phenotype (gene expression levels), which is

partially stochastic, the parasites' phenotype is determined uniquely by

its genotype. Hence parasites of a given genotype share a common

phenotype and parasite population dynamics can be modelled by

deterministic equations (instead of individual-based simulations).

- 192: the variance ... within a population of initial genotype i in

which each individual receives one random mutation.

- 199: Note that developmental noise (the sigma eta term in equation

1) has no significant effect on host fitness once gene expression

dynamics have reached a stable state with all x_i close to 0 or 1.

This is because host growth depends only on mean gene expression

levels and parasite attack depends on rounded values (whether x_i is

< or > 0.5). However, noise is important, because when amplified in

the GRN, it may cause the initial dynamics to converge to a

different stable state (i.e., noise-induced switches in gene

expression). Regaring phenotypic variances, we discuss two types of

robustness: mutational robustness and developmental robustness.

Mutational robustness describes the "ability" the ability of a GNR

to prevent changes in gene-expression patterns due to mutations in

the network structure, and developmental robustness describes the

ability to prevent changes in expression patterns due to white-noise

fluctuations in gene expression level during the initial network

dynamics.

- 209: We started simulations with N host cells of an initial

gene-regulatory network created by creating random paths with

probability alpha = 0.3.

- 210: What are the 50 samples? 50 individuals from a population, or

50 replicated populations?

- 211: The initial condition for each cell is given by a state in

which only the input genes are expresses: all other genes are

initialized with expression levels drawn uniformly between 0 and

0.05. This same initial condition was also applied to the daughter

cells after each cell division.

- 214: that went extinct

- 216: add "in the absence of noise (sigma = 0)" at the end of the

sentence

- 229: distribution of parasite genotypes

- 231: Note that in the absence of noise...

- 233: Next, we computed

- 233: expression levels

- 236: Independent of the level of developmental noise, the Figure

shows a marked increase...

- 238: This sentence is redundant

- 255: These correspond to two strategies

- Fig. 3, new sentence: Note that the plot starts with the first

rather than te 0th generation, and that host types that went extinct

during the simulation have been exluded. This explains the different

initial variance levels between parameter combinations.

- 261: after evolution of the system reaches a stable state

- 263: These results are independent of the choice of the initial

network

- 287: to both noise and mutation

- 294: whereas V_ip remains small

- 299: Maybe: For moderate c > c_t, the genotypes with the fittest

type 111 coexists with genotypes having one target gene switched off

..., whereas for c = 4, even more host types appear, leading to

coexistnence of all eight target gene expression patterns.

- 307: The sentence starting with "In addition" seems redundant with

respect to the previous one.

- 309: Variances in expression levels are highest for target genes

- 320: Note that, in Fig. 5, V_g decreases slightly with sigma once

sigma exceeds the critical value sigma_c. This reflects previous

results for c = 0 (i.e., in the absence of parasites), which showed

that mutational robustness decreases as a correlated response to the

evolution of developmental robustness against increasing noise

levels (Fig. S8).

- 329: As mentioned above, once expression levels of host target genes

have settled down close to one of the fixed points x = 0 or 1,

continued variation due to the white-noise term in Fig. 1 has no

effect on the interaction with the parasite. The noise is important,

however, for the initial network dynamics that determine which of

these fixed points is approached in the first place. In order to

separate the contributions of noise-induced on/off-switches from

those of the noise itself, we define a variance V_ip' which only

reflects the former. More precisely, V_ip'(i) is the variance

of the on/off variable z_i = ..., computed in the same manner as

V_ip(i).

- 336: After the first sentence, first discuss Fig. 6: "The key

difference to Fig. 5 is that for sigma > 0.02, V_ip' remains

essentially constant for both target and output genes, suggesting

that the variance due to noise-induced switches cannot be increased

arbitrarily, and that the continued increase of V_ip seen in Fig. 5

is merely a reflection of the increase in white noise."

- Then go on to discuss Fig. S9: E.g., "Evolutionary time

courses for both V_ip and V_ip' are shown in Fig. S9. For target

genes, both measures increase from their initial value and then stay

at a high level, with no evidence for evolution of developmental

robustness. For non-target output genes, in contrast, V_ip and V_ip'

decrease over time relative to the values in the initial random

networks. This shows that ..."

- 345-351: I don't understand what you are trying to say here.

- 352-355: Thanks for adding this, but it be better placed earlier in

the paper, for example in the model description, after line 151.

- in the 0th generation?

- 368: closer than what?

- 370: remove "the" before developmental robustness

- 379: after developmental robustness has evolved

- 381: and let the system evolve

- Caption Fig. 7: It seems that there is a problem with the sentence

right after the blue insert.

- 418: within a regulatory network

- 421: the phenotypic variance of the host increases as strains evolve

adaptation to decrease attack by specific parasites.

- 427: which enhances phenotypic variation and helps...

- 428: The new sentence seems redundant.

- 430: both variances

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: No: The program text is promised to be published in the future.

Reviewer #2: No: Simulation code should be made available

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Anton S. Zadorin

Reviewer #2: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

References:

Review your reference list to ensure that it is complete and correct. If you have cited papers that have been retracted, please include the rationale for doing so in the manuscript text, or remove these references and replace them with relevant current references. Any changes to the reference list should be mentioned in the rebuttal letter that accompanies your revised manuscript.

If you need to cite a retracted article, indicate the article’s retracted status in the References list and also include a citation and full reference for the retraction notice.

12 Oct 2021

Submitted filename: responce to reviewers.docx

Click here for additional data file.

19 Oct 2021

Dear Prof. Kaneko,

We are pleased to inform you that your manuscript 'Evolution of Phenotypic Fluctuation under Host-Parasite Interactions' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Jacopo Grilli

Associate Editor

PLOS Computational Biology

Natalia Komarova

Deputy Editor

PLOS Computational Biology

***********************************************************

1 Nov 2021

PCOMPBIOL-D-21-00026R3

Evolution of Phenotypic Fluctuation under Host-Parasite Interactions

Dear Dr Kaneko,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Zsofia Freund

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol