High Dimensional Random Walks Can Appear Low Dimensional: Application to Influenza H3N2 Evolution
James Moore Hasan Ahmed
q-bio.PE

Sectoring in Multi-cell Massive MIMO Systems
    Shahram Shahsavari, Parisa Hassanzadeh, Alexei Ashikhmin, and Elza Erkip
S. Shahsavari, P. Hassanzadeh  and  E. Erkip are with the ECE Department of New York University, Brooklyn, NY. Email: {shahram.shahsavari,ph990, elza}@nyu.eduA. Ashikhmin is with Bell Labs, Nokia, Murray Hill, NJ, USA. Email:  alexei.ashikhmin@nokia.com
    December 30, 2023
==========================================================================================================================================================================================================================================================================================================================================




One important feature of the mammalian immune system is the highly specific binding of antigens to antibodies.
Antibodies generated in response to one infection may also provide some level of cross immunity to other infections. One model to describe this cross immunity is the notion of antigenic space, which assigns each antibody and each virus a point in ℝ^n. Past studies
have suggested the dimensionality of antigenic space, n, may be small.
In this study we show that data from hemagglutination assays suggest a high dimensional random walk (or self avoiding random random walk).
The discrepancy between our result and prior studies is due to the fact that random walks can appear low dimensional according to a variety of analyses.
including principal component analysis (PCA)
and multidimensional scaling (MDS).


1cm


§ INTRODUCTION



 §.§ Antigenic Space

During a viral infection, antibodies bind to viral antigens by recognizing specific epitopes on their surface. The same antibody may provide protection against
other strains of the same virus if the antigens are not too dissimilar. In 1979, Alan Perelson and George Oster defined the idea of antigenic space <cit.>.
They supposed that each antibody and antigen might be described by a vector in ℝ^n. This idea was used as a basis to study the dynamics of the immune response, with n assumed to be a small number <cit.>. Subsequent work to estimate n based on the frequency of cross reactivity resulted in the conclusion that n was around five to eight <cit.>.

The notion of antigenic space has proven particularly popular for understanding the evolution of influenza H3N2 <cit.>. This strain has been circulating in the human population since 1968 and gradually mutating. These mutations can in principle be represented as the movement of the virus through antigenic space. As the antigen moves
it can evade the antibodies elicited by older strains and thus reinfect individuals.
The distance between a viral strain and an antibody can be measured via the hemagglutination inhibition (HAI) assay, in which a viral strain and a serum of antibodies are both added to a culture of red blood cells. If the
antibodies are ineffective against the viral strain then the virions stick to the red blood cells causing them to cluster together (hemagglutinate). However if the antibodies are effective, they will neutralize the virions and inhibit their hemagglutination of the red blood cells. In the former case, the strain and the serum are distant antigenically, whereas in the latter case they are close. By performing serial dilutions of the antibody serum, one can quantify just how close a serum and antibody are.

Points in antigenic space can be inferred from a distance matrix via multidimensional scaling (MDS). 
Low dimensional reconstructions of antigenic space can reproduce the HAI data with high fidelity, and adding new dimensions beyond n=5 does not improve
the quality of the fit <cit.>. Thus it may be tempting to conclude that influenza is evolving in an antigenic space of no more than five dimensions.



 §.§ Outline of results

In this work we will argue that influenza H3N2 is evolving in a very high dimensional space, and that it may appear to be low dimensional due to the nature of random walks. Our argument consists of three parts.


  * High dimensional random walks contain most of their variance in a small number of dimensions. Specifically, one would expect at least 60% of the variance
to occur in a single dimension. We show this via principal component analysis. However, we note that the true dimensionality of random walks can be revealed by out of sample PCA.

  * This apparent low dimensionality also occurs with multidimensional scaling — the method used to analyze HAI titer data. We simulate HAI data generated using a high dimensional random walk and find we can accurately reproduce the data with points taken from a low dimensional space; increasing the dimensionality of our representation
does not improve our fit. We also show that non-metric MDS can be used to reconstruct an infinite dimensional random walk in a single dimension.

  * Finally, we show that H3N2 data has characteristics of a high dimensional random walk, indicating that it was unlikely to result from a random walk of dimensionality less than n=10. This is even the case when we consider that the random walk of H3N2 is likely self avoiding.




 §.§ Why a random walk?

Throughout this paper we argue for a high dimensional random walk as a model for influenza evolution. A random walk may seem a priori to
be a poor model for viral evolution, as immunological memory should prevent a virus from revisiting areas of antigenic space. Therefore we should expect the path of viral evolution to be self avoiding. In high dimensions an unbiased random walk and self avoiding random walk will behave very similarly, because a high dimensional random walk is already extremely unlikely to cross itself.
In an n dimensional random walk the distance between points i and j is a random variable D_ij. Its distribution is

    D_ij^2∼α |i-j|χ^2_n

where α is a constant of proportionality. This means that for large n the distances increase in a very predictable manner as the χ_n^2 distribution narrows. The probability of the random walk approaching a previous point is essentially zero, so we need not include any further tendency for self avoidance. However, in the latter part of the paper we will
address the question as to whether low dimensional self avoiding random walk could also be consistent with the data.



 §.§ True dimensionality vs effective dimensionality

Let x_i ∈ℝ^n represent distinct viral strains and/or antisera.
The antigenic dissimilarity of the two different strains x_i and x_j
is the euclidean distance D_ij=x_i-x_j_2. D_ij represents the true distance between
these two strains and n is the true dimensionality of antigen space.

If we can represent each point x_i with a corresponding point y_i ∈ℝ^k
such that the distances from y_i to y_j is approximately D_ij, then we say that the effective dimensionality of {x_i} is k.

There are several possible ways to obtain the points y_i and to evaluate how closely the reconstructed distances match the true or measured distances. In this paper we shall make use of three: classical MDS, metric MDS and nonmetric MDS. Classical MDS is also commonly known as principal component analysis (PCA). The points y_i are found by orthogonally
projecting the points x_i onto the dimensions that contain the most variance.
Metric MDS simply minimizes the residual between the true distances and the reconstructed distances. Non-metric MDS is similar to metric MDS but the relationship
between the true distances and the reconstructed distances is only assumed to be monotonically increasing as opposed to directly proportional.



§ RESULTS



 §.§ Truly high dimensional random walks have low effective dimensionality according to PCA

We will now show that random walks with high true dimension can have a very low effective dimension. Our analysis will focus
on infinite dimensional random walks, but we will show via simulations that finite random walks have similar behavior.
Let x_i∈ℝ^n be the ith step in our random walk in n dimensions.
Let X be an m by n matrix whose ith row is x_i. The rows of x_i are determined by a seres of random steps b_i ∈ℝ^n,
each entry of which is an independent sample from 𝒩(0,1/√(n))

    x_1   =b_1
    
    x_i-x_i-1   =b_i    i>1

We can state this relationship compactly as LX=B, where the rows of B are b_i and L is an m by m matrix with one on the diagonal
and negative one on the subdiagonal.

To perform PCA, we first center each column of X by premultiplying be the projection matrix P=I-^T. Note that centering each column is just a translation of the data, so all pairwise distances
are preserved. We then compute the eigenvalues of A=PXX^TP, λ_1>λ_2>⋯λ_m. These eigenvalues indicate the variance in each principal
component of X. For this section we shall define effective dimension as

    ED_a=min{d:∑_i=1^k λ_k≥ a tr(A)}

where 0<a<1 is the fraction of variance that we require within the first k dimensions.

The random matrix A=PL^-1BB^TL^-TP is full rank and has a wishart distribution if the number of steps is no greater than the number of dimensions, i.e. m ≤ n. If number of steps is greater than
the number of dimension (n>m), A is singular and has a pseudo-wishart distribution.



  §.§.§ Infinite dimensional random walks

The matrix W=BB^T is an uncorrelated Wishart matrix. The joint probability distribution of the eigenvalues of Wishart matrices is known, but the formula is cumbersome <cit.>. Therefore we will focus our analysis on the limiting behavior as n→∞, i.e. infinite dimensional random walks.
 As n increases, W approaches the m by m identity matrix. (To see this, either note that the i,jth entry of W is <b_i,b_j>). Using this simplification, the matrix A approaches

    A=PL^-1L^-TP

Noting that L^-1 is an upper triangular matrix with all ones, we can calculate tr(A)=(m-1)(m+1)/6.

We can first observe that this matrix has one eigenvalue v= with corresponding eigenvalue λ_1=0. To find the remaining eigenvalues we use
the pseudo-inverse of A, A which has the simple form

    A=
    [  1 -1  0  ⋯  0; -1  2 -1  ⋱  ⋮;  0  ⋱  ⋱  ⋱  0;  ⋮  ⋱ -1  2 -1;  0  ⋯  0 -1  1 ]

The eigenvalues of this matrix are μ_s=2-2cos(π s/m). Therefore the eigenvalues of A are

    λ_k=1/(2-2cos(π s/m))


The fraction of the variance in the first k dimensions is

    ∑_s=1^k λ_s/tr(A)   =∑_s=1^k (1-cos(π s/m))^-1/(m-1)(m+1)   Finite m
       =6/π^2∑_s=1^k 1/s^2    Infinite m

For infinite dimensional random walks, the first principal component contains at least  6/π^2 or roughly 60% of the total variance. The first two components contain 80% of the variance, but the subsequent convergence is slow as 12 components are needed to account for 95% of the variance (Fig. <ref>A,B). These numbers are the limiting behavior as the length of the random walk grows. For shorter walks, more of the variance is in the first few components (Fig. <ref>C, E) and thus the effective dimensionality may be even lower.





  §.§.§ Finite dimensional random walks

In a finite number of dimensions n, the matrix A is a random variable with a correlated Wishart or Pseudo-Wishart distribution. We therefore use simulation to investigate the effective dimensionality. As the true dimensionality increases, the effective dimensionality asymptotically approaches the value for the infinite case, in which 12 dimensions contain 95% of the variance. However, the effective dimensionality for finite random walks is stochastic and in a few cases we find that it may exceed 12, although the mean is lower (Fig <ref>A).

Like the effective dimensionality, the percentage of variance in the first principal component is a random variable, Ψ_1. Regardless of the length of the walk or the true dimensionality, the mean of Ψ_1 is very similar to the infinite dimensional case, and the variance of Ψ_1 decreases with the true dimensionality but is independent on the length of the walk (Fig <ref>C,D). This gives the somewhat counter intuitive result that lower dimensional
random walks may have less variance in their first principal component than higher dimensional walks.
When we consider the percentage of variance in the first 5 components Ψ_5, this pattern disappears (Fig. <ref>E,F). The values of Ψ_5 tend to decrease monotonically
as both the true dimensionality and the number of steps increase, asymptotically approaching the prediction of roughly 89% from the infinite dimensional case.



  §.§.§ True dimensionality of random walks can be revealed by subsampling

We have demonstrated that when we perform PCA on a random walk, we find most of the variance in only one dimensoin.
Therefore, it is be difficult to tell whether the data originate from a high or low dimensional random walk using PCA. In this section we demonstrate a way around this problem
using out of sample PCA.

We consider two scenarios: a uniform random walk in which movement in every direction is equally likely, and a nonuniform random walk in which movement is ten times greater — and thus the variance is 100 times greater — in 14 dimensions compared to the other 86. PCA finds the majority of variance is in the first dimension for both of these scenarios (Fig <ref>).
However if we take the projection matrix derived from performing PCA on the first half of the random walk and apply it to the second half, we see that the bias towards the first dimension largely disappears. The approximation improves further if we divide the second half into subsequences, apply the projection to each subsequence and then average the results (see methods for further details). This latter method shows clearly that in the uniform random walk, each principal component has roughly the same variance. It also captures the fact that
14 dimensions contain significantly more variance in the nonuniform random walk. Therefore out of sample PCA can in principle allow one to recover the true dimensionality of a random walk.





 §.§ Truly high dimensional random walks look low dimensional according to MDS



  §.§.§ H3N2 evolution

Derek Smith, Alan Lapedes and colleagues have used dimensional reduction techniques to reconstruct two-dimensional antigenic maps of the evolution of Influenza H3N2 <cit.>.
In these maps, some points represent viruses whereas others represent antisera. Distances can only be computed between a virus and an antisera as opposed to between different
viruses directly. Thus classic MDS, which requires measured distances between all pairs of points, cannot be used for this problem. Instead, they use metric MDS with a slight modification to handle
sub-threshold values. Using this method, they have reconstructed a two-dimensional map of antigenic space. To validate this map, they challenged it to predict
the antigenic distance between strains and antisera that were not used as inputs to build the map. They report a good correlation (roughly 80%) between the predicted and actual
values, and furthermore that the accuracy of their predictions is similar for a two dimensional map and for higher dimensions. This is consistent with prior work showing that no more than five dimensions is required to satisfactorily reconstruct HAI data <cit.>. 

Although these authors make no claims
about the true dimensionality of the underlying space, it may be tempting to conclude that antigenic space has a relatively low dimensionality.
Our work, however, suggests that the antigenic space of H3N2 may appear to have low dimensionality even though the true dimensionality is very high. We therefore created artificial data to to mimic the evolution of H3N2 via a high dimensional random walk.  We then followed the procedure for evaluating the quality of reconstructions outlined in <cit.>.



  * We partitioned the real and artificial distance matrices into a training set (90% of the data) and a test set (10%)

  * Using only the distances in the training set we reconstructed points in antigenic space corresponding to each virus and serum. We varied the dimensionality of antigenic space (k) between 2 and 20. Note our artificial data was created using points in 100 dimensional space.

  * We then uses the reconstructed points to predict the antigenic distances in the test set. We evaluated the quality of these fits both by the root-mean-square error and the correlation
between the actual and predicted distance. These were the same criteria used in <cit.>.


We find that increasing k beyond 2, i.e. adding more dimensionality to our reconstructed antigen space, did not improve the quality of the distance estimates (Fig <ref>A,B). This suggests that even data originating from a high dimensional space can be approximated in two dimensions.





  §.§.§ An extreme case: High dimensional random walks can look one dimensional in non-metric MDS

The goal of metric MDS is to find points x ∈ℝ^k whose pairwise euclidean distances approximately match a set of given distances. In particular, the points x must minimize the stress

    Metric Stress=√((x_i-x_j-d_ij)^2/∑ d^2_ij)

where d_ij is the distance between point i and j.

In practice, the d_ij used as an input to MDS may not represent physical distances but rather a generic measure of dissimilarity that does not behave like a metric. For example,
consider three points with dissimilarities d_12=1, d_23=2 and d_13=100. The dissimilarity cannot be a metric as it violates the triangle inequality. However, there may be an
underlying metric space that generated the points, and the underlying dissimilarity represents a transformation, f, of that metric.
The goal of non-metric MDS is to find both f and values for x that minimize

    Non-metric Stress=√((f(x_i-x_j)-d_ij)^2/∑ d^2_ij)


The distance-squared between steps in an n dimensional random walk follows a χ^2_n distribution.

    D_ij^2∼α |i-j|χ^2_n

where α is a constant of proportionality. In this case we let the dissimilarity d_ij be the euclidean distance
D_ij.
As n →∞, the dissimilarity becomes exactly proportional to √(|i-j|). Thus we can construct a one dimensional map of our random walk
via

    Dissimilarity between point i and point j     d_ij^2    =α |i-j|
    Location of points     x_i    =α i
    Monotonic Transform   f(·)    =√(·)

Plugging (<ref>) into (<ref>) we find that the stress is identically equal to zero. Therefore an infinite dimensional random walk can
be exactly represented in one dimension using non-metric MDS.




 §.§ Distinguishing between low dimensional and high dimensional random walks using distance measures

So far we have demonstrated that high dimensional random walks can appear low dimensional using methods such as MDS or PCA. We therefore ask whether there is any
good way to distinguish a low dimensional from a high dimensional random walk. One characteristic of high dimensional random walks is the curve that forms when the first two principal components are plotted against each other <cit.>. This pattern emerges as the number of dimensions approaches infinity. Fig. <ref>A shows the first two principal
components by year, whose path is reminiscent of the quadratic shape predicted in <cit.>. This pattern is suggestive of a high (FIg <ref>B), rather than low dimensional (Fig <ref>C) random walk.




We next sought an objective measure to distinguish a low dimensional from high dimensional walk.
Once again we take advantage of the fact that the square of the distance between the ith and jth step in an n dimensional random walk is proportional to a χ^2 distribution with n degrees of freedom.

    D_ij^2∼α |i-j|χ^2_n

Let D̅^2_s denote the average distance between two points separated by s steps, i.e.

    D̅^2_s=1/m-s∑_i=1^i-m+sD^2_i, i+s

where m is the number of steps.
As n→∞, the coefficient of variation of χ^2_n → 0, so we should expect D̅^2_s ≈α s for some α.
In other words, for high dimensional random walks plotting D̅^2_s vs s should give a linear relationship.
We therefore grouped all antigenic map points in the H3N2 data by year, representing both viral strains and anti-sera, and computed the mean of each group.
Computing the distance-squared between the annual means and plotting versus the number of years elapsed does indeed produce a linear pattern. In this regard, the H3N2 data once again qualitatively resembles a high dimensional random walk rather than a low dimensional one.(Fig <ref>D-F).

To quantify this resemblance, let τ be the coefficient of variation of x_s=1/sD̅^2_s, i.e

    τ=√(∑_s=1^m-1 x^2_s/(∑_s=1^m-1 x_s)^2-1)

 For high dimensional random walks,
we expect τ to be close to zero. We simulated τ for large numbers of random walks and compared to the value computed for the H3N2 data (Fig <ref>).
We find that the H3N2 data is unlikely to have been produced by a random walk of less than ten dimensions (n<10). Note also that the value of τ for the data may be inflated due to any number of sources of variability not accounted for in the model, such as measurement error or a heavy tailed distribution in annual step size. Therefore d=10 is the lower bound for the dimensionality antigenic space and dimensionality of 30 or more is quite likely. These findings extend to self avoiding random walks, as decribed in the methods section, which behave similarly in this regard.



















§ DISCUSSION

High dimensional random walks can appear to be low dimensional. We have shown that most of their variance will appear in a single principle component, and that their pointwise distances can be well approximated by embedding in a two dimensional space. Therefore when using principal component analysis (PCA) or multidimensional scaling (MDS) we may erroneously find that only a few dimensions are important. We demonstrate that by splitting the data set into two, we can use out of sample PCA to potentially recover the true dimensionality of the random walk. If only MDS is possible, perhaps due to missing data, we demonstrate that high dimensional and low dimensional random walks can be distinguished by testing the relationship between distance and step number.

We apply these ideas to the evolution of influenza H3N2. Using data from the HAI assay different strains from 1968 to 2003 can be readily represented in a two dimensional
antigenic map. The distances on this map represent the level of cross immunity that an antibody response to one strain may provide against another. Given this it may be tempting to conclude that influenza is constrained to move in a two dimensional space or at least is biased to move in two dimensions. However, we caution that antigenic space may in fact be very high dimensional or even infinite dimensional and we would still expect to be able to reconstruct it in two dimensions. In fact, we find that a low dimensional random walk or self avoiding random walk is very unlikely to have produced the H3N2 HAI titer data.

One interesting aspect of influenza H3N2 evolution is that despite the annual variation in strain, the strains have not been diverging. Different strains circulating in a given year tend to be close together antigenically. There have been many models proposed to explain this phenomenon <cit.>. Here we have only shown that the data is consistent with a high dimensional random walk or self avoiding random walk. We assume that the direction of the next step is equally likely to be in any direction, provided that this does not bring the walk back to an already visited region. Therefore our model is consistent with antigenic drift hypotheses as described in <cit.> but not the antigenic thrift hypothesis described in <cit.>.



§ METHODS



 §.§ Out of sample PCA



  * Start with a length 2m random walk in n dimensions.

  * Split the random walk into two halves.

  * Apply PCA to the first half. This will yield an orthogonal projection matrix P∈ℝ^n× n which projects any point onto the principal component axis.

  * Split the second half of the walk into subsequences with a short length, μ. These subsequences consist of steps m+1 to m+μ, m+μ+1 to m+2μ, etc. Let S_i ∈ℝ^n×μ denote the ith such subsequence.

  * Apply P to each subsequence S_i, i.e. compute B_i=PS_i.

  * Let σ^2_i be the variance of each row of B_i.

  * Let σ^2=1/ν∑_i=1^νσ^2_i be the sum of all i, where ν is the number of subsequences (i.e. νμ=m).




 §.§ Generating simulated data

We want to generate N_v viral strains and N_s antisera. We will represent each strain and antisera as a point in ℝ^n, where d is our chosen dimensionality of antigenic space. The generation of the artificial data makes repeated use of random vectors in ℝ^n. We draw these random vectors from N(0,I_d), i.e. their
entries are i.i.d unit normal random variables.
We first generate the viral strains via a random walk. Let V_i ∈ℝ^n represent the ith viral strain in the walk, then

    V_1   =b_1
    
    V_i-V_i-1   =b_i    i>1

where b_i N(0,I_n)

Next we use these strains to randomly generate antisera. Let S_j ∈ℝ^n be the jth antiserum, then

    S_j   =V_p_j+b

where p_j is a random integer between 1 and N_v and c_j N(0,I_n).




 §.§ Multi Dimensional Procedure

In this paper we use the form of multidimensional scaling described in <cit.>. We perform this scaling in two distinct cases: starting from
a matrix of measured HI and starting from artificial data generated as described in the previous section.

Not every viral strain is measured against each serum. Let N_Meas be the number of measurements and let 0<l≤ N_Meas index those measurements.
We can then define

    i_l      Strain corresponding to the lth measurement
    
    j_l      Serum corresponding to the lth measurement
    
    d_l      Distance between strain i_l and serum j_l


In the case of real HI titer data, we convert the titer to a distance using the method described in <cit.>.
When using artificial data, we simply compute d_l=V_i_l-S_j_l.

Next we seek a set of points V̂_i∈ℝ^k and Ŝ_j∈ℝ^k that minimizes the objective function

    F({V̂_i},{Ŝ_j}|{d_l}, {i_l}, {j_l})=∑_l=0^N_Meas(d_l-V̂_i_l-Ŝ_j_l)^2

We use values of k between 2 and 20.

We minimize F numerically using the optim function in R. We use the BFGS method and provide the analytically derived gradient of F for efficient computation.
This is an iterative algorithm that requires an initial guess. We generate these guesses using classic MDS, which is equivalent to principal component analysis. The classic MDS algorithm requires a complete distance matrix. In the case of the artificial data, we can calculate this distance matrix directly. In the case of the real data, we construct a plausible distance matrix by finding the shortest path between any two points whose distance was not directly measured. Note that although this plausible distance matrix is likely incorrect, it is only used to find a good starting guess for our algorithm. The solution found using this method compared favorably to the solution found in <cit.>



 §.§ Self avoiding random walk

The self avoiding random walk behaves similarly to an unbiased random walk accept that the new steps are excluded from getting to close to previous ones.
Each new step in the walk is determined according to the following algorithm.


  * Given a current step x_i, propose a new step x̂=x_i+b∈ℝ^n, where the entries of b are i.i.d with distribution 𝒩(0,1/√(n).

  * Calculate the probability of acceptance via

    P   =∏_j=0^i g(x̂|x_i)
    
    g(x̂|x_i)   =1-e^-x̂-x_i^2


  * With probability P accept the new step and let x_i+1=x̂. With probability 1-P return to step one and propose a new step.





alpha

