
[Arc Institute](https://arcinstitute.org/) recently unveiled the [Virtual Cell Challenge](https://virtualcellchallenge.org/). Participants are required to train a model capable of predicting the effect of silencing a gene in a (partially) unseen cell type, a task they term _context generalization_.
For ML engineers with little to no biology background, the jargon and required context can seem quite daunting. To encourage participation, we recapitulate the challenge in a form better suited to engineers from other disciplines. 

# Background

> **Goal** \
> Train a model to predict the effect on a cell of silencing a gene using CRISPR.

Doing things in the world of atoms is expensive, laborious and error prone. What if we could test thousands of drug candidates without ever touching a petri dish?
This is the goal of the virtual cell challenge — a model (most likely a neural network) that can simulate exactly what happens
to a cell when we change some parameter. Given that tightening your feedback loop is often the best way to speed up progress,
a model capable of doing this accurately would have significant impact.

To train this neural network, we will need data. For the challenge, Arc has curated a dataset of ~300k single-cell RNA sequencing profiles. It may be worthwhile to revisit the [Central Dogma](https://www.khanacademy.org/science/biology/gene-expression-central-dogma/central-dogma-transcription/v/rna-transcription-and-translation) before continuing. This essay will build off of this to provide the ~minimum biology knowledge you'll need for the challenge.

## Training data

The training set consists of a sparse matrix and some associated metadata. More specifically, we have 220k cells, and
for each cell we have a [transcriptome](https://en.wikipedia.org/wiki/Transcriptome). This transcriptome is a sparse row vector, where each
entry is the **raw count of RNA molecules** (transcripts) that the corresponding gene (our column) encodes for. Of the 220k cells,
~38k are _unperturbed_, meaning no gene has been silenced using CRISPR. These control cells are crucial as we will see shortly. 

To understand the dataset more concretely, let's select a gene, TMSB4X (the most frequently silenced gene in the dataset) and compare the number of RNA molecules detected for a control cell and a
perturbed cell. 

![TMSB4X](/avcc/TMSB4X.png)

We can see that the cell with TMSB4X silenced has a greatly reduced number of transcripts compared with the control
cells.

### Modelling the challenge 

The astute among you may be wondering why you don't just measure the count of the RNA molecules before and after
silencing the gene — why do we need the control cells at all? Unfortunately, **reading the transcriptome destroys the cell**, which is a problem reminiscent of the [observer effect](https://en.wikipedia.org/wiki/Observer_effect_(physics)). 

This inability to measure the cell state before and after introduces many issues, as we are forced to use a population of **basal**
(a.k.a control, unperturbed) cells as a reference point. The control cells and perturbed cells are not entirely
homogeneous even prior to the perturbation. This means that we have to now separate out our true signal, the perturbation, from
noise induced by the heterogeneity.

More formally, we can model observed gene expression in perturbed cells as:

$$
\hat{X}_p \sim \hat{T}_p(\mathcal{D}_{\text{basal}}) + H(\mathcal{D}_{\text{basal}}) + \varepsilon, \quad \varepsilon \sim P_\varepsilon 
$$

where:
- $\hat{X}_p$: The observed gene expression measurements in cells with perturbation $p$
- $\mathcal{D}_{\text{basal}}$: The distribution of the unperturbed, baseline cell population.
- $\hat{T}_p(\mathcal{D}_{\text{basal}})$: True effect caused by perturbation $p$ on the population.
- $H(\mathcal{D}_{\text{basal}})$: Biological heterogeneity of the baseline population.
- $\varepsilon$: Experiment-specific technical noise, assumed independent of the unperturbed cell state and $\mathcal{D}_{\text{basal}}$.

# STATE: The baseline from Arc

Prior to the Virtual Cell Challenge, Arc released [STATE](https://arcinstitute.org/manuscripts/State), their own attempt to solve the challenge
using a pair of transformer based models. This serves as a strong baseline for participants to start with, so we will
explore it in detail.

STATE consists of two models, the **State Transition Model** (ST) and the **State Embedding Model** (SE). SE is designed to produce rich semantic embeddings of cells in an effort to improve cross cell type generalization. ST is the "cell simulator", that takes in either a transcriptome of a control cell, or an embedding of a cell produced by SE, along with a one hot encoded vector representing the perturbation of interest, and outputs the perturbed transcriptome. 

## State Transition Model (ST) 

![ST](/avcc/ST.png)

The State Transition Model is a relatively simple transformer with a Llama backbone that operates upon the following: 
1. A set of transcriptomes (or SE embeddings) for covariate matched basal cells.
2. A set of one hot vectors representing our gene perturbation for each cell.

Using a covariate matched set of control cells with paired target cells should assist the model in discerning the
actual effect of our intended perturbation. Both the control set tensor and the perturbation tensor are fed through independent encoders, which are simply 4 layer MLPs with GELU activations.
If working directly in gene expression space (i.e producing a full transcriptome), they pass the output through a learned
decoder.

ST is trained using [Maximum Mean Discrepancy](https://en.wikipedia.org/wiki/Kernel_embedding_of_distributions). Put simply, the model learns to minimize the difference between the two probability distributions.

## State Embedding Model (SE) 

![SE](/avcc/SE.png)

The State Embedding Model is a BERT-like model trained using a masked prediction task. To understand this more deeply, first we have to
take a little detour for some more biological grounding. 

### A little biological detour 

![Alternative Splicing](/avcc/alt_splicing.gif)

A gene consists of _exons_ (protein coding sections) and _introns_ (non-protein coding sections). DNA is first _transcribed_ into pre-mRNA, as shown above. The cell then performs [Alternative Splicing](https://en.wikipedia.org/wiki/Alternative_splicing). This is basically "pick and choose exons", cut out all introns. You can think of the gene as an IKEA manual for making a table. One could also construct a 3 legged table, perhaps an odd bookshelf with some effort, by leaving out some parts. These different objects are analogous to **protein isoforms**, proteins coded for by the same gene.

### Back to the model

With this basic understanding, we can move on to how the SE model works. Remember, our core goal for SE is to create **meaningful
cell embeddings**. To do this, we must first create meaningful gene embeddings.

To produce a single gene embedding, we first obtain the amino acid sequence (e.g $\texttt{SDKPDMAEI}$... for TMSB4X) of all the different protein isoforms encoded for by the gene in question. We then feed these sequences to [ESM2](https://huggingface.co/facebook/esm2_t48_15B_UR50D), a 15B parameter Protein Language Model from FAIR. ESM produces an embedding _per amino acid_, and we mean pool them together to obtain a "transcript" (a.k.a protein isoform) embedding. 

Now we have all of these protein isoform embeddings, we then just mean pool those to get the gene embedding. Next, we project these gene embeddings to our model dimension using a learned encoder as follows:

$$
\tilde{g}_j = \text{SiLU}(\text{LayerNorm}(g_j \mathbf{W}_g + \mathbf{b}_g))
$$

We've now obtained a gene embedding, but what we really want is a _cell embedding_. To do this, Arc represents each cell
as the top 2048 genes ranked by [log fold expression level](https://en.wikipedia.org/wiki/Fold_change#Fold_changes_in_genomics_and_bioinformatics).

We then construct a "cell sentence" from our 2048 gene embeddings as follows:
$$
\tilde{\mathbf{c}}^{(i)} = \left[\mathbf{z}_{\text{cls}}, \tilde{\mathbf{g}}_1^{(i)}, \tilde{\mathbf{g}}_2^{(i)}, \ldots, \tilde{\mathbf{g}}_L^{(i)}, \mathbf{z}_{\text{ds}}\right] \in \mathbb{R}^{(L+2) \times h}
$$

We add a $\texttt{[CLS]}$ token and $\texttt{[DS]}$ token to our sentence. The $\texttt{[CLS]}$ token ends up being used as our "cell embedding" (very BERT-like)
and the $\texttt{[DS]}$ token is used to "disentangle dataset-specific effects". Although the genes are sorted by log fold
expression level, Arc further enforces the magnitude of each genes expression by incorporating the transcriptome in a
fashion analogous to positional embeddings. Through an odd ["soft binning" algorithm](https://github.com/ArcInstitute/state/blob/main/src/state/emb/nn/model.py#L374) and 2 MLPs, they create some
"expression encodings" which they then add to each gene embedding. This should modulate the magnitude of each gene
embedding by how intensely it is expressed in the transcriptome.

To train the model, they mask 1280 genes per cell, and the model is tasked with predicting them. The 1280 genes are
selected such that they have a wide range of expression intensities. For the graphically inclined, the below
demonstrates the construction of the cell sentence. 

![SE PATH](/avcc/SE_path.png)

# Evaluations 

Understanding how your submission will be evaluated is key to success. The 3 evaluation metrics chosen by Arc are **Perturbation Discrimination**, **Differential Expression** and **Mean Average Error**. Given that Mean Average Error is simple and exactly as it sounds, we will omit it from our analysis.

## Perturbation Discrimination

![Perturbation Discrimination](/avcc/pert_disc.png "Perturbation Discrimination")

Perturbation Discrimination intends to evaluate how well your model can uncover _relative differences_ between
perturbations. To do this, we compute the Manhattan distances for all the measured perturbed transcriptomes in the test set (the ground
truth we are trying to predict, $y_t$ and all other perturbed transcriptomes, $y_p^n$) to our predicted transcriptome $\hat{y}_t$. We then rank where the
ground truth lands with respect to all transcriptomes as follows:

$$
r_t = \sum_{p \neq t} \mathbf{1}\{d(\hat{y}_t, y_p) < d(\hat{y}_t, y_t)\}
$$

After, we normalize by the total number of transcriptomes:

$$
\text{PDisc}_t = \frac{r_t}{T}
$$

Where $0$ would be a perfect match. The overall score for your predictions is the mean of all $$\text{PDisc}_t$$. This is then normalized to: 

$$
\text{PDiscNorm} = 1 - 2\text{PDisc}
$$ 

We multiply by 2 as for a random prediction, ~half of the results would be closer and half would be further away.

## Differential Expression

Differential Expression intends to evaluate what fraction of the truly affected genes did you correctly identify as significantly affected. Firstly, for each gene compute a $p$-value using a [Wilcoxon rank-sum test with tie correction](https://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U_test). We do this for both our predicted perturbation distribution and the ground truth perturbation distribution.

Next, we apply the [Benjamini-Hochberg procedure](https://en.wikipedia.org/wiki/False_discovery_rate#Benjamini%E2%80%93Hochberg_procedure), basically some stats to modulate the $p$-values, as with $20,000$ genes and a $p$-value threshold of $0.05$, you'd expect $1,000$ false positives. We denote our set of predicted differentially expressed genes $G_{p,pred}$, and the ground truth set of differentially expressed genes $G_{p,true}$. 

If the size of our set is less than the ground truth set size, take the intersection of the sets, and divide by the true number of differentially expressed genes as follows:

$$
DE_p = \frac{G_{p,pred} \cap G_{p,true}}{n_{p,true}}
$$

If the size of our set is greater than the ground truth set size, select the subset we predict are most differentially expressed (our "most confident" predictions, denoted $\tilde{G}_{p,pred}$), take the intersection with the ground truth set, and then divide by the true number.

$$
DE_p = \frac{\tilde{G}_{p,pred} \cap G_{p,true}}{n_{p,true}}
$$

Do this for all predicted perturbations and take the mean to obtain the final score.

# Conclusion

If a virtual cell can accurately model the change in a cells state in response to perturbations, we can look forward to a
very interesting time in pharma. I hope this post accelerated your understanding of the challenge. 

May the best team win.


Arc recently unveiled the Virtual Cell Challenge - we explore the challenge from a non-biologist's perspective in the hope of encouraging participation from a range of disciplines.

Arc Virtual Cell Challenge: A Primer

With exploding demand for AI inference, many hardware startups are designing
Domain Specific Architectures. Working backwards from the Transformer workload, we will identify optimal design
choices, promising hardware and how to predict the future of inference.


Domain specific architectures for AI inference

Rotary Position Encoding (RoPE) can seem intimidating at first. This post endeavours to rediscover RoPE by iteratively improving upon proposed positional encoding schemes.

You could have designed state of the art Positional Encoding

Smartglasses are a form of Augmented Reality (AR) providing significant compute power. What are the current fundamental limitations prohibiting the consumer adoption of smartglasses, and when can we expect their widespread adoption?

A first principles analysis of Smartglasses


In order to perform classification on a functional brain scan it first
undergoes many preprocessing steps. One of these steps is the transformation from
the timeseries output of the fMRI scan, transforming a $m{\times}n$ (where $ m $ is
the number of timepoints recorded and $ n $ is the number of brain regions used)
to an $n{\times}n$ matrix of similiarity values (called a connectome). This
similarity value is a measure of the neural synchronization between the 2 regions.

So how do we quantify the similarity of 2 different timeseries? This blog post will
explore the common ways of quantifying time series similarity in a neuroscientific
setting. Before we get into the actual methods used to calculate time series similarity,
we need to cover the corner stone of almost all of the methods we are about to
explore - covariance.

## Covariance

**Covariance** is simply a measure of how two random variables change together.
Below is the formula for calculating covariance:

$$
\operatorname{\Sigma} = \operatorname{E}\left[\left(X-\operatorname{E}[X]\right ) \left (X-\operatorname{E}[X]\right)^\mathrm{T}\right]
$$

In the case of fMRI, we have a multivariate random variable, allowing us to use
Maximum Likelihood Estimation to estimate the covariance matrix. Below is a toy
example of our estimated covariance matrix.

$$
\begin{bmatrix}
var(X) & cov(X,Y) & cov(X,Z) \\
cov(X,Y) & var(Y) & cov(X,Y) \\
cov(X,Z) & cov(X,Y) & var(Z)
\end{bmatrix}
$$

This estimated covariance matrix is called the empirical covariance matrix. In
practice we don't use the covariance matrix for our fMRI analyis.
This for a number of reasons:

-   Covariance is bound between $ -\infty $ and $ \infty $ making it less
    suitable for downstream classifiers.
-   Covariance coefficients are not standardized and cannot be used to quantify the
    strength of the relationship.
-   Symmetric Positive semi-definite (SPD) matricies do not naturally form a Euclidean space (this will be
    important later).
-   When the number of features is large relative to the number of observations,
    the sample/empirical covariance matrix has excessive estimation error.

Many of the following approaches will aim to address some or all of the above
drawbacks with empirical covariance.

To address the excessive estimation error, it is common to perform a
transformation to the covariance coefficients, known as **"shrinkage"**. In
their seminal paper, "Honey, I shrunk the Sample Covariance Matrix" \[[1](http://www.ledoit.net/honey.pdf)\], Ledoit \& Wolf
proposed using shrinkage to regularize the sample covariance matrix.
"Shrinkage" as the name implies, pulls the most extreme covariance coefficients
towards more central values. This not only resolves our excessive estimation
error, it can also make the matrix easily invertable by encouraging numerical
stability. (The interested reader should consult the scikit-learn docs
\[[2](https://scikit-learn.org/stable/auto_examples/covariance/plot_covariance_estimation.html)\] )

Now that we have a well conditioned covariance matrix, we can attempt to address
some of the other identified drawbacks.

## Canonical Approaches

**Pearson's correlation coefficient** (Pearson's $ R $) or simply correlation, is
the most commonly used method to quantify similarity between 2 fMRI timeseries.
Correlation is a linear metric computed from the covariance of the 2 timeseries.
Below is the mathematical formula to compute correlation for a pair of random
variables:

$$
{\displaystyle \rho _{X,Y}={\frac {\operatorname {cov} (X,Y)}{\sigma _{X}\sigma _{Y}}}}
$$

Correlation is widely used in neuroscience as it has long statistical history
and is bound between -1 and 1. However, correlation does have some disadvantages.
The below figure should demonstrate one clearly:

![Correlation Decay](/corr.gif "Correlation Decay")

Due to correlations linear nature, the same timeseries being slightly out of
phase causes a huge decrease in the correlation value. Additionally, correlation
provides no distinction between whether 2 regions are directly
connected or indirectly connected via another region. To account for this, we
can use partial correlation!

**Partial correlation** is a variant of PCC that attempts to address
distinguishing between direct and indirect connections. This is done by
computing correlation between regions after regressing all other timeseries.
Partial correlation is computed from the inverse of the covariance matrix
(this is where the shrinkage comes in handy), also known as the precision matrix.
Below is the mathematical formula for computing partial correlation for a pair
of random variables:

$$
\rho = -{\frac {p_{ij}}{\sqrt {p_{ii} p_{jj}}}}
$$

where $p$ is the respective precision matrix (obtained through inversion of our
well conditioned covariance matrix).

## Novel Approaches

**Tangent space parameterization** is a recently proposed solution to the problematic
nature of SPD matricies, in that they don't naturally form Euclidean spaces.
Varoquax et al. \[[3](https://arxiv.org/abs/1008.5070)\] proposed a solution to this
sticky problem. Whilst I won't profess to be an expert on differential
manifolds or Riemmenian geometry, here is how I understand it.

Any covariance matrix is a symmetric positive semi-definite matrix. In order to
preserve the geometry during our transformations, we need to
perform a projection into a space where our geometry is preserved.
As SPD matricies form a diffentiable manifold, we can use this fact and
Riemannian geometry to keep our geometric intutions in tact.

![Manifold](/manifold.jpg "Manifold")

<a href="https://juliahub.com/docs/PosDefManifold/vb3YZ/0.4.8/introToRiemannianGeometry/#Intro-to-Riemannian-Geometry-1" styles="text-align: center;">
Image source
</a>

The above image makes it clear that if we perform Euclidean operations like
subtraction in our manifold, the distances and therefore the geometry will break
down. By defining a homomorphic (ie structure preserving) tangent space, we can
approximate the distance on the manifold using Euclidean distance in the tangent
space!

But we aren't quite out of the woods yet, since we still need to define the
reference point $ G $ on the manifold. As all of our covariance matricies need
to be projected into the same tangent space, we need to choose an appropriate
reference point. To do this, Varoquax et al. compute the geometric mean of the
matricies to define a reference space.

In their paper, Dadi et al. \[[4](https://hal.inria.fr/hal-01824205v3)\] showed
that across many fMRI classification tasks, tangent space parameterization
should be preferred over both correlation and partial correlation.

### Dynamic Time Warping

**Dynamic Time Warping** (DTW) is not like the previously discussed methods, as it does
not rely on covariance matricies. DTW is one of the canonical algorithms used for measuring
similarity between timeseries. It is commonly used in scenarios where timeseries
may have a different length or be temporally offest. DTW is a non linear metric,
and calculates the optimal match between two given sequences with certain
restrictions.

DTW was first proposed for usage with fMRI timeseries by Meszlényi et al. \[[5](https://www.frontiersin.org/articles/10.3389/fnins.2017.00075/full)\].
They propose its value for fMRI timeseries analysis as it has been noted that
"dynamic switch of brain states" can cause non-stationary time lags. If there
are known or unknown mechanisms in the brain that cause regional time differences
in hemodynamic response - DTW could be a great way to account for this.

Below is a Rust implementation of windowed DTW.
The window refers to the maximal warping distance.

```rust
let m = s.len() + 1;
let n = t.len() + 1;
let mut dtw = Array::from_elem((m, n), f64::MAX);

dtw[[0, 0]] = 0.;

let max_window = i32::max(*window, i32::abs((n - m) as i32));
for si in 1..n {
    let lower_bound = i32::max(1, si as i32 - max_window);
    let upper_bound = i32::min(m as i32, si as i32 + max_window);
    for ti in lower_bound as usize..upper_bound as usize {
        let cost = distance_fn(&s[si - 1], &t[ti - 1]);
        dtw[[si, ti]] = cost + 
            f64::min(
                f64::min(dtw[[si - 1, ti]], dtw[[si, ti - 1]]), 
                dtw[[si - 1, ti - 1]]
            );
    }
}
```

Dynamic Time Warping does have some disadvantages. One particularly
troublesome one being the computational complexity - $ O(n^2) $. The paper from
Meszlényi et al. does come with a DTW implementation for fMRI data on GitHub,
however it is extremely cumbersome and difficult to get running (particularly
on HPC systems where you don't have root access.)

To address this I wrote RustDTW, a Python module backed by Rust code that works
natively within Python on numpy matricies. Check out an example of using DTW
for fMRI classification [here.](https://github.com/FL33TW00D/rustDTW/blob/master/examples/classification/ABIDE_classification.ipynb)
My basic testing shows that it should be nearly 10x faster than a mutlithreaded
Python implementation.

## Classification Performance

Whilst not meant to be an exhaustive exploration of the different connectivity
metrics explored above (and itself is not an exhaustive exploration of available
connectivity metrics, for more please see Pervaiz et al. \[[6](https://www.biorxiv.org/content/10.1101/741595v1.full.pdf)\]
I thought I would demonstrate a comparison between them using a large open source dataset.
You can check out the inner workings and run it for yourself on Colab [here](https://colab.research.google.com/github/FL33TW00D/rustDTW/blob/master/examples/classification/ABIDE_classification.ipynb)!

For those just interested in the results, the below graph summarizes the
performance of each of the similarity metrics.


![Accuracy](/accuracy.png "Accuracy")

The results are pretty interesting! We can see that both non-linear matching and
transforming the covariance matricies into a natural space can provide improvements
over standard correlation. Perhaps a future metric will incorporate
characteristics from both tangent space parameterization and DTW to improve
performance.


A key step in diagnosing pathology using functional magnetic resonance imaging (fMRI) data is computing the similarity of timeseries, as this is what determines if 2 brain regions are "connected" or not. This post explores different methods of computing timeseries similarity

Essays