Antigenic Cartography

In this Science paper by Katzelnik et al. (2015), use multidimensional scaling to analyze their virus neutralization data by ‘antigenic cartography.’ The in-depth methods are described in their supplement. Unfortunately, as of this writing, the link that is provided in the method is not (no longer?) available at the web location specified. It looks like lispmds has become ACMAPS, based on the website they also cited: https://albertine.antigenic-cartography.org:1168. The notation behind the method is given below. As a disclaimer, some of the text is copied verbatim from the supplement.

We start with a neutralization titer for virus $i$ with antiserum $j$, called $N_{i,j}$.

The table antigenic distances, $D_{i,j}$, are calculated as the difference between the log titer for the virus best neutralized by a given serum, $b_j$, and the log titer for each virus $i$ against serum $j$:

\begin{equation} D_{ij} = \log_2(b_j) - \log_2(N_{ij}) \end{equation}

The map antigenic distances, $d_{i,j}$, are found by minimizing the error function:

\begin{equation} E=\Sigma_{ij} e(D_{ij}, d_{ij})^2 \end{equation}

where the errors for any given ${i,j}$ pair are given by:

\begin{equation} e(D_{ij}, d_{ij})^2=(D_{ij}-d_{ij})^2 \end{equation}

The limits of detection for the assay range from the dilutions 1:10 to 1:20,480. When the neturalization was below the 1:10 limit of detection, the following error formula was used:

\begin{equation} e(D_{ij}, d_{ij})^2=\frac{(D_{ij}-1-d_{ij})^2}{1+e^{-10(D_{ij}-1-d_{ij})}} \end{equation}

At this point, we look at the optimization technique used to obtain the estimates for the $d_{i,j}$.

Conjugate Gradient Optimization

The R code, following the basic conjugate gradient algorithm on Wikipedia, is provided below.

conjgrad <- function(A, b, x, ...){
	r <- b - A%*%x
	p <- r
	rsold <- t(r) %*% r

	for(i in 1:length(b)){
		Ap <- A%*%p
		alpha <- rsold / (t(p)%*%Ap)
		x <- x + alpha%*%t(p)
		r <- r - t(alpha%*%t(A%*%p))
		rsnew <- t(r) %*% r
		if(sqrt(rsnew) < 1e-10){
			break;
		}
		p <- r + t((rsnew/rsold)%*%t(p))
		rsold <- rsnew
	}

	return(t(x))
}

An example can be tested with the following values for A, b, and x:

A <- matrix(c(4,1,1,3), nrow=2, ncol=2)
b <- c(1,2)
x <- c(2,1)

conjgrad(A,b,x)

The result is:

           [,1]
[1,] 0.09090909
[2,] 0.63636364

There is also the optim function in R, that can be used with method='cg' in order to specify the conjugate gradient algorithm for optimization. References for the conjugate gradient have been described by Fletcher and Reeves (1964), and modified by Polak-Ribiere, and Beale-Sorenson. There is a pre-condidtioned method of conjugate gradient optimization that improves the algorithm under certain circumstances.

In this case, the error function E is what we would like to minimize (the fn in optim), and the output values for x should be map values that minimize the value of E. The starting values for x would be random starting values for the map distances. We will specify the initial values of the “parameters” to be optimized with the par argument in optim. At this moment, I am not sure what are good starting values for x, but perhaps they should be roughly within the range of differences that we observe for the table distances.

The optimal map distances are generated by generating random starting coordinates and optimizing, then repeating the method for 5,000 independent methods. The minimum error map plots the distances such that one grid unit represents one two-fold dilution difference (an antigenic unit, AU), and moving away from a point by each additional grid unit multiplies that difference by 2 again.

The software for antigenic mapping (on the website https://albertine.antigenic-cartography.org:1168) is available for local use at: https://github.com/acjim/AcmacsDesktop.