Double-replica theory for evolution of genotype-phenotype interrelationship

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.


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 spinspin 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 → 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.

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 ∈ {−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 T .Those spins that do not contribute to the fitness are called non-target.In general, the fitness Ψ is some field that acts on J ij but whose value depends only on s i , i ∈ 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 Ψ 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 = − 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 := var(J ij ) = N −1 .The coupling matrix J includes interactions between target spins J (tt) ij for i ∈ T and j ∈ T ; those between non-target spins J (oo) ij for i ∈ T and j ∈ T ; and those between target spin and non-target spin J (to) ij for i ∈ T and j ∈ T .Let S T and S O denote the subsystem of target spins (with their interactions J (tt) ) and the subsystem of non-target spins (with the couplings J (oo) among them), respectively.The Hamiltonian of the full system denoted by S can be decomposed into where H T and H O are the Hamiltonian of the subsystems S T and S O , respectively, while H T O describes the interactions between these subsystems.
Now we introduce the effective potential to obtain the distribution of J [24].For it, we consider a continuous Langevin-type dynamics for the couplings where V = V J is the potential of all the couplings and ξ ij is the white noise whose intensity equal to the temperature T J .The factors 1/N and 1/ √ 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 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 where N t is the size of T .Under this fitness that favors the alignment of target spins, the couplings are necessarily subjected to a fitness field K [25]: or equivalently, V needs to be modified from V 0 into Without the evolution of genotypes, the Hamiltonians H T , H O and H T O dictate the spins to adapt to a set of fixed couplings 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 (tt) ij and s i s j for i, j ∈ T ; that between the state of J (oo) ij and s i s j for i, j ∈ T and that between the state of J (to) ij and s i s j for i ∈ T , j ∈ 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: This is equivalent to have a potential of the form [26] The stochastic process induced by Eq. ( 2) under this potential admits an equilibrium-like stationary joint distribution P(J (tt) , J (oo) , J (to) ) of Boltzmann-form (with associated temperature T J ), i.e., P(J (tt) , J (oo) , J (to) ) = e −β J Va /Z total , where Z total = {J} e −β J Va .Instead of calculating this distribution, we introduce our approximate approach, in which J (to) ij are assumed to always attain equilibrium well before J (tt) ij and J (oo) ij and hence can be adiabatically eliminated.As a consequence, only the weights of equilibrium configurations of J (to) contribute to the stationary distributions where P T (J (tt) , τ ) and P O (J (oo) ) are the time-dependent solutions of the corresponding Forrker-Planck equations with the effective potentials V tt for J (tt) and V oo for and J (oo) , respectively [27].To obtain the distribution only of J (tt) and J (oo) , we first replace 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 J (to) on J (tt) (J (oo) ) can be taken into account in the effective potential V tt (V oo ).The joint effect of H T O and H T on the dynamics of target spins suggests that the dependence of J (tt) ij and J (to) ik on each other arises from any triad formed between (i, j) ∈ T and k ∈ T , i.e., via jk .This influence is represented by the frustration [16] which implies that an optimal spin configuration (s jk > 0 holds (optimality in this context means that H T O and H T can be lowered simultaneously by (s * i , s * j , s * k )).Since for any given pair of (i, j) ∈ T there are N − N t triads ∆ k formed between it and a non-target spin k ∈ T , the total effect of frustration is given by Likewise, frustration among all N t triads ∆k formed between a given pair of non-target spins i ∈ T and j ∈ T The proposed scheme just allows us to define the effective potential V tt for J (tt) and V oo for and J (oo) as The stationary distributions induced by the diffusion process in Eq.
(2) with these effective potentials have a Boltzmann-form P T (J (tt) ) = e −β J Vtt /Z T and P O (J (oo) ) = e −β J Voo /Z O where Z T and Z O are the partition function of the genotypes J (tt) and that of the genotypes J (oo) , respectively [28].Here we replace J (to) by the replica matrix ik , 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 Z T and Z O as In writing these equations, we assume that the fitness acts only on J (tt) ij [29] and hence in Z O we replace the fitness field K by a constant 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 (oo) ij 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 firstorder approximation of the fitness's effect, while a term that affects the dynamics of J (oo) ij and J (to) ij is considered to be of higher order.
We here propose to interpret σ k i as the k-th replica of another variable σ i ∈ {−1, 1} that is also located at the site i of the graph (generally σ i = s i ).To distinguish these different types of replica from each other, we call s a i spin-replica and σ k i coupling-replicas.Following this interpretation, apart from n that appears as the number of spin-replicas 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 Φ .Let E[•] and Ẽ[•] denote ensemble average over P T (J (tt) ) and P O (J (oo) ), respectively.These order parameters are given by Similarly, for the non-target spins we have In the thermodynamics limit, N → ∞ and N t → ∞, while keeping p = N t /N fixed, using a replica symmetric ansatz for the variables m a = m, q ab = q and Q kk = Q, ma = m; qab = q and Qkk = Q, as well as, M ak = M and Mak = M , where , we obtain the following free energy densities: 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 Ĩk are given in the SM [31].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 T s 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 where r and r are defined and computed in the SM.

PHASE DIAGRAM
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 meanfield equations for m, q, Q, while computing Φ 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 √ 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 = J (to) J (to) starts to become non-zero.As can be anticipated from Eq. (12a), the mean value of J (tt) ij 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 RS T for small m and q, the transition between paramagnet and t-ferro occurs at T P→F s = κ, where κ = 2 −Nt Nt z=0 Nt z |N t − 2z|/N t , while the spin-glass to t-ferro transition occurs at 1 + (n − 1)q(β SP→F s ) • β SP→F s κ = 1.We also check that both the magnetisation m of the non-target spins and the average value Φ of J (oo) ij are always zero as the non-target spin subsystem remains frustrated all the time, while the spin overlap 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.
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 ≥ 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 m = q = Q = 0).Adjacent to this phase on the side of high phenotypic noise Here SP1 denotes the spin-glass phase with m = m = Q = Q = 0 and q, q > 0; SP2 denotes the spin-glass phase with m = Q = 0 and q, m > 0 (within this region, both q and Q can be either zero or nonzero, see SM); P1 denotes the paramagnet phase with m = m = q = q = Q = Q = 0; P2 denotes the paramagnet phase with m = m = q = q = Q = 0 but Q > 0; R denotes the robust fitted phase with m = q = Q = 0 but m, q, Q > 0.
Here Nt = 10, N = 100.Note the y-axis is on logarithmic scale.
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 q 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 [33].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 T and S 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 ferro- magnetic ordering is destroyed by genotypic fluctuations.
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.

STRUCTURE OF THE ROBUST FITTED PHASE
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 c , Q goes to zero, and above T (2) c , 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 ∈ [0.04, 0.5], the temperature T zoom-in of a small dip of m at this point).On the other hand, T (2) c 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 ∼ 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 nontargets.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 , the fitness m does not depend on p.It follows a unique curve, independent of p.Even though the in-crease 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 (Q) J 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 (Q) J 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 (m) J at which the fitness becomes non-zero is plotted versus p. Similar behavior is observed for others c ].The result supports p = 0.5 as the maximal value of T (m) J , suggesting the existence of an optimal fraction of target spins to acquire high fitness in this intermediate range of T s [34].

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 [9,13,[35][36][37][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 δΨ i (δJ jk ) denote the change of the average local magnetisation of a target spin i [43] upon mutating a genotype J jk → J jk + δJ jk .The mutational susceptibility of this target spin w.r.t such a change M i,jk then can be defined as In general, J jk ∈ J = J (tt) ∪ J (oo) ∪ J (to) .However, since fitness is determined solely by the configurations of target spins at equilibrium, we consider only J jk ∈ J (tt) and show that the average of this mutational susceptibility over all triples (i, j, k) ∈ T is equal to where χ m := − lim h→0 ∂ 2 f /∂h 2 is the susceptibility of target spins.The quantities M and χ 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 χ m [44].In fact, let X := lim h→0 ∂ 3 f T /∂h 3 .For L given in the SM Eq. (18a) (L/β J has the meaning of an effective Hamiltonian that is defined in the combined space {s a i , σ k i } of s-replicas and J-replicas), according to the definitions Tr (s a s b s c e L )/Tr (e L ) , Tr (s a e L )/Tr (e L ) , Tr (s a s b e L )/Tr (e L ) , the symmetry between different replicas in the robust fitted phase implies that the third moment X RS is proportional to the product of the first and second moments m RS χ RS m .Therefore, approximately, M ∝ χ RS m .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 χ m .

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 ≤ T (2) c , 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 ≥ T (1) c 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 (1) c and T 2 c .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 fluctuationresponse 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 highdimensional 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 phe-notypic 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 J (to) that couple the subsystem S T to S O on these subsystems' own dynamics by the equilibrium correlations J (to) J (to) .Such correlations are then incorporated separately into each of the dynamics of the J (tt) and 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 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 geneexpression 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].
We acknowledge support from Novo Nordisk Foundation and would like to thank Ayaka Sakata, Koji Hukushima and Yoshiyuki Kabashima for stimulating discussion.
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 Ts.One might want to consider a different form of K that explicitly takes into account the Boltzmann weights corresponding to the Hamiltonian HS (i.e. e −βsH 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.
[27] Requiring that the couplings should remain bounded as τ → ∞, we need to introduce a decay term −λijJij to the couplings dynamics.In order to keep the genotypic selection pressure the same for all the couplings regardless of their types, we choose −λijJ ij , where p = Nt/N is the fraction of target spins.
[28] Note that ZT and ZO are different from the partition functions of the target spins' subsystem Z1(J (tt) ) and that of the non-target spins' subsystem Z1(J (oo) ) given in [26].
[29] This can already be seen directly from Eq. ( 5) as the expression inside the bracket does not depend on J (oo) ij and J (to) ij .
[30] Though n appears as integer number here, we will analytically continue to real positive n in subsequent steps of calculations.
[31] Here I kz = I kz (m, q, Q, r, M ) and Ĩk = Ĩk ( m, q, Q, r, M ) with r, M, r, 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 ( m, q, Q), respectively, they do not provide additional information about the structure of the model phase diagram.

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 ∈ {−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 T .Those spins that do not contribute to the fitness are called non-target.The fitness Ψ at a noise level T s is determined by the spin configurations at equilibrium as where • Ts 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 = − 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 := var(J ij ) = N −1 .The model Hamiltonian of the full system is given by 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 Pr J → J = min 1, e β J ∆Ψ , where ∆Ψ = Ψ( J) − Ψ(J) and β J ≡ 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 (s) according to the energy landscape H S for given genotype (J), whereas the evolution of genotype is given by the stochastic change of (J) according to the fitness Ψ(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 Ĩk The partition functions are given in terms of the target and the non-target free energy densities, f T (m, q, r, Q, M) and f O ( m, q, r, Q, M), respectively, by Z T = Dm Dq Dr DQ DM e −β J N pf T (m,q,r,Q,M) (16a) where

FIG. 1 .
FIG.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 nontarget spin that share a common non-target spin Q (C).Averaged frustration among target spins Φ (D).Here Nt = 10, N = 100.Note the y-axis is on logarithmic scale.

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

FIG. 3 .
FIG. 3. The phase diagrams obtained by combining the behaviour of m as function of Ts and TJ with that of Q for different values of p, at sufficiently low TJ .Here N = 100.

FIG. 4 .
FIG. 4. (A) Magnetisation of target spins m as function of Ts at fixed TJ = 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.Ts ∈ [T (1) c , T (2) c ]) at fixed TJ = 0.005.(D) The highest value of TJ at which Q remains non-zero as function of p. (E) Magnetisation of target spins m as function of TJ at fixed Ts = 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 TJ at which the magnetisation of target spins m drops to zero at fixed Ts = 1.3 for different p from (E).Similar behaviour to (E) and (F) is observed for others Ts ∈ [T (1) c , T (2) c ].Here N = 100.
kk ln Z1, where the partition function of the target spins' subsystem and that of the non-target spins' subsystem are Z1 :=