# Double-replica theory for evolution of genotype-phenotype interrelationship

Tuan Minh Pham  The Niels Bohr Institute, University of Copenhagen,
Blegdamsvej 17, Copenhagen, 2100-DK, Denmark  Kunihiko Kaneko  The Niels Bohr
Institute, University of Copenhagen, Blegdamsvej 17, Copenhagen, 2100-DK,
Denmark  Center for Complex Systems Biology, Universal Biology Institute,
University of Tokyo, Komaba, Tokyo 153-8902, Japan

###### Abstract

The relationship between genotype and phenotype plays a crucial role in
determining the function and robustness of biological systems. Here the
evolution progresses through the change in genotype, whereas the selection is
based on the phenotype, and genotype-phenotype relation also evolves. Theory
for such phenotypic evolution remains poorly-developed, in contrast to
evolution under the fitness landscape determined by genotypes. Here we provide
statistical-physics formulation of this problem by introducing replicas for
genotype and phenotype. We apply it to an evolution model, in which phenotypes
are given by spin configurations; genotypes are interaction matrix for spins
to give the Hamiltonian, and the fitness depends only on the configuration of
a subset of spins called target. We describe the interplay between the genetic
variations and phenotypic variances by noise in this model by our new approach
that extends the replica theory for spin-glasses to include spin-replica for
phenotypes and coupling-replica for genotypes. Within this framework we obtain
a phase diagram of the evolved phenotypes against the noise and selection
pressure, where each phase is distinguished by the fitness and overlaps for
genotypes and phenotypes. Among the phases, robust fitted phase, relevant to
biological evolution, is achieved under the intermediate level of noise
(temperature), where robustness to noise and to genetic mutation are
correlated, as a result of replica symmetry. We also find a trade-off between
maintaining a high fitness level of phenotype and acquiring a robust pattern
of genes as well as the dependence of this trade-off on the ratio between the
size of the functional (target) part to that of the remaining non-functional
(non-target) one. The selection pressure needed to achieve high fitness
increases with the fraction of target spins.

##  I Introduction

Over decades, evolution under given fitness landscape, determined as a
function of genes (genotypes), has been studied extensively [1, 2]. Here
genotypes are changed (mutated) in the reproduction process, and those that
produce higher function, called phenotypes, are selected. These phenotypes
determine fitness, the rate of offsprings that survive. However, the evolution
of phenotypes whose configurations are shaped by the genetic evolution remains
poorly explored. Here, phenotypes are a result of dynamics whose rule is
determined by genotypes. Such dynamics are stochastic in general. Cells
involve stochastic gene expression dynamics, whereas protein folding dynamics
to give the protein shape is under thermal noise [3, 4, 5, 6]. Phenotypes
hence are variable under noise, while fitted phenotypic states are better to
be preserved under noise, i.e., to keep robustness to noise. In addition, they
can also be varied by genetic mutation during the evolution of genotypes, and
the robustness to mutation will also be required. The achievement of
robustness of phenotypes to noise and to mutation is important to the
evolution, as has been discussed recently [7, 8, 9, 10]. Now considering
stochastic dynamics of phenotypes, a general formulation of the evolution of
such genotype-phenotype mapping and phenotypic robustness is hence wanted.

Underlying such stochastic dynamics are the interactions among a vast number
of elements that constitute a biological system. A cell consists of a huge
variety of interacting molecules and its constitute polymers (proteins) are
composed of many monomers (residues). Now the states of such interacting
elements that shape the evolution of phenotypes can be properly described by
statistical physics [11, 12]. To this kind of study, use of spin models is
relevant, where phenotypes are spin configurations that are updated by
Hamiltonian with spin-spin interactions under thermal noise, whereas genotypes
specify such spin-spin interactions, and fitness is given by a function of
configuration of a subset of spins, termed as target spins. It is then
important to identify possible phases with regards to genotypes and
phenotypes, using the set of order parameters, a concept rooted in statistical
physics.

In the present paper, to address these problems in a systematic way, we
formulate double replica methods both for spins (phenotypes) and couplings
(genotypes), as well as both for target and nontarget parts. Even though we
adopted spin-coupling representation here our formulation can generally be
applied to other problems, in which the interactions among many degrees of
freedom are also dynamical variables with their own dynamics. For the sake of
demonstration, however, we here explain this approach by using specifically a
spin-glass Hamiltonian model developed in [13, 14, 15]. These works uncovered
the transitions between different regimes with regards to the fitness and
robustness upon changing the strength of thermal noise for phenotype and
selection pressure for genotype. In particular, within an intermediate range
of phenotypic noise, the evolved spin-systems could attain high fitness and
robustness to noise and mutation. A systematic way to elucidate the condition
to achieve such robust fitted states with regards to the noise strength,
selection pressure, and relative size of target spins, and to understand
possible relationship between robustness to mutation and to noise, however,
has not been developed yet.

One might expect an application of mean-field methods for disordered systems
[16] in this model since a spin-glass Hamiltonian formulation was adopted by
replacing random couplings among spins (phenotypes) by genotypes that are
evolving [13, 14, 15]. Gradual change in the couplings might fit with partial
annealing approach based on a finite number of replicas $n$. However, the
study is restricted to the case in which the coupling dynamics are affected by
a spin-spin correlation term [17, 18, 19, 20, 21], and is not directly
applicable for our purpose, in which the coupling dynamics depend on the
fitness determined by the spin configurations. Another theoretical method
assumes the quench limit ($n\rightarrow 0$) for a replicated spin system [22].
However, this means that the couplings are treated effectively as ‘static’
(but with a modified distribution), and hence is not suitable to investigate
the evolution of both genotypes (couplings) and phenotypes (spin
configurations). Using this approach, it thus remains elusive to see how
robustness emerges from the interplay between the ordering of spins and that
of the couplings.

In this paper we develop a new mean-field approach that we term as double
replica theory. It describes the evolution of both genotypes and phenotypes by
considering spins and links as two different replica species. With this
formulation, we obtain fitness and replica overlaps for spins and couplings,
which work as the order parameters. Using these order parameters we identify
five regions in the temperature vs selection pressure phase diagram: two non-
fitted paramagnetic phases, fitted- and non-fitted spin glass phases, and a
robust fitted phase. The last phase is the most biologically important, which
can only be achieved under intermediate noise level (temperature) and
sufficient selection pressure, whereas robustness can only arrive at the cost
of lowering the fitness from its maximal value. Dependence of the robust
fitted phase on the ratio of functional to non-functional parts has been
analyzed in depth. As the former ratio is increased, the selection pressure to
achieve this phase is drastically increased, whereas, if achieved, it can
persist for slightly higher noise. This suggests the relevance of having
sufficient non-functional parts in biological systems. In addition,
correlation between robustness to noise and to mutation are formulated as a
proportionality between susceptibilities to external field and to coupling
change.

##  II Double replica theory

Following [14, 15] we study the evolution of the relationship between
phenotype and genotype by representing phenotypes as spin configurations, and
genotypes as interaction matrix for spins. In a system of $N$ spins, each spin
$i$ can take values $s_{i}\in\\{-1,1\\}$ and is linked to exactly $N-1$ other
spins, thus forming a fully-connected network. Here the evolution progresses
through the change in genotype, whereas the selection is based on the
phenotype, resulting in an evolution of the genotype-phenotype relation.
Moreover, fitness is determined by a subset of target spins denoted by
$\mathcal{T}$. Those spins that do not contribute to the fitness are called
non-target. In general, the fitness $\Psi$ is some field that acts on $J_{ij}$
but whose value depends only on $s_{i}$, $i\in\mathcal{T}$. How such
dependence is explicitly described is model-specific and will not limit the
use of our approach. See SM Eq. (14) for an example of $\Psi$ given by the
target spin configurations at equilibrium [14, 15].

Stochastic dynamics of phenotypes are considered as the evolution of spin
configurations at a temperature $T_{s}$ according to a Hamiltonian
$H_{S}=-\sum_{i<j}J_{ij}s_{i}s_{j}$ [23]. Here the couplings $J_{ij}$ are
regarded as _fixed_ over the course of the spin evolution that follows a
Glauber update because they are assumed to evolve on much slower timescale
than that of the spins. Furtheremore, the couplings are symmetric, i.e.
$J_{ij}=J_{ji}$, and, initially, are independently and identically distributed
by a Gaussian distribution with zero mean and the variance $J^{2}:={\rm
var}(J_{ij})=N^{-1}$. The coupling matrix $\mathbb{J}$ includes interactions
between target spins $J_{ij}^{(tt)}$ for $i\in\mathcal{T}$ and
$j\in\mathcal{T}$; those between non-target spins $J_{ij}^{(oo)}$ for
$i\centernot\in\mathcal{T}$ and $j\centernot\in\mathcal{T}$; and those between
target spin and non-target spin $J_{ij}^{(to)}$ for $i\in\mathcal{T}$ and
$j\centernot\in\mathcal{T}$. Let $S_{\mathcal{T}}$ and $S_{\mathcal{O}}$
denote the subsystem of target spins (with their interactions
$\mathbb{J}^{(tt)}$) and the subsystem of non-target spins (with the couplings
$\mathbb{J}^{(oo)}$ among them), respectively. The Hamiltonian of the full
system denoted by $S$ can be decomposed into

$H_{S}=\underbrace{-\sum_{i<j\in\mathcal{T}}J^{(tt)}_{ij}s_{i}s_{j}}_{H_{\mathcal{T}}}\underbrace{-\sum_{i<j\centernot\in\mathcal{T}}J^{(oo)}_{ij}s_{i}s_{j}}_{H_{\mathcal{O}}}\underbrace{-\sum_{\begin{subarray}{c}i\in\mathcal{T}\\\
j\centernot\in\mathcal{T}\end{subarray}}J^{(to)}_{ij}s_{i}s_{j}}_{H_{\mathcal{TO}}},$
(1)

where $H_{\mathcal{T}}$ and $H_{\mathcal{O}}$ are the Hamiltonian of the
subsystems $S_{\mathcal{T}}$ and $S_{\mathcal{O}}$, respectively, while
$H_{\mathcal{TO}}$ describes the interactions between these subsystems.

Now we introduce the effective potential to obtain the distribution of
$\mathbb{J}$ 111Note the difference from the model implementation of [14]
where couplings are updated at discrete time steps.. For it, we consider a
continuous Langevin-type dynamics for the couplings

$\frac{dJ_{ij}}{d\tau}=-\frac{1}{N}\frac{\partial\,V}{\partial
J_{ij}}+\frac{1}{\sqrt{N}}\,\xi_{ij}(\tau)\,,$ (2)

where $V=V\big{(}\mathbb{J}\big{)}$ is the potential of all the couplings and
$\xi_{ij}$ is the white noise whose intensity equal to the temperature
$T_{J}$. The factors $1/N$ and $1/\sqrt{N}$ in front of the potential and the
noise term, respectively, ensure a correct relationship between the drift and
diffusive parts of the Langevin equation. If the couplings were independent
from each other, the potential would simply take the form of the potential of
a free Brownian particle

$V_{0}=\frac{N}{2}\cdot\sum_{i<j}J_{ij}^{2}\,.$ (3)

However, in the presence of fitness we need an additional term. Here we assume
that the fitness would be maximised if a global alignment is established among
target spins, so that we introduce

$\Psi=\frac{1}{N_{t}}\,\Big{|}\sum_{i\in\mathcal{T}}s_{i}\Big{|}\,,$ (4)

where $N_{t}$ is the size of $\mathcal{T}$. Under this fitness that favors the
alignment of target spins, the couplings are necessarily subjected to a
fitness field $K$ 222Technically speaking, the field $K$ is introduced to
break the gauge-symmetry between different local Mattis states that have no
frustration among target spins (i.e. triples
$J_{ij}^{(tt)}J_{jk}^{(tt)}J_{ki}^{(tt)}>0$, $\forall\,i<j<k\in\mathcal{T}$).
As this symmetry is broken, the ferromagnetic state (with all links between
target spins being positive) is distinguished from all other unfrustrated
states and, as a consequence, alignment among target spins can be achieved at
low $T_{s}$. One might want to consider a different form of $K$ that
explicitly takes into account the Boltzmann weights corresponding to the
Hamiltonian $H_{S}$ (i.e. $e^{-\beta_{s}H_{S}}$, as defined in Eq. (4)). Such
a term, however, may not guarantee the evolution towards a subgraph of only
ferromagnetic interactions among target spins. :

$K=\frac{1}{\beta_{J}}\frac{\partial}{\partial
J_{ij}}\ln\Big{(}\sum_{\\{s_{i}\\}}\exp\Big{\\{}\frac{\beta_{J}}{N_{t}}\sum_{i<j\in\mathcal{T}}J_{ij}\Big{|}\sum_{i\in\mathcal{T}}s_{i}\Big{|}\Big{\\}}\Big{)}$
(5)

or equivalently, $V$ needs to be modified from $V_{0}$ into

$\displaystyle
V=V_{0}-\frac{1}{\beta_{J}}\,\ln\Big{(}\sum_{\\{s_{i}\\}}\exp\Big{\\{}\frac{\beta_{J}}{N_{t}}\,\sum_{i<j\in\mathcal{T}}J_{ij}\Big{|}\sum_{i\in\mathcal{T}}s_{i}\Big{|}\Big{\\}}\Big{)}\,.$

Without the evolution of genotypes, the Hamiltonians $H_{\mathcal{T}}$,
$H_{\mathcal{O}}$ and $H_{\mathcal{TO}}$ dictate the spins to adapt to a set
of fixed couplings $\mathbb{J}$ in order to minimise each term of Eq. (1)
through the spin dynamics. Such adaptation results in an accordance between
the state of $J_{ij}^{(tt)}$ and $s_{i}s_{j}$ for $i,j\in\mathcal{T}$; that
between the state of $J_{ij}^{(oo)}$ and $s_{i}s_{j}$ for
$i,j\centernot\in\mathcal{T}$ and that between the state of $J_{ij}^{(to)}$
and $s_{i}s_{j}$ for $i\in\mathcal{T},j\centernot\in\mathcal{T}$. As long as
this kind of accordance exits, it is insufficient to consider the evolving
couplings with selection force given in Eq. (5) only. A link between two spins
hence necessarily needs to adapt to the joint state of these spins. Due to the
time scale separation between the phenotype- and the genotype dynamics, the
direction of change of genotypes is determined by the equilibrium correlation
of the phenotypes, and since $J_{ij}$ is symmetric, it needs to be:

$dJ_{ij}/d\tau\propto\langle s_{i}s_{j}\rangle_{T_{s}}\,.$

This is equivalent to have a potential of the form 333Here we used the
identities $\langle
s_{i}s_{j}\rangle_{T_{s}}=\frac{1}{\beta_{s}}\frac{\partial}{\partial
J_{ij}^{(tt)}}\ln Z_{1}$ and $\langle
s_{k}s_{k^{\prime}}\rangle_{T_{s}}=\frac{1}{\beta_{s}}\frac{\partial}{\partial
J_{kk^{\prime}}^{(oo)}}\ln\tilde{Z}_{1}$, where the partition function of the
target spins’ subsystem and that of the non-target spins’ subsystem are
$Z_{1}:=\sum_{\\{s_{i}\\}_{i\in\mathcal{T}}}e^{-\beta_{s}H_{\mathcal{T}}(\mathbb{J}^{(tt)})}$
and
$\tilde{Z}_{1}:=\sum_{\\{s_{i}\\}_{i\centernot\in\mathcal{T}}}e^{-\beta_{s}H_{\mathcal{O}}(\mathbb{J}^{(oo)})}$,
respectively.

$V_{a}=V-\frac{1}{\beta_{s}}\,\ln\Big{(}\sum_{\\{s_{i}\\}}\exp\Big{\\{}\beta_{s}\,\sum_{i<j}J_{ij}s_{i}s_{j}\Big{\\}}\Big{)}$
(6)

The stochastic process induced by Eq. (2) under this potential admits an
equilibrium-like stationary _joint_ distribution
$\mathbb{P}(\mathbb{J}^{(tt)},\mathbb{J}^{(oo)},\mathbb{J}^{(to)})$ of
Boltzmann-form (with associated temperature $T_{J}$), i.e.,

$\mathbb{P}(\mathbb{J}^{(tt)},\mathbb{J}^{(oo)},\mathbb{J}^{(to)})=e^{-\beta_{J}V_{a}}/\mathcal{Z}_{\rm
total}\,,$

where $\mathcal{Z}_{\rm total}=\sum_{\\{\mathbb{J}\\}}e^{-\beta_{J}V_{a}}$.
Instead of calculating this distribution, we introduce our approximate
approach, in which $J_{ij}^{(to)}$ are assumed to always attain equilibrium
well before $J_{ij}^{(tt)}$ and $J_{ij}^{(oo)}$ and hence can be adiabatically
eliminated. As a consequence, only the weights of equilibrium configurations
of $\mathbb{J}^{(to)}$ contribute to the stationary distributions

$\displaystyle\mathbb{P}_{\mathcal{T}}(\mathbb{J}^{(tt)})$ $\displaystyle=$
$\displaystyle\lim_{\tau\rightarrow\infty}\mathbb{P}_{\mathcal{T}}(\mathbb{J}^{(tt)},\tau)$
$\displaystyle\mathbb{P}_{\mathcal{O}}(\mathbb{J}^{(oo)})$ $\displaystyle=$
$\displaystyle\lim_{\tau\rightarrow\infty}\mathbb{P}_{\mathcal{O}}(\mathbb{J}^{(oo)},\tau)$

where $\mathbb{P}_{\mathcal{T}}(\mathbb{J}^{(tt)},\tau)$ and
$\mathbb{P}_{\mathcal{O}}(\mathbb{J}^{(oo)})$ are the time-dependent solutions
of the corresponding Forrker-Planck equations with the effective potentials
$V_{tt}$ for $\mathbb{J}^{(tt)}$ and $V_{oo}$ for and $\mathbb{J}^{(oo)}$,
respectively 444Requiring that the couplings should remain bounded as
$\tau\rightarrow\infty$, we need to introduce a decay term
$-\lambda_{ij}J_{ij}$ to the couplings dynamics. In order to keep the
genotypic selection pressure the same for all the couplings regardless of
their types, we choose $-\lambda_{ij}J_{ij}^{(tt)}=-\sqrt{p}J_{ij}^{(tt)}$ and
$-\lambda_{ij}J_{ij}^{(oo)}=-\sqrt{1-p}J_{ij}^{(oo)}$, where $p=N_{t}/N$ is
the fraction of target spins. To obtain the distribution only of
$\mathbb{J}^{(tt)}$ and $\mathbb{J}^{(oo)}$, we first replace
$\mathbb{J}^{(to)}$ as given matrix by that obtained self-consistenly from
equilibrium distribution. For it we need to modify $V_{a}$ in such a way that
the influence of $\mathbb{J}^{(to)}$ on $\mathbb{J}^{(tt)}$
($\mathbb{J}^{(oo)}$) can be taken into account in the effective potential
$V_{tt}$ ($V_{oo}$). The _joint_ effect of $H_{\mathcal{TO}}$ and
$H_{\mathcal{T}}$ on the dynamics of target spins suggests that the dependence
of $J_{ij}^{(tt)}$ and $J_{ik}^{(to)}$ on each other arises from any triad
formed between $(i,j)\in\mathcal{T}$ and $k\centernot\in\mathcal{T}$, i.e.,
via $J_{ij}^{(tt)}J_{ik}^{(to)}J_{jk}^{(to)}$. This influence is represented
by the frustration [16] which implies that an optimal spin configuration
$(s_{i}^{*},s_{j}^{*},s_{k}^{*})$ can only be established if the relation
$J_{ij}^{(tt)}J_{ik}^{(to)}J_{jk}^{(to)}>0$ holds (optimality in this context
means that $H_{\mathcal{TO}}$ and $H_{\mathcal{T}}$ can be lowered
simultaneously by $(s_{i}^{*},s_{j}^{*},s_{k}^{*})$). Since for any given pair
of $(i,j)\in\mathcal{T}$ there are $N-N_{t}$ triads $\Delta_{k}$ formed
between it and a non-target spin $k\centernot\in\mathcal{T}$, the total effect
of frustration is given by

$F_{ij}^{(t-o-t)}=(N-N_{t})^{-1}\sum_{k=1}^{N-N_{t}}J_{ik}^{(to)}J_{jk}^{(to)}\,.$
(7)

Likewise, frustration among all $N_{t}$ triads $\tilde{\Delta}_{k}$ formed
between a given pair of non-target spins $i\centernot\in\mathcal{T}$ and
$j\centernot\in\mathcal{T}$ with $k\in\mathcal{T}$ induces a force
$F_{(o-t-o)}$ on the state of $J_{ij}^{(oo)}$:

$F_{ij}^{(o-t-o)}=N_{t}^{-1}\sum_{k=1}^{N_{t}}J_{ik}^{(to)}J_{jk}^{(to)}\,.$
(8)

The proposed scheme just allows us to define the effective potential $V_{tt}$
for $\mathbb{J}^{(tt)}$ and $V_{oo}$ for and $\mathbb{J}^{(oo)}$ as

$\displaystyle V_{tt}$ $\displaystyle=$ $\displaystyle
V_{a}-\frac{1}{\beta_{J}}\,\ln\Big{(}\sum_{\\{J_{ij}^{(to)}\\}}\exp\Big{\\{}\beta_{J}\,\sum_{i<j\in\mathcal{T}}J_{ij}^{(tt)}F_{ij}^{(t-o-t)}\Big{\\}}\Big{)}$
$\displaystyle V_{oo}$ $\displaystyle=$ $\displaystyle
V_{a}-\frac{1}{\beta_{J}}\,\ln\Big{(}\sum_{\\{J_{ij}^{(to)}\\}}\exp\Big{\\{}\beta_{J}\,\sum_{i<j\centernot\in\mathcal{T}}J_{ij}^{(oo)}F_{ij}^{(o-t-o)}\Big{\\}}\Big{)}$

The stationary distributions induced by the diffusion process in Eq. (2) with
these effective potentials have a Boltzmann-form
$\mathbb{P}_{\mathcal{T}}(\mathbb{J}^{(tt)})=e^{-\beta_{J}V_{tt}}/\mathcal{Z}_{\mathcal{T}}$
and
$\mathbb{P}_{\mathcal{O}}(\mathbb{J}^{(oo)})=e^{-\beta_{J}V_{oo}}/\mathcal{Z}_{\mathcal{O}}$
where $\mathcal{Z}_{\mathcal{T}}$ and $\mathcal{Z}_{\mathcal{O}}$ are the
partition function of the genotypes $\mathbb{J}^{(tt)}$ and that of the
genotypes $\mathbb{J}^{(oo)}$, respectively 555Note that
$\mathcal{Z}_{\mathcal{T}}$ and $\mathcal{Z}_{\mathcal{O}}$ are different from
the partition functions of the target spins’ subsystem
$Z_{1}(\mathbb{J}^{(tt)})$ and that of the non-target spins’ subsystem
$\tilde{Z}_{1}(\mathbb{J}^{(oo)})$ given in [26].. Here we replace
$\mathbb{J}^{(to)}$ by the replica matrix $\sigma_{i}^{k}:=J_{ik}^{(to)}$, to
be obtained. (This stepwise scheme is valid, as we are concerned with the
equilibrium property). Denoting $n:=T_{s}/T_{J}$, we can compute
$\mathcal{Z}_{\mathcal{T}}$ and $\mathcal{Z}_{\mathcal{O}}$ as

$\mathcal{Z}_{\mathcal{T}}=\int\prod_{i<j\in\mathcal{T}}dJ^{(tt)}_{ij}\sum_{\big{\\{}s_{i};s_{i}^{\alpha};\sigma_{i}^{k}\big{\\}}_{i\in\mathcal{T}}}\exp\left\\{\beta_{J}\sum_{i<j}\Big{[}-\frac{N_{t}}{2}\,\Big{(}J_{ij}^{(tt)}\Big{)}^{2}+J_{ij}^{(tt)}\Big{(}\frac{1}{N_{t}}\,\Big{|}\sum_{i\in\mathcal{T}}s_{i}\Big{|}+\underbrace{\frac{1}{n}\,\sum_{a=1}^{n}s_{i}^{a}s_{j}^{a}}_{s-{\rm
replicas}}+\underbrace{\frac{1}{N-N_{t}}\,\sum_{k=1}^{N-N_{t}}\sigma_{i}^{k}\sigma_{j}^{k}}_{\sigma-{\rm
replicas}}\Big{)}\Big{]}\right\\}$
$\mathcal{Z}_{\mathcal{O}}=\int\prod_{i<j\centernot\in\mathcal{T}}dJ^{(oo)}_{ij}\sum_{\big{\\{}s_{i}^{\alpha};\sigma_{i}^{k}\big{\\}}_{i\centernot\in\mathcal{T}}}\exp\left\\{\beta_{J}\sum_{i<j}\Big{[}-\frac{N-N_{t}}{2}\,\Big{(}J_{ij}^{(oo)}\Big{)}^{2}+J_{ij}^{(oo)}\Big{(}\tilde{K}+\underbrace{\frac{1}{n}\,\sum_{a=1}^{n}s_{i}^{a}s_{j}^{a}}_{s-{\rm
replicas}}+\underbrace{\frac{1}{N_{t}}\,\sum_{k=1}^{N_{t}}\sigma_{i}^{k}\sigma_{j}^{k}}_{\sigma-{\rm
replicas}}\Big{)}\Big{]}\right\\}$

In writing these equations, we assume that the fitness acts only on
$J_{ij}^{(tt)}$ 666This can already be seen directly from Eq. (5) as the
expression inside the bracket does not depend on $J_{ij}^{(oo)}$ and
$J_{ij}^{(to)}$. and hence in $\mathcal{Z}_{\mathcal{O}}$ we replace the
fitness field $K$ by a constant $\tilde{K}$, which eventually will be set to 0
by virtue of calculations of the observables for non-target spins. Although
neglecting the fitness’s effect on $J_{ij}^{(oo)}$ does not follow exactly the
above-mentioned implementation of the model, we expect that this holds true in
the long times limit because otherwise both target- and non-target
configurations at equilibrium would determine the fitness. This restriction
hence corresponds to a first-order approximation of the fitness’s effect,
while a term that affects the dynamics of $J_{ij}^{(oo)}$ and $J_{ij}^{(to)}$
is considered to be of higher order.

We here propose to interpret $\sigma_{i}^{k}$ as the $k$-th replica of another
variable $\sigma_{i}\in\\{-1,1\\}$ that is also located at the site $i$ of the
graph (generally $\sigma_{i}\neq s_{i}$). To distinguish these different types
of replica from each other, we call $s_{i}^{a}$ spin-replica and
$\sigma_{i}^{k}$ coupling-replicas. Following this interpretation, apart from
$n$ that appears as the number of spin-replicas $s_{i}^{a}$,
$a=\\{1,\cdots,n\\}$, in $\mathcal{Z}_{\mathcal{T}}$ and
$\mathcal{Z}_{\mathcal{O}}$ 777Though $n$ appears as integer number here, we
will analytically continue to real positive $n$ in subsequent steps of
calculations., we thus have $N-N_{t}$ coupling-replicas, $\sigma_{i}^{k}$, for
$i\in\mathcal{T}$ and $k=\\{1,\cdots,N-N_{t}\\}$ and $N_{t}$ coupling-
replicas, $\sigma_{i}^{k}$, for $i\centernot\in\mathcal{T}$ and
$k=\\{1,\cdots,N_{t}\\}$ respectively. As in general, none of these numbers
are zero, our double-replica approach does not correspond to the conventional
quenched limit in spin-glass models [16]. Once setting $K=0$ and neglecting
the terms corresponding to $F^{(t-o-t)}$ and $F^{(o-t-o)}$, we recover Coolen
et. al. model for neural systems with dynamic synapses [18]. In contrast to
the use of a Hamiltonian for the couplings adopted in [22], here, we have
introduced the effective potential for couplings that, by using the time-scale
separation between the dynamics of genotypes and that of phenotypes, allows
for the integration of the spin dynamics specified by the Hamiltonian Eq. (1)
into the Langevin dynamics of the couplings through the second term in Eq.
(6).

Here we characterise the equilibrium behaviour of the model by the average
fitness, $m$, the overlap between different spin replicas $a$ and $b$,
$q_{ab}$, and the correlation between adjacent links $Q$. Additionally, we
want to quantify the mean value of the couplings among the target spins only
$\Phi$ . Let $\mathbb{E}[\cdot]$ and $\mathbb{\tilde{E}}[\cdot]$ denote
ensemble average over $\mathbb{P}_{\mathcal{T}}(\mathbb{J}^{(tt)})$ and
$\mathbb{P}_{\mathcal{O}}(\mathbb{J}^{(oo)})$, respectively. These order
parameters are given by

$\displaystyle
m_{a}=\mathbb{E}\Big{[}s_{i}^{a}\Big{]}_{i\in\mathcal{T}}\,,\,\,\,q_{ab}=\mathbb{E}\Big{[}s_{i}^{a}s_{i}^{b}\Big{]}_{i\in\mathcal{T}}$
(9a) $\displaystyle
Q_{kk^{\prime}}=\mathbb{E}\Big{[}J_{ik}^{(to)}J_{ik^{\prime}}^{(to)}\Big{]}_{\begin{subarray}{c}i\in\mathcal{T}\\\
k,k^{\prime}\centernot\in\mathcal{T}\end{subarray}}$ (9b)
$\displaystyle\Phi=\mathbb{E}\Big{[}J_{ij}^{(tt)}\Big{]}_{\begin{subarray}{c}i,j\in\mathcal{T}\end{subarray}}$
(9c)

Similarly, for the non-target spins we have

$\displaystyle\tilde{m}_{a}=\mathbb{\tilde{E}}\Big{[}s_{i}^{a}\Big{]}_{i\centernot\in\mathcal{T}}\,,\,\,\,\tilde{q}_{ab}=\mathbb{\tilde{E}}\Big{[}s_{i}^{a}s_{i}^{b}\Big{]}_{i\centernot\in\mathcal{T}}$
(10a)
$\displaystyle\tilde{Q}_{kk^{\prime}}=\mathbb{\tilde{E}}\Big{[}J_{ik}^{(to)}J_{ik^{\prime}}^{(to)}\Big{]}_{\begin{subarray}{c}i\centernot\in\mathcal{T}\\\
k,k^{\prime}\in\mathcal{T}\end{subarray}}$ (10b)
$\displaystyle\tilde{\Phi}=\mathbb{\tilde{E}}\Big{[}J_{ik}^{(oo)}\Big{]}_{\begin{subarray}{c}i,j\centernot\in\mathcal{T}\end{subarray}}$
(10c)

In the thermodynamics limit, $N\rightarrow\infty$ and
$N_{t}\rightarrow\infty$, while keeping $p=N_{t}/N$ fixed, using a replica
symmetric ansatz for the variables $m_{a}=m$, $q_{ab}=q$ and
$Q_{kk^{\prime}}=Q$, $\tilde{m}_{a}=\tilde{m}$; $\tilde{q}_{ab}=\tilde{q}$ and
$\tilde{Q}_{kk^{\prime}}=\tilde{Q}$, as well as, $M_{ak}=M$ and
$\tilde{M}_{ak}=\tilde{M}$, where
$M_{ak}=\mathbb{E}\big{[}s^{a}_{i}J_{ik}^{(to)}\big{]}_{\begin{subarray}{c}i\in\mathcal{T}\\\
k\centernot\in\mathcal{T}\end{subarray}}$ and
$\tilde{M}_{ak}=\mathbb{\tilde{E}}\big{[}s^{a}_{i}J_{ik}^{(to)}\big{]}_{\begin{subarray}{c}i\centernot\in\mathcal{T}\\\
k\in\mathcal{T}\end{subarray}}$, we obtain the following free energy
densities:

$\displaystyle f^{\rm RS}_{\mathcal{T}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,\left\\{\frac{q}{n}+\frac{Q}{N-N_{t}}+\frac{(n-1)q^{2}}{2n}+\frac{Q^{2}}{2}+M^{2}\right\\}$
(11a) $\displaystyle-$
$\displaystyle\frac{1}{\beta_{J}}\ln\left[\sum_{z=0}^{N_{t}}\binom{N_{t}}{z}\sum_{k=0}^{N-N_{t}}\binom{N-N_{t}}{k}I_{kz}\right]$
$\displaystyle f^{\rm RS}_{\mathcal{O}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,\left\\{\frac{\tilde{q}}{n}+\frac{\tilde{Q}}{N_{t}}+\frac{(n-1)\tilde{q}^{2}}{2n}+\frac{\tilde{Q}^{2}}{2}+\tilde{M}^{2}\right\\}$
(11b) $\displaystyle-$
$\displaystyle\frac{1}{\beta_{J}}\ln\left\\{\sum_{k=0}^{N_{t}}\binom{N_{t}}{k}\tilde{I}_{k}\right\\}$

From the extremum condition of these free energies we can compute all the
model order parameters via a set of self-consistency equations. These
equations as well as the functions $I_{kz}$ and $\tilde{I}_{k}$ are given in
the SM 888Here $I_{kz}=I_{kz}(m,q,Q,r,M)$ and
$\tilde{I}_{k}=\tilde{I}_{k}(\tilde{m},\tilde{q},\tilde{Q},\tilde{r},\tilde{M})$
with $r,M,\tilde{r},\tilde{M}$ are other variables that the free energy
densities depend on. Since these observables for target and non-target spins
have their behaviour correlated to that of $(m,q,Q)$ and
$(\tilde{m},\tilde{q},\tilde{Q})$, respectively, they do not provide
additional information about the structure of the model phase diagram.. The
use of the replica symmetric is justified in most part of the $(T_{s},T_{J})$
parameter space from the stability analysis [32]. At low $Ts$ and $T_{J}$ the
replica symmetry is broken, which, we will not explore fully. Nevertheless we
will discuss later how robustness of phenotypes, postulated for biological
systems that reproduce similar offspring, is lost in that scenario. The
replica-symmetric free energy densities allow us to derive

$\displaystyle\Phi$ $\displaystyle=$
$\displaystyle\frac{1}{N_{t}}\left[\Psi+m^{2}+r^{2}\right]$ (12a)
$\displaystyle\tilde{\Phi}$ $\displaystyle=$
$\displaystyle\frac{\tilde{m}^{2}+\tilde{r}^{2}}{N-N_{t}}\,.$ (12b)

where $r$ and $\tilde{r}$ are defined and computed in the SM.

##  III Phase diagram

Figure 1: Magnetisation for target spins $m$ (A). Overlap between different
replicas for target spins $q$ (B). Averaged correlation of a pair of couplings
between a target and a non-target spin that share a common non-target spin $Q$
(C). Averaged frustration among target spins $\Phi$ (D). Here $N_{t}=10$,
$N=100$. Note the $y$-axis is on logarithmic scale.

In Fig. 1 we depict the order parameters as function of the temperature
$T_{s}$ and $T_{J}$ for a particular choice of $N=100$ and $N_{t}=10$. Here
for each point $(T_{s},T_{J})$ of the phase diagram, we solve numerically the
set of mean-field equations for $m,q,Q$, while computing $\Phi$ from the
knowledge of these quantities. In terms of only the magnetisation $m$ and the
overlap between spin-replicas $q$ for target spins, we observe three distinct
phases that are typical for spin-glass systems, namely, $m=q=0$ (paramagnet
phase); $m=0,q>0$ (spin-glass phase), and $m>0,q>0$ with $\sqrt{q}>m$, (
target-ferromagnet phase, ’t-ferro’ in short). The transitions between the
phases are second-order at small $T_{s}/T_{J}$, but become discontinuous
(first-order) at large $T_{s}/T_{J}$. At a much lower value of $T_{J}$ there
is a region where, apart from having a non-zero magnetisation of target spins,
the order parameter $Q=\langle J^{(to)}J^{(to)}\rangle$ starts to become non-
zero. As can be anticipated from Eq. (12a), the mean value of $J_{ij}^{(tt)}$
also varies from region to region in accordance to the change of $m$ and that
of $Q$. Note that on sharp contrary to the transition between paramagnet and
spin-glass which is similar to that of the Sherrington-Kirpatrick (SK) model,
the phenotype-genotype coupling results in a repositioning of the boundary
between spin-glass and t-ferro. Such difference arises from the non-zero
correlation of the genotypes. Expanding the free energy $f^{\rm
RS}_{\mathcal{T}}$ for small $m$ and $q$, the transition between paramagnet
and t-ferro occurs at $T_{s}^{\mathbb{P}\rightarrow\mathbb{F}}=\kappa$, where
$\kappa=2^{-N_{t}}\sum_{z=0}^{N_{t}}\binom{N_{t}}{z}|N_{t}-2z|/N_{t}$,while
the spin-glass to t-ferro transition occurs at
$\big{[}1+(n-1)q(\beta_{s}^{\mathbb{SP}\rightarrow\mathbb{F}})\big{]}\cdot\big{(}\beta_{s}^{\mathbb{SP}\rightarrow\mathbb{F}}\kappa\big{)}=1$.
We also check that both the magnetisation $\tilde{m}$ of the non-target spins
and the average value $\tilde{\Phi}$ of $J_{ij}^{(oo)}$ are always zero as the
non-target spin subsystem remains frustrated all the time, while the spin
overlap $\tilde{q}$ can undergo a transition from paramagnet to spin-glass, in
the same way as the SK model. The phase diagrams of these quantities are given
in the SM.

Figure 2: The model phase diagram. Here SP1 denotes the spin-glass phase with
$m=\tilde{m}=Q=\tilde{Q}=0$ and $q,\tilde{q}>0$; SP2 denotes the spin-glass
phase with $\tilde{m}=Q=0$ and $q,m>0$ (within this region, both $\tilde{q}$
and $\tilde{Q}$ can be either zero or non-zero, see SM); P1 denotes the
paramagnet phase with $m=\tilde{m}=q=\tilde{q}=Q=\tilde{Q}=0$; P2 denotes the
paramagnet phase with $m=\tilde{m}=q=\tilde{q}=\tilde{Q}=0$ but $Q>0$; R
denotes the robust fitted phase with $\tilde{m}=\tilde{q}=\tilde{Q}=0$ but
$m,q,Q>0$. Here $N_{t}=10$, $N=100$. Note the $y$-axis is on logarithmic
scale.

Combining the behaviour of the order parameters altogether, we obtain the
model phase structure in Fig. 2. It contains five distinct regions. At low
genotypic selection pressure $T_{J}\geq e^{-1}$, only the first spin-glass SP1
and the paramagnet P1 phases with zero fitness are observed. However, as the
genotypic selection pressure increases other phases emerge. At sufficiently
low $T_{J}$, a robust fitted phase denoted by R ($m,q,Q>0$) emerges in an
intermediate range of $T_{s}$ (here $\tilde{m}=\tilde{q}=\tilde{Q}=0$).
Adjacent to this phase on the side of high phenotypic noise is the second
paramagnet phase P2 where the fitness value is low ($m=q=0$) but there exists
some structure in the genotypes such that $Q>0$. On the other hand, for lower
$T_{s}$, the system is in the second spin-glass phase SP2 with high fitness
but non-robust genoptypes ($m\simeq q\simeq 1,Q=0$). In particular, the
transition from R to SP2 is marked by a replica symmetry breaking (RSB) which
indicates the loss of stability of the replica symmetric (RS) solutions 999The
eigenvalue of the Hessian changes its sign on the border between these phases,
see SM. Note that though we do not calculate the full hierarchy of replica
symmetry breaking (RSB) à la Parisi, but as long as RSB happens this non-
robustness is true for the full RSB solution.. The broad distribution of gene-
gene correlations in the RSB phase implies that the genotype of offsprings is
not preserved, in contrast to the RS phase. In the biological context, this
means that replication is no longer stable so that genotypes are not conserved
over generations.

In overall, the phase diagram agrees with what was observed numerically in
[14, 15]. However, thanks to the explicit account of the coupling-replicas, so
that $Q$ can be treated as an order parameter upon which the free energy
density depends, we discover the existence of the second paramagnet phase P2
that was not reported before. This phase can be interpreted as a _precusor_
region, in which genotypes are structured in such a way that supports
ferromagnetic ordering among target spins, and hence have potentiality to
acquire a high fitness, but due to the high fluctuation induced by $T_{s}$,
this fitness cannot be maintained. Furthermore, by considering separately the
effective dynamics of the target and non-target subsystems, $S_{\mathcal{T}}$
and $S_{\mathcal{O}}$, our approach can differentiate the phase SP1 from SP2.
The previous approach [22] only stressed the distinct arrangement of target
spins in the SP2 region, where the subsystem of target spins becomes
ferromagnetic whereas that of non-target ones remains spin-glass. Our present
approach shows that this is no longer true for a high value of $T_{J}$. Upon
increasing $T_{J}$, this ferromagnetic ordering is destroyed by genotypic
fluctuations.

Figure 3: The phase diagrams obtained by combining the behaviour of $m$ as
function of $T_{s}$ and $T_{J}$ with that of $Q$ for different values of $p$,
at sufficiently low $T_{J}$. Here $N=100$.

In addition, the present analysis allows one to obtain quantitative dependence
of genotypic and phenotypic robustness on the fraction of targets. The phase
diagram in Fig. 2 includes global information of the system including weak
selection region without achieving nonzero fitness, $m$ of target, whereas, of
biological interest is if the fitted state is evolved robustly by the
selection. To this end we focus on the low $T_{J}$ region of the phase diagram
to explore the dependence of the system behaviour on the fraction $p=N_{t}/N$
of target spins. While overall, the phase structure is similar for different
$p$ in Fig. 3, in particular, the robust fitted (yellow) region seems to
change slightly with increasing $p$, the relative size and exact location of
all the other phases vary with $p$. This suggests that a more quantitative
analysis is needed to understand the genotype-phenotype relationship as
function of $p$. We carry on this analysis in the next section.

##  IV  Structure of the robust fitted phase

Figure 4:  (A) Magnetisation of target spins $m$ as function of $T_{s}$ at
fixed $T_{J}=0.005$ for various number of target spins
$p=0.04,0.1,0.2,0.3,0.4,0.5$. (B) The same for genetic overlap $Q$. (C)
Frustration defined as $1-Q$ as function of the number of non-target spins in
the robust fitted phase (i.e. $T_{s}\in[T_{c}^{(1)},T_{c}^{(2)}]$) at fixed
$T_{J}=0.005$. (D) The highest value of $T_{J}$ at which $Q$ remains non-zero
as function of $p$. (E) Magnetisation of target spins $m$ as function of
$T_{J}$ at fixed $T_{s}=1.3$ for various fraction of target spins
$p=0.2,0.3,0.4,0.5,0.6,0.7,0.8$. (F) The value of $T_{J}$ at which the
magnetisation of target spins $m$ drops to zero at fixed $T_{s}=1.3$ for
different $p$ from (E). Similar behaviour to (E) and (F) is observed for
others $T_{s}\in[T_{c}^{(1)},T_{c}^{(2)}]$. Here $N=100$.

The most relevant region in the phase diagram is robust fitted phase R, which
is characterized by both the high fitness ($m>0$) and robustness ($Q>0$). For
a sufficiently low given $T_{J}$ (i.e. high selection pressure), the phase is
bounded by $T_{s}\in[T_{c}^{(1)},T_{c}^{(2)}]$. Below $T_{c}^{(1)}$, $Q$ goes
to zero, and above $T_{c}^{(2)}$, $m$ goes to zero, whereas these transition
points depend on $T_{J}$. We first investigate the dependence on $p$ of the R
region by fixing $T_{J}$. In Fig. 4 (A) and (B) we fixed $T_{J}=0.005$. First,
for a wide range of $p\in[0.04,0.5]$, the temperature $T_{c}^{(1)}$ of the
latter transition in (B) does not depend on $p$ (see SM for the zoom-in of a
small dip of $m$ at this point). On the other hand, $T_{c}^{(2)}$ slightly
increases with increasing $p$ in (A), indicating that the fitness of the high
$p$ case is relatively more robust to noise than the low $p$ case.

In Fig. 4 (C) $1-Q$ is almost constant against $T_{s}$ in the R phase. This
constant value was found to increase with the number of non-target spins. Note
that $Q\sim 1$ implies that the offspring of genotypes are preserved. The
increase in $1-Q$, thus means the increase in redundancy of genotypes, as is
supported by a larger number of non-targets. Such a redundancy has another
meaning in the context of spin-glass systems, where it is indeed equal to the
local frustration (in the $(t-o-t)$ triads with positive $J^{(tt)}$).

In contrast, apart from the two critical points $T_{c}^{(1)}$ and
$T_{c}^{(2)}$, the fitness $m$ does not depend on $p$. It follows a unique
curve, independent of $p$. Even though the increase in genetic heterogeneity
$1-Q$ for more non-targets may perturb the target spin configuration, the
fitness $m$ remains unchanged even for smaller $p$.

Then, we estimate the critical value of $T_{J}$ below which the R phase can
exist. While this critical value denoted by $T_{J}^{(Q)}$ can depend on
$T_{s}$, as see in Fig 2 and 3, it can be approximately identified with the
upper part of the P2 phase from the phase diagram. In Fig. 4 (D) $T_{J}^{(Q)}$
is shown to decrease with $p$. This result together with that in Fig. 4 (A)
mean that the higher $p$ is, the higher selection pressure is needed to
achieve robustness, but once it is achieved, a system with larger $p$ is more
robust with respect to phenotypic noise than one with smaller $p$.

Finally, we examine the fitness as function of $T_{J}$ at fixed $T_{s}=1.3$
for various $p$ in Fig. 4 (E). While fitness decreases with $T_{J}$, its
behaviour with the increase in $p$ is non-monotonic. This behavior is further
shown in Fig. 4 (F) where the critical genotypic noise $T_{J}^{(m)}$ at which
the fitness becomes non-zero is plotted versus $p$. Similar behavior is
observed for others $T_{s}\in[T_{c}^{(1)},T_{c}^{(2)}]$. The result supports
$p=0.5$ as the maximal value of $T_{J}^{(m)}$, suggesting the existence of an
optimal fraction of target spins to acquire high fitness in this intermediate
range of $T_{s}$ 101010This optimality, however, does not happen in the
spinglass phase SP2 with low temperature $T_{s}$, since the selection pressure
needed to achieve non-zero fitness crucially depends on $p$. In fact as $p$ is
increased, higher pressure (i.e. lower $T_{J}$) is needed to achieve non-zero
$m$ of targets, i.e, to transition from SP1 to SP2..

##  V Mutational susceptibility and Phenotypic susceptibility in the robust
fitted phase

Correlation between variances of phenotypes due to genetic changes and to
noise has been discussed both in experiments and simulations, and
relationships to robustness have been discussed both theoretically [35, 36, 9,
37, 13, 38] and experimentally [39, 40, 41, 42]. In statistical physics, this
issue can be analyzed in terms of susceptibility, as it is proportional to the
variance. Then, we need to study the susceptibility due to genetic mutation,
in addition to the standard susceptibility.

In the context of this model, mutations are defined as those change of the
genotypes $J_{ij}$ that might happen spontaneously and independently from the
dynamics specified previously. Let $\delta\Psi_{i}(\delta J_{jk})$ denote the
change of the average local magnetisation of a target spin $i$ 111111More
precisely, in [13] $\Psi_{i}=\langle{\rm sign}\,(m)s_{i}\rangle_{T_{s}}$. upon
mutating a genotype $J_{jk}\rightarrow J_{jk}+\delta J_{jk}$. The mutational
susceptibility of this target spin w.r.t such a change $M_{i,jk}$ then can be
defined as

$M_{i,jk}=\lim_{\delta J_{jk}\rightarrow 0}\delta\Psi_{i}(\delta
J_{jk})/\delta
J_{jk}=\big{\langle}s_{i}s_{j}s_{k}\big{\rangle}-\big{\langle}s_{i}\big{\rangle}\big{\langle}s_{j}s_{k}\big{\rangle}\,.$

In general,
$J_{jk}\in\mathbb{J}=\mathbb{J}^{(tt)}\cup\mathbb{J}^{(oo)}\cup\mathbb{J}^{(to)}$.
However, since fitness is determined solely by the configurations of target
spins at equilibrium, we consider only $J_{jk}\in\mathbb{J}^{(tt)}$ and show
that the average of this mutational susceptibility over all triples
$(i,j,k)\in\mathcal{T}$ is equal to

$M:=\frac{1}{\binom{N_{t}}{3}}\sum_{(i,j,k)}M_{i,jk}=2\beta_{J}^{2}m\chi_{m}-\beta_{J}\frac{\partial^{3}f}{\partial
h^{3}}\Big{|}_{h=0}\\\ $ (13)

where $\chi_{m}:=-\lim_{h\rightarrow 0}\partial^{2}f/\partial h^{2}$ is the
susceptibility of target spins. The quantities $M$ and $\chi_{m}$ correspond
to the susceptibility to mutation of the genotypes and susceptibility to
perturbation by an external field, $h$, respectively. We can expect that in
the robust fitted phase there exists a relation between $M$ and $\chi_{m}$
121212From the replica-symmetric free energy density $f_{\mathcal{T}}^{\rm
RS}$ we get $\chi^{\rm RS}_{m}=\beta_{s}\big{[}1+(n-1)q\big{]}$, while $M$ can
not be obtained directly within this ansatz.. In fact, let
$X:=\lim_{h\rightarrow 0}\partial^{3}f_{\mathcal{T}}/\partial h^{3}$. For $L$
given in the SM Eq. (18a) ($L/\beta_{J}$ has the meaning of an effective
Hamiltonian that is defined in the combined space
$\\{s^{a}_{i},\sigma^{k}_{i}\\}$ of $s$-replicas and $J$-replicas), according
to the definitions

$\displaystyle X\propto\sum_{a,b,c=1}^{n}{\rm Tr}\,(s^{a}s^{b}s^{c}e^{L})/{\rm
Tr}\,(e^{L})\,,$ $\displaystyle m:=\langle
s^{a}\rangle=\frac{1}{n}\,\sum_{a=1}^{n}{\rm Tr}\,(s^{a}e^{L})/{\rm
Tr}\,(e^{L})\,,$ $\displaystyle\chi_{m}:=\lim_{h\rightarrow 0}\frac{\partial
m}{\partial h}=\frac{\beta_{s}}{n}\,\sum_{a,b=1}^{n}{\rm
Tr}\,(s^{a}s^{b}e^{L})/{\rm Tr}\,(e^{L})\,,$

the symmetry between different replicas in the robust fitted phase implies
that the third moment $X^{\rm RS}$ is proportional to the product of the first
and second moments $m^{\rm RS}\chi_{m}^{\rm RS}$. Therefore, approximately,
$M\propto\chi_{m}^{\rm RS}$. This proportionality between the two
susceptibilities, implying a correlation between phenotypic changes due to
genetic variation and those in response to environmental perturbations [45],
does not exist in the RSB phase, as the second term in Eq. (13) is no longer
proportional to $\chi_{m}$.

##  VI Discussion

In the paper we propose a new approach towards biological evolution due to the
interrelationship between genotype and phenotype where fitness is determined
solely by the latter but not by the former. Though the emergence of structured
genotypes from initially random couplings under this relation has been
numerically reported, apart from a study which imposed a specific condition on
the couplings [22], this has not been studied analytically yet. We here are
able to tackle this problem thanks to what we termed _double-replica_ theory.
Within this framework we obtain the phase diagram, that is classified not only
by the fitness but also by the overlap in dual replicas. The diagram is not
only in good agreement with previous studies (including paramagnet, t-ferro
and robust fitted phases, all existing at sufficiently low $T_{J}$), but also
contains previously undiscovered phases. These include the first spin-glass
phase SP1 and the second paramagnet phase P2. The former corresponds to a
system with both target and non-target spins residing in a spin-glass phase
(at low selection pressure), while the latter corresponds to a paramagnetic
phase for all spins but with retaining genetic correlations encoding target-
and non-target couplings (at high selection pressure and high $T_{s}$). Here
even though the genotypes favor a high value of fitness, due to large
fluctuations induced by $T_{s}$, such value can not be maintained over
generations. The existence of the phase suggests that even though the average
fitness is zero due to large noise, there exists genetic precursor to generate
individuals with non-zero fitness. The relevance of this scenario to
evolutionary biology, needs to be explored in future, though.

The system can only acquire high fitness at some $T_{s}\leq T_{c}^{(2)}$,
where the fitness increases discontinuously. If $T_{s}$ is too low, then RSB
will happen, to a phase without _genetic_ overlap, where biologically required
robustness of genotypes is lost. Hence a lower bound of $T_{s}\geq
T_{c}^{(1)}$ is necessary to have RS and robustness, accordingly.

From this approach, the target-fraction dependence of genotypic and phenotypic
robustness can also be understood quantitatively. Such dependence is
quantified via the behaviour of the fitness $m$ and genetic redundancy $1-Q$
in the robust region bounded by $T_{c}^{(1)}$ and $T_{c}^{2}$. Here we find
that a genetically homogeneous population can only be robustly reproduced
under a sufficiently high selection pressure and under a sufficient level of
phenotypic noise (temperature). As the fraction of target spins is increased,
the robust fitted phase is slightly expanded to a higher temperature, whereas
higher selection pressure is needed to achieve it. The existence of an optimal
fraction for attaining high fitness under intermediate phenotypic noise is
suggested. This may explain why, in biological systems, such as in proteins,
the fraction of units that are responsible for function is generally limited,
and a sufficient fraction of non-functional units is needed, providing
redundancy.

In the present theory, the proportionality between the standard thermodynamic
susceptibility and mutational susceptibility is derived in the robust fitted
phase. As the susceptibility measures the change of fitness due to varying
conditions, a correlation between responses to environmental perturbations and
that by genetic changes is suggested. Such correlation, or evolutionary
fluctuation-response relationship [9, 13, 35, 36, 37, 38] has been observed in
experimental data from the evolution of protein dynamics and bacterial protein
expressions, whereas we can derive it here under replica symmetry assumption.
As argued in [13], such correlation can only be achieved in the replica
symmetric region where the original high-dimensional dynamics of the
phenotypes are reduced to a low-dimensional manifold due to evolution towards
robustness. The variation of fitness due to noise and that due to mutation
then happen to occur along the same low-dimensional manifold, resulting in a
correlation between them. If RSB occurs, such restriction of the phenotypic
dynamics no longer exits, because in this case, changes of fitness upon
varying the environmental conditions will vary arbitrarily from realisation to
realisation of the $J_{ij}$’s dynamics. As a result, the system will have
random, uncorrelated responses to noise and to mutations.

In our formulation by assuming that the entire system reaches an equilibrium,
we approximate the effect of the slowly-evolving $\mathbb{J}^{(to)}$ that
couple the subsystem $S_{\mathcal{T}}$ to $S_{\mathcal{O}}$ on these
subsystems’ own dynamics by the equilibrium correlations $\langle
J^{(to)}J^{(to)}\rangle$. Such correlations are then incorporated separately
into each of the dynamics of the $\mathbb{J}^{(tt)}$ and $\mathbb{J}^{(oo)}$
couplings by modifying the effective potentials of these dynamics, thus making
the dynamics of these two different sets of coupling _independent_ of each
other. In an equilibrium statistical physics formulation, this leads to the
necessity of introducing another type of replica, so-called coupling-replicas
into the partition functions, besides the first (standard) spin-replicas that
take care of the effect of the phenotypes on the evolution of genotypes. This
scheme hence allows us to treat the model in a standard mean-field manner. On
one hand, being of mean-field nature, our approach can not provide a formal
argument to support the hypothesis of [13] about the emergent dimensional
reduction from phenotype-genotype co-evolution. On the other hand, the
correlation between the mutational- and environmental susceptibility in the
robust fitted phase suggests the existence of a funnel-like landscape [46, 47]
that reinforces the dynamics to reside in a low-dimensional manifold by a
global attraction.

The current choice of fitness for the sake of simplicity, however, limits the
possibility of having different global maxima in the fitness landscape. One
can enrich the model behavior by determining fitness either by a combination
of $N_{\rm fit}$ different target spin configurations or by a set of gauge-
equivalent configurations.

In the present framework, since the couplings are symmetric, we constructed
the effective potential of the coupling dynamics based partly on the existence
of an energy landscape. For those models in other contexts [48, 49] having
such a landscape picture, we expect a straightforward application of our
approach. Furthermore, the present double-replica theory can be extended to
those stochastic dynamical systems that are not governed by Hamiltonian
dynamics as well. In this case, instead of the effective potential and its
associated partition function, one would need to charaterize the ensemble of
trajectories in the combined space of phenotypes and genotypes, using the
moment generating function [50, 51]. While we so far have solely used
phenotypes and genotypes as the main example of our approach, such an
extension would allow for the applications to co-evolution of gene-expression
patterns and the gene-regulatory networks [9], that of species abundances and
their ecological networks [52], and that of neuronal activities and network
shaped by neural dynamics (learning) [53, 54].

###### Acknowledgements.

We acknowledge support from Novo Nordisk Foundation and would like to thank
Ayaka Sakata, Koji Hukushima and Yoshiyuki Kabashima for stimulating
discussion.

## References

  * Wright [1932] S. Wright, The roles of mutation, inbreeding, crossbreeding and selection in evolution, Proceedings of the XI International Congress of Genetics 8, 209 (1932).
  * de Visser and Krug [2014] J. A. G. de Visser and J. Krug, Empirical fitness landscapes and the predictability of evolution, Nature Reviews Genetics 15, 480 (2014).
  * McAdams and Arkin [1997] H. H. McAdams and A. Arkin, Stochastic mechanisms in gene expression, Proceedings of the National Academy of Sciences 94, 814 (1997).
  * Elowitz  _et al._ [2002] M. B. Elowitz, A. J. Levine, E. D. Siggia, and P. S. Swain, Stochastic gene expression in a single cell, Science 297, 1183 (2002).
  * Furusawa  _et al._ [2005] C. Furusawa, T. Suzuki, A. Kashiwagi, T. Yomo, and K. Kaneko, Ubiquity of log-normal distributions in intra-cellular reaction dynamics, Biophysics 1, 25 (2005).
  * Bar-Even  _et al._ [2006] A. Bar-Even, J. Paulsson, N. Maheshri, M. Carmi, E. O’Shea, Y. Pilpel, and N. Barkai, Noise in protein expression scales with natural protein abundance, Nature Genetics 38, 636 (2006).
  * Wagner [2013] A. Wagner,  _Robustness and evolvability in living systems_ , Vol. 24 (Princeton university press, 2013).
  * Ancel and Fontana [2000] L. W. Ancel and W. Fontana, Plasticity, evolvability, and modularity in RNA, Journal of Experimental Zoology 288, 242 (2000).
  * Kaneko [2007] K. Kaneko, Evolution of robustness to noise and mutation in gene expression dynamics, PLOS ONE 2, 1 (2007).
  * Shreif and Periwal [2014] Z. Shreif and V. Periwal, A network characteristic that correlates environmental and genetic robustness, PLOS Computational Biology 10, 1 (2014).
  * Kubo  _et al._ [1985] R. Kubo, M. Toda, and N. Hashitsume,  _Statistical Physics II: Nonequilibrium Statistical Mechanics_, Springer Series in Solid-State Sciences (Springer, 1985).
  * Sethna [2021] J. P. Sethna,  _Statistical mechanics: entropy, order parameters, and complexity_ , Vol. 14 (Oxford University Press, USA, 2021).
  * Sakata and Kaneko [2020] A. Sakata and K. Kaneko, Dimensional reduction in evolving spin-glass model: Correlation of phenotypic responses to environmental and mutational changes, Phys. Rev. Lett. 124, 218101 (2020).
  * Sakata  _et al._ [2009a] A. Sakata, K. Hukushima, and K. Kaneko, Funnel landscape and mutational robustness as a result of evolution under thermal noise, Phys. Rev. Lett. 102, 148101 (2009a).
  * Sakata  _et al._ [2009b] A. Sakata, K. Hukushima, and K. Kaneko, Statistical-mechanical study of evolution of robustness in noisy environments, Phys. Rev. E 80, 051919 (2009b).
  * Mezard  _et al._ [1986] M. Mezard, G. Parisi, and M. Virasoro,  _Spin Glass Theory and Beyond_ (World Scientific, Singapore, 1986).
  * Penney  _et al._ [1993] R. W. Penney, A. C. C. Coolen, and D. Sherrington, Coupled dynamics of fast spins and slow interactions in neural networks and spin systems, Journal of Physics A: Mathematical and General 26, 3681 (1993).
  * Coolen  _et al._ [1993] A. C. C. Coolen, R. W. Penney, and D. Sherrington, Coupled dynamics of fast spins and slow interactions: An alternative perspective on replicas, Phys. Rev. B 48, 16116 (1993).
  * Penney and Sherrington [1994] R. Penney and D. Sherrington, Slow interaction dynamics in spin-glass models, Journal of Physics A: Mathematical and General 27, 4027 (1994).
  * Dotsenko  _et al._ [1994] V. Dotsenko, S. Franz, and M. Mezard, Partial annealing and overfrustration in disordered systems, Journal of Physics A: Mathematical and General 27, 2351 (1994).
  * Uezu  _et al._ [2009] T. Uezu, K. Abe, S. Miyoshi, and M. Okada, Statistical mechanical study of partial annealing of a neural network model, Journal of Physics A: Mathematical and Theoretical 43, 025004 (2009).
  * Sakata  _et al._ [2012] A. Sakata, K. Hukushima, and K. Kaneko, Replica symmetry breaking in an adiabatic spin-glass model of adaptive evolution, EPL (Europhysics Letters) 99, 68004 (2012).
  * Sherrington and Kirkpatrick [1975] D. Sherrington and S. Kirkpatrick, Solvable model of a spin-glass, Phys. Rev. Lett. 35, 1792 (1975).
  * Note [1] Note the difference from the model implementation of [14] where couplings are updated at discrete time steps.
  * Note [2] Technically speaking, the field $K$ is introduced to break the gauge-symmetry between different local Mattis states that have no frustration among target spins (i.e. triples $J_{ij}^{(tt)}J_{jk}^{(tt)}J_{ki}^{(tt)}>0$, $\forall\,i<j<k\in\mathcal{T}$). As this symmetry is broken, the ferromagnetic state (with all links between target spins being positive) is distinguished from all other unfrustrated states and, as a consequence, alignment among target spins can be achieved at low $T_{s}$. One might want to consider a different form of $K$ that explicitly takes into account the Boltzmann weights corresponding to the Hamiltonian $H_{S}$ (i.e. $e^{-\beta_{s}H_{S}}$, as defined in Eq. (4)). Such a term, however, may not guarantee the evolution towards a subgraph of only ferromagnetic interactions among target spins.
  * Note [3] Here we used the identities $\langle s_{i}s_{j}\rangle_{T_{s}}=\frac{1}{\beta_{s}}\frac{\partial}{\partial J_{ij}^{(tt)}}\mathop{ln}\nolimits Z_{1}$ and $\langle s_{k}s_{k^{\prime}}\rangle_{T_{s}}=\frac{1}{\beta_{s}}\frac{\partial}{\partial J_{kk^{\prime}}^{(oo)}}\mathop{ln}\nolimits\tilde{Z}_{1}$, where the partition function of the target spins’ subsystem and that of the non-target spins’ subsystem are $Z_{1}:=\sum@\slimits@_{\\{s_{i}\\}_{i\in\mathcal{T}}}e^{-\beta_{s}H_{\mathcal{T}}(\mathbf{J}^{(tt)})}$ and $\tilde{Z}_{1}:=\sum@\slimits@_{\\{s_{i}\\}_{i\mathchoice{\mathrel{\hbox to0.0pt{\kern 1.66667pt\kern-2.63887pt$\displaystyle\not$\hss}{\in}}}{\mathrel{\hbox to0.0pt{\kern 1.66667pt\kern-2.63887pt$\textstyle\not$\hss}{\in}}}{\mathrel{\hbox to0.0pt{\kern 1.16666pt\kern-2.05554pt$\scriptstyle\not$\hss}{\in}}}{\mathrel{\hbox to0.0pt{\kern 0.83333pt\kern-1.66664pt$\scriptscriptstyle\not$\hss}{\in}}}\mathcal{T}}}e^{-\beta_{s}H_{\mathcal{O}}(\mathbf{J}^{(oo)})}$, respectively.
  * Note [4] Requiring that the couplings should remain bounded as $\tau\rightarrow\infty$, we need to introduce a decay term $-\lambda_{ij}J_{ij}$ to the couplings dynamics. In order to keep the genotypic selection pressure the same for all the couplings regardless of their types, we choose $-\lambda_{ij}J_{ij}^{(tt)}=-\sqrt{p}J_{ij}^{(tt)}$ and $-\lambda_{ij}J_{ij}^{(oo)}=-\sqrt{1-p}J_{ij}^{(oo)}$, where $p=N_{t}/N$ is the fraction of target spins.
  * Note [5] Note that $\mathcal{Z}_{\mathcal{T}}$ and $\mathcal{Z}_{\mathcal{O}}$ are different from the partition functions of the target spins’ subsystem $Z_{1}(\mathbf{J}^{(tt)})$ and that of the non-target spins’ subsystem $\tilde{Z}_{1}(\mathbf{J}^{(oo)})$ given in [26].
  * Note [6] This can already be seen directly from Eq. (5) as the expression inside the bracket does not depend on $J_{ij}^{(oo)}$ and $J_{ij}^{(to)}$.
  * Note [7] Though $n$ appears as integer number here, we will analytically continue to real positive $n$ in subsequent steps of calculations.
  * Note [8] Here $I_{kz}=I_{kz}(m,q,Q,r,M)$ and $\tilde{I}_{k}=\tilde{I}_{k}(\tilde{m},\tilde{q},\tilde{Q},\tilde{r},\tilde{M})$ with $r,M,\tilde{r},\tilde{M}$ are other variables that the free energy densities depend on. Since these observables for target and non-target spins have their behaviour correlated to that of $(m,q,Q)$ and $(\tilde{m},\tilde{q},\tilde{Q})$, respectively, they do not provide additional information about the structure of the model phase diagram.
  * de Almeida and Thouless [1978] J. R. L. de Almeida and D. J. Thouless, Stability of the sherrington-kirkpatrick solution of a spin glass model, Journal of Physics A: Mathematical and General 11, 983 (1978).
  * Note [9] The eigenvalue of the Hessian changes its sign on the border between these phases, see SM. Note that though we do not calculate the full hierarchy of replica symmetry breaking (RSB) à la Parisi, but as long as RSB happens this non-robustness is true for the full RSB solution.
  * Note [10] This optimality, however, does not happen in the spinglass phase SP2 with low temperature $T_{s}$, since the selection pressure needed to achieve non-zero fitness crucially depends on $p$. In fact as $p$ is increased, higher pressure (i.e. lower $T_{J}$) is needed to achieve non-zero $m$ of targets, i.e, to transition from SP1 to SP2.
  * Sato  _et al._ [2003] K. Sato, Y. Ito, T. Yomo, and K. Kaneko, On the relation between fluctuation and response in biological systems, Proceedings of the National Academy of Sciences 100, 14086 (2003).
  * Kaneko and Furusawa [2006] K. Kaneko and C. Furusawa, An evolutionary relationship between genetic variation and phenotypic fluctuation, Journal of Theoretical Biology 240, 78 (2006).
  * Ciliberti  _et al._ [2007] S. Ciliberti, O. C. Martin, and A. Wagner, Robustness can evolve gradually in complex regulatory gene networks with varying topology, PLOS Computational Biology 3, 1 (2007).
  * Tang and Kaneko [2021] Q.-Y. Tang and K. Kaneko, Dynamics-evolution correspondence in protein structures, Phys. Rev. Lett. 127, 098103 (2021).
  * Landry  _et al._ [2007] C. R. Landry, B. Lemos, S. A. Rifkin, W. J. Dickinson, and D. L. Hartl, Genetic properties influencing the evolvability of gene expression, Science 317, 118 (2007).
  * Félix and Barkoulas [2015] M.-A. Félix and M. Barkoulas, Pervasive robustness in biological systems, Nature Reviews Genetics 16, 483 (2015).
  * Silva-Rocha and de Lorenzo [2010] R. Silva-Rocha and V. de Lorenzo, Noise and robustness in prokaryotic regulatory networks, Annual Review of Microbiology 64, 257 (2010).
  * Uchida  _et al._ [2022] Y. Uchida, S. Shigenobu, H. Takeda, C. Furusawa, and N. Irie, Potential contribution of intrinsic developmental stability toward body plan conservation, BMC Biology 20, 82 (2022).
  * Note [11] More precisely, in [13] $\Psi_{i}=\langle{\rm sign}\,(m)s_{i}\rangle_{T_{s}}$.
  * Note [12] From the replica-symmetric free energy density $f_{\mathcal{T}}^{\rm RS}$ we get $\chi^{\rm RS}_{m}=\beta_{s}\leavevmode@ifvmode{\hbox{$\left[\vbox to1.0pt{}\right.$}}1+(n-1)q\leavevmode@ifvmode{\hbox{$\left]\vbox to1.0pt{}\right.$}}$, while $M$ can not be obtained directly within this ansatz.
  * Lehner [2010] B. Lehner, Genes confer similar robustness to environmental, stochastic, and genetic perturbations in yeast, PLOS ONE 5, 1 (2010).
  * Go [1983] N. Go, Theoretical studies of protein folding, Annual Review of Biophysics and Bioengineering 12, 183 (1983).
  * Onuchic and Wolynes [2004] J. N. Onuchic and P. G. Wolynes, Theory of protein folding, Current Opinion in Structural Biology 14, 70 (2004).
  * Poderoso and Fontanari [2007] F. C. Poderoso and J. F. Fontanari, Model ecosystem with variable interspecies interactions, Journal of Physics A: Mathematical and Theoretical 40, 8723 (2007).
  * Pham  _et al._ [2020] M. T. Pham, I. Kondor, R. Hanel, and S. Thurner, The effect of social balance on social fragmentation, Journal of The Royal Society Interface 17, 20200752 (2020).
  * Martin  _et al._ [1973] P. C. Martin, E. D. Siggia, and H. A. Rose, Statistical dynamics of classical systems, Phys. Rev. A 8, 423 (1973).
  * De Dominicis and Peliti [1978] C. De Dominicis and L. Peliti, Field-theory renormalization and critical dynamics above ${T}_{c}$: Helium, antiferromagnets, and liquid-gas systems, Phys. Rev. B 18, 353 (1978).
  * Barbier  _et al._ [2018] M. Barbier, J.-F. Arnoldi, G. Bunin, and M. Loreau, Generic assembly patterns in complex ecological communities, Proceedings of the National Academy of Sciences 115, 2156 (2018).
  * Kadmon and Sompolinsky [2015] J. Kadmon and H. Sompolinsky, Transition to chaos in random neuronal networks, Phys. Rev. X 5, 041030 (2015).
  * Schuecker  _et al._ [2018] J. Schuecker, S. Goedeke, and M. Helias, Optimal sequence memory in driven random networks, Phys. Rev. X 8, 041029 (2018).

## Supplemental material

### A. Details of the SHK model

In the following we call the model originally introduced in [14, 15] as the
SHK model. In this model, phenotypes are spin configurations, and genotypes
are the interaction matrix for spins. In a system of $N$ spins, each spin $i$
can take values $s_{i}\in\\{-1,1\\}$ and is linked to exactly $N-1$ other
spins, thus forming a fully-connected network. Moreover, fitness is determined
by a subset of target spins denoted by $\mathcal{T}$. Those spins that do not
contribute to the fitness are called non-target. The fitness $\Psi$ at a noise
level $T_{s}$ is determined by the spin configurations at equilibrium as

$\Psi(\mathbb{s})=\frac{1}{N_{t}}\,\left\langle\Big{|}\sum_{i\in\mathcal{T}}s_{i}\Big{|}\right\rangle_{T_{s}}\,,$
(14)

where $\langle\cdot\rangle_{T_{s}}$ is the thermal average according to an
equilibrium distribution over spin configurations only. Such distribution is
computed from the partition function of a spin-glass Hamiltonian
$H_{S}=-\sum_{i<j}J_{ij}s_{i}s_{j}$ [23] in which the couplings $J_{ij}$ are
regarded as _fixed_ over the course of the spin dynamics because they are
assumed to evolve on much slower timescale than that of the spins. Here the
couplings are symmetric, i.e. $J_{ij}=J_{ji}$, and are independently and
identically distributed by a Gaussian distribution with zero mean and the
variance $J^{2}:={\rm var}(J_{ij})=N^{-1}$. The model Hamiltonian of the full
system is given by

$H_{S}=-\sum_{i<j}J_{ij}s_{i}s_{j}$ (15)

Once the spins have relaxed to an equilibrium at a temperature $T_{s}$ via a
Glauber update specified by $H_{S}$, the couplings are next updated with
probability ${\rm
Pr}\big{[}\mathbb{J}\rightarrow\mathbb{\tilde{J}}\big{]}={\rm
min}\,\big{\\{}1,e^{\beta_{J}\Delta\Psi}\big{\\}}$, where
$\Delta\Psi=\Psi(\mathbb{\tilde{J}})-\Psi(\mathbb{J})$ and $\beta_{J}\equiv
1/T_{J}$ is the genotypic selection pressure. These two dynamics are
implemented consecutively one after another until the entire system
equilibrates. Implementing this way, the model captures the evolution of
feedback process between the phenotype and genotype, where the phenotype
dynamics are represented by the stochastic dynamics of spins $(\mathbb{s})$
according to the energy landscape $H_{S}$ for given genotype $(\mathbb{J})$,
whereas the evolution of genotype is given by the stochastic change of
$(\mathbb{J})$ according to the fitness $\Psi(\mathbb{s})$ determined by the
phenotype. On the contrary to more common theories of evolution, this model
hence explicitly considers the co-evolution of these coupled landscapes.

### B. Replica symmetric ansatz solution and the expression of $I_{kz}$ and
$\tilde{I}_{k}$

The partition functions are given in terms of the _target_ and the _non-
target_ free energy densities,
$f_{\mathcal{T}}(\mathbb{m},\mathbb{q},\mathbb{r},\mathbb{Q},\mathbb{M})$ and
$f_{\mathcal{O}}(\mathbb{\tilde{m}},\mathbb{\tilde{q}},\mathbb{\tilde{r}},\mathbb{\tilde{Q}},\mathbb{\tilde{M}})$,
respectively, by

$\displaystyle\mathcal{Z}_{\mathcal{T}}=\int
D\mathbb{m}\,D\mathbb{q}\,D\mathbb{r}\,D\mathbb{Q}\,D\mathbb{M}\,e^{-\beta_{J}Npf_{\mathcal{T}}(\mathbb{m},\mathbb{q},\mathbb{r},\mathbb{Q},\mathbb{M})}$
(16a) $\displaystyle\mathcal{Z}_{\mathcal{O}}=\int
D\mathbb{\tilde{m}}\,D\mathbb{\tilde{q}}\,D\mathbb{\tilde{r}}\,D\mathbb{\tilde{Q}}\,D\mathbb{\tilde{M}}\,e^{-\beta_{J}N(1-p)f_{\mathcal{O}}(\mathbb{\tilde{m}},\mathbb{\tilde{q}},\mathbb{\tilde{r}},\mathbb{\tilde{Q}},\mathbb{\tilde{M}})}$
(16b)

where

$\displaystyle
f_{\mathcal{T}}=\frac{1}{2}\,\Big{\\{}\sum_{a<b}\frac{q_{ab}^{2}}{n^{2}}+\sum_{k<k^{\prime}}\frac{Q_{k,k^{\prime}}^{2}}{(N-t)^{2}}+\frac{1}{n(N-t)}\sum_{a,k}M_{ak}^{2}\Big{\\}}-\frac{1}{\beta_{J}}\ln\sum_{\big{\\{}s_{i};s^{\alpha};\sigma^{k}\big{\\}}_{i\in\mathcal{T}}}e^{L}$
(17a) $\displaystyle
f_{\mathcal{O}}=\frac{1}{2}\,\Big{\\{}\tilde{K}\sum_{a}\frac{\tilde{m}^{2}_{a}}{n}+\tilde{K}\sum_{k}\frac{\tilde{r}^{2}_{k}}{t}+\sum_{a<b}\frac{\tilde{q}_{ab}^{2}}{n^{2}}+\sum_{k<k^{\prime}}\frac{\tilde{Q}_{k,k^{\prime}}^{2}}{t^{2}}+\frac{1}{nt}\sum_{a,k}\tilde{M}_{ak}^{2}\Big{\\}}-\frac{1}{\beta_{J}}\ln\sum_{\big{\\{}s^{\alpha};\sigma^{k}\big{\\}}}e^{\tilde{L}}$
(17b) $\displaystyle L$ $\displaystyle=$
$\displaystyle\frac{\beta_{J}}{4N_{t}}\,\Big{|}\sum_{i\in\mathcal{T}}s_{i}\Big{|}^{2}-\frac{\beta_{J}}{2}\,\Big{|}\sum_{i\in\mathcal{T}}s_{i}\Big{|}\left[\sum_{a=1}^{n}\frac{m^{2}_{a}}{n}+\sum_{k}^{N-N_{t}}\frac{r^{2}_{k}}{N-N_{t}}\right]$
(18a) $\displaystyle+$
$\displaystyle\beta_{J}\left(\frac{1}{N_{t}}\,\Big{|}\sum_{i\in\mathcal{T}}s_{i}\Big{|}\left[\sum_{a=1}^{n}\frac{m_{a}s^{a}}{n}+\sum_{k}^{N-N_{t}}\frac{r_{k}\sigma^{k}}{N-N_{t}}\right]+\sum_{a<b}^{n}\,\frac{q_{ab}s^{a}s^{b}}{n^{2}}+\sum_{k<k^{\prime}}^{N-N_{t}}\,\frac{Q_{kk^{\prime}}\sigma^{k}\sigma^{k^{\prime}}}{(N-N_{t})^{2}}\,+\sum_{a,k}\,\frac{M_{ak}s^{a}\sigma^{k}}{n(N-N_{t})}\right)$
$\displaystyle\tilde{L}$ $\displaystyle=$
$\displaystyle\beta_{J}\left(\frac{\tilde{K}}{n}\sum_{a=1}^{n}\tilde{m}_{a}s^{a}+\frac{\tilde{K}}{N_{t}}\sum_{k}^{N_{t}}\tilde{r}_{k}\sigma^{k}+\frac{1}{n^{2}}\,\sum_{a<b}^{n}\,\tilde{q}_{ab}s^{a}s^{b}+\frac{1}{N_{t}^{2}}\,\sum_{k<k^{\prime}}^{N_{t}}\,\tilde{Q}_{kk^{\prime}}\sigma^{k}\sigma^{k^{\prime}}+\frac{1}{nN_{t}}\,\sum_{a,k}\,\tilde{M}_{ak}s^{a}\sigma^{k}\right)$
(18b)

Denoting $\mathcal{D}x\mathcal{D}y=\frac{e^{-(x^{2}+y^{2})/2}}{2\pi}\,dxdy$
and
$A_{z}(m,r)=\frac{\beta_{J}}{4}\frac{(N_{t}-2z)^{2}}{N_{t}^{2}}-\frac{\beta_{J}}{2}\frac{|N_{t}-2z|}{N_{t}}\big{(}m^{2}+r^{2}\big{)}$
we have

$\displaystyle I_{k,z}$ $\displaystyle=$
$\displaystyle\int\mathcal{D}x\mathcal{D}y\,\exp\left\\{A_{z}+\frac{N-N_{t}-2k}{N-N_{t}}\,\Big{(}y\sqrt{\beta_{J}Q}+\beta_{J}r\frac{|1-2z|}{N_{t}}\Big{)}\right\\}$
(19a) $\displaystyle\times$ $\displaystyle\left[{\rm
cosh}\Big{(}\beta_{s}m\,\frac{|N_{t}-2z|}{N_{t}}+x\frac{\sqrt{\beta_{J}q}}{n}+\frac{\beta_{J}M}{n}\,\frac{N-N_{t}-2k}{N-N_{t}}\Big{)}\right]^{n}$
$\displaystyle\tilde{I}_{k}$ $\displaystyle=$
$\displaystyle\int\mathcal{D}x\mathcal{D}y\,\exp\left\\{\frac{N_{t}-2k}{N_{t}}\,\Big{(}y\sqrt{\beta_{J}\tilde{Q}}+\beta_{J}\tilde{K}\tilde{r}\Big{)}\right\\}\,\left[{\rm
cosh}\Big{(}\beta_{s}\tilde{K}\tilde{m}+x\frac{\sqrt{\beta_{J}\tilde{q}}}{n}+\frac{\beta_{J}\tilde{M}}{n}\,\frac{N_{t}-2k}{N_{t}}\Big{)}\right]^{n}$
(19b)

The argument of the cosh$(\cdot)$ function will be denoted by

$\Omega=\beta_{s}m\,\frac{|N_{t}-2z|}{N_{t}}+x\frac{\sqrt{\beta_{J}q}}{n}+\frac{\beta_{J}M}{n}\,\frac{N-N_{t}-2k}{N-N_{t}}\,,\qquad\tilde{\Omega}=\beta_{s}\tilde{K}\tilde{m}+x\frac{\sqrt{\beta_{J}\tilde{q}}}{n}+\frac{\beta_{J}\tilde{M}}{n}\,\frac{N_{t}-2k}{N_{t}}$
(20)

The replica symmetric free energy densities given in the main text yield the
extremum condition after setting $\tilde{K}=0$

$\displaystyle m$ $\displaystyle=$
$\displaystyle\frac{\displaystyle\sum_{z=0}^{N_{t}}\binom{N_{t}}{z}\sum_{k=0}^{N-N_{t}}\binom{N-N_{t}}{k}\int\mathcal{D}x\mathcal{D}y\exp\left\\{A_{z}+\frac{N-N_{t}-2k}{N-N_{t}}\,\Big{(}y\sqrt{\beta_{J}Q}+\beta_{J}r\theta_{z}\Big{)}\right\\}\,\big{[}{\rm
cosh}(\Omega)\big{]}^{n}{\rm
tanh}(\Omega)}{\displaystyle\sum_{z=0}^{N_{t}}\binom{N_{t}}{z}\sum_{k=0}^{N-N_{t}}\binom{N-N_{t}}{k}\int\mathcal{D}x\mathcal{D}y\,\exp\left\\{A_{z}+\frac{N-N_{t}-2k}{N-N_{t}}\,\Big{(}y\sqrt{\beta_{J}Q}+\beta_{J}r\theta_{z}\Big{)}\right\\}\,\big{[}{\rm
cosh}(\Omega(x,q,M))\big{]}^{n}}$ $\displaystyle q$ $\displaystyle=$
$\displaystyle\frac{\displaystyle\sum_{z=0}^{N_{t}}\binom{N_{t}}{z}\sum_{k=0}^{N-N_{t}}\binom{N-N_{t}}{k}\int\mathcal{D}x\mathcal{D}y\,\exp\left\\{A_{z}+\frac{N-N_{t}-2k}{N-N_{t}}\,\Big{(}y\sqrt{\beta_{J}Q}+\beta_{J}r\theta_{z}\Big{)}\right\\}\,\big{[}{\rm
cosh}(\Omega)\big{]}^{n}\big{[}{\rm
tanh}(\Omega)\big{]}^{2}}{\displaystyle\sum_{z=0}^{N_{t}}\binom{N_{t}}{z}\sum_{k=0}^{N-N_{t}}\binom{N-N_{t}}{k}\int\mathcal{D}x\mathcal{D}y\,\exp\left\\{A_{z}+\frac{N-N_{t}-2k}{N-N_{t}}\,\Big{(}y\sqrt{\beta_{J}Q}+\beta_{J}r\theta_{z}\Big{)}\right\\}\,\big{[}{\rm
cosh}(\Omega)\big{]}^{n}}$ $\displaystyle Q$ $\displaystyle=$
$\displaystyle-\frac{1}{N-N_{t}}+\frac{\displaystyle\sum_{z=0}^{N_{t}}\binom{N_{t}}{z}\sum_{k=0}^{N-N_{t}}\binom{N-N_{t}}{k}\left(\frac{N-N_{t}-2k}{N-N_{t}}\right)^{2}I_{k,z}}{\displaystyle\sum_{z=0}^{N_{t}}\binom{N_{t}}{z}\sum_{k=0}^{N-N_{t}}\binom{N-N_{t}}{k}I_{k,z}}\,,\qquad
r=\frac{\displaystyle\sum_{z=0}^{N_{t}}\binom{N_{t}}{z}\sum_{k=0}^{N-N_{t}}\binom{N-N_{t}}{k}\frac{N-N_{t}-2k}{N-N_{t}}\,I_{k,z}}{\displaystyle\sum_{z=0}^{N_{t}}\binom{N_{t}}{z}\sum_{k=0}^{N-N_{t}}\binom{N-N_{t}}{k}I_{k,z}}$
$\displaystyle M$ $\displaystyle=$
$\displaystyle\frac{\displaystyle\sum_{z=0}^{N_{t}}\binom{N_{t}}{z}\sum_{k=0}^{N-N_{t}}\binom{N-N_{t}}{k}\frac{N-N_{t}-2k}{N-N_{t}}\int\mathcal{D}x\mathcal{D}y\exp\left\\{A_{z}+\frac{N-N_{t}-2k}{N-N_{t}}\,\Big{(}y\sqrt{\beta_{J}Q}+\beta_{J}r\theta_{z}\Big{)}\right\\}\,\big{[}{\rm
cosh}(\Omega\big{]}^{n}{\rm
tanh}(\Omega)}{\displaystyle\sum_{z=0}^{N_{t}}\binom{N_{t}}{z}\sum_{k=0}^{N-N_{t}}\binom{N-N_{t}}{k}\int\mathcal{D}x\mathcal{D}y\,\exp\left\\{A_{z}+\frac{N-N_{t}-2k}{N-N_{t}}\,\Big{(}y\sqrt{\beta_{J}Q}+\beta_{J}r\theta_{z}\Big{)}\right\\}\,\big{[}{\rm
cosh}(\Omega)\big{]}^{n}}$ $\displaystyle\tilde{m}$ $\displaystyle=$
$\displaystyle\frac{\displaystyle\sum_{k=0}^{N_{t}}\binom{N_{t}}{k}\int\mathcal{D}x\mathcal{D}y\exp\left\\{\frac{N_{t}-2k}{N_{t}}\,y\sqrt{\beta_{J}Q}\right\\}\,\big{[}{\rm
cosh}(\Omega(x,\tilde{m},\tilde{q},\tilde{M}))\big{]}^{n}{\rm
tanh}(\Omega(x,\tilde{m},\tilde{q},\tilde{M}))}{\displaystyle\sum_{k=0}^{N_{t}}\binom{N_{t}}{k}\int\mathcal{D}x\mathcal{D}y\,\exp\left\\{\frac{N_{t}-2k}{N_{t}}\,y\sqrt{\beta_{J}Q}\right\\}\,\big{[}{\rm
cosh}(\Omega(x,\tilde{m},\tilde{q},\tilde{M}))\big{]}^{n}}$
$\displaystyle\tilde{q}$ $\displaystyle=$
$\displaystyle\frac{\displaystyle\sum_{k=0}^{N_{t}}\binom{N_{t}}{k}\int\mathcal{D}x\mathcal{D}y\,\exp\left\\{\frac{N_{t}-2k}{N_{t}}\,y\sqrt{\beta_{J}Q}\right\\}\,\big{[}{\rm
cosh}(\Omega(x,\tilde{m},\tilde{q},\tilde{M}))\big{]}^{n}\big{[}{\rm
tanh}(\Omega(x,\tilde{m},\tilde{q},\tilde{M}))\big{]}^{2}}{\displaystyle\sum_{k=0}^{N_{t}}\binom{N_{t}}{k}\int\mathcal{D}x\mathcal{D}y\,\exp\left\\{\frac{N_{t}-2k}{N_{t}}\,y\sqrt{\beta_{J}Q}\right\\}\,\big{[}{\rm
cosh}(\Omega(x,\tilde{m},\tilde{q},\tilde{M}))\big{]}^{n}}$
$\displaystyle\tilde{Q}$ $\displaystyle=$
$\displaystyle-\frac{1}{N_{t}}+\frac{\displaystyle\sum_{k=0}^{N_{t}}\binom{N_{t}}{k}\left(\frac{N_{t}-2k}{N_{t}}\right)^{2}\tilde{I}_{k}}{\displaystyle\sum_{k=0}^{N_{t}}\binom{N_{t}}{k}\tilde{I}_{k}}\,,\qquad\tilde{r}=\frac{\displaystyle\sum_{k=0}^{N_{t}}\binom{N_{t}}{k}\frac{N_{t}-2k}{N_{t}}\,\tilde{I}_{k}}{\displaystyle\sum_{k=0}^{N_{t}}\binom{N_{t}}{k}\tilde{I}_{k}}$
$\displaystyle\tilde{M}$ $\displaystyle=$
$\displaystyle\frac{\displaystyle\sum_{k=0}^{N_{t}}\binom{N_{t}}{k}\frac{N_{t}-2k}{N_{t}}\int\mathcal{D}x\mathcal{D}y\exp\left\\{\frac{N_{t}-2k}{N_{t}}\,y\sqrt{\beta_{J}Q}\right\\}\,\big{[}{\rm
cosh}(\Omega(x,\tilde{m},\tilde{q},\tilde{M}))\big{]}^{n}{\rm
tanh}(\Omega(x,\tilde{m},\tilde{q},\tilde{M}))}{\displaystyle\sum_{k=0}^{N_{t}}\binom{N_{t}}{k}\int\mathcal{D}x\mathcal{D}y\,\exp\left\\{\frac{N_{t}-2k}{N_{t}}\,y\sqrt{\beta_{J}Q}\right\\}\,\big{[}{\rm
cosh}(\Omega(x,\tilde{m},\tilde{q},\tilde{M}))\big{]}^{n}}$

### C. Phase diagram of the order parameters of the non-target spins and that
of the eigenvalue of the Hessian

Here we support the main text with the phase diagram for the order parameters
of the non-target spins in Fig. 5 and that for the third-largest eigenvalue of
the Hessian $\Lambda_{3}$ in Fig. 6. We also depict the dependence of the drop
in fitness at $T_{c}^{(1)}$ denoted by $\Delta m=1-m$ on the fraction of non-
target spins in Fig. 7. At the critical number of target spins $N_{t}^{(c)}$
this change has a minimal value $\Delta m_{\rm min}$. Once subtracted $\Delta
m$ from $\Delta m_{\rm min}$, we find a linear relationship exists between
$\Delta m-\Delta m_{\rm min}$ and $1-Q$. Such relationship is demonstrated
Fig. 7 (C). We finally measure how much fitness changes under a transition
from robust to paramagnet phase at $T_{c}^{(2)}$. We denote this kind of drop
by $m^{*}$ in Fig. 7 (D), where we find a decrease of $m^{*}$ with increasing
$p$, implying that fitness is more robust at higher $p$.

Figure 5: Magnetisation for non-target spins $\tilde{m}$ (A). Overlap between
different replicas for non-target spins $\tilde{q}$ (B). Averaged correlation
of a pair of couplings between a target and a non-target spin that share a
common target spin $\tilde{Q}$ (C). Averaged value of the link $J^{(oo)}$
among non-target spins $\tilde{\Phi}$ (D). Here $N_{t}=10$, $N=100$. Note the
$y$-axis is on logarithmic scale. Figure 6: The third-largest eigenvalue
$\Lambda_{3}$ of the Hessian matrix for the coupling-replicas as function of
the model parameter. Here we find the region with broken replica symmetry
where $\Lambda_{3}<0$. The change in the sign of $\Lambda_{3}$ happens to
coincide with the transition between $Q=0$ and $Q>0$. Here $N_{t}=10$ and
$N=100$. Note the $y$-axis is on logarithmic scale.

Figure 7: (A). A zoom-in of Fig. 4 (A) in the region nearby $T_{c}^{(1)}$.
(B). The drop of fitness subtracted from the minimal value $\Delta m_{\rm
min}$ as a function of the fraction of non-target spins. Here $\Delta m_{\rm
min}$ is defined as the drop of magnetisation at the critical number of target
spins $N_{t}^{(c)}=60$ at $T_{J}=0.005$, while $\Delta
m=1-m\big{(}T_{c}^{(1)}\big{)}$ is the drop of of fitness at $T_{c}^{(1)}$.
(C). Linear relationship between the frustration $1-Q$ and $\Delta m-\Delta
m_{\rm min}$ at $T_{c}^{(1)}$. (D) The drop $m^{*}$ of fitness at
$T_{c}^{(2)}$. Here $N=100$.