ignore ignorerecommendation

††margin:  RESEARCH ARTICLE

Open Access

Open Peer-Review

Open Code

Cite as: Le Rouzic A (2022) Gene network robustness as a multivariate
character. arXiv:2101.01564, ver. 5 peer-reviewed and recommended by Peer
Community in Evolutionary Biology. https://arxiv.org/abs/2101.01564
Posted: 30th March 2022
Recommender:
Frédéric Guillaume
Reviewers:
Charles Rocabert, Diogo Melo and Charles Mullon
Correspondence:
<EMAIL_ADDRESS>

Gene network robustness as a multivariate character

Arnaud Le Rouzic1

1 Laboratoire Évolution: Génomes, Comportement, Écologie; Université Paris-
Saclay, CNRS, IRD – Gif-sur-Yvette, France

This article has been peer-reviewed and recommended by
_Peer Community In Evolutionary Biology_
(https://doi.org/10.24072/pci.evolbiol.100138)

Abstract
Robustness to genetic or environmental disturbances is often considered as a
key property of living systems. Yet, in spite of being discussed since the
1950s, how robustness emerges from the complexity of genetic architectures and
how it evolves still remains unclear. In particular, whether or not robustness
is independent to various sources of perturbations conditions the range of
adaptive scenarios that can be considered. For instance, selection for
robustness to heritable mutations is likely to be modest and indirect, and its
evolution might result from indirect selection on a pleiotropically-related
character (e.g., homeostasis). Here, I propose to treat various robustness
measurements as quantitative characters, and study theoretically, by
individual-based simulations, their propensity to evolve independently. Based
on a simple evolutionary model of a gene regulatory network, I showed that
five measurements of the robustness of gene expression to genetic or non-
genetic disturbances were substantially correlated. Yet, robustness was
mutationally variable in several dimensions, and robustness components could
evolve differentially under direct selection pressure. Therefore, the fact
that the sensitivity of gene expression to mutations and environmental factors
rely on the same gene networks does not preclude distinct evolutionary
histories of robustness components. _Keywords:_ Gene regulatory network;
Transcription regulation; Wagner model; Individual-based simulations;
Canalization

## Introduction

Robustness is the capacity of living organisms to buffer internal or
environmental disturbances. Robustness encompasses, for instance, the ability
to maintain physiological equilibria (homeostasis), to ensure developmental
stability, or to repair and mitigate DNA damage in both soma and germline.
Although robustness is virtually intermingled with the definition of life
itself, its underlying mechanisms and its evolutionary origins remain far from
being clearly understood [58, 41, 66, 22].

Robustness evolves as a consequence of non-linearities in the developmental or
physiological mechanisms, i.e. changes in the magnitude of the effect of some
genetic or environmental factor on the phenotype of interest [44]. The study
of the evolutionary processes leading to robustness roots into the conceptual
and empirical work by C.H. Waddington and the concept of canalization [62, 54,
63, 39]. Canalization is a property of complex developmental systems that
buffers environmental and genetic variation, and maintains actively the
organism in an optimal developmental path. Although the scope and the
definition of canalization varies substantially among authors, canalization is
generally expected to evolve as an adaptation to "canalizing" selection for an
optimal phenotype [13, 9, 19, 32]. However, formal population genetic models
have questioned the unicity of the canalization process. In particular,
robustness to environmental factors appears more likely to evolve as an
adaptation than robustness to genetic (mutational) disturbances, on which
selection seems to be rather weak and indirect even in optimal theoretical
conditions [68, 25, 35].

In this context, the evolution of robustness as a general property of
organisms heavily depends on the genetic and physiological integration of the
different robustness dimensions [16, 17]. If the robustness to environmental
factors and to genetic mutations share the same physiological bases, the
adaptive evolution of environmental canalization can generate a correlated
response of genetic canalization; this hypothesis has been referred to as
"congruent evolution" [8], and have recieved some empirical support [37, 61].
In contrast, if genetic and environmental robustness had independent
biological bases, they would be featured by independent evolutionary
mechanisms, and possibly independent evolutionary histories.

Although this issue would benefit from a better theoretical framework,
modeling the evolution of robustness is not straightforward. The simplest
approach relies on modifiers, i.e. genes that can influence the robustness of
the organism without affecting the phenotype. However, in the case of genetic
robustness, modifier-based models either rely on tricky rescaling or cannot
dissociate the phenotype and the robustness to the phenotype [68, 30, 51]. In
addition, in models where the genotype-phenotype association is arbitrary
(such as the NK model, [29], or the multilinear model, [24]), any correlation
between environmental and genetic robustness would be a modeling choice, and
not an output of the model. More promising to address the congruent evolution
issue are models in which the phenotype is a result of an integrated process
mimicking some developmental or physiological mechanism (referred to as
_causally cohesive genotype phenotype_ models in [50]). In such dynamic
models, robustness to various disturbances appear as an emergent property of
the model complexity, caused by regulatory feedbacks, that cannot be easily
deduced from the model parameters. Although the potential palette of relevant
dynamic models is large and could include morphological development models
[42], RNA folding models [67], or metabolic models [45], evolutionary
biologists have often considered gene regulatory network models as a good
compromise between complexity and numerical tractability for studying the
evolution of canalization and robustness [28, 64, 56, 33].

Such theoretical gene networks have been shown to display enough non-
linearity, leading to epistasis and pleiotropy, to evolve enhanced or reduced
sensitivity to environmental [40, 14, 15] and genetic [65, 3, 10, 2, 53]
perturbations. Interesting observations suggest that environmental or genetic
canalization could be correlated to other robustness properties in such
models. For instance, [7, 27] noticed that robustness to mutations and
robustness to noise was correlated in gene networks — a similar result was
obtained earlier for RNA-folding structures [20]. Furthermore, it has been
shown that network stability, the propensity of the network to maintain stable
(non-cyclic) gene expressions, was correlated to robustness, as selection on
stability alone could drive an indirect response of genetic [55] and
environmental [40, 43] canalization. In contrast, [46] showed that networks
selected to maintain (but not converge to) an equilibrium became both
environmentally sensitive and genetically canalized, suggesting that
environmental and genetic robustness could be theoretically decoupled.
However, no systematic quantitative description of the pleiotropic pattern
underlying different robustness components has ever been attempted.

Here, I aim at extending the study of canalization in theoretical gene
networks to address the multidimensional nature of robustness, by estimating
the evolutionary independence of various robustness components. Four
robustness-related measurements were considered, two of them corresponding to
environmental robustness (early vs. late disturbances), two corresponding to
genetic robustness (early — inherited — or late — acquired — mutations). Gene
expression instability was also included in the set of robustness-related
traits, as it is related to the intrinsic stability of the expression
phenotype. The first part of this study focuses on the multidimensional
patterns of robustness in small and random networks, and the second part on
the evolutionary consequences of the pleiotropic nature of robustness, based
on individual-based simulations.

## Model and Methods

### Gene regulatory network

The network model belongs to the family of gene regulatory network models
sometimes referred to as "Wagner model" (after [64, 65]; see [18] for a
historical record). Two variants of the model were proposed in [64]; the
second one, involving discrete gene expressions scaled between $-1$ and $1$,
has often been reused in the literature [65, 55, 7]. The model described below
is closer to the first model by [64], featuring a continuous gene expression
$\bm{\mathrm{P}}$ between 0 and 1, and a constitutive expression level $0<a<1$
that can be lower than the mid-expression point.

More specifically, the structure of a $n$-gene network is encoded as a
$n\times n$ matrix $\bm{\mathrm{W}}$, while the state of the network is stored
into a vector of size $n$, $\bm{\mathrm{P}}$. In this setting, $W_{ij}$
encodes the influence of gene $j$ on the expression of gene $i$, $W_{ij}<0$
represents a negative interaction (inhibition), $W_{ij}>0$ a positive
interaction (activation), and $W_{ij}=0$ denotes the absence of regulatory
interaction. $P_{i}$ is the expression of gene $i$, ranging between 0 (no
expression) and 1 (maximum expression).

The properties of these gene networks are explored in a discrete dynamic
system:

$~\bm{\mathrm{P}}_{t+1}=F(\bm{\mathrm{W}}\bm{\mathrm{P}}_{t}),$ (1)

where the function $F$ is a vectorized version of a sigmoid scaling function:
$F(x_{1},x_{2},\dots,x_{n})=[f(x_{1}),f(x_{2}),\dots,f(x_{n})]$;

$f(x)=\frac{1}{1+\lambda_{a}e^{-\mu_{a}x}},$ (2)

with $\lambda_{a}=(1-a)/a$ and $\mu_{a}=1/a(1-a)$ [21]. The function $f$ is
scaled such that $f(0)=a$ and $df/dx|_{x=0}=1$; the parameter $a$ thus stands
for the constitutive gene expression (the expression of a gene in absence of
regulators), and this function defines the scale of the matrix
$\bm{\mathrm{W}}$: $W_{ij}=\delta$ ($\delta\ll 1$) means that the expression
of gene $i$ at the next time step will tend to $P_{i,t+1}=a+\delta$ if $i$ is
regulated by a single, fully expressed transcription factor $j$ ($P_{j,t}=1$).
This setting, extensively described in [53], differs mathematically from the
constitutive expression model in [64] that shifts the sigmoid as
$\bm{\mathrm{P}}_{t+1}=F(\bm{\mathrm{W}}\bm{\mathrm{P}}_{t}+a)$.

Gene networks dynamics start from an initial expression $\bm{\mathrm{P}}_{0}$,
and gene expression was updated for $T$ time steps. By default,
$\bm{\mathrm{P}}_{0}=(a,a,...,a)$, since this step immediately follows a
virtual initial state with no expression. The expression phenotype
corresponding to a gene network was determined by averaging gene expressions
during the last $\tau$ time steps for each gene $i$:
$P_{i}^{*}=(1/\tau)\sum_{t=T-\tau}^{T}P_{it}$.

### Robustness indicators

Five robustness indicators were calculated, corresponding to five different
aspects of genetic or environmental robustness in a gene network: robustness
to early (${\rho_{E}}$) and late (${\rho_{e}}$) environmental disturbance, and
robustness to early (${\rho_{M}}$) and late (${\rho_{m}}$) genetic
disturbance, and network stability ${\rho_{S}}$. All indicators were expressed
on a scale homogeneous to log variances in gene expressions; the mode of
calculation is summarized in Table 1, robustness is maximal when the index
$\rho$ is small.

The robustness to early environmental disturbance ${\rho_{E}}$ measures the
capacity of a network to reach a consistent final state starting from
different initial gene expressions. In practice, $R$ replicates of the network
dynamics were run, in which the initial gene expressions
($\bm{\mathrm{P}}_{0}$) were drawn into Gaussian ($\mu=a,\sigma=\sigma_{E}$)
distributions (expression values $<0$ and $>1$ were set to $0$ and $1$,
respectively). The environmental robustness ${\rho_{E}}_{i}$ for each gene $i$
was measured as the log variance in the final gene expression across these
replicates.

The robustness to late environmental disturbance ${\rho_{e}}$ measures the
capacity of a network to recover its equilibrium state after having being
disturbed. Gene expressions after $T$ time steps were disturbed by adding a
random Gaussian noise of standard deviation $\sigma_{e}$ to each gene of the
network, and ${\rho_{e}}_{i}$ was computed for each gene $i$ as the log
variance in gene expression at time step $T+1$ over $R$ replicates.

The robustness to early mutations ${\rho_{M}}$ measures the system robustness
to inherited genetic mutations (modifications of the $\bm{\mathrm{W}}$
matrix). A random non-zero element of the $\bm{\mathrm{W}}$ matrix was shifted
by a random Gaussian number of standard deviation $\sigma_{M}$, and its
consequences on the mean expression of all network genes was recorded. The
procedure was replicated $R$ times, and the robustness score ${\rho_{M}}_{i}$
for each gene $i$ was calculated as the log variance of gene expression across
$R$ replicates.

The robustness to late mutations ${\rho_{m}}$ measured the effect of mutations
in the gene network $\bm{\mathrm{W}}$ after having reached the final state. In
practice, the $\bm{\mathrm{W}}$ matrix was mutated in the same way as for
${\rho_{M}}$ with a standard deviation $\sigma_{m}$, but its consequences on
gene expression were calculated for only one time step, starting from the last
state of the network. The robustness score was calculated as for other
indicators (log variance over $R$ replicates).

Finally, dynamic systems based on the Wagner model often tend to generate
limit cycles and never converge to a stable equilibrium. Network stability
${\rho_{S}}$ quantifies the capacity for a specific network to lead to stable
gene expressions. For consistency with other indicators, this instability was
measured as the log squared difference between the average expression during
the last $\tau$ time steps, and an extra time step.

Indicator |  Robustness component | Computation | Disturbance std. dev.
---|---|---|---
${\rho_{E}}$ |  Early noise in gene expression | ${\rho_{E}}_{i}=\log[\frac{1}{R-1}\sum_{r=1}^{R}(P_{i,r}^{*}-\overline{P_{i}^{*}})^{2}]$ | $\sigma_{E}=0.1$
${\rho_{e}}$ |  Late noise in gene expression | ${\rho_{e}}_{i}=\log[\frac{1}{R-1}\sum_{r=1}^{R}(P_{i,T+1,r}-\overline{P_{i,T+1}})^{2}]$ | $\sigma_{e}=0.1$
${\rho_{M}}$ |  Early (inherited) mutations | ${\rho_{M}}_{i}=\log[\frac{1}{R-1}\sum_{r=1}^{R}(P_{i,r}^{*}-\overline{P_{i}^{*}})^{2}]$ | $\sigma_{M}=0.1$
${\rho_{m}}$ |  Late (aquired) mutations | ${\rho_{m}}_{i}=\log[\frac{1}{R-1}\sum_{r=1}^{R}(P_{i,T+1,r}-\overline{P_{i,T+1}})^{2}]$ | $\sigma_{m}=0.1$
${\rho_{S}}$ |  Expression stability | ${\rho_{S}}_{i}=\log[(P_{i}^{*}-P_{T+1})^{2}]$ |
Table 1:  Summarized calculation of all five robustness indicators. Index $i$
stands for the gene ($1\leq i\leq n$), and $r$ for the replicate ($1\leq r\leq
R$), since all indicators except ${\rho_{S}}$ were estimated by a resampling
procedure. $P_{i}^{*}$ stands for the equilibrium gene expression of gene $i$
(mean expression from the last $\tau$ time steps), and
$\overline{P_{i}^{*}}=(1/R)\sum_{r=1}^{R}P_{i,r}^{*}$ represents the mean over
replicates. Noise in gene expression was simulated by adding a random Gaussian
deviation to the initial state $\bm{\mathrm{P}}_{0}$ of the network (for
${\rho_{E}}$) or to the last state $\bm{\mathrm{P}}_{T}$ of the network (for
${\rho_{e}}$). Mutations were simulated by adding a random deviation to a
random interaction in the network $\bm{\mathrm{W}}$, either before starting
the network dynamics (${\rho_{M}}$) or after the last time step
(${\rho_{m}}$). All robustness indicators are homogeneous to a log variance in
gene expression; robustness increases when the indicator gets smaller, and
sensitivity increases when the indicator increases. The last column indicates
the standard deviation of the corresponding Gaussian disturbance.

All these scores were calculated for every gene $i$ of a given network, and
then averaged over all genes in order to get a series of summary network
descriptors. The magnitude of the score itself is arbitrary, as it depends on
the size of the disturbance. However, indicators happen to increase
approximately linearly with the size of the disturbance (Appendix 1), the
results were thus largely unaffected by a change in the variance of mutational
effects and environmental noise.

### Random networks

Random networks were generated as $n\times n$ $\bm{\mathrm{W}}$ matrices
filled with independent identically-distributed random numbers drawn into a
Gaussian (by default: $\mu_{0}=0,\sigma_{0}=1$) distribution. A density
parameter $1/n\leq d\leq 1$ could be specified, corresponding to the frequency
of non-zero slots in the $\bm{\mathrm{W}}$ matrix. Zeros were placed randomly,
with the constraint that all genes should be regulated by at least another
one.

### Exhaustive exploration of two-gene networks

The main interest of gene-network models is the complexity and the richness of
the underlying genotype-phenotype relationship. As a side effect, such models
are in general difficult to handle mathematically ([6, 33]). Excluding the
one-gene self-regulating case (which already has non-trivial mathematical
properties, [21]), the simplest network (2-by-2 matrix) has four genetic
parameters, which makes the exploration of the parameter set tedious. Here,
the number of dimensions was restricted by considering the set of networks
that lead to a predefined arbitrary equilibrium,
$\bm{\bm{\mathrm{P}}}^{\theta}_{\infty}=(P^{\theta}_{1},P^{\theta}_{2})$. As
$F(\bm{\mathrm{W}}\bm{\mathrm{P}}^{\theta}_{\infty})=\bm{\mathrm{P}}^{\theta}_{\infty}$,
the $\bm{\mathrm{W}}$ matrix can be reduced to two independent parameters,
$W_{11}$ and $W_{21}$:

$\bm{\mathrm{W}}=F\left[\begin{pmatrix}W_{11}&A\\\
W_{21}&B\end{pmatrix}~\begin{pmatrix}P^{\theta}_{1}\\\
P^{\theta}_{2}\end{pmatrix}\right]=\begin{pmatrix}P^{\theta}_{1}\\\
P^{\theta}_{2}\end{pmatrix},$ (3)

with

$\begin{split}A=\frac{1}{P^{\theta}_{2}}[f^{-1}(P^{\theta}_{1})-W_{11}P^{\theta}_{1}],\\\
B=\frac{1}{P^{\theta}_{2}}[f^{-1}(P^{\theta}_{2})-W_{21}P^{\theta}_{1}],\end{split}$
(4)

$f^{-1}(y)=-\frac{1}{\mu_{a}}\log\left(\frac{1-y}{\lambda_{a}y}\right)$ being
the inverse of $f(x)$ (equation 2). This equation can be extended to any
network size, provided that a single element $W_{ij}$ is unknown for each line
$i$ of the matrix:

$W_{ij}=\frac{1}{P^{\theta}_{j}}[f^{-1}(P^{\theta}_{i})-\sum_{j^{\prime}\neq
j}W_{ij^{\prime}}P^{\theta}_{j^{\prime}}].$ (5)

Among the $n^{2}$ elements of a $n$-gene network, there are thus $n(n-1)$
neutral dimensions that can be explored without modifying equilibrium gene
expressions. Large gene networks are thus characterized by a proportionally
larger neutral space.

The $\bm{\mathrm{W}}$ matrix achieving the desired
$\bm{\mathrm{P}}^{\theta\ast}_{\infty}$ equilibrium from a specific pair
$W_{11},W_{21}$ always exists (and is unique), but the stability of the
equilibrium is not guaranteed. Networks which final gene expression
$\bm{\mathrm{P}}^{\ast}=(P^{\ast}_{1},P^{\ast}_{2})$ differed substantially
from the target (in practice, when
$|P^{\ast}_{1}-P^{\theta}_{1}|~+|P^{\ast}_{2}-P^{\theta}_{2}|>0.15$) were
excluded from the analysis. Such discrepancies correspond to either unstable
equilibria (in which case gene expressions were driven away from the
equilibrium) or extreme oscillatory behaviors (large oscillations may hit
expression limits 0 or 1, which drives the average expression away from the
target equilibrium).

### Evolutionary simulations

The evolution of gene networks under various evolutionary constraints was
studied by individual-based simulations. Each individual was featured by its
genotype (a $n\times n$ $\bm{\mathrm{W}}$ matrix, by default $n=6$ to limit
the computational burden), its expression phenotype $\bm{\mathrm{P}}^{*}$, and
the five robustness scores ${\rho_{S}}$, ${\rho_{E}}$, ${\rho_{e}}$,
${\rho_{M}}$, and ${\rho_{m}}$. Individuals were haploid and reproduced
clonally. Mutations consisted in adding a random Gaussian deviate of variance
$\sigma_{\nu}^{2}$ to a random regulatory interaction of the $\bm{\mathrm{W}}$
matrix, with a rate $\nu$ per individual and per generation. Mutational
parameters $\nu$ and $\sigma_{\nu}$ were kept reasonably low to limit the
strength of indirect selection for genetic robustness [68, 53]. Generations
were non-overlapping, and population size $N$ was constant. A generation
consists in sampling $N$ new individuals among the $N$ parents, with a
probability proportional to the individual fitness. Fitness was computed
assuming stabilizing selection around a target (optimal) expression level for
$n^{\prime}\leq n$ genes of the network (by default $n^{\prime}=3$), as
$w=\exp(-\sum_{i=1}^{n^{\prime}}s_{i}(P_{i}^{*}-\theta_{i})^{2})$, where
$s_{i}$ was the strength of stabilizing selection on gene $i$ ($s_{i}=0$
standing for no selection), and $\theta_{i}$ was the optimal expression
phenotype. The $\theta_{i}$ were drawn in a uniform (0,1) distribution at the
beginning of each replicated simulation, and the initial gene network was
empty ($W_{ij}=0$) except for one random element per line, which was
initialized to match the optimal expression using equation (5).

The evolution of robustness components was tracked by estimating ${\rho_{S}}$,
${\rho_{E}}$, ${\rho_{e}}$, ${\rho_{M}}$, and ${\rho_{m}}$ at regular time
points. Components were estimated for each individual, and averaged out over
the population. The response to direct or indirect selection was computed as
the average change from generation 0; the multivariate response was stored as
a 5-dimension vector $\bm{\mathrm{R}}$. Simulation runs were replicated 100
times and the results were averaged out, default parameter values are provided
in Table 2.

Parameter | Symbol | Value
---|---|---
Population size | $N$ | $1000$
Gene network size | $n$ | $6$
Constitutive expression | $a$ | $0.2$
Network time steps | $T$ | $16$
Network measurement steps | $\tau$ | $4$
Network density | $d$ | $1.0$
Simulation replicates |  | $100$
Mutation rate per individual | $\nu$ | $0.01$
Size of mutational effects | $\sigma_{\nu}$ | $0.1$
Number of selected genes | $n^{\prime}$ | $3$
Stabilizing selection coefficient | $s$ | $10$
Directional selection coefficient | $\beta$ | $0$
Number of robustness tests | $R$ | $100$
Size of early environmental noise | $\sigma_{E}$ | $0.1$
Size of late environmental noise | $\sigma_{e}$ | $0.1$
Size of early genetic mutations | $\sigma_{M}$ | $0.1$
Size of late genetic mutations | $\sigma_{m}$ | $0.1$
Table 2:  Default parameter values in the evolutionary simulations.

Directional selection on robustness indicators was also performed in some
simulations, consisting in multiplying individual fitness by
$\exp(\sum_{x\in(S,E,e,M,m)}\beta_{x}\rho_{x})$, where $\beta_{x}$ was the
strength of directional (positive or negative) selection on robustness index
$x$ (in practice, $\beta_{x}=\pm 0.01$). The vector $\bm{\beta}$ is thus
proportional to the multivariate selection gradient on robustness components.
There was no correlated selection (the fitness function is the product of
independent marginal functions applied on gene expressions and robustness
components).

Estimating genetic covariance matrices $\bm{\mathrm{G}}$ was computationally
intractable in simulations (it would require a heavy resampling procedure in
each individual), mutational covariances $\bm{\mathrm{M}}$ from the average
genotype in the population ($\overline{\bm{\mathrm{W}}}$) were used instead to
derive multivariate evolutionary predictions. Mutational covariance matrices
$\bm{\mathrm{M}}=\nu\bm{\mathrm{C}}/5$ were estimated from covariances
$\bm{\mathrm{C}}$ in gene expressions and robustness coefficients among 100
gene networks differing from $\overline{\bm{\mathrm{W}}}$ by 5 mutations
(drawn from the same algorithm as during the simulations). In order to control
for the influence of stabilizing selection on gene expression on the evolution
of robustness, conditional mutational matrices (equivalent to conditional
evolvabilities of $\bm{\mathrm{G}}$ matrices in [23]) were computed as
$\bm{\mathrm{M}}_{c(y|x)}=\bm{\mathrm{M}}_{y}-\bm{\mathrm{M}}_{yx}\bm{\mathrm{M}}_{x}^{-1}\bm{\mathrm{M}}_{xy}$,
where $y$ indicate the $n_{y}$ unconstrained traits and $x$ the $n_{x}$
constrained traits (i.e. the $n^{\prime}=3$ genes under stabilizing
selection). $\bm{\mathrm{M}}_{c(y|x)}$ was thus a $n_{y}\times n_{y}$ matrix
measuring how the unconstrained traits can mutate while traits $x$ remain
constant. Predicted mutational evolvabilities in the direction of selection
$\bm{\beta}$ were calculated as
$e_{\mathrm{pred}}=\bm{\beta}^{\top}\bm{\mathrm{M}}_{c}\bm{\beta}/|\bm{\beta}|^{2}$
[23], and realized (observed) evolvabilities were obtained by projecting the
multivariate response to selection $\bm{\mathrm{R}}$ on the direction of
$\bm{\beta}$: $e_{\mathrm{obs}}=\bm{\mathrm{R}}\bm{\beta}/|\bm{\beta}|$.
Contrary to the genetic covariances $\bm{\mathrm{G}}$, mutational covariances
$\bm{\mathrm{M}}$ cannot be used directly to compute quantitative evolutionary
predictions, as the relationship between $\bm{\mathrm{M}}$ and
$\bm{\mathrm{G}}$ depends on the mutation-selection-drift equilibrium, which
is notoriously difficult to handle mathematically [5]. The following analyses
thus focus on whether mutational evolvabilities are proportional to the
selection responses, assuming that $\bm{\mathrm{G}}$ are proportional to
$\bm{\mathrm{M}}$.

Simulations and data analysis were coded in R [49], except for the core gene
network dynamics that was coded in C++ and embedded in the R code with the
Rcpp package [12]. Scripts to reproduce simulations and figures are available
online (https://doi.org/10.5281/zenodo.6393075), [34].

## Results

### Random networks

Random interaction matrices are regularly used in the literature to study the
general properties of gene networks (e.g. [6, 48]). As such, random networks
are not expected to reflect the properties of biologically-realistic genetic
architectures, as biological networks are far from random. However, such an
approach helps developing a general intuition about the properties of the
underlying model.

Correlations were calculated between all five robustness components over
10,000 random networks (Appendix 2). All robustness components were positively
correlated, correlations ranged from about 0.62 (late genetic vs. early
environmental) to above 0.97 (late environmental vs. late genetic). A
Principal Component Analysis (Figure 1A and B) confirms that robustness
components were partially correlated. The first PC (82% of the total variance)
corresponds to the general robustness of the network, and involves all
robustness indexes. The remaining variance is explained by orthogonal vectors
separating all other robustness components. At least 4 out of 5 PCs,
explaining 10% to 2% of the total variance, did not vanish when increasing the
sample size (Appendix 3). The part of the variance in robustness explained by
the first PC is robust to the network properties, as it remains around 80%
when the mean and the variance in the regulation strengths, the network
density, and the network size vary (Figure 1C, D, E, and F).

Figure 1:  A: Summary of the principal component analysis on the five
robustness indicators over 10,000 random 6-gene networks
($\mu_{0}=0,\sigma_{0}=1$), indicating the position of the five robustness
components on all five (normalized) Principal Components (PC); ${\rho_{S}}$:
Stability, ${\rho_{E}}$: Early environmental, ${\rho_{e}}$: Late
environmental, ${\rho_{M}}$: Early genetic, ${\rho_{m}}$: Late genetic. B:
relative contribution of the five PCs to the total variance. C: Influence of
the average regulation strength ($\mu_{0}$) on the % of the total variance
explained by the first PC (negative values feature inhibitory networks,
positive values activating networks). D: Influence of the standard deviation
of the regulation strength ($\sigma_{0}$). E: Influence of the network
density. F: Influence of the network size.

### Two-gene networks

In the following, I considered an arbitrary case of a two-gene network which
genes are expressed to $\bm{\mathrm{P}}_{\infty}$ = (0.3, 0.6). Equivalent
results could be achieved with a different, arbitrary target. Figure 2
illustrates how the robustness components varied in this constrained 2-gene
network model (red stands for maximum robustness, i.e. minimum scores for
${\rho_{S}}$, ${\rho_{E}}$, ${\rho_{e}}$, ${\rho_{M}}$, and ${\rho_{m}}$). All
the networks considered here converge to the same gene expression, and can
thus be considered as phenotypically equivalent ; the colored space in Figure
2 thus represents a connected neutral network in which populations can evolve,
and thus change the topology and the robustness of the gene network, while
keeping the expression phenotype constant. In the white regions, the
equilibrium was not achieved in numerical simulations for at least three
different reasons (Appendix 4): (i) fluctuations around the equilibrium were
large enough to hit the edges of the (0,1) interval, shifting the mean
expression; (ii) the expression dynamics was slow and the network was unable
to get close to the equilibrium after 16 time steps; (iii) the equilibrium was
not reachable from the default starting point.

Figure 2:  Robustness indicators (${\rho_{E}}$, ${\rho_{e}}$, ${\rho_{M}}$,
${\rho_{m}}$, and ${\rho_{S}}$) estimated for an exhaustive continuum of two-
gene networks with an arbitrary expression equilibrium at
$\bm{\mathrm{P}}_{\infty}=(0.3,0.6)$. Although two-gene networks have four
independent genetic parameters, only two were represented here, the two others
being computed to ensure the desired equilibrium. Red stands for the maximum
robustness (lowest robustness scores); yellow for minimum robustness (highest
scores). For readability, color scales are different across panels. Letters A
to E stand for five example networks illustrated in Appendix 6.

The different robustness components were correlated, but did not overlap
perfectly. In order to assess the variation of the robustness properties, five
networks of contrasted robustness, labeled from A to E, were tracked more
specifically (Figure 2; the corresponding $\bm{\mathrm{W}}$ matrices are
provided in Appendix 5). Appendix 6 illustrates the effect of various sources
of disturbance on each network dynamics. The network denoted as B was robust
to most sources of disturbance, while network E was sensitive to all
components except stability. Network C was unstable, but remained relatively
buffered. Networks A and D illustrate intermediate loss-of-robustness
behaviors, through different mechanisms (instability for network D, and weak
buffering for network A).

This 2-gene network analysis thus confirms the results obtained for large
random networks: robustness components are only partially correlated.
Robustness is not a feature of large and intricate genetic architectures, as
it is already present (and multidimensional) in the simplest gene networks.

### Evolution and evolvability of robustness

The evolution of robustness was studied by individual-based simulations, in
which all individuals were characterized by their genotype (a 6-gene network)
and a set of phenotypes (gene expressions and network robustness). Gene
expressions for 3 out of 6 genes were under stabilizing selection. In addition
to stabilizing selection on gene expression (forcing the network to maintain a
functional role), robustness indicators were directly selected towards more or
less sensitivity. Such direct, artificial selection pressures on robustness
are not designed to reflect realistic selection on gene networks, but they
might reveal evolutionary limits to the evolution of robustness due to
internal constraints. Stabilizing selection on gene expression is expected to
generate a slight selection pressure on the robustness, but this effect was
apparent only for larger or more frequent mutations (Appendix 7).

Direct selection on all robustness components lead to a response, showing that
robustness is evolvable (diagonal panels in Figure 3). Yet, the evolutionary
potential differed substantially among robustness indicators, as indicated by
the differences in the Y-scales. Robustness indicators being all homogeneous
to a sum of squared difference in gene expression (i.e., the variance in gene
expression induced by various disturbances), they could be compared directly.
The most evolvable robustness components were early environmental disturbances
(${\rho_{E}}$) and stability (${\rho_{S}}$), which can differ by up to 25 log
units (11 orders of magnitude) after 10,000 generations of bidirectional
selection. In contrast, robustness to late environmental noise ${\rho_{e}}$
and genetic changes (${\rho_{M}}$ and ${\rho_{m}}$) only differed by 3 to 4
log units (i.e. a factor 10 to 100). For these three robustness components
(${\rho_{e}}$, ${\rho_{M}}$, and ${\rho_{m}}$), the response was clearly
asymmetric (the response towards more robustness was slower). Although the
average response supports a clear evolutionary trend, response to selection
was variable across simulation replicates, as distributions of up and down
responses generally overlap. The selection response was still ongoing after
10,000 generations.

Selection on robustness components also lead to an indirect response of all
other components, which confirms a general genetic correlation. The magnitude
of the correlated response (from 10% to 100% of the direct response) depended
on the correlation across robustness components. Simulations were run to test
the long-term effect of synergistic and antagonistic selection on all pairs of
robustness indicators (Figure 4), and selection responses were compared to the
mutational evolvabilities computed at the beginning of the simulations. There
was a convincing proportional relationship between predicted and observed
evolvabilities on all directions of selection. Selection response was fast in
directions that were mutationally evolvable, and slow in directions that were
not evolvable. Yet, in spite of the variation of evolvability across
directions in the multivariate robustness space, evolution was always
possible, even if reduced proportionally to the mutational variance,
confirming the absence of absolute constraints.

The proportionality between realized and predicted evolvabilities tends to
fade out for long-term selection responses (Appendix 8), which can be due to
the evolution of mutational constraints (the $\bm{\mathrm{M}}$ matrix evolves
compared to the initial network). This was confirmed by tracking the evolution
of mutational correlations across robustness traits through time (Figure 5).
Average correlations did not evolve substantially in control simulations, but
direct selection on robustness components did trigger systematic change in
some (but not all) mutational correlations. For instance, the correlation
between ${\rho_{M}}$ and ${\rho_{m}}$ does not seem to be evolvable, while the
correlation between ${\rho_{M}}$ and ${\rho_{E}}$ changed from $\simeq 0.3$ to
about $0.6$ or $0.15$ depending on the selection regime. All correlations
remained positive. The evolution of correlations was partially driven by the
direction of selection (more or less robustness). Within each specific pair of
robustness components, the evolution of correlation was rather consistent: for
instance, selecting to decrease ${\rho_{E}}$ or ${\rho_{M}}$ (i.e. making the
network more robust) always decreased the correlation between ${\rho_{E}}$ and
${\rho_{M}}$. Yet, there was no general pattern associating the evolution of
robustness and the evolution of correlation; depending on the robustness
component, selecting for more or less robust networks may increase or decrease
the correlations (colored inset in Figure 5). There was no effect of joint
selection; selecting together two robustness components did not make them more
(or less) correlated (Appendix 9).

Figure 3:  The evolution of all five robustness indicators under direct and
indirect selection was recorded for 10,000 generations in individual-based
simulations. The figures show the average and standard deviation of robustness
over 100 replicated simulations. The control simulations (black circles)
correspond to stabilizing selection on gene expression only (no direct
selection on robustness). Colored symbols correspond to simulations in which
robustness indicators were selected up or down (upward or downward triangles),
colors correspond to the observed indicator, columns indicate which indicator
was selected (diagonal panels: direct selection, off-diagonal panels: indirect
selection). Figure 4:  Trajectories of the bivariate response to selection
over 5000 generations (average over 100 simulation replicates) for all
combinations of robustness indicators. Each panel displays the selection
response in eight directions, as illustrated in the legend (four univariate —
colored arrows — and four bivariate — gray arrows — gradients of selection,
same color code as in Figure 3). Mutational and conditional mutational
matrices, estimated from the initial genotypes, are illustrated as ellipses in
each panel (95% ellipses assuming a multivariate Gaussian mutational
distribution). For conditional $\bm{\mathrm{M}}_{c}$ matrices, the
constraining traits were the three gene which expression was under stabilizing
selection. X and Y axes were adjusted so that their scale matches for each
trait comparison (correlational ellipses were not distorted). The colored
inset illustrates the proportionality between the predicted mutational
evolvability (calculated from $\bm{\mathrm{M}}_{c}$) and the observed
evolvability in the direction of selection after 1000 generations (same
color/symbol code as in the rest of the figure, hyphenated line: linear
regression with no intercept).  Figure 5:  Evolution of the average mutational
correlation under univariate selection on robustness components. Panels,
symbols, and colors are the same as in Figure 4. Correlations were estimated
for each $\bm{\mathrm{M}}$ matrix and averaged over 100 simulation replicates.
The color insets shows the (lack of) consistency between the evolution of
correlations when selecting for higher (x axis) and lower (y axis) robustness.
$\Delta r$ stands for the difference between correlation at generation 10,000
and at generation 1. In the inset, the color of the symbol filling corresponds
to the selected robustness component, the color of the symbol border to the
correlated robutsness component.

## Discussion

Whether or not various robustness components of genetic architectures are
independent is central to understand why organisms are robust or sensitive to
genetic or environmental disturbances. Independent genetic bases of robustness
components would call for independent evolutionary histories, while a
pleiotropic genetic architecture could explain the evolution of nonadaptive
robustness components as a result of indirect selection. The analysis of the
genetic correlations between five robustness components, based on a simple
gene network model, results in a balanced answer: robustness components are
largely correlated, but pleiotropy is not an absolute constraint, and pairs of
robustness components evolved in divergent directions under direct, artificial
bivariate selection. Such a quantitative answer to the so-called ’congruence’
hypothesis [8] would explain both how unselected robustness components could
be partly driven by indirect selection and why various robustness-related
features seem to have their own evolutionary history.

### Model limits

Gene regulation networks are popular candidates when attempting to model
complex biological processes: they are at least partly built on solid and
realistic principles (transcription factors can enhance or repress the
expression of other genes), gene regulation plays a crucial role in most
biological, physiological, and developmental mechanisms, and even modest size
regulation networks display a wide diversity of behavior, including
homeostasis (stable equilibrium of gene expressions) [60], cyclic dynamics
[38, 1], or amplification of a weak signal [26]. Conveniently, the phenotypic
level considered as the output of a gene network (the expression level of all
network genes) can be assimilated to a partial transcriptome, which opens the
possibility for confrontation with empirical data.

The gene network model proposed by [64] is particularly popular in
evolutionary biology to model gene network evolution due to its computational
simplicity and efficiency, combined with a direct biological interpretation
(each line of the regulation matrix is the set of transcription factor
fixation sites in the promoter of a gene) (see [57, 18] for review and
alternative models). In practice, multiple variants based on this original
model have been derived, either to address specific questions, or to correct
for unrealistic features. Here, I used a quantitative version of the model, in
which gene expressions were scaled between 0 (no expression) and 1 (maximum
expression), which was first proposed in [64], although later work have often
preferred binary networks (in which genes can be on/off, e.g. [65, 7]), and a
gene expression scaling between -1 and 1. Unlike in [65, 55], mutations had
cumulative effects (the value of the mutant allele was drawn in a Gaussian
centered around the value of the parental allele), which allows for gradual
evolution. Finally, the sigmoid response function was made asymmetrical by
introducing a constitutive expression parameter (as in e.g. [53]) in order to
avoid the unrealistically high expression of unregulated genes (half the
maximum expression) from the default setting. This constitutive expression was
not evolvable in the model, but simulations (Appendix 7) show that two
robustness components (${\rho_{E}}$ and ${\rho_{S}}$) were very sensitive to
this parameter (larger constitutive expression was associated with more robust
networks). It is thus not unlikely that real systems may evolve towards more
robustness by increasing the constitutive expression of key genes, as already
suggested (for different reasons) by [11].

Discrete time and simple matrix algebra made it possible to run evolutionary
individual-based computer simulations, in which the network output needs to be
calculated for thousands of individuals and thousands of generations. Using
more realistic models based on continuous time and differential equations,
non-linear regulation effects, and independent degradation and transcription
rates would make the simulations less practical, with little benefit in terms
of explanatory power. Computational constraints also limit the network size to
a few dozen genes, which was not enough to generate realistic levels of
sparsity — simulated gene networks were too dense to be realistic. Decreasing
network density and smaller network sizes made robustness components slightly
less correlated (Figure 1E and F), suggesting that the integration of
robustness components increases with network complexity (size and number of
connections). The simulated phenotypic target (maintaining a constant set of
gene expressions) were also extremely simple compared to what gene networks
are theoretically able to do (e.g. converging to different equilibria in
different cell types, or controlling a complex dynamic of gene expression
during the development). However, the results are robust to most simulation
parameters (Appendix 7), suggesting that they reflect general properties of
the underlying genetic architecture. In particular, the network size $n$ and
the number of selected genes $n^{\prime}$ do not alter drastically robustness
components, showing that small regulatory motifs are not qualitatively
different from large gene networks in terms of robustness.

In spite of the simplicity of the network model, it appeared that connecting
network features (for instance, the strength of a specific regulation) and
robustness was not trivial, even in very small networks. For instance, in the
$n=2$ gene-network analysis, most robustness components were complex functions
of all four regulation strengths. Throughout this work, robustness was thus
treated as an emergent property of the underlying network, which cannot be
easily deduced from a reductionist approach. Yet, it is possible to interpret
the correlation patterns in terms of network dynamics. Two of the most
correlated components are the robustness to early environmental variation
${\rho_{E}}$ and network stability ${\rho_{S}}$, which both measure the
ability of the network to converge to a given gene expression equilibrium.
Conversely, the correlation between late mutational ${\rho_{m}}$ and
environmental ${\rho_{e}}$ robustnesses can be attributed to the consequences
of such disturbances over a single time step: for a single target gene,
decreasing the concentration of a transcription factor and decreasing the
sensitivity of the promoter to the same transcription factor have very similar
immediate consequences on gene expression. Yet, even if these measurements
happen to be correlated by construction, their partial evolutionary
independence highlights their potential for independent evolvability in real
gene network architectures, which are substantially more complex and subtle
than our gene network model.

In the simulations, selection on robustness components was direct and constant
both in up and down directions (i.e. towards more or less robust genetic
architectures). This setting was not expected to reflect realistic
evolutionary pressures on robustness, which might be more complex,
overlapping, and asymmetric. Stabilizing selection, for instance, selects both
directly for robustness to environment, and indirectly for robustness to
mutations [68]; selection for stability also promotes indirectly robustness to
mutations [55]. Conversely, selecting for lower robustness through the
phenotype may be difficult or even impossible: fluctuating selection does not
promote decanalized genetic architectures [35], and selection for
environmental sensitivity is limited by the inaccuracy of the perception of
the envrionmental signal [52]. Simulation results thus illustrate how
robustness components may evolve independently when individually selected;
whether or not there exists realistic conditions for such selection pressures
is a different — and more complicated — issue.

### Measuring robustness

There are potentially many ways to measure the robustness of a phenotypic
trait. Here, five indicators were proposed to capture various (and potentially
independent) aspects of what is generally defined as robustness. The
sensitivity to inherited mutations (${\rho_{M}}$) is probably the most popular
one, as it is central to the discussion around the evolution of canalization
[63, 65, 16]. The sensitivity to environmental perturbations is also
unavoidable, although its implementation in a gene network model is less
straightforward. Here, it was calculated as both the sensitivity of the
network to disturbance in the initial expression state (${\rho_{E}}$), which
measures the size of the basin of attraction of the optimal expression
pattern, and as the strength of the stability of the equilibrium when
disturbed (${\rho_{e}}$). These two measurements can be interpreted as
developmental robustness and physiological homeostasis, respectively, as they
quantify the response of the network to disturbances in the expression levels
at different time scales. The robustness to mutations occurring after the
network convergence (${\rho_{m}}$) was considered because it sets up an
alternative to the genetic vs. environmental congruence hypothesis: in long-
lived organisms, non-heritable (somatic) mutations participate to the ageing
process [31], ageing being to some extent under direct selection. Thus, the
robustness to somatic mutations could also drive indirectly the evolution of
genetic canalization. Although not strictly a robustness component, the gene
network stability (${\rho_{S}}$, amplitude of the fluctuations of gene
expressions) was also considered because it has been proven to drive an
indirect response of genetic canalization, based on very similar model [55].
Its correlation with other robustness indicators confirms the tight link
between robustness and stability in gene networks.

These indicators were chosen based on the possibility to measure them in
numerical simulations. Although the empirical assessment of the correlation
between robustness components would be way more convincing than a theoretical
study, defining similar measurements from experimental datasets can be
challenging. For instance, ${\rho_{M}}$ and ${\rho_{E}}$ could, at least in
theory, be estimated as the variance in gene expression across genetic
backgrounds or across environmental conditions, respectively. Measuring
${\rho_{m}}$ environmentally is more complicated, as it would likely be
confounded with other ageing mechanisms. In contrast, the empirical
distinction between e.g. ${\rho_{e}}$ and ${\rho_{S}}$ relies on
discriminating internal vs. external sources of noise, and might be in
practice impossible. In all cases, gene expression data are generally quite
noisy and their analysis necessitates heavy corrections to prevent multiple
testing issues. Studying empirically the robustness and evolvability of
molecular and morphological traits has long been considered as a challenging
task, but methodological and technological progress has recently brought new
concrete perspectives [47].

Some popular measurements of developmental robustness were not considered here
for technical reasons. For instance, fluctuating asymmetry (the variance
between the same phenotypic trait measured in the right and the left body
parts of symmetric organisms) is a convenient measurement of
microenvironmental effects on the development [9, 36], but it has no
equivalent at the level of gene expression in a regulation network. The
deterministic sensitivity to a directional environmental gradient could also
be used to measure phenotypic plasticity, which is central to the question of
phenotypic robustness. Yet, there are several ways to model phenotypic
plasticity in a gene network [40, 4], and it requires a specific selection
setup (different expression optima as a function of the environment). Because
of this additional complexity, adaptive phenotypic plasticity was excluded
from the focus of this work, although the evolution of plasticity of gene
expression remains an intriguing and fundamental question. In particular,
phenotypic plasticity (i.e. an adaptive lack of robustness to some
environmental signal) may itself be canalized to genetic or other
environmental disturbances [59]; considering reaction norms (a measurement of
plasticity) as quantitative traits thus opens challenging questions about the
adaptive evolution of the canalization of robustness traits.

## Data accessibility

Scripts to reproduce simulations and figures are available online
(https://doi.org/10.5281/zenodo.6393075), [34].

## Acknowledgements

Version 5 of this preprint has been peer-reviewed and recommended by Peer
Community In Evolutionary Biology
(https://doi.org/10.24072/pci.evolbiol.100138). I thank the recommender and
the three reviewers for their constructing comments which have substantially
improved the manuscript. Many thanks to Laurent Loison for insightful
discussions. Simulations were partly performed on the Core Cluster of the
Institut Français de Bioinformatique (IFB) (ANR-11-INBS-0013) .

## Funding

No specific funding beyond the CNRS basic support to the author.

## Conflict of interest disclosure

The author of this preprint declare that they have no financial conflict of
interest related to the content of this article.

## References

  * [1] Ozgur E Akman, David A Rand, Paul E Brown and Andrew J Millar  “Robustness from flexibility in the fungal circadian clock”  In _BMC Systems Biology_ 4.1 Springer, 2010, pp. 88  DOI: 10.1186/1752-0509-4-88
  * [2] Ricardo BR Azevedo et al.  “Sexual reproduction selects for robustness and negative epistasis in artificial gene networks”  In _Nature_ 440.7080 Nature Publishing Group, 2006, pp. 87–90  DOI: 10.1038/nature05189
  * [3] Aviv Bergman and Mark L Siegal  “Evolutionary capacitance as a general feature of complex gene networks”  In _Nature_ 424.6948 Nature Publishing Group, 2003, pp. 549–552  DOI: 10.1038/nature01765
  * [4] Ewen Burban, Maud Irene Tenaillon and Arnaud Le Rouzic  “Gene network simulations provide testable predictions for the molecular domestication syndrome”  In _bioRxiv_ Cold Spring Harbor Laboratory, 2021  DOI: 10.1093/genetics/iyab214
  * [5] Reinhard Bürger and Russell Lande  “On the distribution of the mean and variance of a quantitative trait under mutation-selection-drift balance.”  In _Genetics_ 138.3 Oxford University Press, 1994, pp. 901–912  DOI: 10.1093/genetics/138.3.901
  * [6] Mauricio O Carneiro, Clifford H Taubes and Daniel L Hartl  “Model transcriptional networks with continuously varying expression levels”  In _BMC Evolutionary Biology_ 11.1 BioMed Central, 2011, pp. 363  DOI: 10.1186/1471-2148-11-363
  * [7] Stefano Ciliberti, Olivier C Martin and Andreas Wagner  “Innovation and robustness in complex regulatory gene networks”  In _Proceedings of the National Academy of Sciences_ 104.34 National Acad Sciences, 2007, pp. 13591–13596  DOI: 10.1073/pnas.0705396104
  * [8] J Arjan GM de Visser et al.  “Perspective: evolution and detection of genetic robustness”  In _Evolution_ 57.9 Wiley Online Library, 2003, pp. 1959–1972  DOI: 10.1111/j.0014-3820.2003.tb00377.x
  * [9] Vincent Debat and Patrice David  “Mapping phenotypes: canalization, plasticity and developmental stability”  In _Trends in Ecology & Evolution_ 16.10 Elsevier, 2001, pp. 555–561  DOI: 10.1016/S0169-5347(01)02266-2
  * [10] Jeremy Draghi and Gunter P Wagner  “The evolutionary dynamics of evolvability in a gene network model”  In _Journal of Evolutionary Biology_ 22.3 Wiley Online Library, 2009, pp. 599–611  DOI: 10.1111/j.1420-9101.2008.01663.x
  * [11] Jeremy Draghi and Michael Whitlock  “Robustness to noise in gene expression evolves despite epistatic constraints in a model of gene networks”  In _Evolution_ 69.9 Wiley Online Library, 2015, pp. 2345–2358  DOI: 10.1111/evo.12732
  * [12] Dirk Eddelbuettel and James Joseph Balamuta  “Extending extitR with extitC++: A Brief Introduction to extitRcpp”  In _PeerJ Preprints_ 5, 2017, pp. e3188v1  DOI: 10.7287/peerj.preprints.3188v1
  * [13] Ilan Eshel and Carlo Matessi  “Canalization, genetic assimilation and preadaptation: a quantitative genetic model”  In _Genetics_ 149.4 Genetics Soc America, 1998, pp. 2119–2133  DOI: 10.1093/genetics/149.4.2119
  * [14] Carlos Espinosa-Soto, Olivier C Martin and Andreas Wagner  “Phenotypic plasticity can facilitate adaptive evolution in gene regulatory circuits”  In _BMC Evolutionary Biology_ 11.1 BioMed Central Ltd, 2011, pp. 5  DOI: 10.1186/1471-2148-11-5
  * [15] Carlos Espinoza-Soto, Olivier C Martin and Andreas Wagner  “Phenotypic robustness can increase phenotypic variability after nongenetic perturbations in gene regulatory circuits”  In _Journal of Evolutionary Biology_ 24.6 Wiley Online Library, 2011, pp. 1284–1297  DOI: 10.1111/j.1420-9101.2011.02261.x
  * [16] Mario A Fares  “The origins of mutational robustness”  In _Trends in Genetics_ 31.7 Elsevier, 2015, pp. 373–381  DOI: 10.1016/j.tig.2015.04.008
  * [17] Marie-Anne Félix and Michalis Barkoulas  “Pervasive robustness in biological systems”  In _Nature Reviews Genetics_ 16.8 Nature Publishing Group, 2015, pp. 483–496  DOI: 10.1038/nrg3949
  * [18] Janna L Fierst and Patrick C Phillips  “Modeling the evolution of complex genetic systems: The gene network family tree”  In _Journal of Experimental Zoology Part B: Molecular and Developmental Evolution_ 324.1 Wiley Online Library, 2015, pp. 1–12  DOI: 10.1002/jez.b.22597
  * [19] Thomas Flatt  “The evolutionary genetics of canalization”  In _The Quarterly Review of Biology_ 80.3 The University of Chicago Press, 2005, pp. 287–316  DOI: 10.1086/432265
  * [20] Walter Fontana  “Modelling ‘evo-devo’with RNA”  In _BioEssays_ 24.12 Wiley Online Library, 2002, pp. 1164–1177  DOI: 10.1002/bies.10190
  * [21] Christophe Guyeux et al.  “Theoretical study of the one self-regulating gene in the modified wagner model”  In _Mathematics_ 6.4 Multidisciplinary Digital Publishing Institute, 2018, pp. 58  DOI: 10.3390/math6040058
  * [22] Benedikt Hallgrimsson et al.  “The developmental-genetics of canalization”  In _Seminars in Cell & Developmental Biology_ 88, 2019, pp. 67–79  Elsevier  DOI: 10.1016/j.semcdb.2018.05.019
  * [23] Thomas F Hansen and David Houle  “Measuring and comparing evolvability and constraint in multivariate characters”  In _Journal of Evolutionary Biology_ 21.5 Wiley Online Library, 2008, pp. 1201–1219  DOI: 10.1111/j.1420-9101.2008.01573.x
  * [24] Thomas F Hansen and Günter P Wagner  “Modeling genetic architecture: a multilinear theory of gene interaction”  In _Theoretical population biology_ 59.1 Elsevier, 2001, pp. 61–86  DOI: 10.1006/tpbi.2000.1508
  * [25] Joachim Hermisson, Thomas F Hansen and Günter P Wagner  “Epistasis in polygenic traits and the evolution of genetic architecture under stabilizing selection”  In _The American Naturalist_ 161.5 The University of Chicago Press, 2003, pp. 708–734  DOI: 10.1086/374204
  * [26] Gil Hornung and Naama Barkai  “Noise propagation and signaling sensitivity in biological networks: a role for positive feedback”  In _PLoS Computational Biology_ 4.1 Public Library of Science, 2008, pp. e8  DOI: 10.1371/journal.pcbi.0040008
  * [27] Kunihiko Kaneko  “Evolution of robustness to noise and mutation in gene expression dynamics”  In _PLoS one_ 2.5 Public Library of Science San Francisco, USA, 2007, pp. e434  DOI: 10.1371/journal.pone.0000434
  * [28] Stuart Kauffman  “Homeostasis and differentiation in random genetic control networks”  In _Nature_ 224.5215 Springer, 1969, pp. 177–178  DOI: 10.1038/224177a0
  * [29] Stuart Kauffman and Simon Levin  “Towards a general theory of adaptive walks on rugged landscapes”  In _Journal of theoretical Biology_ 128.1 Elsevier, 1987, pp. 11–45  DOI: 10.1016/S0022-5193(87)80029-2
  * [30] T.. Kawecki  “The evolution of genetic canalization under fluctuating selection.”  In _Evolution_ 54.1, 2000, pp. 1–12  DOI: 10.1111/j.0014-3820.2000.tb00001.x
  * [31] Scott R Kennedy, Lawrence A Loeb and Alan J Herr  “Somatic mutations in aging, cancer and neurodegeneration”  In _Mechanisms of Ageing and Development_ 133.4 Elsevier, 2012, pp. 118–126  DOI: 10.1016/j.mad.2011.10.009
  * [32] Christian Peter Klingenberg  “Phenotypic plasticity, developmental instability, and robustness: The concepts and how they are connected”  In _Frontiers in Ecology and Evolution_ 7 Frontiers, 2019, pp. 56  DOI: 10.3389/fevo.2019.00056
  * [33] Yann Le Cunff and Khashayar Pakdaman  “Phenotype-genotype relation in Wagner’s canalization model.”  In _Journal of Theoretical Biology_ 314, 2012, pp. 69–83  DOI: 10.1016/j.jtbi.2012.08.020
  * [34] A. Le Rouzic  “Github Repository for the manuscript "Gene network robustness as a multivariate character"”, 2022  DOI: 10.5281/zenodo.6393075
  * [35] Arnaud Le Rouzic, José M Álvarez-Castro and Thomas F Hansen  “The evolution of canalization and evolvability in stable and fluctuating environments”  In _Evolutionary Biology_ 40.3 Springer, 2013, pp. 317–340  DOI: 10.1007/s11692-012-9218-z
  * [36] Larry J Leamy and Christian Peter Klingenberg  “The genetics and evolution of fluctuating asymmetry”  In _Annual Reviews in Ecology, Evolution, and Systematics_ 36 Annual Reviews, 2005, pp. 1–21  DOI: 10.1146/annurev.ecolsys.36.102003.152640
  * [37] Ben Lehner  “Genes confer similar robustness to environmental, stochastic, and genetic perturbations in yeast”  In _PloS One_ 5.2 Public Library of Science San Francisco, USA, 2010, pp. e9035  DOI: 10.1371/journal.pone.0009035
  * [38] Jean-Christophe Leloup and Albert Goldbeter  “Toward a detailed computational model for the mammalian circadian clock”  In _Proceedings of the National Academy of Sciences_ 100.12 National Acad Sciences, 2003, pp. 7051–7056  DOI: 10.1073/pnas.1132112100
  * [39] Laurent Loison  “Canalization and genetic assimilation: Reassessing the radicality of the Waddingtonian concept of inheritance of acquired characters”  In _Seminars in Cell & Developmental Biology_ 88, 2019, pp. 4–13  Elsevier  DOI: 10.1016/j.semcdb.2018.05.009
  * [40] J Masel  “Genetic assimilation can occur in the absence of selection for the assimilating phenotype, suggesting a role for the canalization heuristic”  In _Journal of Evolutionary Biology_ 17.5 Wiley Online Library, 2004, pp. 1106–1110  DOI: 10.1111/j.1420-9101.2004.00739.x
  * [41] Joanna Masel and Mark L Siegal  “Robustness: mechanisms and consequences”  In _Trends in Genetics_ 25.9 Elsevier, 2009, pp. 395–403  DOI: 10.1016/j.tig.2009.07.005
  * [42] Lisandro Milocco and Isaac Salazar-Ciudad  “Is evolution predictable? Quantitative genetics under complex genotype-phenotype maps”  In _Evolution_ 74.2 Wiley Online Library, 2020, pp. 230–244  DOI: 10.1111/evo.13907
  * [43] Shintaro Nagata and Macoto Kikuchi  “Emergence of cooperative bistability and robustness of gene regulatory networks”  In _PLoS Computational Biology_ 16.6 Public Library of Science San Francisco, CA USA, 2020, pp. e1007969  DOI: 10.1371/journal.pcbi.1007969
  * [44] H Frederik Nijhout  “The nature of robustness in development”  In _Bioessays_ 24.6 Wiley Online Library, 2002, pp. 553–563  DOI: 10.1002/bies.10093
  * [45] H Frederik Nijhout, Janet A Best and Michael C Reed  “Systems biology of robustness and homeostatic mechanisms”  In _Wiley Interdisciplinary Reviews: Systems Biology and Medicine_ 11.3 Wiley Online Library, 2019, pp. e1440  DOI: 10.1002/wsbm.1440
  * [46] Andreas Odorico, Estelle Rünneburger and Arnaud Le Rouzic  “Modelling the influence of parental effects on gene-network evolution”  In _Journal of Evolutionary Biology_ 31.5 Wiley Online Library, 2018, pp. 687–700  DOI: 10.1111/jeb.13255
  * [47] Joshua L Payne and Andreas Wagner  “The causes of evolvability and their evolution”  In _Nature Reviews Genetics_ 20.1 Nature Publishing Group, 2019, pp. 24–38  DOI: 10.1038/s41576-018-0069-z
  * [48] Ricardo Pinho, Elhanan Borenstein and Marcus W Feldman  “Most networks in Wagner’s model are cycling”  In _PloS One_ 7.4 Public Library of Science, 2012, pp. e34285  DOI: 10.1371/journal.pone.0034285
  * [49] R Core Team  “R: A Language and Environment for Statistical Computing”, 2020  R Foundation for Statistical Computing  URL: https://www.R-project.org/
  * [50] Hannah Rajasingh, Arne B Gjuvsland, Dag Inge Vaåge and Stig W Omholt  “When parameters in dynamic models become phenotypes: a case study on flesh pigmentation in the Chinook salmon (Oncorhynchus tshawytscha)”  In _Genetics_ 179.2 Oxford University Press, 2008, pp. 1113–1118  DOI: 10.1534/genetics.108.087064
  * [51] Etienne Rajon and Joanna Masel  “Compensatory evolution and the origins of innovations”  In _Genetics_ 193.4 Genetics Soc America, 2013, pp. 1209–1220  DOI: 10.1534/genetics.112.148627
  * [52] Thomas E Reed et al.  “Phenotypic plasticity and population viability: the importance of environmental predictability”  In _Proceedings of the Royal Society B: Biological Sciences_ 277.1699 The Royal Society, 2010, pp. 3391–3400  DOI: 10.1098/rspb.2010.0771
  * [53] Estelle Rünneburger and Arnaud Le Rouzic  “Why and how genetic canalization evolves in gene regulatory networks”  In _BMC Evolutionary Biology_ 16.1 BioMed Central, 2016, pp. 1–11  DOI: 10.1186/s12862-016-0801-2
  * [54] Ivan Ivanovich Schmalhausen  “Factors of evolution: the theory of stabilizing selection.”  Blakiston, 1949
  * [55] Mark L Siegal and Aviv Bergman  “Waddington’s canalization revisited: developmental stability and evolution”  In _Proceedings of the National Academy of Sciences_ 99.16 National Acad Sciences, 2002, pp. 10528–10532  DOI: 10.1073/pnas.102303999
  * [56] Paul Smolen, Douglas A Baxter and John H Byrne  “Modeling transcriptional control in gene networks —methods, recent results, and future directions”  In _Bulletin of Mathematical Biology_ 62.2 Springer, 2000, pp. 247–292  DOI: 10.1006/bulm.1999.0155
  * [57] Alexander Spirov and David Holloway  “Using evolutionary computations to understand the design and evolution of gene and cell regulatory networks”  In _Methods_ 62.1 Elsevier, 2013, pp. 39–55  DOI: 10.1016/j.ymeth.2013.05.013
  * [58] Stephen C Stearns  “Progress on canalization”  In _Proceedings of the National Academy of Sciences_ 99.16 National Acad Sciences, 2002, pp. 10229–10230  DOI: 10.1073/pnas.172388999
  * [59] Stephen C Stearns and Tadeusz J Kawecki  “Fitness sensitivity and the canalization of life-history traits”  In _Evolution_ 48.5 Wiley Online Library, 1994, pp. 1438–1450  DOI: 10.1111/j.1558-5646.1994.tb02186.x
  * [60] Michael D Stern  “Emergence of homeostasis and “noise imprinting” in an evolution model”  In _Proceedings of the National Academy of Sciences_ 96.19 National Acad Sciences, 1999, pp. 10746–10751  DOI: 10.1073/pnas.96.19.10746
  * [61] Stephen J Tonsor, Tarek W Elnaccash and Samuel M Scheiner  “Developmental instability is genetically correlated with phenotypic plasticity, constraining heritability, and fitness”  In _Evolution_ 67.10 Wiley Online Library, 2013, pp. 2923–2935  DOI: 10.1111/evo.12175
  * [62] Conrad H Waddington  “Canalization of development and the inheritance of acquired characters”  In _Nature_ 150.3811 Nature Publishing Group, 1942, pp. 563–565  DOI: 10.1038/150563a0
  * [63] Conrad H Waddington  “Canalization of development and genetic assimilation of acquired characters”  In _Nature_ 183.4676 Nature Publishing Group, 1959, pp. 1654–1655  DOI: 10.1038/1831654a0
  * [64] Andreas Wagner  “Evolution of gene networks by gene duplications: a mathematical model and its implications on genome organization”  In _Proceedings of the National Academy of Sciences_ 91.10 National Acad Sciences, 1994, pp. 4387–4391  DOI: 10.1073/pnas.91.10.4387
  * [65] Andreas Wagner  “Does evolutionary plasticity evolve?”  In _Evolution_ 50.3 Wiley Online Library, 1996, pp. 1008–1023  DOI: 10.1111/j.1558-5646.1996.tb02342.x
  * [66] Andreas Wagner  “Robustness and evolvability in living systems”  Princeton university press, 2013  DOI: 10.1515/9781400849383
  * [67] Andreas Wagner and Peter F Stadler  “Viral RNA and evolved mutational robustness”  In _Journal of Experimental Zoology_ 285.2 Wiley Online Library, 1999, pp. 119–127  DOI: 10.1002/(SICI)1097-010X(19990815)285:2%3C119::AID-JEZ4%3E3.0.CO;2-D
  * [68] Günter P Wagner, Ginger Booth and Homayoun Bagheri-Chaichian  “A population genetic theory of canalization”  In _Evolution_ 51.2 Wiley Online Library, 1997, pp. 329–347  DOI: 10.1111/j.1558-5646.1997.tb02420.x

## Appendix

## Appendix 1

### Sensitivity of the robustness measurements to the magnitude of the
disturbance

Four out of five robustness indicators (${\rho_{E}}$, ${\rho_{e}}$,
${\rho_{M}}$, ${\rho_{m}}$) depend on the magnitude of the disturbance
($\sigma_{E}$, $\sigma_{e}$, $\sigma_{M}$, and $\sigma_{m}$, respectively).
The figure displays the influence of the size of the disturbance on the
robustness measurement (left: 10 random networks, right: 10 evolved networks).
Vertical dotted lines stand for the values used in the simulations. Robustness
scores are not completely consistent for random networks, as some of them can
be differentially robust to large or small disturbances. The consistency is
better in evolved networks (the rank of different genotypes in terms of
robustness rarely depends on the size of the disturbance).

## Appendix 2

### Correlations among robustness indexes among random networks

Correlations between all five robustness components among 10,000 random 6-gene
networks ($\mu_{0}=0,\sigma_{0}=1$).

## Appendix 3

### Sampling effects on Principal Components

Influence of the sampling effect (number of networks and number of replicates
$R$ to estimate robustness) on the relative weight of the principal
components. All PCs except the last one are robust to sampling.

## Appendix 4

### Reasons for not reaching the desired equilibrium

Although equation 4 guarantees that an equilibrium exists at the target
phenotypic expression, the equilibrium might not be reachable in practice when
simulating the gene network dynamics. The colored area in the figure
corresponds to networks that failed to produced the target phenotype, each
color representing a distinct reason; Yellow: network dynamics was slow and
the final gene expression has not been reached yet after 16 time steps; Gray:
an alternative equilibrium was reached (most of the time implying that one or
both genes are either completely silenced to fully expressed). Red: The
network steady state featured oscillations that were so large that they hit
the maximum or minimum expression, shifting the average expression away from
the target expression.

## Appendix 5

### Two-gene example networks

| $W_{11}$ | $W_{21}$ | $W_{12}$ | $W_{22}$
---|---|---|---|---
A | 0.70 | 0.20 | -0.21 | 0.38
B | -0.30 | 0.30 | 0.29 | 0.33
C | -0.40 | 0.80 | 0.34 | 0.08
D | -1.00 | -0.80 | 0.64 | 0.88
E | 1.50 | 3.50 | -0.61 | -1.27

The five two-gene networks detailed in Figure 2 and Appendix 6.

## Appendix 6

### Illustration of the robustness scores

The figure displays a subset of the replicated tests for four robustness
indexes. Rows A to E correspond to the five networks described in Appendix 5.
Four (out of five) robustness measurements rely on a resampling procedure
(corresponding to the four columns of the figure). In each panel, the default
(undisturbed) network kinetics is displayed as plain lines (black for gene 1,
red for gene 2), while 10 disturbed networks are indicated as pale lines. By
construction, all networks have an equilibrium at (0.3, 0.6). The network
stability can be assessed from the amplitude of the cycles in the undisturbed
kinetics (thick lines), and does not rely on a stochastic algorithm. The
network robustness to genetic disturbance was estimated by mutating the gene
network before the first time step (early genetic mutation, first column) or
before the last time step (late genetic mutation, second column).
Environmental robustness was estimated by disturbing the gene expression,
without changing the genotype, before the first time step (early
environmental, third column) and before the last time step (late
environmental, fourth column).

## Appendix 7

### Exploration of the parameter set

Influence of simulation parameters (mutation rate $\nu$, mutation size
$\sigma_{\nu}$, population size $N$, constitutive expression $a$, total number
of genes $n$, number of selected genes $n^{\prime}$, network density $d$, and
strenght of selection $s$) on fitness and robustness indexes after 5000
generations (default settings except for the target parameter). The figure
reports the mean $\pm$ standard deviation across 20 replicated simulations.
Vertical dotted lines stand for the default parameter values.

## Appendix 8

### Accuracy of the prediction vs. simulation time

Effect of the number of generations on the proportionality relationship
between predicted and observed evolvabilities of robustness components. The
figure displays the $r^{2}$ of a linear regression (without intercept) between
the predicted evolvability from the conditional $\bm{\mathrm{M}}_{c}$ matrix
measured at the first generation and the observed evolvability in the
direction of selection for all replicated simulations. The regression at
generation 1,000 is illustrated in the colored inset in Figure 4.

## Appendix 9

### Evolution of correlations

Evolution of the mutational correlation among robustness components after
10,000 generations of evolution ($\Delta r=r_{10,0000}-r_{0}$), averaged over
100 simulation replicates. For each pair of robustness components, nine
selection gradients were simulated (including control simulations without
selection on robustness, central slot).