on Neutral Evolution of Model Proteins Di4usion in Sequence Space and Overdispersion

J.theor.Biol.(1999)200,49}64Article No.jtbi.1999.0975,available online at onNeutral Evolution of Model Proteins:Di4usion in Sequence Space and OverdispersionU GO B ASTOLLA*-,H.E DUARDO R OMAN?AND M ICHELE V ENDRUSCOLO A*H¸RZ,Forschungszentrum Ju lich,D-52425Ju lich,Germany,?Dipartimento di Fisica and INFN,;niversita di Milano,I-20133Milano,ItalyA Department of Physics of Complex Systems,=eizmann Institute of Science,Rehovot76100,Israel(Received on12November1998,Accepted in revised form on21May1999)We stimulate the evolution of model protein sequences subject to mutations.A mutation isconsidered neutral if it conserves(1)the structure of the ground state,(2)its thermodynamicstability and(3)its kinetic accessibility.All other mutations are considered lethal and arerejected.We adopt a lattice model,amenable to a reliable solution of the protein foldingproblem.We prove the existence of extended neutral networks in sequence space*sequencescan evolve until their similarity with the starting point is almost the same as for randomsequences.Furthermore,we"nd that the rate of neutral mutations has a broad distribution insequence space.Due to this fact,the substitution process is overdispersed(the ratio betweenvariance and mean is larger than1).This result is in contrast with the simplest model of neutralevolution,which assumes a Poisson process for substitutions,and in qualitative agreementwith the biological data.1999Academic Press1.IntroductionA recent study on the Protein Data Bank(PDB) shows that the distribution of pairwise sequence identity between structurally homologous pro-teins presents a large peak at8.5%sequence identity,only slightly larger than what expected in the purely random case B(Rost,1998).This is -Author to whom correspondence should be addressed. Present address:Freie Universita t Berlin,FB Chemie, Takustr.6,D-14195Berlin,Germany.E-mail:ugo@chemie. fu-berlin.de.B The number of amino acid matches obtained by pairing two random sequences of the same length is given by the binomial distribution with p" if one assumes that the20 amino acids have the same probability to occur.For se-quences of length N there will be on average pN identical amino acids,with a variance Np(1!p).For random se-quences,95%of pairwise comparisons yield a sequence identity between1and9%.an interesting result which means that the struc-tural similarity does not imply sequence similarity. An intensive computational study on second-ary structures of RNA molecules(Schuster et al., 1994),which is a problem much simpler than protein folding,and can be studied through e$-cient and reliable algorithms,showed that an exponentially large number of sequences corres-ponds in average to a single structure,and the distribution of structures in sequence space is quite inhomogeneous(it follows a Zipf law).Se-quences folding into the most common structures form connected&&neutral networks''that perco-late sequence space.These neutral networks directly arise from the non-uniqueness of the relation between sequence and structure. These results are important to understand how evolution works at the molecular level.Kimura0022}5193/99/017049#16$30.00/0 1999Academic Press(1968,1983)and Jukes&King(1969)proposed a long time ago that most of the evolutionary events at the molecular level consist in the substi-tution of one allele with another one selectively equivalent(neutral evolution).A keystone of neu-tral evolution is the existence of a&&phenotypic threshold'':below a critical value of structural similarity,natural selection cannot avoid the spreading of variants of the protein in a popula-tion through&&di!usion'',called&&genetic drift''by genetists.Neutral evolution in sequence space conserves the structure of the protein(thus taking place below the phenotypic threshold)but can drive to sequence similarities expected for ran-domly chosen sequences.The purpose of this paper is to investigate these subjects in a simpli"ed model of protein folding.Our aim is to show that,as in the case of RNA sequences,protein sequence space is traver-sed by vast neutral networks.The attempt to study this problem requires some methodological choices.The"rst and most severe condition is that,given a sequence of amino acids,we should be able to identify its native conformation.This is tantamount of solv-ing the protein folding problem.We decided to adopt a lattice model of proteins,since,to date, this is the only case in which protein folding is routinely solvable.This choice restricts drasti-cally the possibility to represent active sites and structural motifs,which might play an important role in constraining neutral evolution.The second choice is how to represent muta-tions.We mutate amino acids and not nucleo-tides.In nature,mutations a!ect DNA,so that mutations happen with larger frequency between amino acids whose codes di!er by just one nu-cleotide.We do not take this fact into account and consider an uniform probability of mutation for any nucleotide.This does not change the characteristics of the&&neutral networks'',and it should not change the relevant features of di!u-sion in sequence space.More important is the fact that we work with sequences of"xed length (number of amino acids)N"36,considering only point mutations,without insertion and dele-tions that would make the model much more complicated.The third choice is how to represent natural selection.We introduce a stochastic algorithm which at each time step attempts to change one amino acid.The mutation is accepted if three conditions are met:the ground state should be conserved,thermodynamically stable and easily accessible.Since we are limited by computational re-sources,we decided to study one particular neutral network,namely,we followed the evolution from one single common ancestor.The extent to which this network can be considered paradigmatic will be commented.We believe that the features of neutral evolu-tion which are found in real proteins are universal enough to be captured by the present simpli"ed model.Support for this view comes from a recent study by Babajide and co-workers (Babajide et al.,1997).They also found evidence for the presence of neutral networks in sequence space.Their work is similar in spirit to the pres-ent one,but signi"cantly di!erent methodologi-cally.They used a description of the protein structure based on the C?and C@coordinates as taken from the PDB,and an approximate cri-terion of fold recognition based on the Z score (Bowie et al.,1991;Casari&Sippl,1992).The paper is organized as follows:in Section2, we present an outline of the theory of neutral evolution.In Section3,we describe our model protein and our protocol to simulate neutral evolution.In Section4,we illustrate our results concerning the genetic drift.In Section5,we deal with population genetic considerations and study the dispersion index of the process of substitu-tions.Section6presents an overall discussion, where we discuss the applicability of our results to the interpretation of existing biological data.2.The Neutral Theory of EvolutionThe theory of evolution explains how di!erent living species evolve adapting to an ever chang-ing environment.The key ingredients of the explanation are reproduction,mutation and natural selection,and the complex ecological in-teractions that modify the environment.The pat-tern of change in time at the phenotypic level is quite irregular:long periods of stasis are followed by abrupt bursts of activity,species suddenly appear and disappear and very large extinction events happen,showing complex features(Gould50U.BASTOLLA E¹A¸.&Eldredge,1977;Sole et al.,1997).In contrast to this,evolution at the molecular level seems to be more regular.The"rst studies on this subject showed that the rate of change in time of the amino acid sequence of a given protein,which can be estimated from the di!erence of homolog-ous sequences of di!erent species,does not vary much from one species to another in the same class,even if the species compared have rather di!erent population size and environmental conditions(Zuckerland&Pauling,1962)(on the other hand,the evolution rate changes drastically from one protein to another).This property makes each protein a sort of&&molecular clock'' (Zuckerland&Pauling,1962;Kimura,1983; Gillespie,1991;Ratner et al.,1996),with its own measure of time,and allows to transform a dis-tance in sequence into a distance in time and to reconstruct phylogenetic trees from molecular data.In order to explain this and other features of molecular evolution,as for example the presence of many ploymorphic loci in the genome,30years ago Kimura(1968)and Jukes&King(1969) proposed a new interpretation of molecular evolution,that was named the neutral theory of molecular evolution,reviewed in Kimura(1983). According to this theory,most of the changes in protein sequences happen not because better variants of the protein are found but because many mutations do not modify signi"cantly the e$ciency of the protein,so that natural selection cannot avoid their spreading through the popu-lation by random genetic&&drift''(&&di!usion''in physical language).There have been,however,more recent and accurate studies which question the very exist-ence of molecular clocks,noting that most of the clocks underwent in some instances(i.e.in some periods of the evolution of some lineages)drastic accelerations or decelerations(Ayala,1997)so that the validity of the clock hypothesis is du-bious.In some cases,the adaptive nature of the acceleration of the substitution rate was clearly demonstrated,like in the case of the changes that hemoglobin underwent in the time when verte-brates colonized the earth(Ratner et al.,1996). Some authors believe that these discrepancies of the clocks are such that the neutral hypothesis should be completely rejected.This dispute pro-duced a vast literature.See the book by Gillespie (1991)for a nice and detailed discussion of the selectionist point of view.One decade after Kimura's proposal,the the-ory of the molecular quasi-species(Eigen et al., 1989)showed the possible existence of an evolu-tionary phase,at high mutation rate and low selective pressure,with the di!usive features of neutral evolution:thus this mode of evolution is expected not only for strictly neutral mutations, but whenever the damage brought by frequently occurring mutations is smaller than a threshold depending on the mutation rate.Neutral evolu-tion RNA molecules,whose secondary structure can be easily computed,was studied in detail by the Vienna group(Schuster et al.,1994;Huynen et al.,1996;Fontana&Schuster,1998).The prop-erties of neutral networks in sequence space have been characterized(Schuster et al.,1994).Adap-tive evolution of a population of RNA molecules towards a target structure has been simulated (Huynen et al.,1996).This alternates periods of stasis where evolution happens on neutral net-works of phenotypically equivalent structures and fast transitions towards better adapted struc-tures.The population dynamics in the periods of stasis is equivalent to neutral evolution and has the features of di!usion in sequence space.The transitions between di!erent structures suggest to de"ne and a metric in structure space related to the contiguity of neutral networks of di!erent structures(Fontana&Schuster,1998).With regard to proteins,Shakhnovich and Gutin estimated,in the framework of their ran-dom heteropolymer model,the probability that a mutation does not change the ground state (Shakhnovich&Gutin,1991).Their conclusion is that this probability has a"nite limit for increas-ing system size.A more recent study considering a lattice model and an optimized sequence showed that indeed the probability that a muta-tion is neutral can be as high as70%(Tiana et al., 1998).A numerical study of neutral networks in sequence space was performed by Bornberg-Bauer(Bornberg-Bauer,1997).He considered a two-dimensional lattice model with two amino acid types(HP)and sequences of length up to18, for which exact enumeration allows to deter-mine exactly the correspondence between se-quences and structures.He found out that someNEUTRAL EVOLUTION OF MODEL PROTEINS51structures are very frequent(indeed,the ranking of structures according to the corresponding number of sequences follows a generalized Zipf law)but,contrary to the case of RNA,the subsets of sequences corresponding to the most frequent structures are localized in sequence space.This is probably an artifact of the HP model:our study, considering20amino acid types,reaches the op-posite conclusion.These"ndings suggest the question of how many di!erent species of amino acids are needed to create the possibility of ex-tended neutral networks in sequence space.An interesting conclusion of Bornberg}Bauer is that neutral networks of di!erent structures are un-likely to be connected by very few point muta-tions.This conclusion seems to hold also for the 20amino acid model presented in this work.In this work,we investigate with a minimal model the possibility of neutral evolution of proteins.We"nd that the model supports this possibility**.Interestingly,we"nd also that the rate R between the variance of the number of substitutions and its average value is larger than 1.Originally,Kimura proposed that the substitu-tion of amino acids in proteins follows a Pois-sonian stochastic process,which implies R"1. For most proteins,however,real data give R sig-ni"cantly larger than1,till values between30and 50(Gillespie,1991).This problem is produced by the drastic accelerations of the substitution pro-cess mentioned above,and has been regarded as a severe evidence against the neutral theory. Thus it is remarkable that our&&microscopic'' model leads to similar results considering only neutral evolution.This fact suggests that the neu-tral theory can be reconciled with at least part of the data.**This possibility does not imply that neutral evolution should occur.In particular,our model describes neutral evolution with the additional hypothesis that the environ-ment remains constant.Some selectionist theories,on the other hand,consider a rapidly changing environment,so that amino acids substitutions are needed to"t the protein to the new environment(Gillespie,1991).While it is un-debatable that dramatic climatic and ecological changes happen on the time-scale of molecular evolution,it seems not unreasonable that the homeostatic properties of cells protect the cellular environment against such changes.In-deed,the most irregularly(thus less&&neutrally'')evolving proteins seem to be hormones(Gillespie,1991),which are responsible for intercellular communication,while most enzymes evolve in a more regular fashion.Finally,a remark about terminology:by &&mutation''we mean the modi"cation of a triplet of the genetic code in one individual lineage(by means of environmental damages,errors in a rep-lication process or other causes).By&&substitu-tion''we mean the much more complex process inwhich an allele dominant in a biological popula-tion is replaced by a new one arising from muta-tion.We consider here only individual lineages. Passing from the level of individual lineages(mutation and reproduction)to the level ofa population(substitution)requires the tools ofpopulation genetics.However,Kimura(1968,1983)showed that the passage is very simple inthe case in which the only mutations occurring are either lethal or completely neutral.In thiscase,the substitution rate in the population isQS@Q" x,(1) where is the&&microscopic''mutation rate(at thelevel of the single lineage)and x is the fraction of mutations which are selectively neutral[see alsoBastolla&Peliti(1991)].Note that the rate ofevolution does not depend on the size of thepopulation,whereas it does if positive naturalselection plays a role.Under these assumptions,the model of neutral evolution of single lineages that we present here is equivalent to a model ofevolution at the level of a population,and thus itcan be in principle compared with the biologicaldata.3.Model of Neutral Evolution in ProteinsIn this section,we describe in detail the latticemodel used to represent protein structures andour model of evolution in sequence space.3.1.DEFINITION OF THE MODEL PROTEINSWe represent protein con"gurations as self-avoiding walks on the simple cubic lattice.Eachof the N occupied sites represents an amino acid,chosen among the20possible ones.We associate each con"guration of the protein to its contactmap C.This is a N;N matrix whose element CGH is1if residues i and j are nearest neighbors on thelattice(but not along the chain)and0otherwise.Note that the correspondence between contact52U.BASTOLLA E¹A¸.maps and con "gurations is not unique,and the smaller is the &&density of contacts ''N A/N of a con-tact map,where N A " G H C GHis the number ofcontacts,the larger is the number of con "gura-tions to which it corresponds (Vendruscolo et al .,1999).We measure similarity in structure space through the overlap q (C ,C )between the contact maps C and Cq (C ,C )"1N*A G HC GH C GH ,(2)where N*A is the maximal between N A and N A,thenumber of contacts,respectively,of two contact maps C and C .Two maps are identical if and only if q "1.We study sequences of length N "36,denoted by S "+s ,2,s ,,, G3+1,2,20,.The energyE (S ,C )of con "gurations of sequence S whose contact map is C is de "ned asE (S ,C )" ,G HC GH ;(s G ,s H ),(3)where ;(a ,b )is a 20;20symmetric interaction matrix,which gives the energy gain obtained when amino acids of species a and b are brought into contact.We use a matrix ;(a ,b )derived from the Miyazawa }Jernigan interaction matrix (Miyazawa &Jernigan,1985).We identify the native structure with the ground state of the model.We choose as native state the structure represented in Fig.1.Since its contact map C *has the highest number of con-tacts for this chain length,N A"40,it de "nesuniquely the structure apart from trivial symmet-ries.We note that C *is an atypical structure for the interaction parameters that we choose:since ;(a ,b )has an average value zero and variance 0.3,we would expect open structures to be ener-getically favored.Indeed,typical random se-quences with N "36and contact interactions whose average vanishes have a ground state with approximately 30}32contacts (Bastolla,unpub-lished result),being thus less than maximally compact.The contact map C*was studied by Shak-hnovich and co-workers in a computer experi-ment of inverse folding (Abkevich et al .,1994).F IG .1.The &&native state ''of our model protein.They designed a sequence S *with ground state on C *using the procedure of Shakhnovich &Gutin (1993),showing that S *has good prop-erties of kinetic foldability and thermodynamic stability at the temperature where the folding is the fastest.The lower part of the energy land-scape of this sequence is remarkably smooth:all the structures with low energy have a high over-lap qwith the ground state.The lowest energy ofcon "gurations with a "xed value of qdecreasesregularly as qapproaches one.This correlatedenergy landscape,reminiscent of the &&funnel ''paradigm (Bryngelson &Wolynes,1987),is the reason of the good folding properties of the se-quence,which is very di !erent from a random one.In Tiana et al .(1998)it was shown that the sequence S *is also very stable against mutations.It was estimated that about 70%of the point mutations performed on S *result in sequences with the same ground state and good folding properties.Thus,energy minimization makes C *stable not only in structure space,but also in sequence space.We now de "ne three measures of distance in sequence space.The "rst one,directly related to point mutations,is the Hamming distance,which counts how many amino acids are di !erent between two sequences S and S of the same length N :D (S ,S )",G[1! (s G ,s G)],(4)NEUTRAL EVOLUTION OF MODEL PROTEINS53where is the Kronecker symbol and sG takes20di!erent values.Equivalently,one can use the overlap Q(S,S ),which is equal to one minus the normalized Hamming distance D/N:Q(S,S )"1N,G(s G,s G).(5)If we divide the20amino acid types in just two classes,eight hydrophobic ones and12hydro-philic ones,we obtain as the corresponding Ham-ming distance the&&hydrophobicity distance''D &.(S,S ),de"ned as above but withG3+1,2,.The third distance that we use is D3(S,S ),whichmeasures the di!erence in the native interactions. This is given by the quadratic di!erence of the native interactions in two di!erent sequences:D 3(S,S )"G HC*GH[;(SG,SH)!;(SG,SH)] ,(6)where C*is the native contact map.3.2.MODEL OF EVOLUTIONWe consider a reference structure C*as the biologically active native structure,and we set as starting point of neutral evolution a sequence S* which folds onto C*.The simulated molecular evolution is realized by the following iterative procedure:1.At t"0we start from S(0)"S*.2.At time step t we mutate at random oneamino acid of S(t!1),producing a newsequence S (t).3.We submit the new sequence to selectionaccording to the criteria speci"ed below.Ifthe sequence survives we accept the muta-tion(S(t)"S (t)),otherwise we restoreS(t)"S(t!1).The selection of sequences is governed by the following conditions:Conservation of the phenotype The ground state C of S (t)must have an overlap with C*not smaller than a&&phenotypic threshold''qRFP.q(C*,C)*qRFP ,(7)In our calculations we imposed strict conserva-tion of the phenotype by setting qRFP"1.¹hermodynamic stability We de"ne thermodyn-amic stability through the condition1q(C*,C)2*1q2RFP,(8)where1)2represent a Boltzmann average at thetemperature¹of the simulation and1q2RFPisa"xed parameter.This condition implies that allthe thermodynamically relevant states have anoverlap greater than1q2RFPwith the target state.Kinetic accessibility We use the PERM method(Grassberger,1997;Frauenkron et al.,1998;Bastolla et al.,1998),a new Monte Carlo algo-rithm which is particularly successful in"ndingthe ground state of lattice polymers.For everysequence S to be tested,we run iterations of thealgorithm.If the lowest energy structure found atthis point does not coincide with C*,we discardthe sequence because the kinetic accessibilitycondition is violated.Otherwise,we continue torun the algorithm for a time to check that nolower energy structure is found.If this conditionis met,we start again the Monte Carlo algorithmand run it for a time2 .Only if C*is once againfound as the lowest energy structure we concludethat the ground state of S really coincides withC*.Note that there is no bias towards C*in ourMonte Carlo algorithm.We never found in the second independent runof the MC algorithm a structure with lower en-ergy than the putative ground-state C*.This factencourages us in believing that the algorithm wasreally e!ective in"nding the ground state.An-other support to this conclusion comes from thefact that all of the selected sequences had a re-markably correlated energy landscape,whichmakes the search for the ground state easier.Onthe other hand,when a sequence is rejected,weare not sure whether we were able to determineits ground state.The di!erence is due to tworeasons:"rst,we run rejected sequences fora shorter time( instead of2 is the rejection ismade at the"rst decision stage,1run instead of2independent ones if it is made at the secondstage).Second,rejected sequences have typically54U.BASTOLLA E¹A¸.a much less correlated energy landscape,so that it is expected that the determination of the ground state is more di $cult.Nevertheless,we shall discuss at the end also data about rejected sequences,since they are interesting and refer to a very large number of sequences.The three conditions for the acceptance of a mutation enforce the conservation of biological activity of the new sequence.We have to stress that conservation of the fold is not a necessary condition for selective neutrality in real world,and it is not even su $cient,since the active site has also to be conserved and the environment has to remain reasonably stable.Thus with our model we represent the neutral evolution of the part of the chain not involved in chemical activ-ity,in a stable chemical environment.Never-theless,we think that these limitations do not prevent from the applicability of the model to real molecular evolution.3.3.GENETIC DRIFTFor a given sequence S of N amino acids,wede "ne the rate of neutral mutation x as the frac-tion of acceptable non-synonymous mutationsx (S )"120N, GO s G?G (S ),(9)where ?G(S )"1if assigning the amino acid ofspecies at position i on the sequence S does not change the native state and 0otherwise.In Kimura 's (1983)ordinary neutral theory it is assumed that x (S )"x is indeed independent of the sequence.With this hypothesis,the evolu-tion of the overlap with the starting sequence,Q (t )"Q (S (0),S (t ))(the overbar denotes an average with respect to the mutation process)is given byQ (t )+ #(1!)e \VR ,.(10)Since every mutation is either neutral or lethal,the evolution of the wild-type genome of the population coincides with that of a single repro-ductive line (Kimura,1983),and we may interpret di !erent realizations of our evolutionary process as di !erent &&species ''originating from a commonancestor.Thus QS@Q"x /N is the substitution rateboth of a single lineage and of the population,and the time t represents in our model the num-ber of mutational events.As we shall see,how-ever,the hypothesis of constance of x (S )is not in agreement with the results of the simulations.3.4.NEUTRAL SETWe de "ne the neutral network (C *,S *)as follows: (C *,S *)is formed by all the sequences S that have C *as their ground state and are connected to S *through a neutral path (a path in sequence space which satis "es the three selection criteria given above).We then de "ne the neutral set (C *)as(C *)"8I (C *,S *I)(11)with running over all disconnected neutral net-works.If more than one disconnected componentof the neutral set exists it would be possible to distinguish between convergent (di !erent ances-tors,same "nal structure)and divergent (same ancestor)evolution:if the evolution of two homo-logous proteins takes place in two disjoint neu-tral networks,it must be convergent.We studied only one starting sequence S *,thus this work is concerned only with divergent evolu-tion.However,since the neutral networks that we studied is so spread that typical pairs of se-quences have a Hamming distance comparable with that of random sequences,it turns out that it is not possible to distinguish between divergent and convergent evolution on the basis of the sequence homology alone.Under this respect,our results are equivalent to those of Rost (1997).4.Results4.1.EVOLUTIONARY DRIFTWe generated eight trajectories with the following values of the selection parameters:q RFP "1,1q 2RFP"0.9,¹"0.16, "1.6;10 .In each trajectory,a number of sequences ranging from 875to 2709were tested.The CPU time needed for a single trajectory was about 4}5weeks on Sun Ultra stations.Di !erent realizations di !erentiated greatly from the &&common ancestor ''S (0)and from eachNEUTRAL EVOLUTION OF MODEL PROTEINS55。