Table Of ContentMean-fieldmethodsinevolutionaryduplication-innovation-lossmodelsforthe
genome-levelrepertoireofproteindomains
A. Angelini,1 A. Amato,1 G. Bianconi,2 B. Bassetti,1 and M. Cosentino Lagomarsino1
1Universita` degli Studi di Milano, Dip. Fisica. Via Celoria 16, 20133 Milano, Italy
2Physics Department, Northeastern University, 111 Forsyth St.,
111 Dana Research Center, Boston 02115, MA, USA
(Dated: January22,2010)
Wepresentacombinedmean-fieldandsimulationapproachtodifferentmodelsdescribingthedy-
namicsofclassesformedbyelementsthatcanappear,disappearorcopythemselves. Thesemodels,
0
1 related to a paradigm duplication-innovation model known as Chinese Restaurant Process, are de-
0
visedtoreproducethescalingbehaviorobservedinthegenome-widerepertoireofproteindomains
2
n of all known species. In view of these data, we discuss the qualitative and quantitative differences
a
ofthealternativemodelformulations,focusinginparticularontherolesofelementlossandofthe
J
2 specificityofempiricaldomainclasses.
2
]
N
I. INTRODUCTION
G
.
o
Understanding the complexity of genomes and the drives that shape them is a fundamental problem
i
b
- of contemporary biology [1], which poses a number of challenges to contemporary statistical mechanics.
q
[
Considering this problem from a large-scale viewpoint, the basic observables to account for are the distri-
2
v butions of the different “functional components” (such as genes, introns, non-coding RNA, etc.) encoded
3
5 bysequencedgenomesofvaryingsize.
2
4 Whenthesegenome-widedataareparametrizedbymeasuresof“genomesize”(suchasthenumberof
.
0
basesorthenumberofgenesinagenome),thereareimportantemerging“scalinglaws”bothfortheclasses
1
9
ofevolutionaryrelatedgenes[2–4],thefunctionalcategoriesofgenes[5,6]andsomenon-codingpartsof
0
:
v genomes [7]. These scaling laws are the signs of universal invariants in the processes and constraints that
i
X
gave rise to the genomes as they can be observed today. A current challenge is the understanding of these
r
a
laws using physical modeling concepts and the comparison of the models to the available whole-genome
data. Thiseffortcanhelpdisentangleneutralfromselectiveeffects[8,9].
Here we consider the statistical features of the set of proteins expressed by a genome, or proteome. A
convenient level of analysis is a description of the proteome in terms of structural protein domains [10].
Domainsaremodular“topologies”,orsub-shapes,formingproteins[11]. Adomaindeterminesasetofpo-
tentialbiochemicalorbiophysicalfunctionsandinteractionsforaprotein,suchasbindingtootherproteins
2
or DNA and participation in well-defined classes of biochemical reactions [6, 10]. Despite the practically
unlimitednumberofpossibleproteinsequences,therepertoireofbasictopologiesfordomainsseemstobe
relativelysmall[12]. Withalooseparallel,domainscouldbeseenasan“alphabet”ofbasicelementsofthe
proteinuniverse. Understandingtheusageofdomainsacrossorganismsisasimportantandchallengingas
decodinganunknownlanguage.
Thecontentofagenomeisdeterminedprimarilybyitsevolutionaryhistory,inwhichneutralprocesses
and naturalselection play interdependent roles. In particular, thecoding parts of genomesevolve by some
well-defined basic “moves”: gene loss, gene duplication, horizontal gene transfer (the transfer of genetic
material between unrelated species), and gene genesis (the de novo origin of genes). Since domains are
modular evolutionary building blocks for proteins, they are coupled to the dynamics followed by genes.
In particular, a new domain topology can emerge by genesis or horizontal transfer, and new domains of
existingdomaintopologiescanemergebyduplicationorbelost. Finally,topologiescanbecompletelylost
byagenomeifthelastdomainthatcarriesthemislost(seeFig.1).
Large-scale data concerning structural domains are available from bioinformatic databases, and can be
analyzed at the genome level. These coarse-grained data structures can be represented as sets of “domain
classes” (the sets of all realizations of the same domain topology in proteins), populated by domain re-
alizations. In particular, much attention has been drawn by the intriguing discovery that the population
of domain classes have power-law distributions [13–17]: the number F(j,n) of domain classes having j
members follows the power-law ∼ 1/jz, where the exponent z typically lies between 1 and 2. An inter-
estingthreadofmodelingworkascribestheemergenceofpower-lawstoagenericpreferential-attachment
principle due to gene duplication. Growth models are formulated as nonstationary, duplication-innovation
models[13,18,19]andasstationarybirth-death-innovationmodels[14,20–22].
Intriguingly,aswehaverecentlyshown[23],thedomaincontentofgenomesalsoexhibitsscalinglaws
as a function of the total number of domains n, indicating that even evolutionarily distant genomes show
commontrendswheretherelevantparameteristheirsize.
(i) The number of domain classes (or distinct hits of the same domain) concentrates around a master
curveF(n)thatappearstobemarkedlysublinearwithsize,perhapssaturating.
(ii) Thefittedexponentofthepower-law-likedistributionF(j,n)ofdomainclasseshavingj members,
in a proteome of size n decreases with genome size. In other words, there is evidence for a cutoff
thatincreaseslinearlywithn.
(iii) The occurrence of fold topologies across genomes is highly inhomogeneous - some domain super-
families are found in all genomes, some rare, with a sigmoid-like drop between these two cate-
3
gories[35].
Werecentlyreportedtheabovecollectivetrends,andshowedhowthescalinglawsinthedatacouldbe
reproduced using universal parameters with non-stationary duplication-innovation models. Our results in-
dicatethatthebasicevolutionarymovesthemselvescandeterminetheobservedscalingbehaviorofdomain
content,aprioriofmorespecificbiologicaltrends.
Thismodelingapproach,whilesimilarinformulationtothatofpreviousinvestigators[13]whodidnot
considerthesescalinglaws,hasimportantmodifications,mostlyrelatedtothescalingwithnoftherelative
probabilityofaddingadomainbelongingtoanewclassandduplicatinganexistingone. Toreproducethe
observed trends, newly added domain classes cannot be treated as independent random variables, but are
conditionedbythepreexistingproteomestructure.
In this paper, we give a detailed account of this modeling approach, considering different variants of
duplication-innovation-loss models for protein domains, and relate them to available results in the mathe-
matical and physical literature. In particular, we will focus on mean-field approaches for the models and
comparison with direct simulation, and we will show how they can be generally used to obtain the main
qualitativeandquantitativetrends.
Thefirstpartofthepaperisdevotedtotheminimalmodelformulation,whichonlyincludesduplication
and innovation moves for domains, and relates to the so-called Chinese Restaurant Process (CRP) of the
mathematicalliterature[24,25]. Wewillreviewthemainknownresultsforthismodel,deriveanalytically
solvable mean-field equations, and show how they compare to the available rigorous results and the finite
sizebehavior.
The rest of the paper is devoted to biologically motivated variants of the main model related to two
mainfeatures: includingtheroleoflossofdomains,whichisafrequentlyreportedevent,andbreakingthe
exchange symmetry of domain classes, which is unrealistic, as specific protein domains perform different
biologicalfunctions. Forthesevariants,wewillpresentmean-fieldandsimulationsresults,andcharacterize
their phenomenology in relation with empirical data. In particular, while in general for these variants the
rigorousmathematicalresultsexistingfortheCRPbreakdown,wewillshowhowtheuseofsimplemean-
fieldmethodsprovestobearobusttoolforaccessingthequalitativephenomenology.
4
II. GENERICFEATURESOFTHEMODEL
Themodelrepresentsaproteomethroughitsrepertoireofdomains. Domainshavingthesametopology
are collected in domain classes (Fig. 1). Thus the relevant data structures are partitions of elements (do-
mains) into classes. The basic observables considered are the following: n, the total number of domains,
f(n), a random variable indicating the number of classes (distinct domain topologies) at size n, a random
variable k , the population of class i, n , the size at birth of class i, and f(j,n): the number of domain
i i
classes having j members. We will generally indicate mean values by capitalized letters (e.g. F(n) is the
meanvalueoff(n),K themeanvalueofk ,etc.).
i i
The model is conceived as a stochastic process based on the elementary moves available to a genome
(Fig. 1) of adding and losing domains, associated to relative probabilities: p , the probability to duplicate
O
anolddomain(modelinggeneduplication),p ,theprobabilitytoaddanewdomainclasswithonemember
N
(whichdescribesdomaininnovation,forexamplebyhorizontaltransfer),andp ,alossprobability(which
L
we will initially disregard, and consider in a second step). Iteratively, either a domain is added or it is lost
withtheprescribedprobabilities.
An important feature of the duplication move is the (null) hypothesis that duplication of a domain has
uniformprobabilityalongthegenome,andthusitismoreprobabletopickadomainofalargerclass. This
is a common feature with previous models [13]. This hypothesis creates a “preferential attachment” [26]
principle,statingthefactthatduplicationismorelikelyinalargerdomainclass,which,inthismodelasin
previousones,isresponsiblefortheemergenceofpower-lawdistributions.
In mathematical terms, if the duplication probability is split as the sum of per-class probabilities pi ,
O
thishypothesisrequiresthatpi ∝ k , wherek isthepopulationofclassi, i.e.theprobabilityoffindinga
O i i
domainofaparticularclassandduplicatingitisproportionaltothenumberofmembersofthatclass.
Itisimportanttonoticethatinthismodel,whilencanbeusedasanarbitrarymeasureoftime,theratio
ofthetime-scalesofduplicationandinnovationisnotarbitrary,andissetbytheratiop /p . Inthemodel
N O
ofGersteinandcoworkers,thisistakenasaconstant,astheinnovationmoveconsideredtobestatistically
independentfromthegenomecontent. Inparticular, bothprobabilitiesareconsideredtobeconstant. This
choice has two problems. First, it cannot give the observed sublinear scaling of F(n). Indeed, if the
probability of adding a new domain is constant with n, so will be the rate of addition, implying that this
quantity will increase on average linearly with genome size. Moreover, the same model gives power-law
distributionsfortheclasseswithexponentlargerthantwo,incontrastwithmostoftheavailabledata.
Previous investigators did not consider the fact that genomes cluster around a common curve [23], and
thought of each of them as coming from an independent stochastic process with different parameters [13].
5
Coarse-grained view of aproteome:
n domains
Fdomain
classes
3 2 3 2 1 2
i.Duplication of anexisting domain
ii.Innovation, genesis ortransfer of a domain
Possible moves:
iii.Loss of anexisting domain
FIG.1: (Coloronline)Schemeofthegenericfeaturesoftheduplication-innovation-lossmodel. Top: Thesequences
of proteins coded by a genome can be broken down into domains, represented by colored boxes. All the boxes of
thesamecolor,representingdomainswiththesametopologyarecollectedinthesamedomainclass. Bottom: basic
moves. Elementsofaclasscanbeduplicatedorlost,andnewclassescanbeformedwithprescribedprobabilities.
Furthermore, choosing constant p implies that for larger genomes the influx of new domain families is
N
heavilydominantonthefluxofduplicateddomains.
As noted by Durrett and Schweinsberg [19], constant p makes sense only if one thinks that new fold
N
topologiesemergefromaninternal“nucleation-like”processwithconstantrate,ratherthanfromanexternal
flux. Thisprocesscouldbepicturedasthegenesisofnewtopologiesfromsequencemutation. Empirically,
whilegenesiseventsarereported[27]andmustoccur,itisclearthatdomaintopologiesareverystable,and
theexplorationofsequencespaceisnotfree, butconditionedbyanumberofadditionalimportantfactors,
includingchromosomalpositionandexpressionpatternsofgenes,andtheirroleinbiologicalnetworks[28].
Moreover, in prokaryotes, it is known that a large contribution to the innovation of coding genomes is
provided by horizontal gene transfer [27], the exchange of genetic material between species, which can
be reasonably represented in a model as an external flux, as opposed to the internal nucleation process
representinggenesis.
Foreukaryotes,horizontalgenetransferislessimportant,andtherecanbemultiplerelevantinnovation
6
processes including exonization, loss of exons, alternative start sites changing the protein. We have not
attemptedtomodelthedetailedprocessesleadingtoinnovation,andbecauseoftheirhighercomplexityin
eukaryotes,weprefertocomparethemodeltotheprokaryotedatasetalone. However,wecanpointoutthat
in principle the same model has good agreement with the set of prokaryotes and eukaryotes together [23].
Ineukaryotes,thechangeinnumberofclasseswithrespecttosizechangeisgenerallysmall,butseemsto
haveatrendthat”glues”quitewellwithprokaryotes(andinparticular,innovationdecreaseswithsize).
Motivated by the sublinear scaling of the number of domain classes, and taking into account in an
effective way the role of processes that condition the addition of new domain topologies, we consider
statisticallydependentmoves. Ongeneralgrounds,ifagenomeisacomplexsystemwheresub-components
interactinclustersandnon-locally,domaintopologiesaswellhavetobecoordinatedwithotherpartsofthe
system,sothatitisreasonablethatevolutionarymovesareconditionedbywhatisalreadypresent,andthat
theactualnumberofdomaintopologiesneednottobetriviallyanextensivequantity.
Thesimplestwaytoimplementthischoiceistoconcentrateontheinnovationprocess. Letusconsider
the indicator ξ(n), taking value 1 if a new domain class is born at size n, and value 0 otherwise. The
number of classes at size n will be f(n) = (cid:80) ξ(j). If the random variables ξ(n) are independent and
j≤n
identicallydistributed,i.e.p(ξ(n) = 1) = p =constant,f(n)followsaBernoullidistributionwhosemean
N
valueislinearinn. Moreover, f(n) willbeincreasinglyconcentratedwithincreasingnonthedeterministic
n
valuep . Ifviceversatherandomvariablesξ(n)arestatisticallydependentoralsosimplynotidentically
N
distributed, the mean value may not be linear in n, and the concentration phenomenon may not occur.
Both features, dependence and lack of concentration, are important. The former is necessary to obtain
the observed sublinear behavior, the latter might create an intrinsic “diversity” in the genome ensemble,
independentlyonthefinitesizeofobservedgenomes(however,thecurrentlyavailabledataareinsufficient
toestablishthisempirically).
III. SIMPLESTFORMULATION:CHINESERESTAURANTPROCESS
We investigate this process using analytical asymptotic equations and simulations. We start by consid-
ering only growth moves, by duplication and innovation, postponing the inclusion of domain loss in the
model. We will see that the resulting model contains the basic qualitative phenomenology of the scaling
laws and can thus be regarded as the paradigmatic case. One can arrive at the defining equations with dif-
ferentarguments. Asimplewayistoassumethatdomainduplicationisarareevent,describedbyaPoisson
distribution with characteristic time τ, during which there is a flux of external or new domain topologies
θ. Then p = θ . In this case the variables ξ(n) are independent but not identically distributed. It is
τ N n+θ
7
immediate to verify that f(n) has mean value given by (cid:80)n θ and thus grows as ∼ θlogn. The same
j=1 θ+j
resultcanbeobtainedbythinkingofdomainadditionasadependentmove,conditionedonn,f(n)orboth.
It is possible to consider different intermediate scenarios where the pool of old domain classes is in
competition with the universe explorable by the new classes. The simplest scheme, which turns out to be
quitegeneral,canbeobtainedbychoosingtheconditionalprobabilitythatanewclassisborn(ξ(n+1) = 1)
giventhefactthatf(n) = f atsizen
θ+αf
p = , (1)
N
n+θ
hence
n−αf
p = , (2)
O
n+θ
whereθ ≥ 0andα ∈ [0,1].
Consideringper-classduplicationprobability,onecanchoosethefollowingexpression,thatasymptoti-
callyestablishesthepreferentialattachmentprinciple:
k −α
pi = i . (3)
O n+θ
Here,θrepresentsacharacteristicnumberofdomainclassesneededforthepreferentialattachmentprinciple
tosetin,anddefinesthebehavioroff(n)forsmalln(n → 0). αisthemostimportantparameter,which
sets the scaling of the duplication/innovation ratio. Intuitively, the smaller α, the more the growth of f is
depressedwithgrowingn,andsincep isasymptoticallyproportionaltotheclassdensityf/nitisharder
N
to add a new domain class in a larger, or more heavily populated genome. As we will see, this implies
p /p → 0asn → ∞,correspondingtoanincreasinglysubdominantinfluxofnewfoldclassesatlarger
N O
sizes. This choice reproduces the sublinear behavior for the number of classes and the other scaling laws
describedinproperties(i-iii).
ThiskindofmodelhaspreviouslybeenexploredinstatisticsunderthenameofPitman-Yor,orChinese
RestaurantProcess(CRP)[24,25,29,30],whereitisknownasoneoftheparadigmaticprocessesthatgen-
erate partitions of elements into classes that are symmetric by swapping, or “exchangeable”. This process
is used in Bayesian inference and clustering problems [31]. In the Chinese restaurant (with table sharing)
parallel, individual domains correspond to customers and tables are domain classes. A domain belonging
to a given class is a customer sitting at the corresponding table. In a duplication event, a new customer is
8
seatedatatablewithapreferentialattachmentprinciple,andinaninnovationevent,anewtableisadded.
Mean-fieldtheoryfortheCRP
A simple mean-field treatment of the CRP allows to access its scaling behavior. Rigorous results for
the probability distribution of the fold usage vector {k ,...,k }, for f(n) = F, are in good agreement
1 F
withmean-fieldpredictions. Itisimportanttonotethatforthisstochasticprocess,theusuallarge-deviation
theorems do not hold, so that large-n limit values of quantities such as f(n) do not converge to numbers,
F(n)
but rather to random variables [24]. Despite of this non-self-average property, it is possible to understand
the scaling of the averages K and F (of k and f respectively) at large n, writing simple “mean-field”
i i
equations,forcontinuousn. NotethatrigorouslythemeanvalueK (n)isstillarandomvariable,function
i
ofthe(stochastic)birthtimen ofclassi.
i
Fromthedefinitionofthemodel,weobtain
K (n)−α
i
∂ K (n) = (4)
n i
n+θ
αF(n)+θ
∂ F(n) = . (5)
n
n+θ
These equations have to be solved with initial conditions K (n ) = 1, and F(1) = 1. Hence, for α (cid:54)= 0,
i i
onehas
n+θ
K (n) = (1−α) +θ , (6)
i
n +θ
i
and
1 (cid:20) (cid:18)n+θ(cid:19)α (cid:21)
F(n) = (α+θ) −θ ∼ nα , (7)
α θ
while,forα = 0,
F(n) = θlog(n+θ) ∼ log(n) . (8)
TheseresultsimplythattheexpectedasymptoticscalingofF(n)issublinear,inagreementwithobservation
(i).
Themean-fieldsolutioncanbeusedtocomputetheasymptoticofP(j,n) = F(j,n)/F(n), following
9
the same line of reasoning used by Barabasi and Albert for the preferential attachment model [26]. This
works as follows. From the solution, j > K (n) implies n > n∗, with n∗ = (1−α)n−θ(j−1), so that the
i i j−α
cumulativedistributioncanbeestimatedbytheratioofthe(average)numberofdomainclassesbornbefore
size n∗ and the number of classes born before size n, P(K (n) > j) = F(n∗)/F(n). P(j,n) can be
i
obtainedbyderivationofthisfunction. Forn,j → ∞,andj/nsmall,wefind
P(j,n) ∼ j−(1+α) , (9)
forα (cid:54)= 0,and
θ
P(j,n) ∼ , (10)
j
forα = 0. Theaboveformulasindicatethattheaverageasymptoticbehaviorofthedistributionofdomain
class populations is a power law with exponent between 1 and 2, in agreement with observation (ii). In
contrast, thebehaviorofthemodelofGersteinandcoworkers[13,18,19]canbefoundinthisframework
by taking improperly α = 1, that is for constant p , p . It gives a linearly increasing F(n) and a power-
N O
law distribution with asymptotic exponent 2 for the domain classes. Note that the phenomenology of the
Barabasi-Albertpreferentialattachmentscheme[26]isreproducedbyaCRP-likemodelwhereateachstep
anewdomainclass(correspondingtothenewnetworknode)withonaveragemmembers(theedgesofthe
node)isintroduced,andatthesametimemdomainsareduplicated(theedgesconnectingoldnodestothe
newlyintroducedone).
It is possible to obtain the same results through a different route compared to the above reasoning, by
writingthehierarchyofmean-fieldequationsforF(j,n),usingamasterequation-likeapproach. Similarly
aswhathappensfortheZero-Rangeprocess[32],theseequationscontainsourceandsinktermsgoverning
the population dynamics of classes. Duplications create a flux from classes with j − 1 to classes with j
members,whileonlyF(1,n)hasasourcetermcomingfromtheinnovationmove:
∂∂nnFF(1(n,n))== αFn(+n)θ+θα−Fn((+n1)θ−+θα)Fn(+1,θn) . (11)
∂nF(·2·,·n) = (1−α)Fn(+1,θn)·−··(2−α)Fn(+2,θn)
We consider the limit of large n, and use the ansatz F(j,n) = χ F(n). This ansatz can be justified
j
10
300 j = 1 400 j = 1
j = 2 j = 2
j = 3 j = 3
j = 4 300 j = 4
n)200 n)
(j, (j,
F F200
100
100
0 0
0 200 400 600 800 0 500 1000
F(n) F(n)
(a) (b)
400 j = 1 400 j = 1
j = 2 j = 2
j = 3 j = 3
300 j = 4 300 j = 4
) )
n n
(j, (j,
F200 F200
100 100
0 0
0 500 1000 0 500 1000
F(n) F(n)
(c) (d)
FIG. 2: (Color online) The classes with j members form a finite constant fraction of the total number of classes.
F(j,n)(thenumberofclasseswithj members)isplottedversusF(n)fordifferentvaluesofj inempiricaldata(a)
andsimulationsoftheCRP(α = 0.32,θ = 50)without(b)andwithdomainloss(δ = 0.2)inmodelwithuniform
loss(c)andweightedloss(d). Intheplotsfromsimulations,thelinesinlightercolorsrepresent10differentrealiza-
tions,whilethepointswitherrorbarsarethemeanvalueoverthem. Thefigureshowsthataddingdomainlosstothe
modelhasnoqualitativeconsequencesinthebehaviorofthedomainclasshistograms.
empirically and by simulations, as shown in Fig. 2, which compares this feature in empirical data and in
simulationsofdifferentvariantsofthemodel.
Wehave
∂nαFχ(1n=) = α−α(1Fn(−n)α)χ1
αχ = (1−α)χ −(2−α)χ , (12)
2 1 2
αχ··j·= (j −1−α)χj−1−(j −α)χj
