Modelling the evolution ofprotein coding sequences
sampled from measurably
evolving populations
Matthew Goode, Stephane Guindon andAllen Rodrigo
The Bioinformatics Institute (New Zealand),
Department of Statistics, University of Auckland
andAllan Wilson Centre for Molecular Ecology and Evolution

Outline
 Measurably Evolving Populations
 Modelling codon evolution
 Conclusions

HIV env Evolution

Pt.9

Shankarappa et al. 1999

Pt.2
HIV1U35926

Pt.7

Patient #6 fromWolinsky et al.

HIVU95460
HIV1U36148

Pt.6

HIV1U36073
HIV1U36015
HIV1U35980

Pt.8

Pt.3

Pt.1

HIV env Evolution
Shankarappa et al. 1999Viral Divergence - SC

10%

8%6%
4%2%02468

Years Post Seroconversion

1010%Pt.5

Measurably evolving populations (MEPs)
MEP: Any population evolving fastenough so that a statistically significant
accumulation of substitutions between
serially sampled sequences can bedetected.

Present timepoint(n = 5)

Earlier timepoint(n = 5)

Rapidly evolving pathogens, e.g.,.HIV, FIV, Influenza.
Ancient DNA: so far mostly
mitochondrial, e.g.
• Adelie penguins• Pleistocene bears

Developments in the Analysis of MEPs
1: Uniform Rate (TipDate: Rambaut,2000)
 Estim ates uniform rate2t1

R e c o n s tru c tin gs e ria l g e n e a lo g ie su s in g s U P G M A(D ru m m o n d &R o d rig o , 2 0 0 0 )(ω) over entire sam pling
period.

τ3
ω(τ 3 − τ 1 )h1

 Strict m olecular clock.

τ2h2

τ1

 U se M L to optim izeh4branch lengths, estim ate
param eters h, ω.

h3h5h6

 M axim ise

L(h, ω) = P(D | T, h, ω);Multiple Rates with Dated Tips (M RDT)
(Drum m ond et al. 2001)
τ3
τ∗
τ2

ω 2 (τ 3 − τ ∗ )

rate
ω1

ω 1 (τ ∗− τ 1 )τ1

The s-coalescent in a likelihood framework(Rodrigo & Felsenstein, 1998)

ω2

τ1
M ultiple rates: sam pling tim es
not coincident with rate changes.

τ2 τ∗ τ3tim e

 The likelihood function:P( D | N τ , ω ) = ∑ P( D | G , ω )P (G | N τ )Gwhereω = µ /τn2P(G | N τ ) ≈∏p= n 2 − c + 1

 p( p − 1) 2exp  −ν p


2 Nτ
2 Nτ

 (n − c )(n 2 − c − 1)
× exp  − 2(v c − ν 1 ) 
2 Nτ×n 1 +n 2 − c2q( q − 1)

∏ 2 N τ exp  − 2 N τ
q= 2
νq

Estimating evolutionary rates
SingleRate with Date Tips (SRDT)
2t1
Estimates uniform rate (ω)
over entire sampling period.

τ3
ω(τ3− τ1)

h1

τ2h2

τ1h4

Strict molecular clock.
Use ML to optimize branch
lengths, estimate parametersh, ω.

h3h5
h6

MaximiseL(h, ω) = P(D | T, h, ω);Measurably evolving populations (MEPs)

Modelling Codon Evolution
The ratio of the rate of nonsynonymoussubstitutions (dN) to the rate of synonymous
substitutions (dS)This ratio is symbolised by ωω = dN /dS

Modelling Codon Evolution
Codons evolve under different selective regimes
Positive, diversifying selection
dN > dS

ω>1

Negative, purifying selectiondN < dS

ω<1Neutrality
dN = dS,

ω=1

Codon Evolution
In Nielsen & Yang Codon Model M2, a particularsite is assumed to evolve under one, and only
one, of the three classes
With probabilities p0, p1, p2, for ω=0, ω=1, ω>1
respectively.
Across the tree a site never changes selection
classNielsen & Yang's (1998) model

ω1=0.0ω2=1.0ω3>1.0Time (within Subject)
P9

P2

P8

P3

P1P7

P5

P6Selection
10

Positionalong
sequence

20

30

HIV-1 env gene C3 - V5
region codons 1-118

40

Ross & Rodrigo, 2002

5060

P r o p o r tio n o f S ite sw ith ω > 170< 0 .0 1
800 .0 1 - 0 .0 30 .0 3 - 0 .0 7900 .0 7 - 0 .1 5100

> 0 .1 5
No sequence110

A new model of codon evolution
for serially sampled sequences

A site is allowed to be in different
selection classes before and after thesplit.
Instantaneous rate matrix andtransition probabilities change acrosssplit, but still easy to calculate
likelihood.

TimeSubstitution model changes at an apriori specified timepoint.

Single split

Before Split

After Split
Negative

NeutralPositive

Negative

p1

p2

p3

Neutral

p4

p5

p6

Positive

p7

p8

p9

Before Split

After Split
Negative

NeutralPositiveNegative

p02

p0p1

p0p2

Neutral

p1p0

p12

p1p2

Positive

p2p0

p2p1

p22

Meaneffectp(wb=i) = p(wa=i)p(wb=i, wa=j) = p(wa=i)p(wb=j)

Estimation &HierarchicalLikelihood Ratio
Tests

Before Split (p)After Split (q)Negative

NeutralPositive

Negative

p0q0

p0q1

p0q2

Neutral

p1q0

p1q1

p1q2

Positive

p2q0

p2q1

p2q2

Meaneffectp(wb=i) = p(wa=i)p(wb=i, wa=j) = p(wa=i)p(wb=j)

Maineffectsp(wb=i) ≠ p(wa=i)p(wb=i,wa=j) = p(wa=i)p(wb=j)

Estimation &HierarchicalLikelihood Ratio
Tests

Before Split

After Split
Negative

NeutralPositive

Negative

p1

p2

p3

Neutral

p4

p5

p6

Positive

p7

p8

p9

Meaneffectp(wb=i) = p(wa=i)p(wb=i, wa=j) = p(wa=i)p(wb=j)

Main
effectsp(wb=i) ≠ p(wa=i)p(wb=i,wa=j) = p(wa=i)p(wb=j)Interaction p(wb=i) ≠ p(wa=i)effectsp(wb=i,wa=j) ≠ p(wa=i)p(wb=j)

Estimation &HierarchicalLikelihood RatioTestsBefore Split

After Split
Negative

NeutralPositive

Negative

p1

p2

p3

Neutral

p4

p5

p6

Positive

p7

p8

p9

Meaneffectp(wb=i) = p(wa=i)p(wb=i, wa=j) = p(wa=i)p(wb=j)

Main
effectsp(wb=i) ≠ p(wa=i)p(wb=i,wa=j) = p(wa=i)p(wb=j)

Estimation &HierarchicalLikelihood Ratio
Tests

Not Nielsen & Yang!

Interaction p(wb=i) ≠ p(wa=i)effectsp(wb=i,wa=j) ≠ p(wa=i)p(wb=j)

Before Split

After Split
Negative

NeutralPositive

Negativep000

Neutral0

p10

Positive

00p2

Meaneffect
Maineffectsp(wb=i) = p(wa=i)p(wb=i, wa=j) = p(wa=i)p(wb=j)p(wb=i) ≠ p(wa=i)p(wb=i,wa=j) = p(wa=i)p(wb=j)Interaction p(wb=i) ≠ p(wa=i)effectsp(wb=i,wa=j) ≠ p(wa=i)p(wb=j)

Estimation &HierarchicalLikelihood Ratio
Tests

Not Nielsen & Yang!
p(wb=i, wa=j) = p(w=i), i=j0 otherwiseBefore Split

After Split
Negative

NeutralPositive

Negative

p1p2p3Neutral

p4

p5

p6

Positive

p7

p8

p9

Meaneffectp(wb=i) = p(wa=i)p(wb=i, wa=j) = p(wa=i)p(wb=j)

Maineffectsp(wb=i) ≠ p(wa=i)p(wb=i,wa=j) = p(wa=i)p(wb=j)

Estimation &HierarchicalLikelihood Ratio
Tests

Not Nielsen & Yang!
p(wb=i, wa=j) = p(w=i), i=j
0 otherwise

Interaction p(wb=i) ≠ p(wa=i)effectsp(wb=i,wa=j) ≠ p(wa=i)p(wb=j)

An HIV-1 env example• HIV-1 envelope (env) sequences (60 sequences of 660 bases)from infected patient.

• sampled at days 0, 214, 671, 699 and 1005.
• Monotherapy with zidovudine was initiated after day 409.

Likelihood Ratio Test
Estimate parameters (p’s, ω’s, and κ’s) using
maximum likelihood for NY‐M2 and fully
saturated model.
8 degrees of freedom difference between models

ML EstimatesNielsen‐Yang M2
Log Likelihood ‐2873.4

Saturated Model
Log Likelihood ‐2855.8

Parameter estimatesωafter=0

ωafter=1

ωafter=∞ωbefore=0

0.425

0.065

0.000

Marginal p(before)
0.490

ωbefore=10.368

0.000

0.000

0.368

ωbefore=7.9

0.139

0.000

0.003

0.142

Marginal p(after)

0.932

0.065

0.003

1.000

Future work To date, most methods on changing evolutionary parameters

deal with lineage‐independent changes.
 This is suitable for species (and higher taxa) phylogenies.
 Forces that influence rates of evolution may act differently in different

lineages Not necessarily suitable for intraspecific phylogenies.
 External influences act on the population as a whole.
 Also true for some taxonomic phylogenies

Progression of HIV Infection
CD4+T cells

Viralload

Acuteinfection

Asymptomatic Phase

AIDSTime (t)Changing models as a function of time

Model ofevolution, Q(t) or
mutation rate, µ(t)Changing models of evolution as afunction of time (commutable models)Commutable models of evolution

Q(t )× Q(t' ) = Q(t' )× Q(t )If Q changes as a function of time, we can calculate the transition
probabilities as:

PN (T ) = e ∫Q ( t ) dt

Rodrigo et al. (2008)
Phil Trans Roy Soc B

Conclusions
We have developed a codon model of evolutionthat permits:
Changes to the ratio of non‐synonymous to synonymous
substitution rates over time.
Different proportions of sites in each selective class.

The model is based on a simultaneous change inrate across all lineages.
Consequently, it is better for intraspecific phylogenies thaninterspecific phylogenies.

Acknowledgements
David Bryant
Alexei Drummond
Joseph HeledHoward Ross