INTRODUCTION Darwin meets Perutz and Kendrew In evolutionary dynamics, we consider the manner in which mutations occur at the genotypic level, the phenotypic consequences, how that affects the fitness of the organism in its environment, and how this influences the fate of this mutation. It was clear to Darwin that the environment of any organism was largely determined by the properties, behaviours, and interactions with other organisms, what Darwin referred to as the ‘tangled bank’, composed of ‘many plants of many kinds, with birds singing on the bushes, with various insects flitting about, and with worms crawling through the damp earth … and dependent upon each other in so complex a manner’. This network of interactions extends over the entire range of time and length scales, from the creation of the oxygen-rich atmosphere by bacteria, to the selection for communal behaviour found in termites, mole rats, and humans, to the ‘red queen’ dynamics of hosts and their parasites, to the evolution of genes and their gene products. In particular, as soon as Perutz and Kendrew demonstrated the high-resolution structure of myoglobin, it was obvious that proteins are structured, highly integrated entities formed of linked amino acids. Selection acts on proteins based on their ability to function and (for most proteins) their capacity to fold and be stable under physiological conditions. These selective constraints are holistic constraints, that is, they are properties of the entire protein acting as a unit, explicitly including arrangements of amino acids at different sites. Structure and stability involve the construction of hydrogen bonds, ion pairs, packing interactions, and the formation of hydrophilic surfaces. Function involves forming binding sites and specific geometries of catalytic residues. The changes that occur at a given location in the protein can only be understood in the context of the rest of the protein. Changes in selection can result from shifts in function, structure, or physiological role. They can also result from interactions between sites in the protein; just as the evolutionary dynamics of a pathogen influences the selection acting on the immune system of the host, so does the evolutionary changes that occur in one part of the protein has an effect on the evolutionary changes that occur at another location (as well as in other genes, gene products, and regulatory regions). And as with the tangled bank, there will be a complex network of feedback loops; substitutions at one focal position will affect the selection acting on other sites; the corresponding substitutions at these other sites will, in turn, influence the selection acting on the focal position. One aspect of this is the so-called ‘Stokes Shift’, where the rest of a protein adopts itself to the new amino acid resident at a given site so that, as long as that amino acid remains at the site, it becomes increasingly suitable for that position. These interactions result in great complexity, so it is not surprising that they are ignored, neglected, and disregarded in most modelling of protein evolution. The result was the production of many powerful phenomenological approaches, where each location in a protein was assumed to evolve in a manner independent (and often identical) to the evolution at other locations. Heuristic approaches were brought in that encoded the effect of these interactions, while still preserving the independent-site conditions. Often the heuristic nature of these approaches was forgotten, and a connection was lost to the basic biological phenomena underlying these models. As sequence data and computational methods become more advanced, we have the opportunity to go beyond these simple phenomenological models, to reconnect with the underlying mechanisms of molecular evolution. At the same time, it has become increasingly clear how important these epistatic interactions are between locations in proteins, how much information is available by modelling these interactions, and how perilous it is to ignore them. In this paper, we first attempt to unify some of the current mechanistic and phenomenological models, and to make a connection between these models and the underlying ‘tangled bank’ reality. We when we describe a different approach to conceptualising the process of protein evolution, one explicitly considering the effects of protein epistasis on the substitution process at a particular site. RESULTS Fisher and Dayhoff meet Kimura and Bruno It has been observed many times that mutations with small effect are more likely to be accepted than mutations with large effect. This is observed, for instance, in the standard substitution matrices created by Dayhoff and others, where the rate QIV of the conservative substitution from isoleucine to valine is much higher than the rate QCP for the non-conservative substitution from cysteine to proline. There was an inherent symmetry in these models; both isoleucine to valine and valine to isoleucine transitions would be fast, while both cysteine to proline and proline to cysteine substitutions would be rare. This seems a simple consequence of an argument first made by Fisher, who considered the focusing of a microscope; if the sample in the microscope was more or less in focus, turning the fine adjustment knob is much more likely to yield an improvement than turning the course adjustment knob. If you are near a fitness optimum, large steps will tend to move you away from the optimum, not nearer to the peak. Fisher than predicts that conservative substitutions, as they are more likely to cause a small change in the protein, are more likely to be accepted than non-conservative substitutions. Protein biophysicists and bioinformaticians tend to see changes in amino acid not in terms of relative distance (how similar or dissimilar two amino acids are), but how they compare on an external scale, how appropriate the two different amino acids are for a given site. (This is, for instance, the principle behind hidden Markov models used to classify and cluster protein sequences.) In addition, this appropriateness is specific to each individual site, with locations preferring hydrophobic or helix-forming or aliphatic or specific amino acids depending upon structural and functional constraints. In contrast, standard substitution models (such as the Dayhoff matrix, as well as other popular models such as JTT, WAG, LG) generally assume that the same model applies at all locations in the protein at all points in evolutionary time (or consider a site-specific scaling factor, that does not affect relative rates). Bruno and collaborators developed this alternative conception in his mutation-selection models, by representing the substitution rate in terms of the mutational rate times the fixation probability dependent on the relative ‘fit’ between the amino acids at that position. In this model, conservative changes would be close to neutral, and would be accepted at rates similar to the neutral rate. Non-conservative substitutions that result in more significant changes in the protein would be either slower or faster, depending upon whether the changes result in deleterious or advantageous changes in the protein, using the fixation probability formalism of Kimura. As a result, if a non-conservative change from amino acid X to Y was highly deleterious, and thus was unlikely to be accepted, resulting in a very low substitution rate QXY, we would expect the non-conservative change from Y to X, QYX, to be substantially higher than neutral, not lower. Site-heterogeneity was explicitly included in these models, with the relative fitnesses of the amino acids defined at each location. It is exactly this variation in rates that can reconcile the standard substitution models with the biophysical description. Instead of rates, let us consider the fluxes between the different amino acids, that is the number of times a given transition occurs. For instance, the expected flux of substitutions from cysteine to proline would be equal to πC QCP, where πC is the equilibrium frequency of cysteine. We will, in general, observe few changes from an appropriate amino acid X to an inappropriate amino acid Y, as this would have a small fixation probability and thus a small QXY. We would, also, observe few changes from an inappropriate amino acid Y to an appropriate amino acid X, not because the rate of the transition is slow, but rather because the probability of having inappropriate amino acid Y at that site (πY) would be small. The largest flux would be between relatively appropriate amino acids, where the equilibrium frequencies of the initial amino acid is high but where the rate of substitutions to the new amino acid is sufficiently fast. We would then observe predominantly substitutions between relatively appropriate amino acids at each position. According to the biophysical model, the selective constraints, and thus the relative equilibrium frequencies and substitution rates, will depend on site in the protein and may change over time. The standard substitution models are computed based on the observed number of substitutions observed, that is, on the product of πX and QXY. With selection constraint heterogeneity, the observed number of substitutions would depend on , averaged over the different sites and the different evolutionary times of the dataset, with the estimated rate given by the flux divided by the equilibrium frequency, . We can write this as (1.1) where is the covariation between πX and QXY. Note that the second term on the right can be equal in magnitude (but potentially of opposite sign) to the first term.