-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathflux.tex
834 lines (720 loc) · 80 KB
/
flux.tex
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
\documentclass[11pt,oneside,letterpaper]{article}
% graphicx package, useful for including eps and pdf graphics
\usepackage{graphicx}
\DeclareGraphicsExtensions{.pdf,.png,.jpg}
% basic packages
\usepackage{color}
\usepackage{parskip}
\usepackage{float}
% text layout
\usepackage{geometry}
\geometry{textwidth=15cm} % 15.25cm for single-space, 16.25cm for double-space
\geometry{textheight=22cm} % 22cm for single-space, 22.5cm for double-space
% helps to keep figures from being orphaned on a page by themselves
\renewcommand{\topfraction}{0.85}
\renewcommand{\textfraction}{0.1}
% bold the 'Figure #' in the caption and separate it with a period
% Captions will be left justified
\usepackage[labelfont=bf,labelsep=period,font=small]{caption}
% review layout with double-spacing
%\usepackage{setspace}
%\doublespacing
%\captionsetup{labelfont=bf,labelsep=period,font=doublespacing}
% cite package, to clean up citations in the main text. Do not remove.
\usepackage{cite}
%\renewcommand\citeleft{(}
%\renewcommand\citeright{)}
%\renewcommand\citeform[1]{\textsl{#1}}
% Remove brackets from numbering in list of References
\renewcommand\refname{\large References}
\makeatletter
\renewcommand{\@biblabel}[1]{\quad#1.}
\makeatother
\usepackage{authblk}
\renewcommand\Authands{ \& }
\renewcommand\Authfont{\normalsize \bf}
\renewcommand\Affilfont{\small \normalfont}
\makeatletter
\renewcommand\AB@affilsepx{, \protect\Affilfont}
\makeatother
% notation
\usepackage{amsmath}
\usepackage{amssymb}
\newcommand{\virus}{\mathbf{x}} % virus coordinate
\newcommand{\serum}{\mathbf{y}} % serum coordinate
\newcommand{\viruses}{\mathbf{X}} % set of virus coordinates
\newcommand{\sera}{\mathbf{Y}} % set of serum coordinates
\newcommand{\ve}{v} % virus avidity
\newcommand{\se}{s} % serum potency
\newcommand{\ves}{\mathbf{v}} % set of virus avidities
\newcommand{\ses}{\mathbf{s}} % set of serum potencies
\newcommand{\point}{f_{\scriptscriptstyle \vert}} % point likelihood
\newcommand{\threshold}{f_{\textstyle \lrcorner}} % threshold likelihood
\newcommand{\interval}{f_{\sqcup}} % interval likelihood
\newcommand{\mdssd}{\varphi} % MDS standard deviation
\newcommand{\virussd}{\sigma_x} % virus / diffusion standard deviation
\newcommand{\serumsd}{\sigma_y} % serum standard deviation
\newcommand{\drift}{\mu} % drift / advection
\newcommand{\tree}{\tau} % phylogeny
\newcommand{\vn}{n} % number of viruses
\newcommand{\sn}{k} % number of sera
\newcommand{\normal}{\mathcal{N}} % normal distribution
\newcommand{\bwithin}{\beta_w} % within clade drift coefficient
\newcommand{\bsister}{\beta_s} % sister clade drift coefficient
\newcommand{\bother}{\beta_t} % across clade drift coefficient
\newcommand{\incclade}[1]{y_\mathrm{#1}}
\newcommand{\driftclade}[1]{x_\mathrm{#1}}
\setlength{\arraycolsep}{2pt}
\newcommand{\smalltwomatrix}[2]{\scriptsize \Big( \begin{matrix} #1 \\ #2 \end{matrix} \Big)} % pretty inline matrix
\newcommand{\smallfourmatrix}[4]{\scriptsize \Big( \begin{matrix} #1 & #2 \\ #3 & #4 \end{matrix} \Big)} % pretty inline matrix
\newcommand{\twomatrix}[2]{\left( \begin{matrix} #1 \\ #2 \end{matrix} \right)} % pretty inline matrix
\newcommand{\fourmatrix}[4]{\left( \begin{matrix} #1 & #2 \\ #3 & #4 \end{matrix} \right)} % pretty inline matrix
%%% TITLE %%%
\title{\vspace{1.0cm} \Large \bf
Integrating influenza antigenic dynamics with molecular evolution
}
\author[1]{Trevor Bedford}
\author[2,3,4]{Marc A. Suchard}
\author[5]{Philippe Lemey}
\author[1]{Gytis Dudas}
\author[6]{Victoria Gregory}
\author[6]{Alan J. Hay}
\author[6]{John W. McCauley}
\author[7,8]{Colin A. Russell}
\author[7,8,9]{Derek J. Smith}
\author[1,10,11]{Andrew Rambaut}
\affil[1]{Institute of Evolutionary Biology, University of Edinburgh, Edinburgh, UK}
\affil[2]{Department of Biomathematics, David Geffen School of Medicine at UCLA, University of California, Los Angeles CA, USA}
\affil[3]{Department of Human Genetics, David Geffen School of Medicine at UCLA, University of California, Los Angeles CA, USA}
\affil[4]{Department of Biostatistics, UCLA Fielding School of Public Health, University of California, Los Angeles CA, USA}
\affil[5]{Department of Microbiology and Immunology, Katholieke Universiteit Leuven, Leuven, Belgium}
\affil[6]{Division of Virology, MRC National Institute for Medical Research, Mill Hill, London, UK}
\affil[7]{Centre for Pathogen Evolution, Department of Zoology, University of Cambridge, Cambridge, UK}
\affil[8]{WHO Collaborating Center for Modeling, Evolution, and Control of Emerging Infectious Diseases, University of Cambridge, Cambridge, UK}
\affil[9]{Department of Virology, Erasmus Medical Centre, Rotterdam, Netherlands}
\affil[10]{Fogarty International Center, National Institutes of Health, Bethesda, MD, USA}
\affil[11]{Centre for Immunity, Infection and Evolution, University of Edinburgh, Edinburgh, UK}
\date{}
\begin{document}
\maketitle
%%% ABSTRACT %%%
\begin{abstract}
Influenza viruses undergo continual antigenic evolution allowing mutant viruses to evade host immunity acquired to previous virus strains.
Antigenic phenotype is often assessed through pairwise measurement of cross-reactivity between influenza strains using the hemagglutination inhibition (HI) assay.
Here, we extend previous approaches to antigenic cartography, and simultaneously characterize antigenic and genetic evolution by modeling the diffusion of antigenic phenotype over a shared virus phylogeny.
Using HI data from influenza lineages A/H3N2, A/H1N1, B/Victoria and B/Yamagata, we determine patterns of antigenic drift across viral lineages, showing that A/H3N2 evolves faster and in a more punctuated fashion than other influenza lineages.
We also show that year-to-year antigenic drift appears to drive incidence patterns within each influenza lineage.
This work makes possible substantial future advances in investigating the dynamics of influenza and other antigenically-variable pathogens by providing a model that intimately combines molecular and antigenic evolution.
\end{abstract}
%%% IMPACT %%%
% Combined evolutionary and antigenic analysis shows that human influenza viruses differ dramatically in rates of antigenic drift and these rates significantly impact seasonal incidence patterns.
\pagebreak
%%% INTRODUCTION %%%
\section*{Introduction}
Seasonal influenza infects between 10\% and 20\% of the human population every year, causing an estimated 250,000 to 500,000 deaths annually \cite{flufactsheet}.
While individuals develop long-lasting immunity to particular influenza strains after infection, antigenic mutations to the influenza virus genome result in proteins that are recognized to a lesser degree by the human immune system, leaving individuals susceptible to future infection.
The influenza virus population continually evolves in antigenic phenotype in a process known as antigenic drift.
A large proportion of the disease burden of influenza stems from antigenic drift, which allows individuals to be infected multiple times throughout their lives.
Although influenza vaccines may lack efficacy for a variety of reasons \cite{Osterholm12}, antigenic drift causes efficacy of a fixed vaccine formulation against circulating viruses to decline over time.
A thorough understanding of the process of antigenic drift is essential to public health efforts to control mortality and morbidity through the use of a seasonal influenza vaccine.
Before 2009, there were four major lineages of influenza circulating within the human population: the H3N2 and H1N1 subtypes of influenza A, and the Victoria and Yamagata lineages of influenza B.
In the case of influenza A, subtypes A/H3N2 and A/H1N1 refer to the genes, hemagglutinin (H or HA) and neuraminidase (N or NA), that are primarily responsible for the antigenic character of a strain.
In the case of influenza B, Victoria (B/Vic) and Yamagata (B/Yam) refer to antigenically distinct lineages which diverged from a single lineage prior to 1980 \cite{Rota90}.
Mutations to the HA1 region of the hemagglutinin protein are thought to drive the majority of antigenic drift in the influenza virus \cite{Wiley81, Nelson07NatRevGenet}.
Experimental characterization of antigenic phenotype is possible through the hemagglutination inhibition (HI) assay \cite{Hirst43}, which measures the cross-reactivity of one virus strain to serum raised against another strain through challenge or vaccination.
Sera raised against older strains react poorly to more recent viruses resulting in new strains having a selective advantage over previously established strains.
The results of many HI assays across a multitude of viruses of a single subtype can be combined to yield a two-dimensional map, quantifying antigenic similarity and distance \cite{Smith04}.
The antigenic map of influenza A/H3N2 has shown substantial evolution of the influenza virus population since its emergence in 1968.
Evolution of antigenic phenotype appears punctuated with episodes of more rapid innovation interspersed by periods of relative stasis, while genetic evolution appears more continuous \cite{Smith04}, suggesting that a relatively small number of genetic changes or combinations of genetic changes may drive changes in antigenic phenotype \cite{Koel13}.
The process of antigenic drift results in the rapid turnover of the virus population, so that despite mutation, genetic diversity among contemporaneous viruses remains low.
Such population turnover is supported by phylogenetic analysis that shows a characteristically `spindly' tree with a single predominant trunk lineage and transitory side branches that persist for only 1--5 years \cite{Fitch97}.
Previously, the antigenic and genetic patterns of influenza evolution have been analyzed essentially in isolation.
An antigenic map is constructed from a panel of HI measurements, and a phylogenetic tree is constructed from sequence data.
However, the opportunity for a combined approach exists as both the antigenic map and the phylogenetic tree often contain many of the same isolates.
Here, we implement a flexible Bayesian approach to jointly characterize the antigenic and genetic evolution of the influenza virus population.
We apply this approach to investigate the dynamics of A/H3N2, A/H1N1, B/Vic and B/Yam viruses, and, for the first time, present detailed reconstructions of the antigenic dynamics of all four circulating influenza lineages.
%%% RESULTS %%%
\section*{Results and discussion}
\subsection*{Antigenic and evolutionary cartography}
In order to assess patterns of antigenic evolution among influenza strains, we implemented a Bayesian probabilistic analog of multidimensional scaling (MDS), referred to here as BMDS (see Materials and methods).
In this model, viruses and sera are given $N$-dimensional locations, thus specifying an `antigenic map', such that distances between viruses and sera in this space are inversely proportional to cross-reactivity.
In the BMDS model, a map distance of one antigenic unit translates to an expectation of a 2-fold drop in HI titer between virus and sera.
Maps that produce pairwise distances most congruent with the observed titers will have a high likelihood and will be favored by the BMDS model.
We integrate over sources of uncertainty, such as antigenic locations, in a flexible Bayesian fashion.
We apply this model to HI measurements of virus isolates against post-infection ferret antisera for influenza A/H3N2, A/H1N1, B/Vic and B/Yam.
We begin with Bayesian analogs of the models used by Smith et al.\ \cite{Smith04}, in which viruses and sera are represented as $N$-dimensional locations as described in the `Antigenic cartography' section of Materials and methods.
In this case, `serum potencies' are fixed to the maximum titers exhibited by particular ferret sera and give the baseline expectation for titer when virus and serum are antigenically identical.
Potency differs between serum isolates due to experimental noise (e.g.\ variation in serum concentration), but also due to differential ferret immune responses, causing some serum isolates to inhibit hemagglutination at generally higher titers than other isolates.
In this model, virus and serum locations follow an uninformative diffuse normal prior.
After comparing models of differing dimensions, Smith et al.\ \cite{Smith04} arrive at a 2D model as the preferred model for their data.
Smith et al.\ \cite{Smith04} implement a form of MDS, seeking to optimize virus and serum locations such that the sum of squared errors between expected and observed titers is minimized (eq.\ \ref{mds} of Materials and methods).
Here, in implementing BMDS, we provide a likelihood function for the probability of observing HI data given virus and serum locations (eq.\ \ref{bmds} of Materials and methods) and seek to estimate model parameters through Bayesian inference using Markov chain Monte Carlo (MCMC).
However, the basic antigenic model describing drop in HI titer as proportional to Euclidean distance between virus and serum locations is identical between these methods.
We test model performance by constructing training datasets representing 90\% of the HI measurements for each of the four influenza lineages and test datasets representing the remaining 10\% of the measurements for each lineage.
By fitting the BMDS model to the training dataset, we are able to predict HI titers in the test dataset and compare these predicted titers to observed titers.
We find that a two dimensional model has better predictive power than models of lower or higher dimension in all four influenza lineages (models 1--5; Table \ref{errortable}).
We find that this 2D model performs well, yielding an average absolute predictive error of between 0.78 and 0.91 log$_2$ HI titers across influenza lineages (model 2; Table \ref{errortable}), in line with the results of Smith et al.\ \cite{Smith04}.
Consequently, we specify a two dimensional model in all subsequent analyses.
The finding of a low-dimensional map across influenza lineages extends previous studies in A/H3N2 \cite{Smith04} and remains an interesting and fundamental empirical observation.
%%% errortable %%%
\begin{table}[h]
\centering
\caption{\textbf{Average absolute prediction error of log$_2$ HI titer for test data across models and datasets.}}
\label{errortable}
\makebox[\textwidth][c]{
\small
\begin{tabular}{ c c c c c c c c c c }
\hline
\multicolumn{6}{c}{} & \multicolumn{4}{c}{Test error} \\
Model & Data & Dimen & Location prior & Serum potency & Virus avidity & A/H3N2 & A/H1N1 & B/Vic & B/Yam \\
\hline
1 & HI & 1D & Uninformed & Fixed & None & 1.35 & 0.94 & 0.90 & 1.08 \\
2 & HI & 2D & Uninformed & Fixed & None & 0.91 & 0.78 & 0.82 & 0.90 \\
3 & HI & 3D & Uninformed & Fixed & None & 0.93 & 0.80 & 0.85 & 0.92 \\
4 & HI & 4D & Uninformed & Fixed & None & 0.98 & 0.84 & 0.90 & 0.97 \\
5 & HI & 5D & Uninformed & Fixed & None & 1.04 & 0.89 & 0.98 & 1.04 \\
6 & HI/year & 2D & Drift & Fixed & None & 0.91 & 0.75 & 0.77 & 0.83 \\
7 & HI/year/seq & 2D & Diffusion/Drift & Fixed & None & 0.89 & 0.74 & 0.74 & 0.83 \\
8 & HI/year/seq & 2D & Diffusion/Drift & Estimated & None & 0.77 & 0.73 & 0.66 & 0.75 \\
9 & HI/year/seq & 2D & Diffusion/Drift & Fixed & Estimated & 0.80 & 0.72 & 0.69 & 0.75 \\
10 & HI/year/seq & 2D & Diffusion/Drift & Estimated & Estimated & 0.76 & 0.71 & 0.64 & 0.72 \\
\hline
\end{tabular}
}
\end{table}
Previous work on influenza antigenic and genetic evolution has shown that antigenic distance accumulates with increasing genetic distance \cite{Hay01, Smith04, Russell08}.
Here, we examine pairwise relationships between viruses and observe a correlation between amino acid mutations and antigenic distance (Figure~\ref{seq_grid}) and a similar correlation between phylogenetic distance, measured in years, and antigenic distance (Figure~\ref{seq_grid}).
Thus, genetic relationships between viruses provide some predictive power to estimate antigenic distances in the absence of HI data.
However, the magnitudes of the coefficients of determination $R^2$ are low (Figure~\ref{seq_grid}), suggesting that genetic relationships alone will not completely resolve antigenic distances.
%%% map %%%
\begin{figure}[h]
\centering
\includegraphics[width=1.0\textwidth]{figures/seq_grid}
\caption{\textbf{Pairwise correlations between genetic distance, measured as amino acid mutations or as phylogenetic distance, and antigenic distance for influenza A/H3N2, A/H1N1, B/Vic and B/Yam.}
The top row shows correlations between number of amino acid mutations in HA1 and average antigenic distance between 10,000 random pairs of viruses.
The bottom row shows correlations between average phylogenetic distance, measured in terms of years, and average antigenic distance between 10,000 random pairs of viruses.
Dashed lines show linear model fits, with $R^2$ and slope noted, while solid lines show LOESS fits.
Antigenic distances derive from model 2 of Table \ref{errortable}.
}
\label{seq_grid}
\end{figure}
Consequently, we seek to flexibly incorporate genetic data by modeling antigenic phenotype as an evolutionary diffusion \cite{Lemey10}, wherein a virus's antigenic character state evolves along branches of the phylogenetic tree according to a Brownian motion process (see Materials and methods).
The phylogenetic diffusion process acts as a prior on virus locations, so that genetically similar viruses are expected to share similar antigenic locations.
The antigenic diffusion process includes both systematic drift with time and covariance induced by phylogenetic proximity.
We examine the effects of including only systematic drift (model 6; Table \ref{errortable}) and systematic drift plus phylogenetic diffusion (model 7; Table \ref{errortable}), finding a small increase in predictive accuracy of between 0.02 and 0.08 log$_2$ HI titers when both processes are included.
The systematic drift process informs virus and serum locations by dates of isolation and the phylogenetic diffusion process informs virus locations by genetic sequences.
Thus, in these models, antigenic locations are inferred using both genetic data and HI data and will differ from locations inferred from HI data alone.
If HI data is rich, then we expect minor differences in antigenic locations with the inclusion of genetic data (as may be the case for A/H3N2), while if HI data is spare, then we expect genetic data to play a larger role in determining antigenic locations (as may be the case for B/Vic and B/Yam).
We further extend the model by estimating serum potencies, rather than fixing serum potencies at maximum titers.
Serum potencies differ across isolates due to experimental variation in serum extraction and processing or due to variation in ferret immune response.
Serum potency determines the baseline expectation of titer when virus and serum have identical antigenic locations.
However, if serum potency is fixed to the serum's maximum titer, this will often not be the case, as the virus giving the maximum titer may be antigenically distinct from the serum.
Thus, fixing serum potencies will tend to under-estimate effect size; we observe a mean effect of 10.42 log$_2$ HI titers for A/H3N2 when fixing serum potencies and a mean of 10.94 when estimating serum potencies.
We find that estimating serum potencies improves test error further (model 8, Table \ref{errortable}), with improvements of between 0.01 and 0.12 log$_2$ HI titers.
Additionally, we include and estimate `virus avidity' in an analogous fashion, which is intended to represent differences in overall HI reactivity between viruses.
Experimental work has demonstrated that influenza variants exist that differ in HA binding activity for cell surface glycan receptors, and that these high-avidity variants may arise in the presence of antibody pressure \cite{Hensley09}.
The presence of differential virus avidity has been previously shown to distort antigenic maps constructed from a model that disregards avidity effects \cite{Li13}.
Here, with virus avidities estimated, baseline titer derives from both the virus and the serum used in the HI reaction.
We find that including virus avidities further improves test error, either with fixed serum potencies (model 9, Table \ref{errortable}) or with estimated serum potencies (model 10, Table \ref{errortable}).
With fixed serum potencies, the inclusion of virus avidities results in improvements of between 0.02 and 0.09 log$_2$ HI titers and with estimated serum potencies, the inclusion of virus avidities results in improvements of between 0.01 and 0.05 log$_2$ HI titers.
We find that the average absolute error in predicted log$_2$ HI titer is nearly constant with antigenic distance (Pearson correlation, $r = 0.098$), thus supporting our model assumption that the drop in log$_2$ titer is proportional to the Euclidean distance separating viruses and sera on the antigenic map.
Additionally, we find that the absolute error in predicted titer is nearly constant with time (Pearson correlation, $r = -0.007$).
Antigenic locations inferred by the model are well resolved; estimates of antigenic distance between pairs of viruses show relatively little variation across the posterior.
We estimate that virus distances have, on average, a 50\% credible interval of $\pm0.45$ antigenic units for A/H3N2, $\pm0.57$ units for A/H1N1, $\pm0.76$ units for B/Vic, and $\pm0.65$ units for B/Yam.
We find strong correspondence between our results and previous results by Smith et al.\ \cite{Smith04}, with equivalent models producing globally consistent antigenic maps and other models producing locally consistent maps with a small degree of global inconsistency (Figures~\ref{smith_comparison}--\ref{smith_comparison_effects}).
When implementing the same underlying model, differences in the MDS and BMDS approaches reflect greater philosophical differences between maximum-likelihood and Bayesian statistical approaches, with the former seeking the single most likely explanation for the data, and the latter seeking to fully characterize model uncertainty.
Additionally, the BMDS method improves flexibility, allowing extensions to the basic cartographic model, such as the incorporation of virus avidities and evolutionary priors, that improve fit and add biological interpretability.
\subsection*{Antigenic evolution across influenza lineages}
Through our analysis, we reveal the antigenic, as well as evolutionary, relationships among viruses in influenza A/H3N2, A/H1N1, B/Vic and B/Yam, quantifying both antigenic and evolutionary distances between strains (Figure~\ref{map}).
Over the time period of 1968 to 2011, influenza A/H3N2 shows substantially more antigenic evolution than is exhibited by A/H1N1 over the course of 1977 to 2009 or B/Vic and B/Yam over the course of 1986 to 2011.
We observe prominent antigenic clusters in A/H3N2 and A/H1N1, but less prominent, though still apparent, clustering in B/Vic and B/Yam.
Antigenic clusters show high genetic similarity, so that we observe very few mutation events leading to each cluster, rather than the repeated emergence of clusters.
This analysis makes the fate of antigenic clusters obvious, with two clusters in A/H3N2 (Victoria/75 and Beijing/89) appearing to be evolutionary dead-ends.
Labeling of prominent antigenic clusters in figures \ref{map} and \ref{drift} is intended as a rough guide for orientation and not as exhaustive catalog of antigenic variation.
%%% map %%%
\begin{figure}[h]
\centering
\includegraphics[width=1.0\textwidth]{figures/map}
\caption{\textbf{Antigenic locations of A/H3N2, A/H1N1, B/Vic and B/Yam viruses showing evolutionary relationships between virus samples.}
Circles represent a posterior sample of virus locations and have been shaded based on year of isolation.
Antigenic units represent two-fold dilutions of the HI assay.
Absolute positioning of lineages, e.g.\ A/H3N2 and A/H1N1, is arbitrary.
Lines represent mean posterior diffusion paths when virus locations are fixed.
Prominent antigenic clusters are labeled after vaccine strains present within clusters, and are abbreviated from Hong Kong/68, England/72, Victoria/75, Bangkok/79, Sichuan/87, Beijing/89, Beijing/92, Wuhan/95, Sydney/97, Fujian/02, California/04, Wisconsin/05, Brisbane/07, Perth/09 (A/H3N2), USSR/77, Singapore/86, Beijing/95, New Caledonia/99, Solomon Islands/06 (H1N1), Victoria/87, Hong Kong/01, Malaysia/04, Brisbane/08 (Vic), Yamagata/88, Shanghai/02, Florida/06, Wisconsin/10 (Yam).}
\label{map}
\end{figure}
HI assays lack sensitivity beyond a certain point, so that for A/H3N2, cross-reactive measurements only exist between strains sampled at most 14 years apart, leaving only threshold titers, e.g.\ `$<$40', in more temporally distant comparisons.
Because of the threshold of sensitivity of the HI assay, it's difficult to distinguish a linear trajectory in 2D antigenic space from a slightly curved trajectory (see Materials and Methods).
To solve this problem of identifiability, we assumed a weak prior that favors linear movement in the 2D antigenic space (present in models 6 through 9; Table \ref{errortable}), with the slope of the linear relationship and the precision of the relationship incorporated into the Bayesian model (see Materials and methods).
Because of this, we interpret map locations locally rather than globally, and assess rates of antigenic movement without making strong statements about the larger configuration under which the movement occurs.
We find that influenza A/H3N2 evolved along antigenic dimension 1 at an estimated rate of 1.01 antigenic units per year (Figure~\ref{drift}, Table~\ref{paramtable}).
However, we observe occasional large jumps in antigenic phenotype (Figure~\ref{drift}), corresponding to cluster transitions identified by Smith et al. \cite{Smith04}.
Most variation is contained within the first antigenic dimension, but dimension 2 occasionally shows variation when two antigenically distinct lineages emerge and transiently coexist (Figure~\ref{map}), as is the case with the previously identified Beijing/89 and Beijing/92 clusters.
%%% drift %%%
\begin{figure}[h]
\centering
\includegraphics[width=0.75\textwidth]{figures/drift}
\caption{\textbf{Antigenic drift of A/H3N2, A/H1N1, B/Vic and B/Yam viruses showing evolutionary relationships between virus samples.}
Antigenic drift is shown in terms of change of location in the first antigenic dimension through time.
Circles represent a posterior sample of virus locations and have been shaded based on year of isolation.
Antigenic units represent two-fold dilutions of the HI assay.
Relative positioning of lineages, e.g.\ A/H3N2 and A/H1N1, in the vertical axis is arbitrary.
Lines represent mean posterior diffusion paths when virus locations are fixed.
Prominent antigenic clusters are labeled after vaccine strains present within clusters, and are abbreviated from Hong Kong/68, England/72, Victoria/75, Bangkok/79, Sichuan/87, Beijing/89, Beijing/92, Wuhan/95, Sydney/97, Fujian/02, California/04, Wisconsin/05, Brisbane/07, Perth/09 (A/H3N2), USSR/77, Singapore/86, Beijing/95, New Caledonia/99, Solomon Islands/06 (H1N1), Victoria/87, Hong Kong/01, Malaysia/04, Brisbane/08 (Vic), Yamagata/88, Shanghai/02, Florida/06, Wisconsin/10 (Yam).}
\label{drift}
\end{figure}
We find that other lineages of influenza evolved in antigenic phenotype substantially slower than A/H3N2 (Figure~\ref{drift}, Table~\ref{paramtable}).
Influenza A/H1N1 evolved at a rate of 0.62 units per year, but showed a similar pattern of punctuated antigenic evolution with occasional larger jumps in phenotype, such as the emergence of the Solomon Islands/06 cluster.
Influenza B/Victoria and B/Yamagata evolved slower still, with mean estimated rates 0.42 units per year and 0.32 units per year, respectively.
Punctuated evolution is less obvious in B/Yam and B/Vic compared to A/H3N2 and A/H1N1, but antigenic clusters are still apparent, with recent transitions to the Brisbane/08 cluster in B/Vic \cite{Barr10} and to the Wisconsin/10 cluster in B/Yam \cite{Klimov12}.
Interestingly, a minor lineage of B/Vic, denoted B/Hubei-Songzi/51/2008 \cite{Barr10}, has persisted through 2011, while remaining antigenically distinct from B/Brisbane/60/2008 viruses (Figure~\ref{drift}).
Although we observe significantly different drift rates between lineages, we observe less variation in diffusion volatility (Table~\ref{paramtable}).
This is reflected in Figure \ref{drift}, where all four lineages exhibit similar levels of standing antigenic variation, despite A/H3N2 drifting more quickly in antigenic phenotype.
%%% paramtable %%%
\begin{table}[h]
\centering
\caption{\textbf{Estimates of drift rate $\drift$ (in units per year), diffusion volatility $\virussd^2$ (in units$^2$ per year) and scaled effective population size $N_e\tau$ (in years) for influenza A/H3N2, A/H1N1, B/Vic and B/Yam including posterior means and 95\% highest posterior density intervals.}}
\label{paramtable}
\begin{tabular}{ c c c c}
\hline
Lineage & Drift $\drift$ & Volatility $\virussd^2$ & Effective pop size $N_e\tau$ \\
\hline
A/H3N2 & 1.01 (0.98--1.04) & 1.25 (0.98--2.35) & 5.03 (4.42--5.73) \\
A/H1N1 & 0.62 (0.56--0.67) & 0.92 (0.65--1.56) & 6.38 (4.99--8.12) \\
B/Vic & 0.42 (0.32--0.51) & 1.22 (0.85--2.25) & 10.40 (8.42--12.80) \\
B/Yam & 0.32 (0.25--0.39) & 0.71 (0.46--1.36) & 9.48 (7.76--11.50) \\
\hline
\end{tabular}
\end{table}
These patterns of antigenic drift influence the corresponding virus phylogenies (Figure~\ref{trees}).
Influenza A/H3N2 has a characteristically spindly tree showing rapid turnover of the virus population, while A/H1N1 and B have trees that show greater degrees of viral coexistance (Figure~\ref{trees}).
The scaled effective population size $N_e \tau$ measures the timescale of coalescence of a phylogeny and quantifies the visual distinction between a `spindly' tree and a `bushy' tree \cite{Bedford11}.
In this case, $N_e$ represents the number of concurrent infections in a panmictic population with generation interval $\tau$ (time separating infections up the genealogical tree), so that $N_e \tau$ is measured in terms of years and gives the expected waiting time for two randomly chosen lineages to coalesce in the genealogical tree.
We see that $N_e \tau$ broadly correlates with the rate of antigenic drift (Table~\ref{paramtable}), with A/H3N2 showing fast drift and reduced effective population size as expected from basic epidemiological models \cite{Bedford12}.
Antigenic drift results in the replacement of antigenically primitive lineages by antigenically advanced lineages, thereby reducing genealogical diversity.
Thus, we observe a faster rate of antigenic drift in influenza A/H3N2 than in A/H1N1 or either lineage of influenza B.
Previous work using general epidemiological models has suggested that rates of antigenic drift may be influenced by both the fundamental reproductive number $R_0$ and the rate at which mutation decreases cross-immunity \cite{Gog02,Lin03}.
Correspondingly, models specific to influenza evolution ascribe differences in the rate of antigenic drift of A/H3N2 relative to A/H1N1 and influenza B to either greater $R_0$ or greater mutation rate \cite{Ferguson03,Bedford12}.
Without more detailed epidemiological modeling, the present study cannot conclusively distinguish between these causal possibilities.
%%% trees %%%
\begin{figure}[h]
\centering
\includegraphics[width=0.85\textwidth]{figures/trees}
\caption{\textbf{Time-resolved phylogenetic trees of A/H3N2, A/H1N1, B/Vic and B/Yam viruses.}
The maximium-clade credibility (MCC) tree is shown for each virus.
These trees show genealogical relationships, so that branches are measured in terms of years rather than substitutions.}
\label{trees}
\end{figure}
\subsection*{Punctuated evolution and its epidemiological consequences}
We sought to summarize year-to-year patterns of antigenic drift by calculating the difference in mean virus location between consecutive years (Figure~\ref{jumps}A).
We estimate year-to-year antigenic drift for years 1992 to 2011 by calculating the average location along dimension 1 of phylogenetic lineages present in the tree at year $i$ and comparing this location to the average location of phylogenetic lineages present in the tree at year $i-1$.
There may often be large discontinuities in virus locations across the population; our use of difference in mean location is meant to capture both the distance between antigenic clusters and also the change in cluster frequency over consecutive years.
We observe greater heterogeneity in year-to-year antigenic drift in type A than in type B lineages (Figure~\ref{jumps}B), with standard deviation of year-to-year antigenic drift equal to 0.97 units in A/H3N2, 0.66 units in A/H1N1, 0.46 units in B/Vic and 0.26 units in B/Yam.
This analysis classifies drift only to the level of consecutive years; some coarse-graining of the timings of transition events will necessarily occur.
%%% jumps %%%
\begin{figure}[h]
\centering
\includegraphics[width=1.0\textwidth]{figures/jumps}
\caption{\textbf{Year-to-year antigenic drift between 1992 and 2011 in A/H3N2, A/H1N1, B/Vic and B/Yam viruses.}
(A) Timeseries of year-to-year antigenic drift between 1992 and 2011 in A/H3N2, A/H1N1, B/Vic and B/Yam viruses.
Colored lines represent year-to-year antigenic drift, where drift for year $i$ is measured as the mean of antigenic dimension 1 of phylogenetic lineages in year $i$ compared to the mean of antigenic dimension 1 of phylogenetic lineages from the previous year $i-1$.
For example, 2000 represents difference in antigenic dimension 1 between viruses from 1999 and 2000.
Error bars represent 50\% Bayesian credible intervals of year-to-year drift.
Gray dotted lines represent lineage-specific seasonal incidence in the USA taken as average influenza-like illness (ILI) multiplied by proportion of viruses attributable to a lineage for each season.
Here, 2000 represents the 2000/2001 influenza season.
(B) Distribution of year-to-year antigenic drift between 1992 and 2011 in A/H3N2, A/H1N1, B/Vic and B/Yam viruses.
}
\label{jumps}
\end{figure}
We investigate the relationship between rates of antigenic drift and seasonal incidence in the USA in A/H3N2, A/H1N1, B/Vic and B/Yam.
We measure seasonal incidence from USA CDC influenza surveillance reports for each virus lineage (A/H3N2, A/H1N1, B/Vic, B/Yam) by taking the average influenza-like (ILI) percentage in a season and multiplying this by the relative proportion of virus isolations attributable to a particular influenza lineage in a season (see Materials and methods).
This measure of incidence has previously been shown to have have predictive power in the analysis of seasonal influenza trends \cite{Goldstein11}.
We analyze incidence from the 1998/1999 to the 2008/2009 seasons to avoid possible complications from the 2009 pandemic.
We begin by comparing overall rates of antigenic drift (Figure~\ref{drift}, Table~\ref{paramtable}) to overall levels of seasonal incidence across influenza lineages, finding a significant correlation between rate of antigenic drift and relative incidence across the four lineages (Pearson correlation, $r=0.97$, $p = 0.041$).
We follow-up this analysis with a more detailed analysis of year-to-year variation in antigenic drift and lineage-specific incidence, comparing incidence in a season to antigenic drift of viruses coming into this season (Figure~\ref{jumps}A).
For example, we compare antigenic drift of viruses from 2000 to 2001 to incidence in the 2001/2002 season.
Within each virus lineage, we find that years with pronounced antigenic drift tend to show increased incidence (Figure~\ref{corr}), finding Pearson correlation coefficients of 0.51, 0.29, 0.44 and 0.14 for A/H3N2, A/H1N1, B/Vic and B/Yam respectively.
We calculated significance using bootstrap permutation tests finding $p$-values of 0.056, 0.201, 0.097 and 0.341 respectively.
We applied a similar bootstrap permutation test to calculate the significance of finding the observed degree of correlation across all four lineages, arriving at a $p$-value of 0.018.
The fact that we observe periods of pronounced antigenic drift preceding increased incidence in each of the four influenza lineages suggests a causal relationship, in which antigenic evolution drives increased incidence.
%%% corr %%%
\begin{figure}[h]
\centering
\includegraphics[width=0.7\textwidth]{figures/corr}
\caption{\textbf{Relationship between antigenic drift and seasonal incidence for years 1998 to 2009 in influenza A/H3N2, A/H1N1, B/Vic and B/Yam.}
Antigenic drift from year $i-1$ to year $i$ is compared to incidence in the season $i$ / $i+1$.
For example, year-to-year antigenic drift from 2000 to 2001 is measured against incidence in the 2001/2002 season.
}
\label{corr}
\end{figure}
However, it is possible that if sampling count influences cartographic estimates then a spurious correlation could arise in which years with greater incidence have higher sample counts and artifactually high estimates of drift.
We controlled for this possibility by testing to see if viral isolate count influences estimates of year-to-year drift by correlating drift between years $i$ and $i-1$ against the ratio of the number of isolates from year $i$ to the number of isolates from year $i-1$.
We found little correlation when combining data across lineages (Figure~\ref{counts_corr}, Pearson's $r = -0.01$).
We tested significance following the same bootstrap procedure we used to assess the correlation between drift and incidence, finding a $p$-value of 0.514.
Separating lineages gave $p$-values of 0.717, 0.246, 0.337, 0.504 for A/H3N2, A/H1N1, B/Vic and B/Yam, respectively.
These findings suggest our results to be unbiased with regard to sample count.
%%% counts_corr %%%
\begin{figure}[h]
\centering
\includegraphics[width=0.4\textwidth]{figures/counts_corr}
\caption{\textbf{Relationship between antigenic drift and sample counts for years 1998 to 2011 in influenza A/H3N2, A/H1N1, B/Vic and B/Yam.}
Antigenic drift from year $i-1$ to year $i$ is compared to the ratio of sample counts in year $i$ to counts in year $i-1$.
Only comparisons which had one or more samples in years $i-1$ and $i$ were retained, leaving 11 A/H3N2, 7 A/H1N1, 9 B/Vic and 10 B/Yam comparisons.
Points are colored according to influenza lineage based on the color scheme in Figure \ref{corr}.
}
\label{counts_corr}
\end{figure}
%%% CONCLUSIONS %%%
\section*{Conclusions}
Understanding antigenic evolution in seasonal influenza is crucial to our efforts of surveillance and control.
Cartographic methods allow complex HI datasets to be compressed to more approachable location-based summaries that quantify antigenic relationships between strains, including relationships not directly measured via HI.
In this study we provide a foundation for evolutionary antigenic cartography, which seeks to simultaneously assess antigenic phenotype and antigenic evolution.
We use this approach to characterize competitive dynamics across influenza lineages A/H3N2, A/H1N1, B/Vic and B/Yam and show that antigenic evolution within each lineage drives strain replacement and contributes to seasonal incidence patterns.
We find that influenza A/H3N2 evolves faster in antigenic phenotype than A/H1N1, which in turn evolves faster than B/Vic or B/Yam.
Consequently, the influenza A/H3N2 virus population turns over more quickly than A/H1N1 or influenza B, exhibiting a smaller effective population size and a `spindlier' phylogenetic tree.
Furthermore, we observe a correlation between antigenic drift and viral incidence both across and within influenza lineages.
The finding that antigenic evolution correlates with subsequent increased incidence within a lineage suggests a causal role for antigenic drift driving influenza incidence patterns.
The correlation between incidence and drift further suggests the possibility of using HI data at the start of an influenza season to predict which lineage will subsequently predominate.
The statistical framework presented here represents a baseline to which further advancements in modeling antigenic phenotype and evolution may be made.
For example, our likelihood-based model facilitates the inclusion of possible covariates affecting immunological titer, which could include experimental factors such as red blood cell type used in the HI assay \cite{Lin12} and whether oseltamivir is included in the HI reaction \cite{Lin10}.
Additionally, this framework should be ideally suited to uncovering genetic determinants of antigenic change, as both the sequence state and antigenic location of internal nodes in the phylogeny may be estimated.
In this fashion, it should be possible to correlate sequence substitutions directly to antigenic diffusion.
Identifying viruses that will come to predominate in the global virus population while they are still at low frequency remains an enormous challenge.
However, combining evolutionary and antigenic information may eventually prove useful in identifying low-frequency, but expanding, lineages of antigenically novel viruses that represent ideal targets for vaccine strain selection.
%%% METHODS %%%
\section*{Materials and methods}
\subsection*{Antigenic cartography}
Antigenic characteristics of viral strains are often assessed through immunological assays such as the hemagglutination inhibition (HI) assay \cite{Hirst43}.
At heart, these assays compare the reactivity of one virus strain to antibodies raised against another virus strain via challenge or vaccination.
In the case of HI, the measurement of cross-reactivity takes the form of a titer representing the dilution factor at which serum raised against a particular virus ceases to be effective at inhibiting the binding of another virus to red blood cells.
These factors are commonly assessed by serial dilution, so that HI titers will form a log series, 40, 80, 160, etc \dots.
Because experimental HI titers typically differ by factors of two, we find it convenient to work in log$_2$ space and represent the titer of virus $i$ against serum $j$ as $H_{ij} = \mathrm{log}_2 \mbox{(HI titer)}$, i.e.\ a titer of 160 has $H_{ij} = 7.32$.
Due to experimental constraints, most comparisons cannot be made, leading to a sparse observation matrix $\mathbf{H} = \{H_{ij}\}$.
Further, measurements are usually interval and truncated, e.g.\ inhibition may cease somewhere between the serial titers of 160 and 320, or inhibition may be absent at all titers assayed, suggesting a threshold somewhere between 0 and 40.
Previous work \cite{Smith04, Cai10} has used multidimensional scaling (MDS) to place viruses and sera on an `antigenic map'.
These methods heuristically optimize locations of viruses and sera by seeking to minimize the sum of squared errors between titers predicted by map locations and observed titers.
Antigenic maps produced by these methods have proven useful in categorizing virus phenotypes \cite{Smith04}, but the extension of these methods to integrate genetic data remains notably lacking.
Here, we follow previous models in representing antigenic locations as points in a low $P$-dimensional antigenic map.
One of our initial goals is to find an optimal projection of the high-dimensional distance matrix $\mathbf{H}$ into this lower dimensional space.
We conduct this projection using Bayesian multidimensional scaling (BMDS) \cite{Oh01} in which we construct a probabilistic model to quantify the fit of a particular configuration of cartographic locations to the observed matrix of serological measurements.
Typically, $P = 2$, but higher or lower dimensions may better reflect the data.
Let $\virus_i \in \mathbb{R}^{P}$ represent the cartographic location of virus $i$ for $i = 1,\ldots,\vn$, so that $\virus_i = (x_{i1}, x_{i2})^{\prime}$ for $P=2$.
Similarly, let $\serum_j$ represent the cartographic location of serum $j$ for $j = 1,\ldots,\sn$, so that $\serum_j = (y_{j1},y_{j2})^{\prime}$ for $P=2$.
For notational compactness, we collect together all virus coordinates into an $\vn \times P$ matrix $\viruses = (\virus_1, \ldots, \virus_{\vn})^{\prime}$ and all serum coordinates into an $\sn \times P$ matrix $\sera = (\serum_{1},\ldots,\serum_{\sn})^{\prime}$.
Virus and serum may be isolated from / raised against the same strain and have different cartographic locations, and separate serum isolates raised against the same strain may also have different cartographic locations.
This gives a set of distances between virus and serum cartographic locations
\begin{equation}
\delta_{ij} = || \virus_i - \serum_j ||_2,
\end{equation}
where $|| \cdot ||_2$ is an $L_2$ norm.
Traditional approaches to antigenic cartography \cite{Smith04} begin by defining immunological distance as
\begin{equation}
d_{ij} = \se_j - H_{ij},
\end{equation}
where $H_{ij}$ is the log$_2$ titer of virus $i$ against serum $j$ and serum potency $\se_j = \max ( H_{1j},\ldots,H_{\vn j} )$ is fixed.
In following multidimensional scaling, these approaches attempt to optimize over unknown $\viruses$ and $\sera$ such that
\begin{equation} \label{mds}
\sum_{(i,j) \in \cal I}
\left(
\delta_{ij} - d_{ij}
\right)^2
\end{equation}
is minimized, where $\mathcal{I} = \{ (i,j) : H_{ij} \mbox{ is measured} \}$.
In the case of threshold measurements, this error function is modified slightly; see \cite{Smith04} for further details.
Here, we instead assume a probabilistic interpretation in which an observed titer is normally distributed around its cartographic expectation with variance $\mdssd^2$,
\begin{equation} \label{hij}
H_{ij} \sim \normal( \se_j - \delta_{ij}, \, \mdssd^2 ).
\end{equation}
Consequently, the likelihood of observing an exact titer given the placement of antigenic locations is
\begin{equation}
\point(H_{ij}) = \phi \left( \frac{ H_{ij} + \delta_{ij} - \se_j }{ \mdssd } \right),
\end{equation}
where $\phi(\cdot)$ represents the standard normal probability density function (PDF).
Previous BMDS has employed a sampling density truncated to strictly positive quantities since $d_{ij}$ are directly observed, non-negative quantities.
In the antigenic setting, these remain random and can be negative since neither $\se_j$ is known nor is $H_{ij}$ observed with much precision.
HI assays sometimes show no inhibition at all measured titrations, e.g.\ a measurement can be reported as `$<$40'.
In this case, the likelihood of observing the threshold measurement follows the cumulative density of the lower tail of the normal distribution
\begin{equation}
\threshold(H_{ij}) = \Phi \left( \frac{ H_{ij} + \delta_{ij} - \se_j }{ \mdssd } \right),
\end{equation}
where $\Phi(\cdot)$ represents the standard normal cumulative distribution function (CDF).
Although it is simplest to assume that immunological measurements represent point estimates, it seems more natural to assume that the threshold for inhibition occurs between two titers, e.g.\ we observe inhibition at 1:160 dilution and no inhibition at 1:320 dilution.
Rather than taking the HI titer as 160, we can instead treat this as an interval measurement, assuming that the exact titer for inhibition would occur somewhere between 160 and 320.
HI titers are usually reported as the highest titer that successfully inhibits virus binding, so that in this case, we calculate the likelihood of an interval measurement as
\begin{equation}
\interval(H_{ij}) = \Phi \left( \frac{ H_{ij} + \delta_{ij} - \se_j + 1 }{ \mdssd } \right) - \Phi \left( \frac{ H_{ij} + \delta_{ij} - \se_j }{\mdssd} \right).
\end{equation}
These likelihoods are illustrated in Figure~\ref{hij_likelihood}.
Throughout our analyses, we use interval likelihoods $\interval$ rather than point likelihoods $\point$ unless otherwise noted.
%%% hij_likelihood %%%
\begin{figure}[tb]
\centering
\includegraphics[width=0.6\textwidth]{figures/hij_likelihood}
\caption{\textbf{Likelihood of HI titers in the BMDS model.}
Here we show the likelihoods of observing three different outcomes given $\delta_{ij} = 4$, $\mdssd = 0.95$ and $\se_j = \mathrm{log}_2 (1280) = 10.32$.
The likelihood of observing a threshold titer of `$<$40' is equal to the lower tail of the probability density function $\threshold(5.32) = 0.146$.
The likelihood of observing a point measurement with an exact inhibiting titer of `90.5' is equal to the density function $\point(6.5) = 0.413$.
The likelihood of observing an interval measurement with an inhibiting titer somewhere between `160' and `320' is equal to $\interval(7.32) = 0.129.$
}
\label{hij_likelihood}
\end{figure}
We calculate the overall likelihood by multiplying probabilities of individual measurements
\begin{equation} \label{bmds}
L(\viruses,\sera) = \prod_{(i,j) \in \cal I} f(H_{ij}),
\end{equation}
using probability functions $\point$, $\threshold$ and $\interval$ as appropriate.
We begin by assuming independent, diffuse normal priors on virus and serum locations
\begin{eqnarray}
\virus_i &\sim& \normal(\mathbf{m},\boldsymbol{\Sigma}) \nonumber \\
\serum_j &\sim& \normal(\mathbf{m},\boldsymbol{\Sigma}),
\end{eqnarray}
where $\mathbf{m} = (0,\ldots,0)^{\prime}$ and $\boldsymbol{\Sigma}$ is a diagonal matrix with diagonal elements all equal to $10000$.
\subsection*{Virus avidity and serum potency}
The preceding model represents immunological distance as a drop in titer against the most reactive comparison for a particular serum.
However, this model may be biased in some circumstances.
In one example, if a particular serum $j$ is only measured against distant viruses, its maximum titer will be artificially low and the likelihoods concerning this serum will appear poor.
To address this issue, we relax the assumption of fixed $\se_j$ values and treat the expected log$_2$ titer when $\delta_{ij}=0$ as a random variable.
In this case, $H_{ij}$ still follows equation \ref{hij} with expectation $\se_j - \delta_{ij}$, but the vector of `serum potencies' $\ses = (\se_1,\ldots,\se_{\sn})$ is random and estimated rather than fixed.
We assume that $\se_j$ values are hierarchically distributed according to a normal distribution.
We take an Empirical Bayesian approach in specifying the mean and variance of this distribution, set to the empirical mean and empirical variance of the set of maximum titers across sera $\{ \max ( H_{1j},\ldots,H_{\vn j} ) : j = 1,\ldots,\sn \}$.
This formulation assumes that particular sera are more reactive in general than other sera.
Additionally, we follow the same logic and assume that some virus isolates are more reactive than other virus isolates and include a parameter for `virus avidity' $\ve_i$ representing the general level of reactivity across HI assays.
With virus avidity included, observed titers follow
\begin{equation}
H_{ij} \sim \normal \left( \frac{\ve_i+\se_j}{2} - \delta_{ij}, \, \mdssd^2 \right),
\end{equation}
and the vector of virus avidities $\ve_i$ for $i = 1,\ldots, \vn$ is estimated in an analogous hierarchical fashion, with $\ves$ normally distributed with mean and variance equal to the empirical mean and variance of the set of maximum titers across viruses $\{ \max ( H_{i1},\ldots,H_{i \sn} ) : i = 1,\ldots,\vn \}$.
\subsection*{Drift model of antigenic evolution}
As presented, multiple configurations of virus and serum locations $\viruses$ and $\sera$ will give the same likelihood of an observed data matrix $\mathbf{H}$.
An example of this phenomenon is shown in Figure~\ref{schematic_map}.
In this case, it is impossible to determine from the HI data at hand whether the blue and yellow viruses are antigenically similar (Figure~\ref{schematic_map}A) or antigenically divergent (Figure~\ref{schematic_map}B).
This presents an issue of model identifiability, where absolute, as opposed to relative, antigenic locations cannot be determined from observing the serological data alone.
Thus, in order to achieve more a interpretable model we impose a weak prior on global locations.
In influenza, it's clear that antigenic distance between strains increases with time \cite{Smith04,Cai10}.
To capture this, we replace our previous diffuse prior with an informed prior in which the expected location of viruses and sera increases with date of sampling along dimension one, and each virus and serum location follows an independent normal distribution centered around this temporal expectation, so that
\begin{eqnarray}
x_{i1} &\sim& \drift \, t_i + \normal(0, \virussd^2) \nonumber \\
y_{j1} &\sim& \drift \, t_j + \normal(0, \serumsd^2),
\end{eqnarray}
where $t$ is the difference between the date of the indexed virus or serum and the date of the earliest sampled virus or serum, and other dimensions follow $x_{im} \sim \normal(0, \virussd^2)$ and $y_{jm} \sim \normal(0, \serumsd^2)$ for $m\ge2$.
Thus, this model assumes that virus and serum locations drift in a line across the antigenic map at rate $\drift$.
The parameter $\virussd$ determines the breadth of the cloud of virus locations at each point in time, while $\serumsd$ determines the breadth of the cloud of serum locations.
%%% schematic_map %%%
\begin{figure}[tb]
\centering
\includegraphics[width=0.95\textwidth]{figures/schematic_map}
\caption{\textbf{Schematic antigenic map with three viruses and two sera.}
(A) Map with virus 1 and virus 3 antigenically similar.
(B) Map with virus 1 and virus 3 antigenically divergent.
Virus 1 is shown in blue, virus 2 is shown in red and virus 3 is shown in yellow.
Virus isolates are represented by filled circles, sera raised against viruses are shown as open circles and map distances $\delta_{ij}$ are shown as solid lines connecting viruses and sera.
Sera from virus 1 is compared against viruses 1 and 2, while sera from virus 2 is compared against viruses 2 and 3.
Configurations (A) and (B) represent cartographic models that would give equal likelihoods to a set of serological data $\{H_{11},H_{21},H_{22},H_{32}\}$.
}
\label{schematic_map}
\end{figure}
\subsection*{Phylogenetic diffusion model of antigenic evolution}
We simultaneously model antigenic locations and genetic relatedness by assuming that virus locations are influenced by evolution following a Brownian motion process \cite{Lemey10}.
To do this, we replace the previous prior specifying independent virus locations with a prior that incorporates covariance based on shared evolutionary history
\begin{equation} \label{ebpprior}
\viruses \sim \left( \begin{matrix} \drift \, t_1 & 0 \\ \vdots & \vdots \\ \drift \, t_n & 0 \end{matrix} \right) + \, \mbox{Evolutionary Brownian Process}(\virussd, \tree)
\end{equation}
for $P=2$, where $\virussd$ is the volatility parameter of the Brownian motion over virus locations and $\tree$ is a phylogeny specifying tree topology and branch lengths.
Thus, viruses which are genetically similar are induced to have prior locations close to one another on the antigenic map.
In the evolutionary Brownian process, the tips of the phylogeny $\tree$ correspond to the set of virus locations $(\virus_1, \ldots, \virus_\vn)$, and the probability of observing tip locations depends on the locations of internal nodes $(\virus_{n+1}, \ldots, \virus_{2\vn-2})$ and on the location of the root node $\virus_{2\vn-1}$.
This process assumes that a virus location $\virus_i$ follows from the location of its parent virus $\virus_{f(i)}$, and with the addition of drift along dimension 1, is distributed as
\begin{equation}
\virus_i \sim (\drift \, d_i, \, 0)^{\prime} + \normal(\virus_{f(i)}, d_i \, \boldsymbol{\Sigma})
\end{equation}
for $P=2$, where $f(i)$ is a function that maps nodes to parental nodes, $d_i$ is the length of the branch connecting virus $i$ to parent virus $f(i)$, and $\boldsymbol{\Sigma}$ is a diagonal matrix with diagonal elements all equal to $\virussd^2$.
The root virus location $\virus_{2\vn-1}$ is assumed to follow a normal distribution with expectation $(\drift \, t_{2\vn-1}, 0)^{\prime}$ for $P=2$ and variance determined by the diffusion volatility $\virussd$ \cite{Lemey10}.
The probability of virus locations $p(\viruses | \drift, \virussd, \tree)$ is determined through analytical integration across internal states following the methods introduced in \cite{Lemey10}.
This formulation corresponds to a Wiener process with drift, in which the drift term $\drift$ only influences the expected states of nodes along the phylogeny, but does not influence the covariance structure among these nodes, which remains the same as it does in a standard Wiener process \cite{BrownianMotionHandbook}.
This allows the separation in equation \ref{ebpprior} between drift terms affecting only expectations and the evolutionary Brownian process that includes covariance among virus locations $\virus_1,\ldots,\virus_n$.
Here, the phylogenetic tree $\tree$ is estimated using sequence data for viruses $1,\ldots,\vn$ according to well established methods implemented in the software package BEAST \cite{BEAST17}.
\subsection*{Posterior inference}
Top-level priors for $1/\mdssd^2$, $\drift$, $1/\virussd^2$, and $1/\serumsd^2$ are assumed to follow diffuse $\mbox{Gamma}(a, b)$ distributions with $a=0.001$ and $b=0.001$.
These diffuse priors were chosen to be non-informative and provide little-to-no weight on the resulting posterior distributions.
Under the full model, the posterior probability of observing virus and serum locations given immunological data is factored
\begin{equation}
p(\viruses,\sera | \mathbf{H}) \propto p(\mathbf{H} | \viruses, \sera, \ses, \ves, \mdssd) \;
p(\viruses | \drift, \virussd, \tree) \;
p(\sera | \drift, \serumsd) \;
p(\ses, \ves, \mdssd, \drift, \virussd, \serumsd, \tree).
\end{equation}
We sample from this posterior distribution using the MCMC procedures implemented in the software package BEAST \cite{BEAST17}.
Metropolis-Hastings proposals include transition kernels that translate individual virus and serum locations $\virus_i$ and $\serum_j$ and individual virus avidities $\ve_i$ and serum potencies $\se_j$, and other transition kernels that scale the entire set of virus and serum locations $\viruses$ and $\sera$ and that scale parameters $\mdssd$, $\drift$, $\virussd$ and $\serumsd$.
For the present analysis, a two-step approach was taken to sample phylogenies, where a posterior sample of phylogenies was gathered using sequence data and then, in the cartographic analysis, trees from this set were randomly proposed and accepted following the Metropolis-Hastings algorithm \cite{Pagel04}.
\subsection*{Genetic, antigenic and surveillance data}
We compiled an antigenic dataset of hemagglutination inhibition (HI) measurements of virus isolates against post-infection ferret sera for influenza A/H3N2 by collecting data from previous publications \cite{Hay01,Smith04,Russell08,Barr10}, NIMR vaccine strain selection reports for 2002 and 2008--2012 \cite{NIMR02,NIMRMarch08,NIMRFeb09,NIMRFeb10,NIMRSep10,NIMRSep11,NIMRFeb12} and the Feb 2011 VRBPAC report \cite{Cox11FDA}.
We queried the Influenza Research Database \cite{IRD} and the EpiFlu Database \cite{GISAID} for HA nucleotide sequences by matching strain names, e.g.\ A/HongKong/1/1968, and only strains for which sequence was present was retained.
If a strain had multiple sequences in the databases we preferentially kept the IRD sequence and preferentially kept the longest sequence in IRD.
Many strains had full length HA sequences, while other strains only possessed HA1 sequences.
Sequences were aligned using MUSCLE v3.7 under default parameters \cite{MUSCLE}.
This dataset had 2051 influenza isolates (present as either virus or serum in HI comparisons) dating from 1968 to 2011.
However, the majority of isolates were present from 2002 to 2007.
Because we are interested in longer-term antigenic evolution, we subsampled the data to have at most 20 virus isolates per year, preferentially keeping those isolates with more antigenic comparisons.
We then kept only those serum isolates that are relatively informative to the antigenic placement of viruses, dropping serum isolates that are compared to 4 or fewer different virus isolates.
This censoring left 402 virus isolates, 519 serum isolates and 10,059 HI measurements.
Each virus isolate was compared to an average of 21.9 serum isolates, and each serum isolate was compared to an average of 18.0 virus isolates.
Antigenic data for influenza A/H1N1 was collected from previous publications \cite{Kendal78,Webster79,Nakajima79,Nakajima81,Chakraverty82,Pereira82,Chakraverty86,Cox83,Daniels85,Raymond86,Stevens87,Donatelli93,Hay01,Daum02,McDonald07,Barr10} and NIMR vaccine strain selection reports for 2002--2010 \cite{NIMR02,NIMR03,NIMR04,NIMRFeb05,NIMRSep05,NIMRMarch06,NIMRSep06,NIMRMarch07,NIMRSep07,NIMRMarch08,NIMRSep08,NIMRFeb09,NIMRFeb10}.
The same procedure that was followed for A/H3N2 was repeated to match sequence data and to subsample antigenic comparisons.
This procedure yielded 115 virus isolates, 77 serum isolates and 1882 HI measurements over the course of 1977 to 2009.
Each virus isolate was compared to an average of 10.0 serum isolates, and each serum isolate was compared to an average of 16.2 virus isolates.
Antigenic comparisons for influenza B/Victoria were collated from previous publications \cite{Rota90, Hay01, Muyanga01, Shaw02, Ansaldi04, Puzelli04, Xu04, Barr06, Daum06, Lin07} and vaccine strain selection reports for 2002--2012 \cite{AusWHO06, NIMR02, NIMR03, NIMR04, NIMRFeb05, NIMRSep05, NIMRMarch06, NIMRSep06, NIMRMarch07, NIMRSep07, NIMRMarch08, NIMRFeb09, NIMRSep09, NIMRFeb10, NIMRSep10, NIMRFeb11, NIMRSep11, NIMRFeb12}.
Here, the sequence matching and subsampling procedure yielded 179 virus isolates, 70 serum isolates and 2003 HI measurements over the course of 1986 to 2011.
Each virus isolate was compared to an average of 6.5 serum isolates, and each serum isolate was compared to an average of 16.7 virus isolates.
Antigenic comparisons for influenza B/Yamagata were collected from previous publications \cite{Rota90, Kanegae90, Nakajima92, Nerome98, Hay01, Muyanga01, Nakagawa02, Abed03, Ansaldi03, Ansaldi04, Matsuzaki04, Puzelli04, Shaw02, Xu04, Barr06, Daum06, Lin07} and vaccine strain selection reports for 2002--2012 \cite{AusWHO06, NIMR02, NIMR03, NIMR04, NIMRFeb05, NIMRSep05, NIMRMarch06, NIMRSep06, NIMRMarch07, NIMRSep07, NIMRMarch08, NIMRFeb09, NIMRSep09, NIMRFeb10, NIMRSep10, NIMRFeb11, NIMRSep11, NIMRFeb12}.
For B/Yamagata, the matching and subsampling procedure resulted in 174 virus isolates, 69 serum isolates and 1962 HI measurements over the course of 1987 to 2011.
Each virus isolate was compared to an average of 6.9 serum isolates, and each serum isolate was compared to an average of 17.3 virus isolates.
Surveillance data was obtained from the Centers of Disease Control and Prevention FluView Influenza Reports from the yearly summaries of influenza seasons 1997--1998 to 2010--2011 \cite{CDCReports}.
As an example, one report states ``collaborating laboratories in the United States tested 195,744 respiratory specimens for influenza viruses, 27,682 (14\%) of which were positive. Of these, 18,175 (66\%) were positive for influenza A viruses, and 9,507 (34\%) were positive for influenza B viruses. Of the 18,175 specimens positive for influenza A viruses, 7,631 (42\%) were subtyped; 6,762 (87\%) of these were seasonal influenza A (H1N1) viruses, and 869 (13\%) were influenza A (H3N2) viruses.''
In this case, we estimate the relative proportion of A/H3N2 of the four lineages as $0.66 \times 0.13 = 0.09$.
Similar calculations were performed for A/H1N1, B/Vic and B/Yam.
\subsection*{Implementation}
Phylogenetic trees were estimated for A/H3N2, A/H1N1, B/Vic and B/Yam using BEAST \cite{BEAST17} and incorporated the SRD06 nucleotide substitution model \cite{Shapiro06}, a coalescent demographic model with constant effective population size and a strict molecular clock across branches.
MCMC was run for 60 million steps and trees were sampled every 50,000 steps after allowing a burn-in of 10 million steps, yielding a total sample of 2000 trees.
These trees were treated as a discrete set of possibilities when subsequently sampled in the BMDS analysis \cite{Pagel04}.
However, it would be possible to jointly sample from sequence data and serological data using these methods.
MCMC was used to sample virus locations $\viruses$, serum locations $\sera$, virus avidities $\ves$, serum potencies $\ses$, MDS precision $1/\mdssd^2$, antigenic drift rate $\drift$, virus location precision $1/\virussd^2$, serum location precision $1/\serumsd^2$ and phylogenetic tree $\tree$.
MCMC chains were run for 500 million steps and parameter values sampled every 200,000 steps after a burn-in of 100 million steps, yielding a total of 2000 MCMC samples.
In all cases when drift parameter $\drift$ was included the MCMC chain mixed well and arrived at the same estimated posterior distribution from different starting points.
However, without drift parameter $\drift$, maps for A/H3N2 showed some degree of metastability, where some chains would converge on one solution and other chains would converge on a different solution.
We favor models that include $\drift$, because its inclusion, in addition to correcting most identifiability issues, yields much improved mixing of antigenic locations.
There is some difficulty summarizing posterior cartographic samples, as sampled virus and serum locations represent only relative quantities, and because of this, over the course of the MCMC, virus locations may shift.
Our prior on virus and serum locations removes much of this issue, orienting the antigenic map along dimension 1 and fixing it to begin at the origin.
However, local isometries are often still a problem.
For example, in A/H3N2 the HK/68, EN/72 and VI/75 clusters may rotate in relation to other clusters.
Consequently, it may be difficult to fully align MCMC samples using Procrustes analysis.
For the present study, we take a simple approach and sample a single MCMC step and visualize the antigenic locations at this state (Figures~\ref{map}, Figure~\ref{drift}).
Then, for specific quantities of interest, like rate of antigenic drift and rate of diffusion at different points along the phylogeny, we calculate the quantity across MCMC samples to yield an expectation and a credible interval.
This approach accurately characterizes uncertainty that may be hidden in an analysis of a single antigenic map.
We summarize diffusion paths of viral lineages (Figure~\ref{map}, Figure~\ref{drift}) by taking each virus and reconstructing $x$ and $y$ locations along antigenic dimensions 1 and 2 backward through time.
We use MCMC to sample tip locations, but when outputting trees sample internal node locations using a peeling algorithm as described in \cite{Pybus12}.
Thus, after the MCMC is finished we have a posterior sample of 2000 trees each tagged with estimated tip locations and internal node locations.
We post-processed each posterior tree by conducting a linear interpolation between parent-child node locations to arrive at $x$ and $y$ values at intervals of 0.05 years for each virus.
Then, for each interval, $x$ and $y$ values are averaged across the sample of posterior trees.
We draw lines between these locations to approximate mean posterior diffusion paths.
As virus lineages coalesce backwards through time down the phylogeny these diffusion paths will also coalesce.
\subsection*{Comparison with previous results}
Here, we attempt to compare antigenic locations inferred by our BMDS model to antigenic locations previously inferred by the error minimization methods of Smith et al.\ \cite{Smith04}, referred to here as antigenic cartography by MDS.
For this comparison, we use exactly the same HI data used to produce the results in \cite{Smith04}, consisting of 273 virus isolates, 79 serum isolates and a total of 4252 HI measurements taken between 1968 and 2003.
We begin with a BMDS analog of the antigenic model used in \cite{Smith04}, where serum potencies are taken as the maximum titer of a particular ferret serum and the expected log$_2$ drop in HI titer is proportional to Euclidean distance between virus and serum locations.
To bring models into further alignment, we use a Uniform$(-100,100)$ distribution over virus locations and serum locations.
Unsurprisingly, we find that this BMDS model produces results that are strongly congruent with MDS results (Figure~\ref{smith_comparison}).
Antigenic cluster locations are consistent between methods (Figure~\ref{smith_comparison}A--B) and antigenic distances between pairs of viruses are consistent between temporally similar and temporally divergent viruses (Figure~\ref{smith_comparison}C--E), suggesting that the resulting maps are consistent at both local and global scales.
Credible intervals of antigenic distances for the BMDS model remain narrow across the temporal spectrum (Figure~\ref{smith_comparison}C--E), implying a fair degree of rigidity to the map.
%%% smith_comparison %%%
\begin{figure}[h]
\centering
\includegraphics[width=1.0\textwidth]{figures/smith_comparison}
\caption{\textbf{Comparison of A/H3N2 antigenic locations estimated by Smith et al.\ \cite{Smith04} using MDS and an equivalent BMDS model.}
(A) MDS antigenic locations, reoriented so that the primary dimension lies on the $x$-axis rather than on the $y$-axis as in Figure 1 of \cite{Smith04}.
(B) A posterior sample of antigenic locations from an equivalent BMDS model.
In (A) and (B), viruses are shown as colored circles, with color denoting antigenic cluster inferred by \cite{Smith04}, and sera are shown as gray points.
(C) Antigenic distances between A/Bilthoven/16190/1968 and all other viruses determined for both methods.
(D) Antigenic distances between A/Fujian/411/2002 and all other viruses determined for both methods.
(E) Antigenic distances between 750 random pairs of viruses determined for both methods.
In (C), (D) and (E) red points show distances for the MDS model and gray bars show the 95\% credible interval of distances for the BMDS model, while the red dashed line shows a LOESS regression to MDS distances, and the black dashed line shows a LOESS regression to the BMDS distances.
The BMDS model has a Uniform$(-100,100)$ prior on antigenic locations and serum potencies fixed at maximum titer values.
}
\label{smith_comparison}
\end{figure}
Smith et al.\ \cite{Smith04} show that there exist at least two solutions in their assignment of antigenic locations, involving the rotation of clusters HK/68, EN/72 and VI/75 (shown in Figure S2 of \cite{Smith04}).
We observe the same metastable behavior in our analysis; some MCMC chains converge on the solution shown in Figure \ref{smith_comparison}B, while other MCMC chains converge on the alternative solution shown in Figure \ref{smith_comparison_rotation}B.
The distribution of likelihood values appear highly similar between these two solutions, suggesting that they represent global optima.
The rotation of the HK/68, EN/72 and VI/75 clusters creates a map that bends slightly, so that temporally distant viruses appear closer in the rotated solution than in the original solution (Figure~\ref{smith_comparison_rotation}C--E).
In this case, it is clear that the solutions are locally consistent between viruses up to $\sim$15 years divergent, even if there some degree of global flexibility.
%%% smith_comparison_rotation %%%
\begin{figure}[h]
\centering
\includegraphics[width=1.0\textwidth]{figures/smith_comparison_rotation}
\caption{\textbf{Comparison of A/H3N2 antigenic locations estimated by Smith et al.\ \cite{Smith04} using MDS and an equivalent BMDS model under an alternative solution.}
(A) MDS antigenic locations, reoriented so that the primary dimension lies on the $x$-axis rather than on the $y$-axis as in Figure 1 of \cite{Smith04}.
(B) A posterior sample of antigenic locations from an equivalent BMDS model that has converged on the alternative solution.
In (A) and (B), viruses are shown as colored circles, with color denoting antigenic cluster inferred by \cite{Smith04}, and sera are shown as gray points.
(C) Antigenic distances between A/Bilthoven/16190/1968 and all other viruses determined for both methods.
(D) Antigenic distances between A/Fujian/411/2002 and all other viruses determined for both methods.
(E) Antigenic distances between 750 random pairs of viruses determined for both methods.
In (C), (D) and (E) red points show distances for the MDS model and gray bars show the 95\% credible interval of distances for the BMDS model, while the red dashed line shows a LOESS regression to MDS distances, and the black dashed line shows a LOESS regression to the BMDS distances.
The BMDS model has a Uniform$(-100,100)$ prior on antigenic locations and serum potencies fixed at maximum titer values.
}
\label{smith_comparison_rotation}
\end{figure}
As discussed in the main text, the presence of multiple optima with different degrees of 2D curvature implies an issue of identifiability; the HI likelihood model alone cannot distinguish between these possibilities.
Because of this issue, and to more easily estimate rates of antigenic drift, we include a model of systematic drift in antigenic location that favors linear movement in the antigenic map.
We find that including this drift prior on antigenic locations removes the problem of identifiability.
Antigenic locations produced by this model remain locally consistent with MDS results between viruses $\sim$15 years divergent, but global comparisons show that this BMDS model has partitioned more variance to the first antigenic dimension (Figure~\ref{smith_comparison_drift}).
We additionally find that including the drift prior on antigenic locations often results in greater predictive power, with a slight improvement of test error for the A/H1N1, B/Vic and B/Yam datasets (Table~\ref{errortable}).
%%% smith_comparison_drift %%%
\begin{figure}[h]
\centering
\includegraphics[width=1.0\textwidth]{figures/smith_comparison_drift}
\caption{\textbf{Comparison of A/H3N2 antigenic locations estimated by Smith et al.\ \cite{Smith04} using MDS and an extended BMDS model that includes date-informed priors on antigenic locations.}
(A) MDS antigenic locations, reoriented so that the primary dimension lies on the $x$-axis rather than on the $y$-axis as in Figure 1 of \cite{Smith04}.
(B) A posterior sample of antigenic locations from a BMDS model that includes date-informed priors on antigenic locations.
In (A) and (B), viruses are shown as colored circles, with color denoting antigenic cluster inferred by \cite{Smith04}, and sera are shown as gray points.
(C) Antigenic distances between A/Bilthoven/16190/1968 and all other viruses determined for both methods.
(D) Antigenic distances between A/Fujian/411/2002 and all other viruses determined for both methods.
(E) Antigenic distances between 750 random pairs of viruses determined for both methods.
In (C), (D) and (E) red points show distances for the MDS model and gray bars show the 95\% credible interval of distances for the BMDS model, while the red dashed line shows a LOESS regression to MDS distances, and the black dashed line shows a LOESS regression to the BMDS distances.
The BMDS model has a date-informed prior on antigenic locations and serum potencies fixed at maximum titer values.
}
\label{smith_comparison_drift}
\end{figure}
Our final BMDS model (model 9, Table~\ref{errortable}) differs from antigenic model used by Smith et al.\ \cite{Smith04} in including temporally- and phylogenetically-informed priors on antigenic locations and also in estimating serum and virus avidities.
Here, we investigate the impact on antigenic locations of estimating virus avidity and serum potency in the BMDS model.
To isolate this difference, we use a Uniform$(-100,100)$ prior on antigenic locations.
Surprisingly, estimating virus avidity and serum potency results in a more linear antigenic map (Figure~\ref{smith_comparison_effects}), resembling the appearance of the map incorporating the antigenic drift prior, while preserving local consistency.
We generally observe congruence between MDS and BMDS antigenic locations for viruses less than $\sim$10 years divergent (Figure~\ref{smith_comparison_effects}E).
However, specific viruses may be affected, for instance A/Bilthoven/16190/1968 (Figure~\ref{smith_comparison_effects}C), which appears more distant from all other viruses when serum and virus avidities are included.
%%% smith_comparison_effects %%%
\begin{figure}[h]
\centering
\includegraphics[width=1.0\textwidth]{figures/smith_comparison_effects}
\caption{\textbf{Comparison of A/H3N2 antigenic locations estimated by Smith et al.\ \cite{Smith04} using MDS and an extended BMDS model that estimates serum and virus avidities.}
(A) MDS antigenic locations, reoriented so that the primary dimension lies on the $x$-axis rather than on the $y$-axis as in Figure 1 of \cite{Smith04}.
(B) A posterior sample of antigenic locations from a BMDS model that estimates virus avidity and serum potency.
In (A) and (B), viruses are shown as colored circles, with color denoting antigenic cluster inferred by \cite{Smith04}, and sera are shown as gray points.
(C) Antigenic distances between A/Bilthoven/16190/1968 and all other viruses determined for both methods.
(D) Antigenic distances between A/Fujian/411/2002 and all other viruses determined for both methods.
(E) Antigenic distances between 750 random pairs of viruses determined for both methods.
In (C), (D) and (E) red points show distances for the MDS model and gray bars show the 95\% credible interval of distances for the BMDS model, while the red dashed line shows a LOESS regression to MDS distances, and the black dashed line shows a LOESS regression to the BMDS distances.
The BMDS model has a Uniform$(-100,100)$ prior on antigenic locations and virus avidities and serum potencies estimated in a hierarchical Bayesian fashion.
}
\label{smith_comparison_effects}
\end{figure}
In this dataset, viruses 15 or more years divergent always yield threshold titers, and hence, their relative locations must be indirectly inferred rather than through direct comparison.
This may explain why we observe local consistency between models at scales less than $\sim$15 years, but some degree of global inconsistency.
Still, these results suggest that, when making local comparisons, such as those used to calculate year-to-year antigenic drift (Figure~\ref{drift}), outcomes are expected to be robust to many model particulars.
\subsection*{Availability}
Source code implementing the cartographic models has been made fully available as part of the software package BEAST \cite{BEAST17}, and can be downloaded from its Google code repository (http://code.google.com/p/beast-mcmc/).
More details on implementing these models can be found at https://github.com/trvrb/flux/tree/master/example-xmls.
Incidence data and HI data used in this analysis is archived with Dryad (doi:10.5061/dryad.rc515).
%%% ACKNOWLEDGMENTS %%%
\subsection*{Acknowledgments}
We thank Richard Reeve, Dan Haydon, and Simon Frost for insights on antigenic modeling and MDS.
We acknowledge the laboratories that provided sequences to EpiFlu database: Centers for Disease Control and Prevention (USA), Chinese Center of Disease Prevention and Control, Hospital Clinic of Barcelona, National Institute of Hygiene of Morocco, National Institute of Infectious Diseases (Japan), National Institute for Medical Research (UK), Norwegian Institute of Public Health, Swedish Institute for Infectious Disease Control, Victorian Infectious Diseases Reference Laboratory (Australia).
%%% AUTHOR CONTRIBUTIONS %%%
\subsection*{Author contributions}
AR, MAS, TB conceived the study.
AR, MAS, TB, DJS, PL designed the antigenic model.
TB, GD, VG, JWM, AJH, CAR, DJS, AR gathered antigenic and genetic data.
AR, MAS, TB, PL implemented statistical phylogenetic procedures.
TB, AR performed the analysis.
TB, AR, DJS, JWM, PL, CAR, AJH, GD, MAS interpreted the results.
TB, MAS, AR, JWM, AJH wrote the paper.
%%% FUNDING %%%
\subsection*{Funding}
TB was supported by a Newton International Fellowship from the Royal Society.
MAS was supported by National Institutes of Health grants R01 HG006139 and R01 AI107034 and National Science Foundation grant DMS126153 and IIS1251151.
The research leading to these results has received funding from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013) under Grant Agreement no. 278433-PREDEMICS and ERC Grant agreement no. 260864.
VG, AJH and JWM were funded by MRC reference no U117512723.
AR acknowledges the support of the Wellcome Trust (grant no. 092807).
Collaboration between MAS, AR and PL was supported by the National Evolutionary Synthesis Center (NESCent), NSF EF-0423641.
DJS and CAR acknowledge EU FP7 programs EMPERIE (223498) and ANTIGONE (278976), Human Frontier Science Program (HFSP) program grant P0050/2008, Wellcome Trust 087982AIA, NIH Director's Pioneer Award DP1-OD000490-01, and National Institute of Allergy and Infectious Diseases NIH CEIRS contract HHSN266200700010C.
CAR was supported by a University Research Fellowship from the Royal Society.
%%% REFERENCES %%%
\bibliographystyle{plos}
\bibliography{flux}
\end{document}