Studying viral populations with tools from quantum spin chains
2National Institute of Standards and Technology, Gaithersburg, MD 20899, USA
3Joint Quantum Institute, University of Maryland, College Park, MD 20742, USA
4Department of Physics and Astronomy, University of California, Riverside, CA 92521, USA
5Department of Physics, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA
)
Abstract
We study Eigen’s model of quasi-species [1], characterized by sequences that replicate with a specified fitness and mutate independently at single sites. The evolution of the population vector in time is then closely related to that of quantum spins in imaginary time. We employ multiple perspectives and tools from interacting quantum systems to examine growth and collapse of realistic viral populations, specifically certain HIV proteins. All approaches used, including the simplest perturbation theory, give consistent results.
1 Introduction
A central concept in population genetics is an idealized fitness , which (in the absence of competition, resource limitation, etc.) governs the exponential growth in the number of individuals according to [2]. However, mutations diversify the genetic make-up (genotype) of the population, and modify its overall fitness. Eigen introduced the concept of quasi-species to describe the ‘cloud’ of closely related genotypes. Eigen’s model considers a population composed of a set of sub-populations (labelled ), each contributing individuals that coexist in a rapidly evolving larger community of (varying) total size [1]. Mutations between quasi-species further diversify the composition of the population. If the sub-populations are assigned fitness values , and mutations from sub-population to occur at rates , the makeup of the population changes over time according to
(1) |
In terms of the fractions , the total population size grows as . The mean fitness , must be non-negative for the quasi-population to survive at long times.
At the molecular level, genetic information is maintained in the sequence of bases of nucleotides (e.g. in the sequence of the DNA or RNA of a virus). For modeling purposes, we assume that this information appears as a sequence of characters (e.g. one of 4 nucleotides of DNA, or 20 amino acids of a protein). For simplicity of presentation (as well as practicality, as in the case of HIV proteins described below) we adapt a binary representation, or , for , to indicate whether the site is in its wild-type (consensus) or a mutated state. In this representation of a genotype, its fitness is some function of its sequence . Note that by our definition of a consensus sequence, it is possible for a sequence to have a higher fitness than the consensus sequence, since the consensus is based only on single-site frequencies.
In a simple model of mutations, the state of each site (independently of other sites) changes away from consensus at the rate , and reverts back to wild-type at the rate . Equation (1) is linear, and its right hand side can be regarded as representing the action of a matrix on a (column) vector containing the sub-population sizes . For , the action of mutations on site can be represented by the Pauli matrix , such that
(2) |
The above equation provides an analogy to the imaginary-time evolution of spins in a quantum chain, governed by the Hamiltonian . The fitness function corresponds to interactions between spins, while mutations are implemented through the transverse magnetic field .
The analogy between Eigen’s quasi-species model and quantum chains has been noted and explored in a number of references [3, 4, 5, 6, 7, 8]. The behaviour of genetic structure of a quasi-species model and magnetization of the corresponding quantum spin model for some idealized fitness landscapes (like an Ising chain or a mean field Hamiltonian) have been studied in Refs. [3, 4, 5]. Hermisson et al. [7] generalize the binary representation of DNA sequences to the four nucleotides, leading to a corresponding four state quantum chain. Saakian & Hu [8] consider mutations between sequences separated by more than one Hamming distance (see also Ref. [9]). As is well known in physics literature, a quantum system in dimensions is equivalent to a classical system in dimensions with discretized time. Leuthäusser [10] explored a similar analogy between the Eigen model and a two dimensional Ising system.
There are, however, important differences between quasi-species evolution, and the dynamics of a quantum chain (in real time):
While the Hamiltonian is real and symmetric, the evolution of is not unitary, in that the overall population size changes as a function of time. By contrast, time variation according to preserves the norm of the vector , such that can be regarded as probabilities. The natural set of probabilities for the evolving population are the proportions .
Even the symmetry of is an artifact of the simplification . In the biological context, it is more likely that forward and mutation rates are not the same, leading to the replacement of for with the asymmetric matrix [8]. Despite these differences, there are aspects of the dynamics that are common to the two systems, in certain special cases. For quantum systems at low temperature, one is often interested in the low energy properties, which are determined by the eigenvectors of corresponding to the lowest eigenvalues. Similarly, the long-time behavior in the quasi-species evolution is governed by the few largest eigenvectors of the matrix , i.e., the exact same eigenstates. In particular, the ground state energy is (minus) the mean fitness, while the ground state vector characterizes the prevalence of mutations in the quasi-species population. Together, they determine the eventual fate of the quasi-species as described below:
Error threshold: Eigen introduced the concept of an error catastrophe [1] by considering a fitness function with a single peak on a uniform background (for one or more mutations), which can be described by
(3) |
In words, the wild type has fitness while any mutation reduces the fitness to . Upon increasing the mutation rate , there is a transition when the fraction of population in the fit state decreases dramatically. This is a genuine singularity (phase transition) in the limit , although the threshold scales as . More precisely, one can show that for , the fraction of the population having mutations falls off exponentially as
(4) |
whereas for , nearly the entire population has roughly mutations.
The fitness function in Eigen’s case is of course highly artificial, and questions remain as to whether the concept of error threshold is applicable to more realistic landscapes (see also Ref. [11]). Nonetheless, sharp transitions like this, in which the relevant quantity is zero on one side and non-zero on the other, are in fact quite common in the context of quantum spin systems. Examples include the transition between ferromagnetic and paramagnetic equilibrium phases [12, 13], frozen versus thermalizing spin glasses [14, 15, 16], and many-body localization-delocalization transitions [17, 18, 19].
Fatal mutation load: In Eigen’s model, the quasi-species population can still grow if , as each sequence is viable. However, the population can also disappear upon increased mutation because of the reduction in overall fitness. As an example, consider the so-called Mount Fuji Landscape [5]:
(5) |
specifically with The most fit sequence has a fitness , with any mutation (independently) carrying a load of . With increasing mutation rate, the quasi-species cloud acquires more mutations, but (unlike the case of error threshold) remains anchored to the wild-type state. Since the problem is equivalent to independent quantum spins, it is easy to identify the ground state as corresponding to all spins tilted along , leading to the mean fitness of . For any mutation rate larger than , is negative and the population will die off, despite the lack of any error threshold.
2 Viral prevalence landscape
Several other forms of fitness landscape have been proposed and explored in the literature. We would like to explore the quantum analogy, identifying mechanisms that lead to population collapse (e.g. by error catastrophe or high mutation load). To this end, we will focus on a particular class of landscapes constructed from the prevalence of observed sequences for a virus.
The massive size of the space of possible protein sequences makes it impossible to estimate virus prevalence simply by counting sequences. Each site in a protein sequence can be one of 20 amino acids, and typical protein sequences have lengths on the order of hundreds of sites. However, the available data typically consists of only a few thousand sequences.
A common approach to this problem is to search for the simplest probability distribution (by which we mean the one with the largest entropy) over the space of protein sequences that is capable of reproducing the single site and pairwise frequencies of amino acids in the data [20]. To further reduce complexity, one can represent protein sequences with a reduced or even binary (zero for ‘wild-type’ and one for ‘mutant’) amino acid alphabet [21, 22, 23]. In the binary case, the maximum entropy model capable of reproducing the empirical correlations is an Ising model with local fields and pairwise couplings . Finding the appropriate fields and couplings to reproduce the correlations in the data is a challenging statistical problem [20]. For the models presented here, we applied the adaptive cluster expansion method [24, 22] to infer field and coupling parameters that reproduced mutant frequencies and correlations observed in sequence data from a variety of HIV proteins [25].
Following the assumption that the most prevalent viral sequences are likely to also have the highest fitnesses, the prevalence landscape can be used as a proxy for fitness. The simple connection between prevalence and fitness could be obscured due to strong interactions between HIV and the immune system, which drives the virus to accumulate mutations rapidly. Careful modeling suggests, however, that prevalence is a good proxy for fitness for viruses such as HIV, which produces chronic infections and stimulates a diverse range of immune responses [26]. This assumption is also well-supported by experiments that tested the effects of different mutations on HIV replication [21, 27, 28].
3 Estimating mean fitness
We study population collapse for a few HIV proteins for which the fitness has been estimated to have the form [25, 29]
(6) |
where are the Pauli matrices and is a constant. Including the effect of mutations as described in Eq. 2 leads to a dynamics governed by the Hamiltonian
(7) | ||||
We apply to this Hamiltonian three simple techniques, all common for studying quantum spin systems: first-order perturbation theory in , second-order perturbation theory, and a (non-perturbative) mean-field approximation.
As will be shown, all three give consistent results, and whatever discrepancies exist can be easily understood.
3.1 First-order perturbation theory
While sophisticated methods from quantum spin chains can certainly be applied to such a system, in practice we find that standard perturbation theory provides a more than adequate tool for identifying population collapse.
Let us denote the energy of the consensus state at by . Recall that the corresponding state is defined to have all . To first order in perturbation theory, the change in energy is the expectation value of in this state. The expectation value of being zero in any eigenstate, the perturbed energy is simply
(8) |
which becomes positive at the threshold value (remember that must be negative)
(9) |
3.2 Second-order perturbation theory
The second-order corrections to the energy, assuming non-degenerate levels, give an expression
(10) |
where the sum is over configurations with one spin flipped. The threshold can then be found by solving the quadratic equation
(11) |
Note from Eq. (10) that if is less than all (as is often the case, and must be if the consensus is the true ground state), then the second-order terms are necessarily negative. The energy is lower than the first-order estimate, and thus the threshold mutation rate is larger.
3.3 Mean-field approximation
As an alternative to the perturbative calculations, we also employ a mean field approach to directly find an approximation to the ground state of the quantum Hamiltonian. In this procedure the interacting Hamiltonian
(12) |
is approximated by that describing non-interacting spins, as
(13) |
The effective field at each site is expressed in terms of the average component of spins at connected sites as
(14) |
To solve the set of coupled equations for , we use an iterative procedure, starting with the , alternately computing the expectation values of spins and the corresponding effective fields until they converge. The mean-fitness is finally computed as
(15) |
from which we obtain the mutation threshold by setting .
3.4 Results on mutation threshold for HIV proteins
We considered five HIV proteins of variable lengths:
p24 (), gag (), integrase (), nef (), and protease ().
The Ising on-site and exchange fields for each protein were estimated previously [25, 29]
from prevalence of different sequences in collected data.
We then multiplied the prevalence landscape by a factor of /day, and added an overall constant of /day
to construct a fitness landscape. These choices were dictated by typical viral loads in patients, as discussed in Ref. [29].
In all cases, we both performed perturbation theory around the consensus state and applied the mean-field approximation. It will be important to note that the consensus is not the fittest state for the gag, nef, and protease proteins, whereas it is the fittest for p24 and integrase. The resulting threshold values are shown in Fig. 1.

All estimates for are quite consistent.
There are some discrepancies, but they are all as expected.
First, it is indeed the case that the first-order estimates are typically lower than the second-order values (the gag protein is the sole exception out of those considered).
For p24 (/site/day) and integrase (/site/day), the consensus sequence corresponds to the ground state and thus we find quite good agreement between the second-order and mean-field results.
Yet for gag, nef, and protease, the perturbation theory values only inform us of the mutation rate beyond which the consensus and all similar sequences die out.
Alternate quasi-species having lower energies will survive for slightly larger .
Since the mean field theory is able to (approximately) identify the true ground state, regardless of whether it agrees with the consensus or not, we expect it to predict higher values for .
This is indeed what we find for the three proteins in question.
We also note that the threshold mutation rates obtained for those are higher as compared to some obtained for other HIV proteins via alternate methods (/site/day) [30, 31].
3.5 Effect of mutations on comparison with experimental fitness
As discussed in Sec. 2, the prevalence landscape is able to provide us an approximation of the fitness landscape, and is well supported by experiments. One quantity often measured by experiments is the replication capacity (RC) of a sequence, which is an empirical estimate of its growth rate and thus strongly associated with fitness [27]. Comparing the energy and the RC of some mutatant viruses relative to the consensus sequence (for the gag protein, where multiple measurements are available) reveals a good correlation between the two. This leads us to conclude that indeed the Ising landscape derived for the proteins are a good measure of viral fitness (Figure 2).
In comparisons with experimental growth data (via RC), it may be necessary to include the effect of mutations on the Ising landscape since the observed growth includes mutations as well [27]. We use our perturbation theory approach to compare the corrected energies (Eq. 10), at a mutation rate typically expected in HIV proteins /site/day, with the RC data observed for some mutations in the gag protein in [27]. Figure 2 shows that the effect of mutations is negligible when comparing to experimental fitness. However, to see any generic effect of mutations, we considered a higher mutation rate of 0.001/site/day (which is unrealistic but still below ) and we find that the corrections can vary both in the sign and magnitude for the mutants considered here.

4 Discussion
Here we applied ideas from interacting quantum systems to study Eigen’s quasispecies model. Using past estimates for HIV fitness landscapes, which we verified to be in good agreement with experimental data, we computed the threshold mutation rate beyond which viral populations would be expected to decline due to insufficient fitness. We have found that first-order perturbation theory works remarkably well in predicting . Even though we do not have an exact result against which to compare, neither higher-order terms nor a non-perturbative mean-field approximation gives significantly different results.
Note that the perturbative expression in Eq. (9) is very similar to Eigen’s original result [1], even though the population collapse identified here is not due to an error catastrophe but rather due to negative mean fitness.
Hart and Ferguson recently investigated the possibility of an error catastrophe for HIV using estimated fitness landscapes [31]. They found evidence for a phase transition to a high energy, or low fitness, state in the HIV protein p6, which forms a part of gag. The transition occurs at a temperature larger than one, which is interpreted as a signal of a mutation rate that is higher than what is observed in nature. However, their approach does not estimate the critical mutation rate precisely in terms of a probability of mutation per replication cycle.
In fact, the immune system has special defenses that can force viruses to undergo an “error catastrophe.”
APOBEC proteins cause viral hypermutation, which nearly always results in the production of defective viruses [32].
However, APOBEC proteins can be countered by the HIV protein vif, thus allowing replication of the virus to continue unchecked.
It is clear that a detailed understanding of the “phase diagrams” for HIV proteins will be a useful tool for combating the virus.
Acknowledgements: The authors would like to thank the Galileo Galilei Institute in Florence, Italy where a part of this work was performed. This research was performed while CLB held an NRC Research Associateship award at the National Institute of Standards and Technology. MK is supported by NSF through grant No. DMR-1708280.
References
- [1] Manfred Eigen. Selforganization of matter and the evolution of biological macromolecules. Naturwissenschaften, 58(10):465–523, 1971.
- [2] John Burdon Sanderson Haldane. Mathematical proceedings of the cambridge philosophical society. In A mathematical theory of natural and artificial selection, part V: selection and mutation, volume 23, pages 838–844. Cambridge University Press, 1927.
- [3] Holger Wagner, Ellen Baake, and Thomas Gerisch. Ising quantum chain and sequence evolution. Journal of statistical physics, 92(5-6):1017–1052, 1998.
- [4] Ellen Baake and Holger Wagner. Mutation–selection models solved exactly with methods of statistical mechanics. Genetics Research, 78(1):93–117, 2001.
- [5] Ellen Baake, Michael Baake, and Holger Wagner. Ising quantum chain is equivalent to a model of biological evolution. Physical Review Letters, 78(3):559, 1997.
- [6] Ellen Baake, Michael Baake, and Holger Wagner. Quantum mechanics versus classical probability in biological evolution. Physical Review E, 57(1):1191, 1998.
- [7] Joachim Hermisson, Holger Wagner, and Michael Baake. Four-state quantum chain as a model of sequence evolution. Journal of Statistical Physics, 102(1-2):315–343, 2001.
- [8] David Saakian and Chin-Kun Hu. Eigen model as a quantum spin chain: Exact dynamics. Physical Review E, 69(2):021913, 2004.
- [9] David S Rumschitzki. Spectral properties of eigen evolution matrices. Journal of Mathematical Biology, 24(6):667–680, 1987.
- [10] Ira Leuthäusser. Statistical mechanics of eigen’s evolution model. Journal of statistical physics, 48(1-2):343–360, 1987.
- [11] Stefano Galluccio. Exact solution of the quasispecies model in a sharply peaked fitness landscape. Physical Review E, 56(4):4526, 1997.
- [12] Subir Sachdev. Quantum Phase Transitions. Cambridge University Press, 2011.
- [13] Sei Suzuki, Jun-ichi Inoue, and Bikas K Chakrabarti. Quantum Ising phases and transitions in transverse Ising models, volume 862. Springer, 2012.
- [14] Thomas Jörg, Florent Krzakala, Jorge Kurchan, and AC Maggs. Simple glass models and their quantum annealing. Physical review letters, 101(14):147204, 2008.
- [15] Victor Bapst, Laura Foini, Florent Krzakala, Guilhem Semerjian, and Francesco Zamponi. The quantum adiabatic algorithm applied to random optimization problems: The quantum spin glass perspective. Physics Reports, 523(3):127–205, 2013.
- [16] CL Baldwin, CR Laumann, A Pal, and A Scardicchio. Clustering of nonergodic eigenstates in quantum spin glasses. Physical review letters, 118(12):127201, 2017.
- [17] Arijeet Pal and David A Huse. Many-body localization phase transition. Physical review b, 82(17):174411, 2010.
- [18] Rahul Nandkishore and David A Huse. Many-body localization and thermalization in quantum statistical mechanics. Annu. Rev. Condens. Matter Phys., 6(1):15–38, 2015.
- [19] Dmitry A Abanin, Ehud Altman, Immanuel Bloch, and Maksym Serbyn. Colloquium: Many-body localization, thermalization, and entanglement. Reviews of Modern Physics, 91(2):021001, 2019.
- [20] Simona Cocco, Christoph Feinauer, Matteo Figliuzzi, Rémi Monasson, and Martin Weigt. Inverse statistical physics of protein sequences: a key issues review. Reports on Progress in Physics, 81(3):032601, 2018.
- [21] Andrew L Ferguson, Jaclyn K Mann, Saleha Omarjee, Thumbi Ndung’u, Bruce D Walker, and Arup K Chakraborty. Translating hiv sequences into quantitative fitness landscapes predicts viral vulnerabilities for rational immunogen design. Immunity, 38(3):606–617, 2013.
- [22] John P Barton, Eleonora De Leonardis, Alice Coucke, and Simona Cocco. Ace: adaptive cluster expansion for maximum entropy graphical model inference. Bioinformatics, 32(20):3089–3097, 2016.
- [23] Francesca Rizzato, Alice Coucke, Eleonora de Leonardis, John P Barton, Jérôme Tubiana, Remi Monasson, and Simona Cocco. Inference of compressed potts graphical models. Physical Review E, 101(1):012309, 2020.
- [24] Simona Cocco and Rémi Monasson. Adaptive cluster expansion for inferring boltzmann machines with noisy data. Physical review letters, 106(9):090601, 2011.
- [25] John P Barton, Mehran Kardar, and Arup K Chakraborty. Scaling laws describe memories of host–pathogen riposte in the hiv population. Proceedings of the National Academy of Sciences, 112(7):1965–1970, 2015.
- [26] Karthik Shekhar, Claire F Ruberman, Andrew L Ferguson, John P Barton, Mehran Kardar, and Arup K Chakraborty. Spin models inferred from patient-derived viral sequence data faithfully describe hiv fitness landscapes. Physical review E, 88(6):062705, 2013.
- [27] Jaclyn K Mann, John P Barton, Andrew L Ferguson, Saleha Omarjee, Bruce D Walker, Arup Chakraborty, and Thumbi Ndung’u. The fitness landscape of hiv-1 gag: advanced modeling approaches and validation of model predictions by in vitro testing. PLoS computational biology, 10(8), 2014.
- [28] Raymond HY Louie, Kevin J Kaczorowski, John P Barton, Arup K Chakraborty, and Matthew R McKay. Fitness landscape of the human immunodeficiency virus envelope protein that is targeted by antibodies. Proceedings of the National Academy of Sciences, 115(4):E564–E573, 2018.
- [29] H. Chen and M. Kardar. bioRxiv:518704, 2019.
- [30] Vipul Gupta and Narendra M Dixit. Scaling law characterizing the dynamics of the transition of hiv-1 to error catastrophe. Physical biology, 12(5):054001, 2015.
- [31] Gregory R Hart and Andrew L Ferguson. Error catastrophe and phase transition in the empirical fitness landscape of hiv. Physical Review E, 91(3):032705, 2015.
- [32] Viviana Simon, Nicolin Bloch, and Nathaniel R Landau. Intrinsic host restrictions to hiv-1 and mechanisms of viral escape. Nature immunology, 16(6):546, 2015.