New methods for analyzing serological data with applications to influenza surveillance.

BACKGROUND: Two important challenges to the use of serological assays for influenza surveillance include the substantial amount of experimental effort involved and the inherent noisiness of serological data.
RESULTS: I show that log-transformed serological data exist in an effectively one-dimensional space. I use this result, together with new mechanistic insights into serological assays, to develop computational methods for accurately and efficiently recovering unmeasured serological data from a sample of measured data, for systematically minimizing noise and other types of non-antigenic variation found in the data, and for quantifying and visualizing antigenic variation. The methods can also be applied to data with effective dimensionality greater than one, under certain conditions.
CONCLUSION: Careful application of the methods developed here would enable the collection of better-quality serological data on a greater number of circulating influenza viruses than is currently possible and improve the ability to identify potential epidemic and pandemic viruses before they become widespread. Although the focus here is on influenza surveillance, the described methods are more widely applicable.

Introduction

Serological data on circulating influenza viruses generally contain evolutionarily important information about functional (antigenic) variation in the B‐cell antigens of those viruses. , Because natural selection acts on antigenic variation, serological data can provide important insights into patterns, causes, and epidemiological consequences of influenza viral evolution. , , , , Nevertheless, there are important challenges to the use of serological data. In particular, serological assays require considerable amounts of time and effort to perform, and this limits the number of potential epidemic and pandemic viruses on which serological data can be routinely collected. In addition, serological data are often contaminated by experimental noise (e.g. resulting from the serial dilution of sera), which may cause independently measured data for the same virus and serum sample to vary greatly. Furthermore, the data depend on non‐antigenic variables such as the red cell avidity and the antibody‐inducing capacity of viruses, , which can make it difficult to obtain accurate information about the antigenic variation of influenza viruses.

Smith et al. recently made great progress toward addressing the above challenges, focusing on data obtained from the widely used hemagglutination‐inhibition assay. Those data are typically reported in a tabular format, with rows corresponding to viruses and columns to sera. The entry found in the i’th row and j’th column of the table (called the titer of virus i relative to serum j) represents the reciprocal of the maximum dilution of serum j that can effectively neutralize virus i. Smith et al. used titers from multiple tables to construct an ‘antigenic’ map of viruses and sera in which the Euclidean distance between virus i and serum j was approximated by the negative log (base 2) of the titer of virus i relative to serum j normalized by the maximum titer obtained for that serum. (Use of base 2 in the log‐transformation of titers reflects the fact that titers are typically measured using twofold dilutions of sera. Base 2 is also used here when log‐transforming titers). The constructed antigenic map allowed unavailable normalized titers to be accurately predicted from distances measured on the map, increasing the amount of available data. Because the map was constructed using data from multiple tables, data recovered from it are, in principle, less noisy than data obtained from individual tables.

Here, I introduce computational methods that extend the above antigenic‐map approach into novel directions. First, I show that on average, log‐transformed titers obtained in the same experiment (i.e. reported in the same table) have only one effective dimension. Then, I introduce a computationally efficient method that exploits this low dimensionality to accurately recover a complete table of titers using only a subset of those titers; various quantities derived from titers (e.g. the normalized titers of Ref. 2) can thus be recovered. I also show how ‘noise‐free’ estimates and 95% confidence intervals (CIs) can be computed for titers and, hence, for related quantities. In addition, I derive a mathematical relation that sheds light on the mechanistic basis of titers and allows non‐antigenic variation to be filtered out from these titers [Parts 1 and 2 of Supporting Information (SI)]. I introduce a method for quantifying and visualizing antigenic variation using such filtered titers, and I illustrate this method using titers for novel H1N1 viruses. The above methods were partly motivated by remarkable recent work on the recoverability of data that exist in a low‐dimensional space.

Methods

Computing the effective rank and the coherence of a table of titers

Let be the singular value decomposition (SVD) of an m × n table H of log‐transformed titers. (‘T’ denotes matrix transpose and ‘*’ denotes matrix multiplication). The columns (called eigenvectors) of U (V) are orthonormal bases for the column (row) space of H, whereas the diagonal entries of S are the n singular values of H, , sorted in decreasing order of magnitude. The square of the i’th singular value (i.e. the i’th eigenvalue of ) quantifies the variation explained by the i’th eigenvector; therefore, the fraction of the variation explained by the first r eigenvectors is given by:

The rank of H is defined as the number of its non‐zero singular values. Because some eigenvectors associated with non‐zero singular values may explain negligible amounts of the variation in titers (e.g. because they contain noise), the rank of H may be greater than its effective rank, defined here as the smallest value of r for which F r  = F r+1, r = 1,…,n‐1 (the effective rank is set to n if this condition is not satisfied for any r < n).

Alternatively, assuming that the noise found in log‐transformed titers obtained in the same experiment is normally distributed, the eigenvectors containing noise would be associated with eigenvalues that are significantly smaller than the other eigenvalues. In this case, the effective rank can be estimated by finding the smallest eigenvalue that is significantly greater than all the eigenvalues smaller than it, using Fisher’s variance ratio test. This is equivalent to finding the maximum value of r for which , : where denotes the inverse cumulative function of the F‐distribution with 1 and n‐r degrees of freedom, and α is the level of statistical significance of the test. Here, I used , where the factor (n‐r) represents a Holm‐Bonferroni correction for multiple comparisons. Both of the above methods for determining the effective rank gave similar results.

Let P U (P V) be the orthogonal projection of H onto the first r eigenvectors found in U (V). The coherence of H (more precisely, the maximum coherence of U and V) with respect to the standard bases in R m and R n is given by :

The lower the coherence, the greater the probability that measured titers, selected uniformly at random from H, will contain enough information about the unmeasured titers to enable their recovery by the method described elsewhere.

Recovering titers

Let H denote an m × n table of log‐transformed titers of rank r, m ≥ n, and let Ω denote a subset of ∼rm 1·2log (m) titers randomly selected from H with uniform probability. (r = 1 is used in this study). The unselected (or ‘unmeasured’) titers are recovered by finding an m × n matrix X that minimizes : where μ and η are Lagrange multipliers. denotes the nuclear norm of X, that is, the sum of the singular values of X. When r is known (it is not approximated by the effective rank of H), is replaced by the sum of the r largest singular values of X. Note that (Eqn. 4) was solved using the fixed point continuation algorithm ; for the experiments reported here, it took only a few seconds, on average, to recover titers found in each analyzed table using a Pentium IV machine running Windows XP.

Because titers recovered by the above method are theoretically exactly equal to their ‘noise‐free’ values when the selected titers are free of noise, discrepancies between the recovered and the noise‐free titers are necessarily due to noise found in the selected titers. If there is no systematic bias in the way that noise found in the selected titers induces variation in the recovered titers, then the recovered titers would be randomly distributed about the corresponding noise‐free titers. The distribution of the recovered titers can therefore be used to compute CIs for the noise‐free titers. This is the rationale for the following procedure for computing CIs for titers found in H: 1. Randomly select m 1·2log (m) titers from H. 2. Recover the unselected titers (see above). 3. Repeat steps 1 and 2 until each titer found in H is recovered at least N (=1000) times (on average it will take nN/[n‐m 0·2log(m)] repetitions of steps 1 and 2 for this to happen). Let L i be a list of the k ≥ N values recovered for the ith titer, which are sorted in increasing order. Then, the lower (upper) 95% CI for the ith titer is given by the ‘th (‘th) element of L i, where () denotes the largest (smallest) integer smaller (larger) than x. The mean of the recovered values for each titer is used as a noise‐free estimate for that titer. Note that if H is incomplete, then 95% CIs can still be computed for titers found in any m 1 × n 1 complete sub‐table of H, provided m 1 1·2log(m 1) < m 1 n 1.

In practice, when designing an experiment to measure titers for m viruses relative to n sera, the experimentalist can measure only a random subset of at least m 1·2log(m) titers and subsequently recover the unmeasured titers as well as missing titers. In addition, noise‐free estimates and 95% CIs for all titers (respectively subsets of titers) found in a complete (respectively incomplete) table can be computed and used in place of the original titers, which may be unreliable.

Quantifying and visualizing antigenic differences between viruses

Let be the SVD of a table H of log‐transformed titers. To quantify antigenic differences between viruses found in H, H is projected onto an r‐dimensional subspace of its row space: , where V r denotes the first r eigenvectors found in V. The antigenic difference between viruses i and j is defined as the Euclidean norm of the difference between the i’th and j’th rows of W. If r ≤ 3, then antigenic differences can be visualized by plotting the rows of W. Antigenic differences computed using titers from different tables can be embedded in a common r‐dimensional space by means of probabilistic multidimensional scaling (Part 2 of SI). Note that the above SVD decomposition of H is only feasible when no titers are missing from H. In practice, missing titers should be recovered before H is subjected to SVD (see Part 2 of SI for an SVD approach that is applicable to tables with missing titers).

Results and discussion

Let H denote an m × n table of titers, with m ≥ n. Because it is tedious to measure all mn titers found in H, it will be very helpful if fewer titers can be measured and used to recover the unmeasured titers. Recent theoretical work suggests that if both the rank r and the coherence of H (see Methods for definitions) are not much greater than one, then unmeasured titers can be recovered exactly if only ∼rm 1·2log (m) titers (selected randomly with uniform probability) are measured. I investigated the applicability of this method for recovering unmeasured titers by computing both the rank and the coherence of published tables of empirical titers for influenza A and B viruses (Methods). (Because missing titers can artificially lower the rank of a table, only complete tables were analyzed; see Table S1). The analyzed tables, which typically consist of viruses sampled non‐randomly from multiple influenza seasons, should contain more variation on average than tables consisting of viruses of the same subtype sampled randomly from only one season. Nevertheless, when the titers found in each table are log‐transformed, on average ∼99% of the variation in those titers is explained by only one eigenvector, which is associated with the largest singular value of the table (Figure 1A). This suggests that the analyzed tables have an average effective rank of ∼1. I obtained a similar result by using a statistical F‐test for the effective dimensionality of each table (Methods). In addition, the tables have an average coherence of 1·36 ± 0·12. Note that deviations from an effective rank of 1 may be caused by noise and other factors, including the normalization of titers (Part 3 of SI). Also, based on well‐known results on the ranks of partitioned matrices, when multiple tables of titers are combined the effective rank of the resulting composite table can be >1, especially if the rank of any one of the original tables is also >1 (Part 3 of SI). The above factors may partly explain the higher dimensionality of normalized titers reported previously.

Figure 1

 Dimensionality and recoverability of noise‐free titers. (A) The dimensionality of each of 23 tables of empirical titers (Table S1) was investigated by log‐transforming the titers found in each table and subsequently determining the fraction of the variation in titers (denoted F r) that is explained by the eigenvectors associated with the r largest singular values of each table (Methods). F r (averaged over all 23 tables) is plotted against r, for r = 1,…,5. (B) A table consisting of simulated titers was contaminated by noise (independently drawn from a normal distribution with mean 0 and standard deviation δ), and the resulting noisy titers were used to compute estimates of and 95% CIs for the uncontaminated (or noise‐free) titers (Methods). The fraction of noise‐free titers that occurred within their corresponding CIs is plotted against δ. (C) The mean ratio of the width of the 95% CI for a particular noise‐free titer to the absolute value of that titer is plotted against δ. (D) The mean absolute difference between the estimated and noise‐free titers, between the estimated and noisy titers, and between the noisy and noise‐free titers, is plotted against δ. Bars denote standard deviations.

The fact that both the effective rank and the coherence of the analyzed tables are not much greater than one suggests unmeasured titers can be recovered exactly by the above‐mentioned method. Note, however, that because measured titers are likely to contain noise, exact recovery of unmeasured titers may not be possible. Also, recovery is not possible for rows/columns of H that do not contain any measured titer. Furthermore, accurate (not necessarily exact) recovery is only possible if the number of measured titers is not smaller than the number of degrees of freedom of H, which is given by r (m + n − r), a condition that is satisfied by all the tables analyzed earlier. In practice, each row/column must contain more than r titers in order for titers missing from that row/column to be accurately recovered. It is important to keep in mind that exact recovery of unmeasured titers is still theoretically possible if the effective rank of H is >1, provided conditions (i) and (ii) above are satisfied. In addition to unmeasured titers, the above conditions are also applicable to the recovery of ‘missing’ titers – titers that are presumed to be either too high or too low to be measured experimentally. (Note that some missing titers may actually be within the limits of experimental resolution, but they could not be measured because of noise).

I applied the above method for recovering unmeasured titers to tables of empirical titers for influenza A and B viruses (Table S1). The titers found in each table were log‐transformed, and 90% of those titers were randomly selected and used to recover the unselected (‘unmeasured’) titers (Methods). This procedure was repeated 100 times for each table. [Note that it was necessary to limit this analysis to m × n tables for which 0·9mn ≥ m 1·2log (m).] Both the mean absolute difference (0·66 ± 0·68) and the relative mean absolute difference (0·08 ± 0·10) between the recovered and the unmeasured titers were small, suggesting that the recovery of unmeasured titers was accurate, albeit not exact. A closer look at the above results reveals that the mean absolute difference varies across the analyzed tables, ranging from 0·15 to 2·50. This may reflect the differing amounts of noise found in titers from different tables.

To minimize noise found in empirical titers, I developed a novel method for computing noise‐free estimates and 95% CIs for titers (Methods). The method was tested using 20 × 10 rank‐1 tables consisting of simulated, log‐transformed titers. Each table was constructed as the product of two matrices (of dimensions 20 × 1 and 1 × 10, respectively), the entries of which were independently drawn from a truncated normal distribution with support (1·5, 4); this support was chosen to simulate the dynamic range of log‐transformed, empirical titers (3·32, 13·32). (Note that because empirical titers are measured on a geometric scale, it is reasonable to approximate their distribution by a lognormal distribution. ) To simulate noise found in empirical titers, for each constructed table M, another table M* was obtained by adding a perturbation (independently drawn from a normal distribution with mean 0 and standard deviation δ) to each entry of M. M* was then used to compute estimates of and 95% CIs for the noise‐free titers in M (Methods), for δ = 0, 0·1, 0·2,…, 1. [Note that when δ = 1 there is a 32% chance that, on a geometric scale, a titer will be either multiplied or divided by at least 2 (the dilution factor often used in serological assays), representing a high amount of noise.] Representative results obtained using one of the constructed tables are shown in Figure 1.

The results show that the computed 95% CIs have excellent coverage properties; they contain their corresponding noise‐free titers in >99% of cases and they also have small relative widths (Figure 1B,C). Remarkably, although estimates for the noise‐free titers were computed using noisy titers, the mean absolute difference between those estimates and the noise‐free titers is much smaller than between the estimated and noisy titers and between the noisy and noise‐free titers (Figure 1D). Indeed, the mean absolute difference between the estimated and noise‐free titers grows much more slowly with δ than does the mean absolute difference between the estimated and noisy titers (Figure 1D). These results suggest that the developed method can, in principle, be used to systematically minimize noise found in tables of empirical titers. As an illustration, I computed noise‐free estimates and 95% CIs for the empirical titers analyzed earlier. The average absolute difference between each titer and its noise‐free estimate is small (0·66 ± 0·67). However, 60% of the titers do not occur in their associated 95% CI, suggesting that those titers may be unreliable. This is consistent with the high variation in independently measured normalized titers for the same virus and antiserum that is observed in many cases (Table S1).

In addition to noise, it is also important to minimize other types of non‐antigenic variation found in titers. This is currently not possible because of limited understanding of the fundamental nature of titers. To shed light on the nature of titers, I developed and used a mechanistic model of hemagglutination inhibition to derive the first explicit mathematical equation for the titer of virus i relative to serum j [Eqn. (S23) of SI]: where A j denotes the concentration of antibodies found in serum j, K ij the average affinity of those antibodies for virus i, and J i a dimensionless quantity that depends on such non‐antigenic variables as the avidity of virus i for red cell, the concentration of virus i, etc. (see Part 1 of SI for additional details). A j depends on non‐antigenic variables, including the antibody‐inducing capacity of the virus against which serum j was raised and the immune status of the organisms in which that serum was raised. The derived equation predicts that the normalized titer H  jj/H  ij, a commonly used measure of antigenic difference, is approximately independent of A j, but it depends on both J i and J j. In contrast, a measure of antigenic difference introduced by Archetti and Horsfall – [H  ii H  jj/(H  ij H  ji)]1/2– is predicted to be approximately independent of the non‐antigenic variables A j, J i, and J j, suggesting that it may be more accurate. This is consistent with previous empirical results.

Importantly, the derived equation predicts that by mean‐centering each row (resp. each row and column) of a table of log‐transformed, un‐normalized (resp. normalized) titers the dependence of those titers on non‐antigenic variables would be minimized. I developed a method for quantifying and visualizing antigenic differences between viruses using such mean‐centered tables (Methods). As an illustration, I used the method to construct a one‐dimensional map showing the serological relationships among 123 H1N1 viruses, 118 of which are novel H1N1 viruses (Figure 2). The map shows that the serological responses of novel H1N1 viruses are on average much more similar to each other than to the corresponding responses of earlier H1N1 viruses, particularly a virus isolated in 1930 and another isolated in 1976 (Figure 2). In addition, the map suggests the existence of a well‐defined direction of antigenic evolution of the viruses. Because ‘serological’ maps, such as the one described here, have only one dimension, they may facilitate computational prediction of short‐term changes in the serological attributes of epidemic viruses.

Figure 2

 One‐dimensional serological map of 123 H1N1 viruses. Four tables of titers (Supplementary Data S1) that were recently collected by the US. Centers for Disease Control and Prevention were analyzed. The tables had dimensions 26 by 10, 36 by 10, 40 by 12, and 41 by 11, respectively. Between 10% and 26% of titers were missing from each table. Measured titers found in each table were log‐transformed and used to recover the missing titers, except for titers missing from rows/columns in which fewer than two titers were measured. (In particular, titers for the Brisbane/59/07 virus were not recovered). Each completed table was row‐centered and used to quantify differences between the serological responses of pairs of viruses (Methods). These serological (or ‘antigenic’) differences were subsequently mapped onto a one‐dimensional space by means of probabilistic multidimensional scaling (Part 2 of SI). The viral coordinates thus obtained are shown in the plot. The width of each displayed rectangle represents a unit of serological difference. The novel H1N1 vaccine virus (CA09) is labeled, as are five other viruses that were isolated before 2009. The arrow indicates the apparent direction of temporal changes in the serological properties of the analyzed viruses. Virus labels: CA09 California 07/2009, IL07 Illinois/9(33304)/2007, NJ76 New Jersey/8/1976, SW30 Swine/Iowa/1930, TX08 Texas/14/2008, WI98 Wisconsin/10/1998.

In summary, the computational methods presented in this paper suggest new possibilities for improving the use of hemagglutination‐inhibition titers (and serological data in general) for influenza surveillance. In particular, the method for computing estimates and 95% CIs for noise‐free titers would allow uncertainties associated with titers to be taken into account, for example, when selecting viruses for use in influenza vaccines. Also, the method for minimizing non‐antigenic variation found in titers may help to improve the estimation of antigenic differences between viruses. In addition, the method for recovering unmeasured titers may help to reduce the experimental effort required to collect titers; for example, only 201·2log(20) = 109 titers may need to be measured in order to accurately determine all 200 possible titers for 20 viruses relative to 10 sera. This method can be applied to other types of serological data, especially when the number of viruses under consideration is greater than or equal to the fifth power of the dimensionality of the data. The additional experimental capacity made available by the method would increase the number of viruses, circulating in humans and other organisms, for which serological data can be routinely collected and thereby improve the likelihood that potential epidemic and pandemic viruses will be identified before they become widespread. This is important in light of the recent pandemic spread of an influenza virus whose initial circulation in humans may have gone undetected for several months.

Note that the observation, reported here, that tables of empirical titers have an average effective rank of ∼1 – titers exist in a space that has only one effective dimension – suggests that the effective number of independent variables that determine titers is very small. In other words, while titers depend on many variables [see Equation (5) and Part 1 of SI] some of which may be mutually independent, variation in titers may be dominated by one or more co‐dependent variables. In so far as titers are determined by and contain information about viral phenotypes, the demonstrated low dimensionality of the space of titers and the consequent high recoverability of unknown titers suggests that it may be possible to predict biophysically accessible viral phenotypes from limited information about extant viral phenotypes. Well‐designed theoretical and experimental tests of this idea may yield important new insights and also shed light on questions concerning the evolutionary accessibility of epidemic variants of influenza viruses. In addition to influenza surveillance, the new methods presented in this paper can be used in the surveillance of other pathogens and in the investigation of basic questions concerning pathogen evolution and dynamics. As should be the case for existing methods (e.g. recently developed methods for elucidating the antigenic structure , and the adaptation potential of pathogen populations), continued empirical/experimental assessment of the new methods is necessary to ensure that they will remain useful.

Data S1. List of tables of empirical titers for influenza A and B viruses.

Table S1. Supporting text and list of tables of empirical titers.

Supporting info item

Click here for additional data file.

Supporting info item

Click here for additional data file.