Two fitness inference schemes compared using allele frequencies from 1068 391 sequences sampled in the UK during the COVID-19 pandemic

The data utilized here was sourced from the GISAID repository, spanning from the beginning of the COVID-19 pandemic until 12 August 2023. Subsequently, a remarkable decline in the number of genomes uploaded to the database was observed. The inclusion criteria for our analysis involved selecting only high-quality and full-length genomes, as per the defined standards outlined on the GISAID website. All retained genomes possess a length not exceeding 29 903 base pairs. Each genome is labeled based on its sample collection date, a metadata parameter provided by GISAID. Given the observable bias in alternative submission times to GISAID [11], sequences are systematically stratified weekly, resulting in a total of 179 datasets encompassing 5644 661 sequences. Due to a significant geographical imbalance in the collected samples of SARS-CoV-2, the main analysis reported here is focused on the UK region. For robustness, we show in the appendix also corresponding results for three geographic regions (Colorado, Florida, and Japan) of similar size as the UK, and for most of the time also having sufficiently many sequences per week [12].

The data for the regions in figure 1 satisfy the following minimal criteria:

• In any period of 5 days within the time series, there are at least 20 total samples.
• The number of days in the time series is greater than 20.

Zoom In Zoom Out Reset image size

Applying these criteria, our dataset for the UK region spans 170 weeks, ranging from 2020-03-14 to 2023-06-04. The dataset for the UK comprises 1068 391 sequences in total.

For the sequences in each week, a Multiple Sequence Alignment was constructed through the MAFFT software [13, 14]. Sequences from each week are aligned separately to the reference sequence 'Wuhan-Hu-1', with GISAID accession number EPI_ISL_40 2125 [15]. Note that this is different from the procedure in [16], where pre-aligned MSAs were used. The total number of sequences in that study was much less than that used here.

Each MSA is a matrix $\boldsymbol{\sigma} = \{\sigma_i^n|i = 1,\ldots,L, n = 1,\ldots,N\}$, where N represents the number of genomic sequences in a week while L represents the number of loci in an aligned sequence [17, 18]. Thus, all aligned sequences have a length of $L =$ 29 903, the same as that of the reference sequence, while N by construction varies from week to week, see figure 1. The loci between 256 and 29 674 are referred to as coding region, since they code for the protein-coding genes in the SARS-CoV-2 genome. Each entry $\sigma_i^n$ of the MSA σ is either one of the 4 nucleotides (A,C,G,T), or the alignment gap '-', the minorities like 'KYF...' are changed to the sign of '-' for the sake of simplicity of the following allele frequency analysis.

In figure 2 we show the allele frequency time series for all loci in the UK data, and in figure 3 only for loci in the coding region. These two figures show that at a majority of loci all sequences contain the same symbol, and most of the remaining variation is in the non-coding regions. In a first filtering step we retain in the analysis loci where the most frequent mutation away from wild-type is classified as 'Non-synonymous mutation' as defined in [6] and where the largest mutant frequency is at least 1%. In a second step we retain only those loci which meet the criteria for all weeks. For the UK data used in this study there remains 209 loci.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the other data sets shown in figure 1(lower panel) there would remain respectively 1063 loci in 'Global', 173 loci in 'EU', 328 loci in 'NA' (North America) and 225 loci in 'Asia'. For consistency, the average Hamming distances are however for all computed from the variability at the same 209 loci as in the UK data set.

In both approaches to be considered, the driving forces of evolution are assumed (Darwinian) selection, mutation, recombination, and genetic drift (finite population effects). Effects excluded from consideration are hence e.g. spatial barriers (island models). The genome $\mathbf{x} = \left(x_1,x_2,\ldots,x_L\right)$ where each x_i is an indicator variable of the allele (nucleotide in the set $\{-,N,A,C,G,T\}$ at locus i (position i in the MSA).

• (a)  
  Fitness is assumed to be a function

  $$\begin{align} F\left(\mathbf{x}\right) = \sum_{i = 1}^L\sum_a f^{\,\left(1\right)}_{i,a} \mathbf{1}_{x_i,a} + \sum_{i = 1\lt j}^L\sum_{a,b} f^{\,\left(2\right)}_{ij,ab} \mathbf{1}_{x_i,a} \mathbf{1}_{x_j,b}.\end{align} \tag{ 1 }$$

  The coefficients $f\,^{(1)}_{i,a}$ are called additive fitness and parameterize the selective advantage of allele a at locus i with respect to wild-type. The coefficients $f\,^{(2)}_{ij,ab}$ are called (pair-wise) epistatic fitness and parameterize the selective advantage of alleles a and b at loci i and j beyond what they contribute separately. In tQLE $f\,^{(2)}_{ij,ab}$ are adjustable parameters inferred from the data, which for self-consistency however cannot be too large. In MPL $f\,^{(2)}_{ij,ab}$ are absent.The evolution of genotypes in the population over a short time Δt due to fitness is on the level of a normalized distribution given by

  $$\begin{align} P\left(\mathbf{x},t+\Delta t\right)|_{\mathrm{fit.}} = \frac{e^{\Delta t F\left(\mathbf{x}\right)}} {\sum_{\mathbf{x}} e^{\Delta t F\left(\mathbf{x}\right)} P\left(\mathbf{x},t\right)} P\left(\mathbf{x},t\right).\end{align} \tag{ 2 }$$

  In both approaches our goal is to estimate the coefficients $f^{(1)}_{i,a}$.
• (b)  
  Mutations are random changes of single alleles. In general, they could be parametrized acting as

  $$\begin{align} P\left(\mathbf{x},t+\Delta t\right)|_{\mathrm{mut.}} & = P\left(\mathbf{x},t\right) + \Delta t \sum_{i}\sum_{ab} \mathbf{1}_{x_i,a} \nonumber\\ & \quad \times \left( \mu_i^{ba} P\left(M_i^{ab}\mathbf{x},t\right) - \mu_i^{ab} P\left(\mathbf{x},t\right) \right)\end{align} \tag{ 3 }$$

  where $M_i^{ab}$ is the flip operator which changes allele a at locus i to b, and $\mu_i^{ab}$ rate of this process. These rates are only partially known, and only parametrise a fraction of naturally occurring mutations. As our focus is here on fitness we will follow the original theoretical literature [5, 6] and take them all equal to one overall mutation rate µ.
• (c)  
  Recombination is modelled as a process whereby two genomes combine and give rise to a third. On the level of distributions that is given by

  $$\begin{align} P\left(\mathbf{x},t+\Delta t\right)|_{\mathrm{recomb}} & = P\left(\mathbf{x},t\right)\left(1-r \Delta t\right) \nonumber \\ & \quad + r\Delta t \sum_{\mathbf{x}_m,\mathbf{x}_f} C\left(\mathbf{x};\mathbf{x}_m,\mathbf{x}_f\right) \nonumber \\ & \quad\times P\left(\mathbf{x}_m,t\right) P\left(\mathbf{x}_f,t\right).\end{align} \tag{ 4 }$$

  In above r is an overall recombination rate with dimension inverse time. The function $C(\mathbf{x};\mathbf{x}_m,\mathbf{x}_f)$ is the specific rate at which genomes x_m and x_f combine to yield x. In tQLE it only enters through the derived quantity c_ij which is the probability that alleles at loci i and j are inherited from different parents [5, 10, 19]. In MPL recombination drives the evolution and influences the distribution of evolutionary trajectories, but does not directly enter in the inference formulae. The factor $\left(1-r \Delta t\right)$ serves to normalize the distribution. More general models of recombination in the same context are discussed in [19].
• (d)  
  Genetic drift is the term of stochastic effects due to a finite population. All three evolution equations (2)–(4) are valid on the ensemble level, and can be simulated by evolving several populations in parallel, and then averaging. Single-locus frequencies and two-loci pair frequencies (and other characteristics) will evolve due to both deterministic drift and random noise. In QLE the corresponding stochastic differential equations for single- and two-loci frequencies are derived and discussed [5]. In MPL the stochastic differential equations for single-locus frequencies are central in deriving the path probabilities as discussed below.

The phase of QLE was discovered by Kimura in the study of a two-locus biallelic model [3]. The extension to many loci was investigated by Neher and Shraiman [4, 5]. The generalization to more than two alleles per locus was given in [19]. The two defining properties of QLE (formalized in [10]) are

• 1.  
  Multi-genome probability distributions factorize such that $P_n(\mathbf{x}^{(1)},\mathbf{x}^{(2)},\ldots,\mathbf{x}^{(n)}) = P(\mathbf{x}^{(1)})P(\mathbf{x}^{(2)})\cdots P(\mathbf{x}^{(n)})$. This property is especially important for n = 2 as it allows to model the effects of recombination as a molecular collision in kinetic gas theory (Boltzmann's Stosszahlansatz).
• 2.  
  The single-genome probability distributions are Gibbs distribution with terms no higher than in fitness. For (1) this means the Ising-Potts distributions of equilibrium statistical mechanics

  $$\begin{align} P\left(\mathbf{x}\right)& = \frac{1}{Z\left(\left\{h_t\right\},\left\{J_t\right\}\right)} \exp\left(\sum_{i,a} h_{i,t}\left(a\right)\mathbf{1}_{x_i,a} \right. \nonumber \\ & \quad \left. +\sum_{i\lt j,ab} J_{ij,t}\left(a,b\right) \mathbf{1}_{x_i,a}\mathbf{1}_{x_j,b}\right).\end{align} \tag{ 5 }$$

In (5) we have included a time argument t of the Ising-Potts parameters $h_{i,t}$ and $J_{ij,t}$. This is first to emphasize that (5) is not based on assumptions that lead to thermal equilibrium where h_i and J_ij are constant. Second, and more importantly, when (5) is fit to time-ordered data, as will be described, the inferred parameters $h_{i,t}, J_{ij,t}$ do depend on time t. For simplicity of notation we will from now on however not explicitly write this argument t. Relations between fitness parameters $\{f\,^{(1)}\}$ and $\{f\,^{(2)}\}$ in (1) and stationary Ising/Potts model parameters $(\{h\},\{J\})$ (5) are key quantitative results of QLE theory. In [3] and [5] they were derived in the limit where overall recombination rate r is larger than both overall mutation rate µ and variations in fitness. Direct tests using scatter plots were first given in [9]. Several alternative relations were derived for larger mutation rates in [20], of which one was tested in [20], see also [10].

As already stated above, QLE is a dynamic theory where parameters $(\{h\},\{J\})$ in general change in time. We here introduce the derived abbreviation tQLE to emphasize that we use the formulas for inference on such in principle (and generally in practice) time-changing data.

The equation for $(\{J\})$ is of the relaxation type, and in the theory of [5]

$$\begin{equation} \dot{J}_{ij}\left(a,b\right) = f\,^{\,\left(2\right)}_{ij}\left(a,b\right) - r c_{ij} J_{ij}\left(a,b\right).\end{equation} \tag{ 6 }$$

For large enough r the Potts parameters will hence relax to a stable fixed point, i.e. to $J^*_{ij}(a,b) = \frac{f^{(2)}_{ij}(a,b)}{r c_{ij}}$, which allows to infer epistatic fitness parameters from Potts parameters computed from the data through the formula

$$\begin{equation} f^{\left(2\right),*}_{ij}\left(a,b\right) = r c_{ij} J^*_{ij}\left(a,b\right).\end{equation} \tag{ 7 }$$

This relation was derived in [5], and tested (in stationary state) in [9]. As discussed in [21] since (7) only relates pair-wise quantities, it can also work when the single-nucleotide frequencies change. This could for instance be the case of additive fitness changes in time, say by a change of the fitness landscape of which one example could be the introduction of widespread vaccination against SARS-CoV-2 in the COVID19 pandemic.

The equation for $(\{h\})$ is on the other hand not of the relaxation type [5], equation (24)

$$\begin{equation} \dot{h}_{i}\left(a\right) = f^{\,\left(1\right)}_{i}\left(a\right) + r \sum_{j,b} c_{ij} J_{ij}\left(a,b\right) m_j\left(b\right)\end{equation} \tag{ 8 }$$

where $m_j(b) = \sum_{\mathbf{x}} P(\mathbf{x}) \mathbf{1}_{x_j,b}$ is the frequency of allele b at locus j. Combining (7) and inferred values of $\{h\}$ at two consecutive time intervals lead to the inference formula

$$\begin{align} f^{\left(1\right),*}_{i}\left(a, \Gamma\right) & = \frac{1}{\Delta t}\left[h_{i}\left(a, \Gamma + \Delta t\right) - h_{i}\left(a,\Gamma\right)\right] \nonumber \\ &\quad - \sum_{j,b} f^{\left(2\right),*}_{ij}\left(a,b, \Gamma\right)\, m_j\left(b,\Gamma\right)\,,\end{align} \tag{ 9 }$$

with Γ indicating the discrete time.

The QLE state is lost when the distribution no longer fulfills the two listed criteria. A well-studied loss channel at very low mutation rate and sufficiently low recombination rate is through the emergence of clones [5, 22]. These are groups of identical, high-fitness genomes related by common descent. Instead of an exponential model (as in QLE), the distribution of genotypes is instead a mixture of clones, with one separate distribution for each clone. There may also be only a single clone, in which case all (or most) of the genotypes are the same. Quantitative predictions on the threshold between a QLE phase and a clone-dominated phase were derived in [22]. One aspect of the transition between QLE and a clone-dominated phase is that QLE can only exist as a transient state in a finite population with strictly no mutations. The reason is that in this setting sooner or later the most fit genome takes over as a single dominating clone, see [19] and [9] for a discussion. Hence QLE loss at non-zero mutation rate is relevant. A second loss channel observed at a higher mutation rate leads to a phase of 'noisy clones' coexisting with a QLE-like state. In [10] this new phase was named Non-Random Coexistence (NRC). The transition from QLE to NRC goes through an intermittent phase where the state of the population jumps between QLE and NRC, illustrating the complexity of phenomena not yet completely mapped out even in theoretical models. The dependence of the jump rates on population size N was investigated in [10], and leads to a qualitatively similar behaviour as [22], i.e. for a sufficiently large population only the NRC phase is stable.

The marginal path likelihood (MPL) method [6] is based on the evolution of nucleotide frequencies in Kimura's diffusion approximation [23, 24]. The starting point is thus the joint probability $P(\{m\}^{(1)},\{m\}^{(2)},\ldots,\{m\}^{(L)})$ where $m^{(i)}_a$ is the frequency of allele a on locus i, normalized as $\sum_a m^{(i)}_a = 1$. In the diffusion approximation, this probability satisfies a Fokker-Planck equation

$$\begin{align} \partial_t P & = - \sum_{i,a} \frac{\partial}{\partial m^{\left(i\right)}_a} \left(u^{\left(i\right)}_a P\right)\nonumber\\ & \quad + \sum_{ij,ab} \frac{\partial^2}{\partial m^{\left(i\right)}_a\partial m^{\left(j\right)}_b} \left(D^{\left(ij\right)}_{ab} P\right)\end{align} \tag{ 10 }$$

where the drift vector and diffusion matrix are given by ([6], equations (6) and (S9) and following, notation aligned with the present presentation)

$$\begin{align} u^{\left(i\right)}_a & = m^{\left(i\right)}_a \left(1-m^{\left(i\right)}_a\right)f^{\,\left(1\right)}_i\left(a\right) + \mu\left(1- 2 m^{\left(i\right)}_a\right) \nonumber \\ &\quad + \sum_{j,b} \left( m^{\left(ij\right)}_{ab} - m^{\left(i\right)}_a m^{\left(j\right)}_b\right) f^{\,\left(1\right)}_j\left(b\right)\end{align} \tag{ 11 }$$

$$\begin{eqnarray} D^{\left(ij\right)}_{ab} = \left\{\begin{array}{lr} m^{\left(i\right)}_a m^{\left(i\right)}_b & i = j \\ m^{\left(ij\right)}_{ab} - m^{\left(i\right)}_b m^{\left(j\right)}_b & i\neq j \end{array}\right..\end{eqnarray} \tag{ 12 }$$

The Fokker–Planck equation (10) corresponds to a multidimensional Langevin equation for which the probability of a path sampled at discrete times can be estimated by standard arguments. Maximizing this path probability with a Gaussian prior leads to the central inference formula in MPL

$$\begin{align} f^{\left(1\right),*}_{i}\left(a\right) & = \sum_{j,b} \left[\sum_{k = 1}^K \Delta t_k D^{\left(ij\right)}_{ab}\left(t_k\right)+\gamma\mathbf{1}_{ia,jb}\right]^{-1}_{ia,jb} \nonumber \\ & \quad \times \left[m^{\left(j\right)}_b \left(t_K\right)-m^{\left(j\right)}_b \left(t_0\right) \right. \nonumber \\ & \quad \left. -\mu \sum_{k = 1}^{K-1} \Delta t_k \left(1-2 m^{\left(j\right)}_b\left(t_k\right)\right) \vphantom{m^{\left(j\right)}_b \left(t_K\right)-m^{\left(j\right)}_b \left(t_0\right) -\mu \sum_{k = 1}^{K-1} \Delta t_k}\right] .\end{align} \tag{ 13 }$$

In above a time interval $[t_0,t_K]$ has been divided up in K sampling intervals and the allele frequencies ($m^{(i)}_a$) and drift and diffusion terms (from (11) and (12)) estimated for each. The sampling interval times are defined as $\Delta t_k = t_{k+1}-t_k$. γ is the width of the Gaussian prior, and acts as a regularizer.

In (9) the inferred additive fitness depends on time and is linear in the time derivative of one inferred Potts parameters $h_{i}(a)$. This parameter is in itself a (complicated) function of the single-nucleotide and pair-wise frequencies at that time. In (13) the inferred additive fitness also depends on time, more specifically on a time interval, and is linear in $m^{(i)}_a (t_K)-m^{(i)}_a (t_0)$, the change in all the single-nucleotide frequencies over that interval. Pair-wise frequencies also enter in (13), through the dependence of $D^{(ij)}_{ab}$ as in (12).

SARS-CoV-2 is in classifications such as Pango [25, 26], assumed to evolve by descent, and the growth and subsequent decay of SARS-CoV-2 variants [27–29] is well-known. This can be taken to be the standard view of SARS-CoV-2 evolution, analogous to the evolution of other viruses such as influenza. A complicating factor is that many viral variants seem themselves to evolve, to split into sub-variants and perhaps to recombine [30, 31]. On the other hand, coronaviruses in general exhibit recombination [32–35], a process which has also been observed to occur in SARS-CoV-2 [36–40].

Whether a propensity for recombination between pairs of viral genomes results in observed recombination in a viral population is not a simple question; two viruses must meet, typically in the same host, to recombine. The effectiveness and importance of recombination in the general SARS-CoV-2 viral population has accordingly been questioned [41]. The classic QLE phase (Kimura and Neher & Shraiman) assumes a fully mixing population, which the global SARS-CoV-2 viral population during the pandemic clearly was not. If a QLE-like phase exists in a not fully mixing population, and at which parameter values, has not, as far as we are aware of, been systematically studied in the literature. On the other hand, as found by two of us in previous theoretical work [20], large recombination rate is not the only parameter range that leads to a QLE state in a fully mixing population; a large mutation rate with some recombination rate is also sufficient. An interesting avenue for further theoretical studies would be to try to find out if this mechanism also operates in not fully mixing populations.

The recombination rate r sets a scale for the epistatic contributions to fitness which does not change the order of fitness of sequences as long as inferred epistatic fitness dominates over inferred additive fitness. The MPL inference scheme includes a regularization parameter which plays a similar role; this seems at present to be unavoidable in these studies. Here the value of r is set to 0.027 as indicated in [42].

Given two lists of sequences ordered as to fitness, we compare the rankings by the Spearman correlation coefficient. The Spearman correlation coefficient is a measure of rank correlation obtained by computing the Pearson correlation between the rankings of the values in each vector. Thus, it can robustly measure nonlinear associations between fitness values as well as linear ones. This is computed with the Matlab function 'corr()' with type argument 'Spearman'.

The MSAs over weeks are prepared and filtered according to the processes described in section 2. For tQLE, the maximization of the pseudo-likelihood (PLM) [43, 44] is used on each MSAs to learn the Ising-Potts parameters $h_{i,t}(a)$ and $J_{ij,t}(a,b)$ in equation (5). Then the epistatic fitness parameters $f_{ij}^2(a,b)$ are computed according to equation (7). In which the recombination rate r = 0.027 [42] and the probability that alleles at loci i and j are inherited from different parents $c_{ij} = 0.5$ [9] respectively. With the inferred $h_{i,t}(a)$s and $J_{ij,t}(a,b)$s by PLM, the additive fitness parameters for each week are obtained according to equation (9).

For MPL, the set of all weekly data defines the path. Thus we use the entire path (data set) to compute the additive fitness parameters [7]. The MPL library is available at [45].

The fitness values for each sequence in every week are computed using equation (1), using the week additive fitness parameters from tQLE and the whole set from MPL respectively. Then the rank orders for the fitness of sequences in each week are compared through their Spearman correlation to check the inference similarities of these two schemes.

The weekly stratified datasets are geographically segmented for a detailed analysis. As illustrated in figure 1(upper panel), the number of sequences undergoes dynamic changes across weeks for different regions, including Europe (dashed orange lines), North America (dashed yellow lines), Asia (dashed purple line), three distinct regions of the UK (green line), and the global dataset (blue line). To account for the impact of geographic separation, data from North Ireland is intentionally excluded. Moreover, figure 1(upper panel) reveals a significant downturn in the number of collected samples in 2022, with Europe emerging as the primary source of sequences. The diversity of collected sequences is in figure 1(lower panel) quantified by average Hamming distance. One observes increased diversity after the emergence of Alpha, Delta and Omicron (marked in figure).

We computed allele frequencies over time from SARS-CoV-2 multiple sequence alignments from the UK. Subsequently, the allele frequencies for the nucleotides in the 'Wuhan-Hu-1' reference sequence are selected. Figure 2 illustrates the allele frequencies over all loci ($L = 29,903$ base pairs in total), while figure 3 specifically shows the allele frequencies within the coding region, spanning from the 256th to the 29 674th sites in the sequence. There are sustained fluctuations at the bottom of figure 2, which hence mainly originate from the non-coding region (3'-UTR and 5'-UTR parts) of SARS-CoV-2. Both plots display variations at specific loci within the coding region, of which many of them can be related to listed mutations in known Variants of Concern, compare monthly data and a relation to Omicron reported in [46] and [21] (figure 5). Furthermore, the increased variability in allele frequencies after Omicron took over (after 1Q22) matches the broad peak in sequence diversity shown in figure 1(lower panel).

Epistatic fitness or covariation selection coefficients f_ijs in QLE follow from the theory developed in [3, 5] (in the version for bi-allelic loci). This theory generalizes directly to multi-allelic loci [9, 19, 21], and leads to the inference formulae equation (7). The additive fitness or selection coefficients f_is are in tQLE analogously obtained as the differences between the time derivative of an additive term ($\dot{h}_i$) and a combination of the epistatic terms and allele frequencies ($\sum_j f_{ij} m_j$). In the generalization to multi-allelic loci, and when time derivatives are approximated as discrete time differences, this leads to inference formula equation (9) where t and $t+\Delta$ stand for two different weeks.

As shown in figure 4, the additive fitness by tQLE with the time interval of $\Delta t = 4$ (green) and 8 (red) weeks are consistent with each other. This indicates that equation (9) tQLE fitness values do not strongly depend on the choice of Δt values, which is further shown in figure 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 5 shows the epistatic terms versus f_is inferred by the tQLE for $\Delta t = 4$ weeks. The red squares lay closely around the diagonal. Such high consistency demonstrates that the epistatic term dominates over the time derivative $\dot{h}_i$ in the total inference formula for f_is. This pattern is consistent across all times in the present data set.

The additive fitness inferred by MPL validates its stability and reliability for the SARS-CoV-2 datasets. Here, even with a different stratification strategy for the UK datasets, similar results are obtained with those in [7]. While the great majority of fitness effects of mutations are inferred to be nearly neutral, MPL infers both strongly beneficial and deleterious mutations, as illustrated by the conspicuous deviation of the blue bars in figure 4 from the neutral zero point.

Figure 4 shows the histograms of the additive fitness f_is by the tQLE (green and red bars) and MPL (blue bars) models, respectively. The f_is anticipated by the tQLE model exhibit a distribution proximate to the zero point, constraining in a relatively narrow range. In contrast, fitness effects inferred by MPL are mostly concentrated around zero, but with large deviations for a small number of mutations.

To assess the fitness estimates derived from the tQLE and MPL methods, we used both methods to rank sequences within a time window according to their fitness, and then compared these rankings. The fitness score of a sequence is defined in equation (1), in which the f_is from the first term are the fitness effect while the f_ijs from the second term are the epistatic term. The total fitness of a sequence from the MPL corresponds to the first term of equation (1) while that for the tQLE method is given by both terms of equation (1). For each time window (one week) the tQLE and MPL fitness scores for each sequence are computed and subsequently arranged in descending order.

In figure 6 we show the Spearman rank correlation coefficient c between the two rankings, stratified as to time (mean value per week, upper panel) and as to overall distribution (histogram, lower panel). The agreement is generally good (c > 0.8 for 105 out of 166 weeks). MPL and tQLE hence order sequences as to fitness in closely similar ways, for these UK-sampled SARS-CoV-2 sequences obtained from GISAID.

Zoom In Zoom Out Reset image size

Figure 6 also relates agreement/disagreement in rankings to diversity (or lack thereof) in the sequences obtained in one week. For easier visualization we have plotted Spearman correlation c(t) together with the relative Hamming distance deficit $rH^{(-)}$ (defined in caption to figure 6). The two curves move in concert. This indicates a contravariation of c(t) with average Hamming distance H(t), i.e. that periods of high (low) Spearman correlation are related to periods of low (high) sequence diversity. A possible explanation is that in this data set periods of high sequence variability (large average Hamming distance) appeared when one VoC was in the process of taking over the population, see figure 1(lower panel). At these times the UK viral population resembled a mixture of two clones, different from the one Gibbs-Boltzmann distribution posited in QLE theory from which the tQLE inference method has been derived, see section 3.