Fundamental restriction on epistasis detection and fitness valleys in virus evolution

Epistasis detection and fitness valleys in viral evolution | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Epistasis detection and fitness valleys in viral evolution View ORCID Profile I.V. Likhachev , View ORCID Profile I.M. Rouzine doi: https://doi.org/10.1101/2025.09.05.674565 I.V. Likhachev 1 Sechenov Institute of Evolutionary Physiology and Biochemistry, Russian Academy of Sciences , St Petersburg, Russia Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for I.V. Likhachev I.M. Rouzine 1 Sechenov Institute of Evolutionary Physiology and Biochemistry, Russian Academy of Sciences , St Petersburg, Russia Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for I.M. Rouzine For correspondence: igor.rouzine{at}iephb.ru Abstract Full Text Info/History Metrics Supplementary material Data/Code Preview PDF ABSTRACT Probabilistic prognosis of virus evolution, vital for the design of effective vaccines and antiviral drugs, requires the knowledge of adaptive landscape including epistatic interactions. Although epistatic interactions can, in principle, be inferred from abundant sequencing data by various methods, it has been shown previously that genetic linkage between evolving sites obscures their signature and requires averaging over several independent populations. Here we probe the limits of epistasis detection based on pairwise correlations conditioned on the state of a third site on synthetic sequences evolved in a Monte Carlo algorithm with known epistatic pairs in two scenarios: an initially diverse population in the presence of recombination and mutation, with a uniform epistatic network, and an initially–uniform population in the presence of mutation and a fitness valley. In either scenario, the detection error decreases with the number of independent populations and increases with the sequence length. In the first scenario, the accuracy is enhanced by moderate recombination and is maximal, when epistasis magnitude approaches the point of full compensation. In the second scenario, simulation predicts the spontaneous crossing of a fitness valley and forming a new large strain.. The method is then applied to several thousands of sequences of SARS-CoV-2. Results obtained under equal sampling from world regions are consistent with the existence of fitness valleys connecting groups of viral variants. SIGNIFICANCE The few epistatic pairs of genomic sites hide in genomic data among numerous random correlations caused by common phylogenetic history. We test a method of epistasis detection designed to compensate for this noise. The accuracy is tested using synthetic sequences generated by a Monte Carlo algorithm with known epistatic pairs. The method is applied to several thousands of sequences of SARS-CoV-2 sampled in two different ways. Results obtained under equal sampling from world regions imply the existence of fitness valleys connecting groups of viral variants. INTRODUCTION Elucidating the connection between genotype and phenotype remains a challenging problem. Most traits in organisms are coded by several genomic sites, whose combined effect on the phenotype is not a sum of their effects, but also includes site-site interaction termed “epistasis”. Epistasis is responsible for a significant fraction of genetic inheritance [ 1 , 2 ]. In pathogens, epistasis facilitates the development of drug resistance and immune escape and impedes reversion of drug-resistant mutations [ 3 - 7 ]. Mutations in epitopes assisted by epistasis help a pathogen to escape the cytotoxic immune response in the host, when viral replication capacity lowered by an escape mutation is gradually restored later by compensatory mutations [ 8 - 11 ]. This process, termed fitness valley, is thought to be responsible for the high diversity and rapid HIV evolution in untreated patients [ 12 ]. Epistasis is also implicated in the evolution of variants of SARS-CoV-2 [ 13 ]. Some mutations in neutralizing antibody-binding regions that help SARS-CoV-2 to escape immune recognition during transmission [ 14 - 18 ] negatively affect the function of receptor binding, and have to be compensated elsewhere [ 19 ]. The tradeoff between antibody escape and transmissibility has been inferred in immunocompromised patients [ 20 ]. The same is likely true for mutations conferring escape from innate responses to SARS-CoV-2 [ 21 , 22 ] While pairwise epistasis can be measured from the free energy of binding [ 23 ] and in culture [ 24 ], these results cannot be transferred directly to the biological levels of an organism or population, where epistasis is defined with respect to fitness, i.e. the reproduction number at the within-host or between-host levels. Therefore, a large number of approaches has been proposed to measure epistasis directly from genomic data [ 25 - 27 ]. A popular method of “quasi-linkage equilibrium” assumes that recombination is rapid compared to allele frequency changes, and that population is large, so that recombination and the natural selection with epistasis are the dominant evolutionary factors, while the effects of genetic linkage are relatively small [ 28 - 31 ]. However, none of these methods enable reliable measurement of epistasis in asexual populations or in sexual populations at close loci [ 32 ]. About a century ago, it was realized that the evolution of a population is strongly affected by the fact that the fates of alleles at different loci are linked unless separated by recombination [ 33 ]. Empirically, the effect of strong linkage between genomic sites has been discovered by Thomas Hunt Morgan who used the inter-dependent segregation of alleles in fruit flies to understand the role of chromosomes and sex in inheritance. He showed that some pairs of traits of fruit flies segregate in a correlated fashion, and some in independent fashion. He postulated that the first are close on the same chain of genes, so that recombination cannot split them. Later, the nature of this chain was shown to be DNA in chromosomes. Since then, a large body of experiment and modeling demonstrated that linkage changes evolutionary properties drastically. In the absence of recombination, common phylogeny makes all sequences very similar to their MRCA [ 34 ]. Moreover, even very frequent recombination does not break up linkage effects in a sufficiently long genome [ 35 ]. The increasingly dominant effect of linkage over epistasis in allelic associations originates from the genetic divergence of independent populations in time. All sequences in a population are similar to their most recent common ancestor, which moves away from the origin along a stochastic trajectory [ 34 ]. As a result, any measure of co-variance produces strong noise of random sign. An attempt at detecting epistasis, whether by using covariance measures, such as Pearson coefficient, mutual entropy, universal footprint [ 36 ], or by tree-based methods [ 37 , 38 ] faces the same problem, the overwhelming genetic linkage effects. The only way to resolve this issue is to average the haplotype frequencies over many independent populations evolving under similar conditions. To decrease the required number of populations, it was proposed to use only sequences with a majority allele at a neighbor site of the measured site pair to interrupt a path along sites that creates a false-positive interaction [ 39 ]. In the present work, we test the fidelity of this detection technique using Monte-Carlo simulation in a broad range of parameters. Then, we apply the method to real virus sequences from the adapting world population of SARS–CoV-2 sampled in two different ways to predict the strongest epistatic links. This approach differs from the alternative methods in the following respects. Firstly, it does not make unrealistic assumptions regarding the evolutionary factors at play. Secondly, its accuracy is validated in a broad range of model parameters on simulated sequences, using a generic model of population genetics of the Wright Fisher type for two scenarios of evolution. We demonstrate the accuracy on simulated sequences under conditions resembling the evolution of SARS CoV-2 including the emergence of variants of concern. Thirdly, the application to data is focused on well-balanced sampling between sequences related to major variants. The fourth feature is the use of the conditional correlator. RESULTS Detection method The method is based on the calculation of Universal Footprint correlator between pairs of loci [ 36 ]. This correlator, expressed in terms of haplotype frequencies by applying the second law of thermodynamics to a multi-locus evolutionary model, depends only on epistasis magnitude and its topology, but not on any other evolutionary parameters. The method takes a set of binarized sequences as initial data. We employ parallel computing on more than 10,000 processors of GPU to increase performance. To improve detection accuracy near the point of full compensation, an additional correlator conditioned on an allele at a third locus is calculated. Based on the relationship between the two correlators, an epistatic interaction between the two loci is inferred. The details are described in Methods . Accuracy test 1. The case of standing variation To test the accuracy of the detection method, artificial genomic sequences are obtained from the simulation of the multi-site evolution of a genetically-diverse population by a Monte-Carlo algorithm with a set epistasis network ( Methods ). The input parameters are chosen in a range representative for several viruses ( Table 1 in Methods ). The predicted links are compared with the true links set in the algorithm. The error of detection is defined as the number of false-positive and false-negative interactions among predicted interactions divided by the number of true interactions. The error is calculated as a function of several parameters: genome length (total locus number), the number of independent populations used for averaging haplotype frequencies, the magnitude of epistasis, recombination probability, and the average crossover number ( Fig. 1 ). View this table: View inline View popup Download powerpoint Table 1. Parameters of the model and data processing. Download figure Open in new tab Figure 1. The error of detection of epistatic interactions for the case of standing variation (Test1). The length of the sequence is plotted on the X–axis, and the number of runs is plotted on the Y-axis of each panel. The values inside the squares defined by color map on the right correspond to the sum of the false positive and false negative interactions divided by the number of true interactions. Detection is carried out using the universal footprint of epistasis (UFE) [ 36 ] and three-site haplotype method [ 39 ] ( Methods ). All values greater than 1 are set to have the same red color as values equal to 1. Testing of the method was conducted on the sequences generated by Monte Carlo simulation ( Methods ) with variable parameters (shown): probability of recombination r , number of crossovers M , strength of the epistatic interaction E . The topology of the epistatic links (interactions) comprises isolated double arches, for which E = 0.75 is the threshold value of full compensation [ 36 ]. Fixed parameters: initial allele frequency f 0 = 0.45, population size N = 1000, mutation rate per genome μL = 0.07, simulation time in generations t = 30, selection coefficient s = 0.1. Results for M = 1 are shown in Figure S1 . Variant without recombination is shown in Figure S2 . Results using Pearson coefficient, with and without recombination, are in Figs. S3 and S4 . As expected, the accuracy of epistasis detection increases with the number of independent populations used for averaging and the magnitude of epistasis, and decreases with the genome length, as long as recombination is infrequent or absent. While modest recombination significantly improves the accuracy of detection by compensating linkage effects eclipsing epistasis, very frequent recombination counteracts associations of alleles caused by epistasis and so impairs detection ( Fig. 1 ). Analogous procedure was carried out using Pearson coefficient and UFE [ 36 ] in the absence of recombination ( Supplement , Figs. S2 - S4 ). Accuracy test 2. Initially-monomorphous populations with a fitness valley In real world, not all populations are polymorphous in the beginning of evolution. Therefore, we carried out an additional simulation test with initial allelic frequency at each locus f 0 = 0. The epistatic network was changed from the isolated double arcs to a single star comprising one deleterious allele and 10 rays with compensatory alleles. Recombination is absent. Initially, the population is monomorphous. Then a gradual accumulation of beneficial alleles occurs, which are assumed to take a 10% of loci. Deleterious alleles also accumulate at the remaining 90% of loci due to the hitchhiking effect, i.e., genetic linkage with the beneficial alleles, although more slowly than beneficial alleles. Among deleterious alleles, one is set to be the primary mutation of fitness valley and 10 are set to be compensatory mutations with the degrees of compensation ranging from 80 to 800%. These mutations are deleterious separately, but the combination of the primary allele with the compensatory alleles is advantageous. Sooner or later, deleterious primary mutation becomes abundant in a population due to the hitchhiking effect. Then compensatory alleles become advantageous, and start rapid accumulation, leading to a jump in the curve ( Fig. 2 ). This marks passing the fitness valley and the emergence of a new group of virus variants on the end of fitness valley. Download figure Open in new tab Figure 2. Monte Carlo simulation of the evolutionary dynamics of a population with a fitness valley. X-axis: Time in generations. Y axis: Average numbers of the alleles of different classes (shown) per genome. Parameters: initial allele frequency f 0 = 0, population size N = 1000, mutation rate per genome μL = 0.1,10% of the alleles are beneficial, 90% are deleterious with selection coefficient s = 0.1. Epistatic network represents a star with the primary allele in the center and ten compensatory mutations as the rays, all of which are deleterious as single alleles. Epistatic parameters are given in Methods Recombination is absent. As a result of virus passing the fitness valley, the mutational spectrum averaged over 6 populations at 150 generational steps had two maxima ( Fig. 3 ). The sequences were divided between the two peaks with borders chosen at 6 (the local minimum) and 22 new alleles (upper cutoff). Then we conducted bootstrapping procedure, in which sequences from each peak were oversampled many times with replacement and repetition. The resulting samples were used for the epistasis detection method described above. Download figure Open in new tab Figure 3. Mutational spectrum for six simulated populations, initial state is monomorphous. Y–axis: number of sequences with a number of replacements shown in X–axis. Gray lines at 6 and 22 replacements divide the two groups of sequences numbered by Roman numerals. Simulation parameters are as in Figure 2 . Six Monte Carlo runs at timepoints multiple of 10 until time 150 are used for sampling. The resulting detection error is calculated as the sum of false positive and false negative epistatic interactions out of 10 possible ( Fig. 4 ), as a function of numbers of sequences sampled from peak 1 and 2, normalized to the sequence number of each peak. Thus, about a half of weaker interactions is not detected. Download figure Open in new tab Figure 4. Error of epistatic detection for six simulated populations that are initially monomorphous. Axes X and Y show the rescaled numbers of sequences sampled from parts I and II of the dataset in Fig. 3 , respectively, relative to the total number of sequences in each part. Numbers in cells show the number of false-negative epistatic interactions out of 10 possible. The measured number of false–positives is always zero. Parameters are as in Figure 2 . The number of false positive interactions remained zero, so that the total error was given by false negatives only. Therefore this method is highly accurate, but only in the conservative direction. Our results suggest the need to take either equal sampling from the two peaks or somewhat smaller sampling from the earlier peak (Peak I), as long as they are approximately equal in size, as in our example. Results on epistasis detection in genomic data for SARS CoV-2 are consistent with fitness valleys After testing the method on simulated sequences, we applied it to sequence data on SARS-CoV-2 for two different data sets, as follows. Dataset 1 includes 6000 sequences from 6 regions, 1000 from each: Africa, Asia, North America, Oceania, South America, for the period of one year, 09/24/2020 to 09/25/2021. When each sequence is compared to Wuhan variant, a mutation spectrum is obtained, where sequences segregate naturally into three groups, with gaps between them ( Fig. 5 ).. Download figure Open in new tab Figure 5. Mutational spectrum for Dataset 1. Y–axis: number of sequences with the number of differences shown in X–axis from the consensus sequence, defined as Wuhan variant. Gray lines at 20, 440 and 60 divide the three groups of sequences numbered by the Roman numerals. Fig. 6 shows the inferred epistatic network from groups I and II. To obtain this diagram groups I and II were oversampled 10 times with replacement. As it was shown in Fig. 4 equal sampling is essential for obtaining reliable results. This procedure was also repeated for sequences mixed from groups I and III, where groups I and III were resampled 10 and 20 times respectively since group III almost two times smaller than group I ( Fig. S5 ). Download figure Open in new tab Figure 6. Epistatic interactions for SARS-CoV-2 genome predicted for Dataset 1 for sequence groups I and II ( Fig. 3 ). To obtain the circular diagram using software Circos [ 40 ], each group was oversampled 10 times with replacement. The interacting pairs are shown. The detected network ( Fig. 6 ) is consistent with the presence of fitness valleys, although does not prove their existence directly, presumably, due to the absence of weaker links not captured by the method. The method is very conservative, and its accuracy at small epistatic magnitudes is low ( Figs 1 and 4 and the legend to Fig. 4 ). Hence the epistastic interactions in Figs. 6 are the strongest in possible fitness valleys. For sequences sampled from groups II and III ( Fig. 5 ), no significant epistatic links are found, indicating the absence of a direct passage between them. The results we obtain are consistent with the existence of several fitness valleys connecting large groups of virus variants. We repeated the same procedure for Dataset 2 comprising 20,000 random sequences from the NCBI Virus database for the same period of time as Dataset 1, from 09/24/2020 to 09/25/2021, but without geographic classification nor equal sampling from continents. Results differ strongly from those for Dataset 1, as follows. The number of interactions is smaller, and they are detected only for groups I and II, but not any other pairs of the four groups ( Fig S6 ). Thus, equal sampling from world regions is critical for detection. Biological function of the predicted interacting sites Mutations E484K and N501Y are well-known mutations of the receptor-binding domain of the coronavirus in Spike. The first is responsible for the viral escape from the immune response and is capable of reducing the affinity of antibodies produced to earlier strains by 1-2 orders of magnitude, while the second improves the binding of coronavirus to the human ACE2 receptor on cell surface and increases viral replication in upper respiratory tract cells. [ 41 ] Mutations P681H and P681R, most often found in Alpha (B.1.1.7) and Delta (B.1.617.2) variants of the virus, respectively, are located in the furin cleavage region of Spike. They both confer a change in the amino acids from nonpolar hydrophobic to positively charge. They affect tropism (repertoire of infectable cells) by increasing the cleavage of S1/S2 proteins in human airway epithelial cells [ 42 ]. Mutation D138Y in Spike is believed to affect the protein’s secondary structure thus reducing antibody -binding efficiency [ 43 ]. The R203K mutation in the N protein was shown to increase transmission and virulence of SARS-COV-2 [ 44 ], while the R203M mutation, which creates a start codon for a new open reading frame, was reported to suppress interferon induction, thereby increasing viral infectivity [ 45 ]. V82A and T120I in ORF7a are most frequently found in variant Delta of the coronavirus [ 46 ]. They are located within the di-lysine motif of ORF7a, which is important for the localization of viral proteins in endoplasmic reticulum. Existing data on mutation Q7036H in protein ORF1ab and mutation V71I in ORF7a do not permit to determine their functions definitively; however, their location within viral non-structural proteins hint at their interaction with the immune system of the host. DISCUSSION Detection and measurement of epistasis from genetic data represents a challenge. In our work, we have upgraded a linkage-resistant method of detection and applied it to sequences from SARS-CoV-2. Results are consistent with the existence of fitness valleys that provide passages between large groups of virus variants. As mentioned before, a fitness valley comprises a transient decrease in fitness due to initial mutation(s) forced by changes in the environment, which is followed by compensatory mutations. The fitness valley effect was previously observed or inferred for various viruses, including HIV [ 47 , 48 ] and influenza [ 37 - 39 ], and studied using mathematical models [ 8 , 39 ]. For example, a primary mutation can help a virus to escape the immune response, adapt to a species or a tissue, or expand the pH range in which it can replicate. Mutations conferring these benefits to a virus, but reducing its overall fitness, are followed by compensatory mutations that elevate its fitness above the initial level. The loss of fitness due to escape mutations, and the partial loss of recognition by the immune system, together determine the order, the number, and the evolution rates of escape mutations in HIV infection [ 49 , 50 ]. Compensatory mutations prevent the reversion of escape mutations after the decay of CD8 T cell clones against escaped epitopes. Bacteria developing antibiotic resistance also pass through such fitness valleys [ 5 ]. A population-level simulation study proposed this effect in conjunction with the evolution in immunocompromised patients infected with SARS-CoV-2 [ 51 ]. Another modeling work proposed crossing a deep fitness valley using recombination [ 52 ]. The fates of alleles at different loci are linked due to common phylogeny unless separated by recombination [ 33 ]. For this reason, one-locus models usually do not apply in real populations, and multi-locus models must be used instead [ 53 , 54 ]. These models, based on traveling-wave technique, have been successfully applied to rapidly-evolving viruses, such as HIV [ 55 ] and influenza [ 56 - 58 ]. Linkage effects include clonal interference [ 33 , 59 ], background selection, genetic hitchhiking [ 60 ], enhanced accumulation of deleterious alleles (Muller’s ratchet) [ 61 , 62 ], and the increase of genetic drift at one locus due to selection at another [ 63 ]. Multi-site adaptation has been analyzed analytically in various asexual models, with the result that linkage effects strongly decrease the adaptation rate compared to the limit of independently-evolving sites [ 62 , 64 - 66 ]. Importantly, genetic linkage also creates random associations between pairs of mutations occurring on the same branch of the ancestral tree, which is a serious obstacle to the detection of epistatic interactions in genetic data [ 32 ]. In a single asexual population, stochastic linkage completely overshadows the epistatic footprint, except in a narrow range of times and parameters [ 34 ]. The same limitation exists for the tree-based methods of detection [ 37 , 38 ]. For testing the accuracy, we used a model with directional selection. In fact, SARS-CoV–2 evolves in antibody–binding regions due to collective immunity of population, which is a case of the selection for change (Red Queen effect). Nevertheless, as it has been shown by several groups [ 56 - 58 ], the multi-locus epidemiological model of respiratory viruses is effectively reduced to a Wright-Fisher model with constant directional selection. The selection coefficient is expressed in terms of basic reproduction ratio in a naïve population and the breadth of antibody neutralization, calculated in these works using different approximations. The most popular method for searching for epistasis signature in genomic data is the “quasi linkage equilibrium” approach [ 28 , 30 ], which assumes frequent recombination compared to changes in allelic frequency and a large population size. In this limit, genetic linkage effects are negligible compared to epistasis and recombination. In the present work, we do not rely on QLE approximation, so that the method works for finite population size and for modest recombination or even in its absence. The linkage effects are mostly cancelled out by averaging over well-balanced independent populations and the use of 3-site haplotypes. Such cancellation was not carried out in [ 68 ], hence the abundance of false negative links ( Fig. 7C ) compared to our results ( Fig. 7A ). Download figure Open in new tab Figure 7. Epistatic interactions for SARS-CoV-2 genome: comparison of the present results with two previous studies. A) Results of the present work for Dataset 1 from Fig. 6 obtained from mixing sequence groups I and II in Fig. 5. B, C ) Results from studies [ 67 ] and [ 68 ], respectively. Blue stripes in panel B represent protein domains reported to interact with each other in [ 67 ]. To address the problem of false-negative bonds, another team [ 67 ] compared the sequences of SARS CoV-2 not between themselves, but with other coronaviruses. They used “direct coupling analysis”, which postulates a statistical distribution of sequences (usually, of a Gaussian form) including individual effects of alleles and their pairwise interaction. A physicist would recognize here the Gibbs distribution for the Onsager model of ferromagnetism. In population genetics, this Gaussian form has a justification in the limit of QLE [ 28 , 30 ]. QLE at the base does not work in SARS CoV-2 evolution, because linkage effects, negligible in QLE limit, are quite strong in real life. Parameters of the Gaussian are fit to sequence data using various approximations [ 67 ]. DCA also includes down-weighing sequences with over 80% identity, which is a phenomenological trick to deal with the effects of random linkage. As result, the number of false-positive epistatic bonds due to genetic linkage is sharply decreased ( Fig. 7B ). The problem here is that the nature of such correction is not derived from a model of population genetics. Note that the use of tree-based methods does not automatically compensate for strong linkage effects. This is because, in the absence of sufficiently strong recombination, all sequences in one single population are very similar to their MRCA, and MRCAs for different independent populations are diverging from each other in genetic space. This process was demonstrated to be responsible for strong linkage effects [ 34 ]. To decrease the effects of linkage, one group used synonymous sites as markers of linkage [ 37 , 38 ]. To conclude, the present study confirms the existence of fundamental limitations on epistasis detection due to genetic linkage. Within this limitation, a conditional correlation method combined with balanced sampling from quasi-independent populations and testing the accuracy in simulation allows the testable detection of strong epistatic links. METHODS Generating sequences for method testing Evolution of a haploid asexual population of N binary sequences is simulated by a Monte Carlo code [ 36 , 69 ], as follows. In an individual genome, each locus (site, nucleotide position, amino acid position) numbered i = 1, 2, …, L is occupied by one of two alleles, either the wild-type allele, denoted K i = 0, or the mutant allele, K i = 1.We use a discrete generation scheme in the absence of generation overlap (Wright-Fisher model). The evolutionary factors included in the model are random mutation with rate μL per genome, constant directional selection, random genetic drift due to random sampling of progeny and recombination with rate r . Selection includes an epistatic network with a set strength and topology. The logarithm of the average progeny number of a genome, W , depends on sequence [ K i ] as given by [ 36 ] where T ij = 0 or 1 is the binary matrix that shows epistatically-interacting pairs. Here the selection coefficients s i and s j denote the individual fitness costs of two deleterious mutations, E ij is the compensation degree at sites i and j . In our simulation example in Fig. 1 , we consider a haploid population with the initial frequency of deleterious alleles f 0 = 0.45, mutation rate per genome μL = 0.07, fixed selection coefficient s i ≡ s = 0.1, population size N = 1000 individuals, genome length L = 40 sites, and epistatic network is double arches In our simulation example in Figs. 2 to 4, we consider a haploid population with the initial frequency of deleterious alleles f 0 = 0, mutation rate per genome μL = 0.1, selection coefficients with a fixed absolute value s i ≡ s = 0.1, where 10% of potential mutations are advantageous and 90% are deleterious, population size N = 1000 individuals and genome length L = 100 sites, epistatic network – star with deleterious mutation in the middle and 10 compensatory mutations as rays, degree of compensation gradually increases from 80% to 800% The Monte-Carlo simulation code is written in MATLAB™ and deposited at https://github.com/FireFrog17/Fundamental-restriction-on-epistasis-detection-and-fitness-valleys-in-virus-evolution-Data Processing genomic data of SARS-CoV-2 In Dataset 1, 6000 sequences sampled from 6 regions (Africa, Asia, North America, Oceania, South America) in the period from 09/24/2020 to 09/25/2021 were downloaded from database NCBI Virus. In Dataset 2, we used 20,000 random sequences from the NCBI Virus database for the same period of time as in Dataset 1. The sequences were not stratified by sampling date for two reasons, as follows. Firstly, the necessary size of statistics would be missing. Also, previous experience with hemagglutinin data [ 39 ] demonstrates that combining sequences from different dates is an advantage, because one compares earlier and later strains and thus finds a fitness valley connecting them. In contrast, estimation of selection coefficients, indeed, requires classification by time [ 70 ]. Downloaded data were processed by a custom-made sequence-processing program ProcessGenomes . py . Multiple alignment was performed in MAFT online, using the “Fast calculation of full-length MSA of closely related viral genomes” [ 71 ]. From the total number of downloaded and aligned sequences, sequences in which the percentage of symbols “–” (gap) or “N” (nonread) is greater than Sequence threshold are removed. From the remaining sequences, sites in which the percentage of gaps or letters N in across sequences is higher than Gap threshold are removed. Then, sites with a polymorphism degree (Hamming distance) lower than Monomorphism threshold are removed. We put the monomorphism threshold at two values ( Table 1 ) and make sure that the choice does not affect final results, i.e., results are in a plateau in this parameter. The sequence data are binarized: the variant similar to consensus defined as Wuhan variant is declared 0, the others are 1, with the resulting population matrix K ni . After binarization, the sequence data of virus undergo an additional processing stage, as follows. Based on the number of differences from sample consensus, a mutational spectrum is constructed using a custom-made sequence processing program ProcessGenomes . py ( Figs. 3 and 5 ). Several groups of sequences are identified from this diagram. The new matrix K ni is formed by random sampling from two chosen groups, with the possibility of repetition.. The dependence of the predicted number of epistatic links on the sample size for simulated data is shown in Figs. 1 and 4 . Epistasis detection and measurement Based on matrix K ni obtained as described above, the universal footprint of epistasis [ 36 ] denoted UFE ij is calculated for each possible pair of sites i, j as . Then, the conditional coefficient UFE ij 0 is re-calculated for each locus pair [ 39 ], as follows. Only sequences where at least one of the neighboring loci of the pair has allele 0 are used. The minimum of all UFE over such neighbors is denoted UFE ij 0 . Because the array of data used is large, these coefficients are calculated using parallel computing on the video card with more than 1000 parallel processors with the help of CUDA language. Pairs of loci, where the ratio of UFE ij to UFE ij 0 is equal to or greater than 0.5, are considered “true pairs”, while the remaining pairs are excluded. Results are displayed in the form of circular diagrams ( Figs. 6 and S5). Error calculation . The estimated epistatic pairs are compared with the matrix of true pairs T ij set in the Monte Carlo simulation code that generates genomic sequences for testing the accuracy. The relative error of the method is calculated as the sum of the number of false positive and false negative pairs for an initially–uniform population ( Fig. 4 ). For the case of standing variation, the sum is divided by the total number of true pairs ( Fig. 1 ). Assumptions and approximations The method implies the weak selection limit, defined as max(1/ N, μ ) ≪ s ≪ 1, and is appropriate in the travelling wave regime, when many mutations at different loci are being fixed gradually and compete with each other for the space in population (clonal interference). The traveling wave regime takes place at NμL ≫ 1 and on time scales longer than 1/ s . Next, the fluctuations of mutation rate along the genome, which are typically strong for RNA viruses, do not matter much in the multi-site evolution regime, where the mutation rate enters all the observables in the argument of a large logarithm [ 53 ]. Recombination, although observed for example for Omicron variants, is relatively rare in SARS-CoV-2, because it requires coinfection of an individual with two genetically remote variants, as estimated before for another virus [ 69 , 72 ]. Multiplicity of infection is incorporated into an effective recombination rate measurable directly from genomic data [ 69 , 73 ]. The effect of recombination on detection accuracy was tested an initially-diverse population and found to be modest ( Fig. 1 ). Also, here we consider the population level evolution only, because SARS–CoV-2 is acute and represented by a single variant in 99,5% patients, similar to influenza. In principle, the 0.5% chronic patients could seed new variants, in which case inter-scale evolutionary conflicts are possible [ 74 ]. Funding The study was funded by Russian Science Foundation, grant 24-24-00529. Conflict of interests Authors declare no conflict of interest. Ethical statement No ethical statement is required. Availability of data and software Codes are available at https://github.com/FireFrog17/Fundamental-restriction-on-epistasis-detection-and-fitness-valleys-in-virus-evolution-Data Supplementary material Download figure Open in new tab Figure S1. Accuracy of detection of epistatic pairs from UFE correlator for M = 1. Simulation parameters and notation are as in Figure 1 . Download figure Open in new tab Figure S2. The error of detection of epistatic interactions from UFE correlator in the absence of recombination. Simulation parameters and notation as in Figure 1 . Download figure Open in new tab Figure S3. The error of detection of epistatic interactions from Pearson coefficient. Simulation parameters and notation as in Figure 1 . Download figure Open in new tab Figure S4. The error of detection of epistatic interactions in the absence of recombination from Pearson coefficient. Simulation parameters and notation are as in Figure 1 . Download figure Open in new tab Figure S5. Epistatic interactions for SARS-CoV-2 genome predicted for Dataset 1 between groups I and III in Fig. 5 . Notation as in Fig. 6 . To obtain the circular diagram using software Circos [ 40 ], group I was oversampled 10 times and group III was oversampled 20 times, both with replacement. Download figure Open in new tab Figure S6. Epistatic interactions for SARS-CoV-2 genome predicted for Dataset 2 from UFE correlator. Dataset 2 was constructed using 20,000 random sequences from the NCBI without geographic classification. Virus database for the period from 09/24/2020 to 09/25/2021. Calculations were based on equal mixture of sequences from peak I and II. Acknowledgements We express our gratitude to Richard Neher for stimulating comments. We thank A.A. Makashov and S.V. Lagunov for technical help and discussion. We are grateful to D. A. Gorshkova, A. I. Pikhulya, E. I. Stroganova, and E. A. Vigovskaya for technical help. Funder Information Declared Russian Science Foundation, https://ror.org/03y2gwe85 , 24-24-00529 Footnotes Manuscript content is the same as in previous revision. In this version we added reply to the reviewers as a separate file. https://github.com/FireFrog17/Fundamental-restriction-on-epistasis-detection-and-fitness-valleys-in-virus-evolution-Data REFERENCES 1. ↵ Carlborg O , Jacobsson L , Ahgren P , Siegel P , Andersson L. Epistasis and the release of genetic variation during long-term selection . Nat Genet . 2006 ; 38 ( 4 ): 418 – 20 . Epub 2006/03/15. doi: 10.1038/ng1761 . PubMed PMID: 16532011 . OpenUrl CrossRef PubMed Web of Science 2. ↵ Zuk O , Hechter E , Sunyaev SR , Lander ES . The mystery of missing heritability: Genetic interactions create phantom heritability . Proc Natl Acad Sci U S A . 2012 ; 109 ( 4 ): 1193 – 8 . Epub 2012/01/10. doi: 10.1073/pnas.1119675109 . PubMed PMID: 22223662 ; PubMed Central PMCID: PMCPMC3268279 . OpenUrl Abstract / FREE Full Text 3. ↵ Gonzalez-Ortega E , Ballana E , Badia R , Clotet B , Este JA . Compensatory mutations rescue the virus replicative capacity of VIRIP-resistant HIV-1 . Antiviral Res . 2011 ; 92 ( 3 ): 479 – 83 . Epub 2011/10/27. doi: 10.1016/j.antiviral.2011.10.010 . PubMed PMID: 22027647 . OpenUrl CrossRef PubMed 4. Handel A , Regoes RR , Antia R. The role of compensatory mutations in the emergence of drug resistance . PLoS Comput Biol . 2006 ; 2 ( 10 ): e137 . Epub 2006/10/17. doi: 10.1371/journal.pcbi.0020137 . PubMed PMID: 17040124 ; PubMed Central PMCID: PMCPMC1599768 . OpenUrl CrossRef PubMed 5. ↵ Levin BR , Perrot V , Walker N. Compensatory mutations, antibiotic resistance and the population genetics of adaptive evolution in bacteria . Genetics . 2000 ; 154 ( 3 ): 985 – 97 . Epub 2000/04/11. PubMed PMID: 10757748 ; PubMed Central PMCID: PMCPMC1460977 . OpenUrl Abstract / FREE Full Text 6. Nijhuis M , Schuurman R , de Jong D , Erickson J , Gustchina E , Albert J , et al. Increased fitness of drug resistant HIV-1 protease as a result of acquisition of compensatory mutations during suboptimal therapy . AIDS . 1999 ; 13 ( 17 ): 2349 – 59 . Epub 1999/12/22. PubMed PMID: 10597776 . OpenUrl CrossRef PubMed Web of Science 7. ↵ Piana S , Carloni P , Rothlisberger U. Drug resistance in HIV-1 protease: Flexibility-assisted mechanism of compensatory mutations . Protein Sci . 2002 ; 11 ( 10 ): 2393 – 402 . Epub 2002/09/19. doi: 10.1110/ps.0206702 . PubMed PMID: 12237461 ; PubMed Central PMCID: PMCPMC2384161 . OpenUrl CrossRef PubMed Web of Science 8. ↵ Weissman DB , Desai MM , Fisher DS , Feldman MW . The rate at which asexual populations cross fitness valleys . Theor Popul Biol . 2009 ; 75 ( 4 ): 286 – 300 . Epub 2009/03/17. doi: 10.1016/j.tpb.2009.02.006 . PubMed PMID: 19285994 ; PubMed Central PMCID: PMCPMC2992471 . OpenUrl CrossRef PubMed Web of Science 9. Liu D , Zuo T , Hora B , Song H , Kong W , Yu X , et al. Preexisting compensatory amino acids compromise fitness costs of a HIV-1 T cell escape mutation . Retrovirology . 2014 ; 11 : 101 . Epub 20141119. doi: 10.1186/s12977-014-0101-0 . PubMed PMID: 25407514 ; PubMed Central PMCID: PMCPMC4264250 . OpenUrl CrossRef PubMed 10. Crawford H , Prado JG , Leslie A , Hue S , Honeyborne I , Reddy S , et al. Compensatory mutation partially restores fitness and delays reversion of escape mutation within the immunodominant HLA-B*5703-restricted Gag epitope in chronic human immunodeficiency virus type 1 infection . J Virol . 2007 ; 81 ( 15 ): 8346 – 51 . Epub 20070516. doi: 10.1128/JVI.00465-07 . PubMed PMID: 17507468 ; PubMed Central PMCID: PMCPMC1951305 . OpenUrl Abstract / FREE Full Text 11. ↵ Song H , Pavlicek JW , Cai F , Bhattacharya T , Li H , Iyer SS , et al. Impact of immune escape mutations on HIV-1 fitness in the context of the cognate transmitted/founder genome . Retrovirology . 2012 ; 9 : 89 . Epub 20121030. doi: 10.1186/1742-4690-9-89 . PubMed PMID: 23110705 ; PubMed Central PMCID: PMCPMC3496648 . OpenUrl CrossRef PubMed 12. ↵ Rouzine IM , Coffin JM . Search for the mechanism of genetic variation in the pro gene of human immunodeficiency virus . J Virol . 1999 ; 73 ( 10 ): 8167 – 78 . Epub 1999/09/11. doi: 10.1128/JVI.73.10.8167-8178.1999 . PubMed PMID: 10482567 ; PubMed Central PMCID: PMCPMC112834 . OpenUrl Abstract / FREE Full Text 13. ↵ Rouzine IM . Evolutionary Mechanisms of the Emergence of the Variants of Concern of SARS-CoV-2 . Viruses . 2025 ; 17 ( 2 ). Epub 20250130. doi: 10.3390/v17020197 . PubMed PMID: 40006952 ; PubMed Central PMCID: PMCPMC11861269 . OpenUrl CrossRef PubMed 14. ↵ Carabelli AM , Peacock TP , Thorne LG , Harvey WT , Hughes J , Consortium C-GU , et al. SARS-CoV-2 variant biology: immune escape, transmission and fitness . Nat Rev Microbiol . 2023 ; 21 ( 3 ): 162 – 77 . Epub 20230118. doi: 10.1038/s41579-022-00841-7 . PubMed PMID: 36653446 ; PubMed Central PMCID: PMCPMC9847462 . OpenUrl CrossRef PubMed 15. Shrestha LB , Foster C , Rawlinson W , Tedla N , Bull RA . Evolution of the SARS-CoV-2 omicron variants BA.1 to BA.5: Implications for immune escape and transmission . Rev Med Virol . 2022 ; 32 ( 5 ): e2381 . Epub 20220720. doi: 10.1002/rmv.2381 . PubMed PMID: 35856385 ; PubMed Central PMCID: PMCPMC9349777 . OpenUrl CrossRef PubMed 16. Cobey S , Larremore DB , Grad YH , Lipsitch M. Concerns about SARS-CoV-2 evolution should not hold back efforts to expand vaccination . Nat Rev Immunol . 2021 ; 21 ( 5 ): 330 – 5 . Epub 2021/04/03. doi: 10.1038/s41577-021-00544-9 . PubMed PMID: 33795856 ; PubMed Central PMCID: PMCPMC8014893 . OpenUrl CrossRef PubMed 17. Yewdell JW . Antigenic drift: Understanding COVID-19 . Immunity . 2021 ; 54 ( 12 ): 2681 – 7 . Epub 2021/12/16. doi: 10.1016/j.immuni.2021.11.016 . PubMed PMID: 34910934 ; PubMed Central PMCID: PMCPMC8669911 . OpenUrl CrossRef PubMed 18. ↵ Rouzine IM , Rozhnova G. Evolutionary implications of SARS-CoV-2 vaccination for the future design of vaccination strategies . Commun Med (Lond) . 2023 ; 3 ( 1 ): 86 . Epub 20230619. doi: 10.1038/s43856-023-00320-x . PubMed PMID: 37336956 ; PubMed Central PMCID: PMCPMC10279745 . OpenUrl CrossRef PubMed 19. ↵ Moulana A , Dupic T , Phillips AM , Chang J , Nieves S , Roffler AA , et al. Compensatory epistasis maintains ACE2 affinity in SARS-CoV-2 Omicron BA.1 . Nat Commun . 2022 ; 13 ( 1 ): 7011 . Epub 20221116. doi: 10.1038/s41467-022-34506-z . PubMed PMID: 36384919 ; PubMed Central PMCID: PMCPMC9668218 . OpenUrl CrossRef PubMed 20. ↵ Harari S , Tahor M , Rutsinsky N , Meijer S , Miller D , Henig O , et al. Drivers of adaptive evolution during chronic SARS-CoV-2 infections . Nat Med . 2022 ; 28 ( 7 ): 1501 – 8 . Epub 20220620. doi: 10.1038/s41591-022-01882-4 . PubMed PMID: 35725921 ; PubMed Central PMCID: PMCPMC9307477 . OpenUrl CrossRef PubMed 21. ↵ Reuschl AK , Thorne LG , Whelan MVX , Ragazzini R , Furnon W , Cowton VM , et al. Evolution of enhanced innate immune suppression by SARS-CoV-2 Omicron subvariants . Nat Microbiol . 2024 ; 9 ( 2 ): 451 – 63 . Epub 20240116. doi: 10.1038/s41564-023-01588-4 . PubMed PMID: 38228858 ; PubMed Central PMCID: PMCPMC10847042 . OpenUrl CrossRef PubMed 22. ↵ Thorne LG , Bouhaddou M , Reuschl AK , Zuliani-Alvarez L , Polacco B , Pelin A , et al. Evolution of enhanced innate immune evasion by SARS-CoV-2 . Nature . 2022 ; 602 ( 7897 ): 487 – 95 . Epub 20211223. doi: 10.1038/s41586-021-04352-y . PubMed PMID: 34942634 ; PubMed Central PMCID: PMCPMC8850198 . OpenUrl CrossRef PubMed 23. ↵ Adams RM , Kinney JB , Walczak AM , Mora T. Epistasis in a Fitness Landscape Defined by Antibody-Antigen Binding Free Energy . Cell Syst . 2019 ; 8 ( 1 ): 86 – 93 e3 . Epub 20190102. doi: 10.1016/j.cels.2018.12.004 . PubMed PMID: 30611676 ; PubMed Central PMCID: PMCPMC6487650 . OpenUrl CrossRef PubMed 24. ↵ Lyons DM , Zou Z , Xu H , Zhang J. Idiosyncratic epistasis creates universals in mutational effects and evolutionary trajectories . Nat Ecol Evol . 2020 ; 4 ( 12 ): 1685 – 93 . Epub 20200907. doi: 10.1038/s41559-020-01286-y . PubMed PMID: 32895516 ; PubMed Central PMCID: PMCPMC7710555 . OpenUrl CrossRef PubMed 25. ↵ Chen CC , Schwender H , Keith J , Nunkesser R , Mengersen K , Macrossan P. Methods for identifying SNP interactions: a review on variations of Logic Regression, Random Forest and Bayesian logistic regression . IEEE/ACM Trans Comput Biol Bioinform . 2011 ; 8 ( 6 ): 1580 – 91 . Epub 2011/03/09. doi: 10.1109/TCBB.2011.46 . PubMed PMID: 21383421 . OpenUrl CrossRef PubMed 26. Cordell HJ . Detecting gene-gene interactions that underlie human diseases . Nature reviews Genetics . 2009 ; 10 ( 6 ): 392 – 404 . doi: 10.1038/nrg2579 . PubMed PMID : pmcid: PMC2872761 . OpenUrl CrossRef PubMed Web of Science 27. ↵ Ueki M , Cordell HJ . Improved statistics for genome-wide interaction analysis . PLoS Genet . 2012 ; 8 ( 4 ): e1002625 . Epub 2012/04/13. doi: 10.1371/journal.pgen.1002625 . PubMed PMID: 22496670 ; PubMed Central PMCID: PMCPMC3320596 . OpenUrl CrossRef PubMed 28. ↵ Kimura M. Attainment of Quasi Linkage Equilibrium When Gene Frequencies Are Changing by Natural Selection . Genetics . 1965 ; 52 ( 5 ): 875 – 90 . doi: 10.1093/genetics/52.5.875 . PubMed PMID: 17248281 ; PubMed Central PMCID: PMCPMC1210959 . OpenUrl FREE Full Text 29. Weigt M , White RA , Szurmant H , Hoch JA , Hwa T. Identification of direct residue contacts in protein-protein interaction by message passing . Proc Natl Acad Sci U S A . 2009 ; 106 ( 1 ): 67 – 72 . Epub 2009/01/01. doi: 10.1073/pnas.0805923106 . PubMed PMID: 19116270 ; PubMed Central PMCID: PMCPMC2629192 . OpenUrl Abstract / FREE Full Text 30. ↵ Neher RA , Shraiman BI . Statistical genetics and evolution of quantitative traits . Reviews of Modern Physics . 2011 ; 83 : 1283 . OpenUrl CrossRef 31. ↵ Louie RHY , Kaczorowski KJ , Barton JP , Chakraborty AK , McKay MR . Fitness landscape of the human immunodeficiency virus envelope protein that is targeted by antibodies . Proc Natl Acad Sci U S A . 2018 ; 115 ( 4 ): E564 – E73 . Epub 20180108. doi: 10.1073/pnas.1717765115 . PubMed PMID: 29311326 ; PubMed Central PMCID: PMCPMC5789945 . OpenUrl Abstract / FREE Full Text 32. ↵ Wei W-H , Hemani G , Haley CS . Detecting epistasis in human complex traits . Nature Reviews Genetics . 2014 ; 15 : 722 . doi: 10.1038/nrg3747 https://www.nature.com/articles/nrg3747#supplementary-information . OpenUrl CrossRef PubMed 33. ↵ Fisher RA . The genetical theory of natural selection . Oxford, United Kingdom : Clarendon Press , 1958 ; 1930. 34. ↵ Pedruzzi G , Rouzine IM . Epistasis detectably alters correlations between genomic sites in a narrow parameter window . PLoS One . 2019 ; 14 ( 5 ): e0214036 . Epub 2019/06/01. doi: 10.1371/journal.pone.0214036 . PubMed PMID: 31150393 ; PubMed Central PMCID: PMCPMC6544209 . OpenUrl CrossRef PubMed 35. ↵ Rouzine IM . Long-range linkage effects in adapting sexual populations . Sci Rep . 2023 ; 13 ( 1 ): 12492 . Epub 20230801. doi: 10.1038/s41598-023-39392-z . PubMed PMID: 37528175 ; PubMed Central PMCID: PMCPMC10393966 . OpenUrl CrossRef PubMed 36. ↵ Pedruzzi G , Barlukova A , Rouzine IM . Evolutionary footprint of epistasis . PLoS Comput Biol . 2018 ; 14 ( 9 ): e1006426 . Epub 2018/09/18. doi: 10.1371/journal.pcbi.1006426 . PubMed PMID: 30222748 ; PubMed Central PMCID: PMCPMC6177197 . OpenUrl CrossRef PubMed 37. ↵ Kryazhimskiy S , Dushoff J , Bazykin GA , Plotkin JB . Prevalence of epistasis in the evolution of influenza A surface proteins . PLoS Genet . 2011 ; 7 ( 2 ): e1001301 . Epub 20110217. doi: 10.1371/journal.pgen.1001301 . PubMed PMID: 21390205 ; PubMed Central PMCID: PMCPMC3040651 . OpenUrl CrossRef PubMed 38. ↵ Neverov AD , Kryazhimskiy S , Plotkin JB , Bazykin GA . Coordinated Evolution of Influenza A Surface Proteins . PLoS Genet . 2015 ; 11 ( 8 ): e1005404 . Epub 20150806. doi: 10.1371/journal.pgen.1005404 . PubMed PMID: 26247472 ; PubMed Central PMCID: PMCPMC4527594 . OpenUrl CrossRef PubMed 39. ↵ Pedruzzi G , Rouzine IM . An evolution-based high-fidelity method of epistasis measurement: Theory and application to influenza . PLoS Pathog . 2021 ; 17 ( 6 ): e1009669 . Epub 2021/06/22. doi: 10.1371/journal.ppat.1009669 . PubMed PMID: 34153082 ; PubMed Central PMCID: PMCPMC8248644 . OpenUrl CrossRef PubMed 40. ↵ Krzywinski M , Schein J , Birol I , Connors J , Gascoyne R , Horsman D , et al. Circos: an information aesthetic for comparative genomics . Genome Res . 2009 ; 19 ( 9 ): 1639 – 45 . Epub 20090618. doi: 10.1101/gr.092759.109 . PubMed PMID: 19541911 ; PubMed Central PMCID: PMCPMC2752132 . OpenUrl Abstract / FREE Full Text 41. ↵ Tao K , Tzou PL , Nouhin J , Gupta RK , de Oliveira T , Kosakovsky Pond SL , et al. The biological and clinical significance of emerging SARS-CoV-2 variants . Nat Rev Genet . 2021 ; 22 ( 12 ): 757 – 73 . Epub 20210917. doi: 10.1038/s41576-021-00408-x . PubMed PMID: 34535792 ; PubMed Central PMCID: PMCPMC8447121 . OpenUrl CrossRef PubMed 42. ↵ Johnson BA , Xie X , Bailey AL , Kalveram B , Lokugamage KG , Muruato A , et al. Loss of furin cleavage site attenuates SARS-CoV-2 pathogenesis . Nature . 2021 ; 591 ( 7849 ): 293 – 9 . Epub 20210125. doi: 10.1038/s41586-021-03237-4 . PubMed PMID: 33494095 ; PubMed Central PMCID: PMCPMC8175039 . OpenUrl CrossRef PubMed 43. ↵ Kanakan A , Mehta P , Devi P , Saifi S , Swaminathan A , Maurya R , et al. Clinico-Genomic Analysis Reiterates Mild Symptoms Post-vaccination Breakthrough: Should We Focus on Low-Frequency Mutations? Front Microbiol . 2022 ; 13 : 763169 . Epub 20220303. doi: 10.3389/fmicb.2022.763169 . PubMed PMID: 35308382 ; PubMed Central PMCID: PMCPMC8927057 . OpenUrl CrossRef PubMed 44. ↵ Wu H , Xing N , Meng K , Fu B , Xue W , Dong P , et al. Nucleocapsid mutations R203K/G204R increase the infectivity, fitness, and virulence of SARS-CoV-2 . Cell Host Microbe . 2021 ; 29 ( 12 ): 1788 – 801 e6 . Epub 20211113. doi: 10.1016/j.chom.2021.11.005 . PubMed PMID: 34822776 ; PubMed Central PMCID: PMCPMC8590493 . OpenUrl CrossRef PubMed 45. ↵ Li Y , Li M , Xiao H , Liao F , Shen M , Ge W , et al. The R203M and D377Y mutations of the nucleocapsid protein promote SARS-CoV-2 infectivity by impairing RIG-I-mediated antiviral signaling . PLoS Pathog . 2025 ; 21 ( 1 ): e1012886 . Epub 20250122. doi: 10.1371/journal.ppat.1012886 . PubMed PMID: 39841800 ; PubMed Central PMCID: PMCPMC11771877 . OpenUrl CrossRef PubMed 46. ↵ Cruz CAK , Medina PMB . Temporal changes in the accessory protein mutations of SARS-CoV-2 variants and their predicted structural and functional effects . J Med Virol . 2022 ; 94 ( 11 ): 5189 – 200 . Epub 20220712. doi: 10.1002/jmv.27964 . PubMed PMID: 35764775 ; PubMed Central PMCID: PMCPMC9349927 . OpenUrl CrossRef PubMed 47. ↵ Bonhoeffer S , Chappey C , Parkin NT , Whitcomb JM , Petropoulos CJ . Evidence for positive epistasis in HIV-1 . Science . 2004 ; 306 ( 5701 ): 1547 – 50 . doi: 10.1126/science.1101786 . PubMed PMID: 15567861 . OpenUrl Abstract / FREE Full Text 48. ↵ Zhang TH , Dai L , Barton JP , D. Y , Tan Y , Pang W , et al. Predominance of positive epistasis among drug resistance-associated mutations in HIV-1 protease . PLoS Genet . 2020 ; 16 ( 10 ): e1009009 . Epub 20201021. doi: 10.1371/journal.pgen.1009009 . PubMed PMID: 33085662 ; PubMed Central PMCID: PMCPMC7605711 . OpenUrl CrossRef PubMed 49. ↵ Althaus CL , De Boer RJ . Dynamics of immune escape during HIV/SIV infection . PLoS Comput Biol . 2008 ; 4 ( 7 ): e1000103 . Epub 20080718. doi: 10.1371/journal.pcbi.1000103 . PubMed PMID: 18636096 ; PubMed Central PMCID: PMCPMC2423483 . OpenUrl CrossRef PubMed 50. ↵ Batorsky R , Sergeev RA , Rouzine IM . The route of HIV escape from immune response targeting multiple sites is determined by the cost-benefit tradeoff of escape mutations . PLoS Comput Biol . 2014 ; 10 ( 10 ): e1003878 . Epub 2014/10/31. doi: 10.1371/journal.pcbi.1003878 . PubMed PMID: 25356981 ; PubMed Central PMCID: PMCPMC4214571 . OpenUrl CrossRef PubMed 51. ↵ Smith CA , Ashby B. Antigenic evolution of SARS-CoV-2 in immunocompromised hosts . Evol Med Public Health . 2023 ; 11 ( 1 ): 90 – 100 . Epub 20221111. doi: 10.1093/emph/eoac037 . PubMed PMID: 37007166 ; PubMed Central PMCID: PMCPMC10061940 . OpenUrl CrossRef PubMed 52. ↵ Crona K. Recombination and peak jumping . PLoS One . 2018 ; 13 ( 3 ): e0193123 . Epub 20180301. doi: 10.1371/journal.pone.0193123 . PubMed PMID: 29494618 ; PubMed Central PMCID: PMCPMC5832230 . OpenUrl CrossRef PubMed 53. ↵ Rouzine IM . Mathematical Modeling of Evolution. Volume 1. One-Locus and Multi-Locus Theory and Recombination . Berlin : De Gruyter ; 2020 . 54. ↵ Rouzine IM . Mathematical Models of Evolution. Volume 2: Fitness Landscape, Red Queen, Evolutionary Enigmas, and Applications to Virology . Berlin/Boston : De Gruyter ; 2023 . 305 p. 55. ↵ Rouzine IM , Weinberger LS . The quantitative theory of within-host viral evolution . J Stat Mech: Theory and Experiment . 2013 ; 2013 . doi: 10.1088/1742-5468/2013/01/P01009 . OpenUrl CrossRef 56. ↵ Rouzine IM , Rozhnova G. Antigenic evolution of viruses in host populations . PLoS Pathog . 2018 ; 14 ( 9 ): e1007291 . Epub 2018/09/13. doi: 10.1371/journal.ppat.1007291 . PubMed PMID: 30208108 ; PubMed Central PMCID: PMCPMC6173453 . OpenUrl CrossRef PubMed 57. Yan L , Neher RA , Shraiman BI . Phylodynamic theory of persistence, extinction and speciation of rapidly adapting pathogens . Elife . 2019 ; 8 . Epub 2019/09/19. doi: 10.7554/eLife.44205 . PubMed PMID: 31532393 ; PubMed Central PMCID: PMCPMC6809594 . OpenUrl CrossRef PubMed 58. ↵ Marchi J , Lassig M , Walczak AM , Mora T. Antigenic waves of virus-immune coevolution . Proc Natl Acad Sci U S A . 2021 ; 118 ( 27 ). Epub 2021/06/30. doi: 10.1073/pnas.2103398118 . PubMed PMID: 34183397 ; PubMed Central PMCID: PMCPMC8271616 . OpenUrl Abstract / FREE Full Text 59. ↵ Muller HG . Some genetic aspects of sex . American Naturalist . 1932 ; 66 : 118 – 28 . OpenUrl CrossRef Web of Science 60. ↵ Rice S. Evolutionary Theory: Mathematical and Conceptual Foundations . : Sinauer Associated ; 2004 . 61. ↵ Felsenstein J. The evolutionary advantage of recombination . Genetics . 1974 ; 78 ( 2 ): 737 – 56 . doi: 10.1093/genetics/78.2.737 . PubMed PMID: 4448362 ; PubMed Central PMCID: PMCPMC1213231 . OpenUrl Abstract / FREE Full Text 62. ↵ Rouzine IM , Wakeley J , Coffin JM . The solitary wave of asexual evolution . Proc Natl Acad Sci U S A . 2003 ; 100 ( 2 ): 587 – 92 . Epub 2003/01/15. doi: 10.1073/pnas.242719299 . PubMed PMID: 12525686 ; PubMed Central PMCID: PMC141040 . OpenUrl Abstract / FREE Full Text 63. ↵ Hill WG , Robertson A. The effect of linkage on limits to artificial selection . Genet Res . 1966 ; 8 ( 3 ): 269 – 94 . Epub 1966/12/01. PubMed PMID: 5980116 . OpenUrl CrossRef PubMed Web of Science 64. ↵ Gerrish PJ , Lenski RE . The fate of competing beneficial mutations in an asexual population . Genetica . 1998 ; 102-103 ( 1-6 ): 127 – 44 . Epub 1998/08/28. PubMed PMID: 9720276 . OpenUrl CrossRef PubMed 65. Desai MM , Fisher DS . Beneficial mutation selection balance and the effect of linkage on positive selection . Genetics . 2007 ; 176 ( 3 ): 1759 – 98 . Epub 2007/05/08. doi: 10.1534/genetics.106.067678 . PubMed PMID: 17483432 ; PubMed Central PMCID: PMCPMC1931526 . OpenUrl Abstract / FREE Full Text 66. ↵ Hallatschek O. The noisy edge of traveling waves . Proc Natl Acad Sci U S A . 2010 ; 108 ( 5 ): 1783 – 7 . PubMed PMID: 21187435 . OpenUrl Abstract / FREE Full Text 67. ↵ Rodriguez-Rivas J , Croce G , Muscat M , Weigt M. Epistatic models predict mutable sites in SARS-CoV-2 proteins and epitopes . Proc Natl Acad Sci U S A . 2022 ; 119 ( 4 ). Epub 2022/01/14. doi: 10.1073/pnas.2113118119 . PubMed PMID: 35022216 ; PubMed Central PMCID: PMCPMC8795541 . OpenUrl Abstract / FREE Full Text 68. ↵ Zeng H-L , Dichio V , Rodríguez Horta E , Thorell K , Aurell E. Global analysis of more than 50,000 SARS-CoV-2 genomes reveals epistasis between eight viral genes . Proceedings of the National Academy of Sciences . 2020 ; 117 ( 49 ): 31519 – 26 . doi : doi: 10.1073/pnas.2012331117 . OpenUrl Abstract / FREE Full Text 69. ↵ Batorsky R , Kearney MF , Palmer SE , Maldarelli F , Rouzine IM , Coffin JM . Estimate of effective recombination rate and average selection coefficient for HIV in chronic infection . Proc Natl Acad Sci U S A . 2011 ; 108 ( 14 ): 5661 – 6 . Epub 2011/03/26. doi: 1102036108 [pii] 10.1073/pnas.1102036108 . PubMed PMID: 21436045 ; PubMed Central PMCID: PMC3078368 . OpenUrl Abstract / FREE Full Text 70. ↵ Estimating fitness effects of mutations in the presence of genetic linkage [Internet] . 2025 . Available from: https://www.biorxiv.org/content/10.1101/2025.11.26.690712v1 . 71. ↵ Katoh K , Rozewicki J , Yamada KD . MAFFT online service: multiple sequence alignment, interactive sequence choice and visualization . Brief Bioinform . 2019 ; 20 ( 4 ): 1160 – 6 . doi: 10.1093/bib/bbx108 . PubMed PMID: 28968734 ; PubMed Central PMCID: PMCPMC6781576 . OpenUrl CrossRef PubMed 72. ↵ Rouzine IM , Coffin JM . Evolution of human immunodeficiency virus under selection and weak recombination . Genetics . 2005 ; 170 ( 1 ): 7 – 18 . Epub 2005/03/04. doi: 10.1534/genetics.104.029926 . PubMed PMID: 15744057 ; PubMed Central PMCID: PMCPMC1449738 . OpenUrl Abstract / FREE Full Text 73. ↵ Neher RA , Leitner T. Recombination rate and selection strength in HIV intra-patient evolution . PLoS Comput Biol . 2010 ; 6 ( 1 ): e1000660 . Epub 2010/02/04. doi: 10.1371/journal.pcbi.1000660 . PubMed PMID: 20126527 ; PubMed Central PMCID: PMC2813257 . OpenUrl CrossRef PubMed 74. ↵ Rast LI , Rouzine IM , Rozhnova G , Bishop L , Weinberger AD , Weinberger LS . Conflicting Selection Pressures Will Constrain Viral Escape from Interfering Particles: Principles for Designing Resistance-Proof Antivirals . PLoS Comput Biol . 2016 ; 12 ( 5 ): e1004799 . Epub 20160506. doi: 10.1371/journal.pcbi.1004799 . PubMed PMID: 27152856 ; PubMed Central PMCID: PMCPMC4859541 . OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted May 16, 2026. Download PDF Supplementary Material Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Epistasis detection and fitness valleys in viral evolution Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Epistasis detection and fitness valleys in viral evolution I.V. Likhachev , I.M. Rouzine bioRxiv 2025.09.05.674565; doi: https://doi.org/10.1101/2025.09.05.674565 Share This Article: Copy Citation Tools Epistasis detection and fitness valleys in viral evolution I.V. Likhachev , I.M. Rouzine bioRxiv 2025.09.05.674565; doi: https://doi.org/10.1101/2025.09.05.674565 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Evolutionary Biology Subject Areas All Articles Animal Behavior and Cognition (7642) Biochemistry (17715) Bioengineering (13907) Bioinformatics (42005) Biophysics (21472) Cancer Biology (18624) Cell Biology (25534) Clinical Trials (138) Developmental Biology (13391) Ecology (19935) Epidemiology (2067) Evolutionary Biology (24356) Genetics (15617) Genomics (22529) Immunology (17753) Microbiology (40437) Molecular Biology (17200) Neuroscience (88697) Paleontology (667) Pathology (2840) Pharmacology and Toxicology (4829) Physiology (7653) Plant Biology (15171) Scientific Communication and Education (2046) Synthetic Biology (4304) Systems Biology (9827) Zoology (2272)