# Unitary Howe dualities for fermionic and bosonic algebras and related Dirac
operators

Guner Muarem  Campus Middelheim
Middelheimlaan 1
M.G.221
2020 Antwerpen
België<EMAIL_ADDRESS>
###### Abstract.

In this paper we use the canonical complex structure $\mathbb{J}$ on
$\mathbb{R}^{2n}$ to introduce a twist of the symplectic Dirac operator. As a
matter of fact, these operators can be interpreted as the bosonic analogues of
the Dirac operators on a Hermitian manifold. Moreover, we prove that the
algebra of these symplectic Dirac operators is isomorphic to the Lie algebra
$\mathfrak{su}(1,2)$ which leads to the Howe dual pair
$(\mathsf{U}(n),\mathfrak{su}(1,2))$.

##  1\. Introduction

The CCR (canonical commuting relation) and CAR algebras (canonical
anticommuting relation) are fundamental algebras in theoretical physics used
for the study of bosons and fermions. From a mathematical viewpoint, these
algebras are named the Weyl algebra (or symplectic Clifford algebra) and
Clifford algebra. These algebras can be constructed in a very analogous way.
The Clifford algebra is constructed on a vector space $V$ equipped with a
symmetric bilinear form $B$, whereas the Weyl algebra requires an even
dimensional vector space equipped with a skew-symmetric bilinear form (or
symplectic form) $\omega$. In both cases, one then constructs the tensor
algebra $T(V)$ where an ideal $I(V)$ is divided out. In the orthogonal
setting, this is the ideal $I_{B}(V)$ with elements subject to the relation
$\\{u,v\\}=2B(u,v)$. In the symplectic setting, this is the ideal
$I_{\omega}(V)$ generated by $[u,v]=-\omega(u,v)$.

There is, however, a fundamental difference: the Clifford algebra is finite-
dimensional, whereas the Weyl algebra is infinite-dimensional. For the spinors
(orthogonal versus symplectic) the same infinite-dimensional principle holds
as for the Clifford algebras. As a matter of fact, the symplectic spinors are
the smooth vectors in the metaplectic representation [2]. Using the generators
of the Clifford (resp. Weyl) algebra, one can associate a natural first order
spin (resp. metaplectic) invariant differential operator by contracting the
Clifford algebra elements using the bilinear form $B$ (resp. the symplectic
form $\omega)$ with derivatives. This gives rise to the Dirac operator
$\partial_{x}=\sum_{k=1}^{n}e_{k}\partial_{x_{k}}$ where
$\\{e_{j},e_{k}\\}=-2\delta_{ij}$ and the symplectic Dirac operator
$\sum_{k=1}^{n}\left(iq_{k}\partial_{y_{k}}-\partial_{q_{k}}\partial_{x_{k}}\right)$
where $[\partial_{q_{j}},iq_{k}]=i\delta_{jk}$ are the Heisenberg relations.

The theory which studies the solutions of the Dirac operator is known as
Clifford analysis and can be seen as a hypercomplex function theory. Moreover,
quite some generalisations have occurred in the last two decades. This
involves e.g. Clifford analysis on superspace and Clifford analysis on
(hyper)Kähler spaces (see for instance [3]). It is in the latter framework in
which this paper is situated, but then from a symplectic point of view. More
precisely, we provide the foundations of what we will call a hermitian variant
of symplectic Clifford analysis, where we incorporate the additional datum of
a compatible complex structure $\mathbb{J}$ on the flat symplectic space
$\mathbb{R}^{2n}$. This leads to the study of the symplectic Dirac operator on
a Kähler manifold as was already initiated in [2, 4]. However, the underlying
invariance symmetry and the algebra generated by these type of operators (and
their duals) was never investigated.

##  2\. Rudiments of symplectic Clifford analysis

Let us consider the canonical symplectic space $\mathbb{R}^{2n}$ with
coordinates $({x},{y})$ and the usual symplectic form
$\omega_{0}=\sum_{j=1}^{n}dx_{j}\wedge dy_{j}$ which has the matrix
representation

$\Omega_{0}=\begin{pmatrix}0&I_{n}\\\ -I_{n}&0\end{pmatrix}.$

Recall that the symplectic group $\mathsf{Sp}(2n,\mathbb{R})$ is the group
given by invertible linear transformations preserving the non-degenerate skew-
symmetric bilinear form from above and is given (in terms of matrices) by

$\mathsf{Sp}(2n,\mathbb{R})=\\{M\in\mathsf{GL}(2n,\mathbb{R})\mid
M^{T}\Omega_{0}M=\Omega_{0}\\}.$

The group is non-compact and has dimension $2n^{2}+n$. Moreover, the
corresponding Lie algebra is denoted by $\mathfrak{sp}(2n,\mathbb{R})$. The
main difference with the orthogonal case, lies in the fact that the
metaplectic group (the double cover of the symplectic group) does not admit a
finite dimensional representation (it is not a matrix group). This is a strong
contrast with the spin representation in the orthogonal case. Moreover, the
orthogonal spinors $\mathbb{S}$ are realised as a idempotent left ideal in the
Clifford algebra, which is not the case for the symplectic spinors. As
mentioned, the symplectic equivalent of the spin representations are infinite
dimensional, which means that one needs to work with the theory of unitary
representations.

###  2.1. The Schwartz space and metaplectic representation

For further convenience, we fix notation and define the Schwartz space, which
plays a crucial role in the construction of the metaplectic representation. On
the space $\mathcal{C}^{\infty}(\mathbb{R}^{n},\mathbb{C})$ we define (using
the multi-index notation) the norm

$||f||_{\alpha,\beta}:=\sup_{{q}\in\mathbb{R}^{n}}|{q}^{\alpha}(D^{\beta}f)({q})|$

for all $\alpha,\beta\in\mathbb{N}^{n}$. The Schwartz space
$\mathcal{S}(\mathbb{R}^{n})$ is the subspace of $L^{p}(\mathbb{R}^{n})$ (for
$1\leq p\leq\infty$) consisting of rapidly decreasing functions and is given
by

$\mathcal{S}(\mathbb{R}^{n},\mathbb{C}):=\\{f\in\mathcal{C}^{\infty}(\mathbb{R}^{n},\mathbb{C}):||f||_{\alpha,\beta}<\infty\text{
for all }\alpha,\beta\in\mathbb{N}^{n}\\}.$

We now describe (following [2]) the infinite-dimensional Segal-Shale-Weil
representation (also oscillator or metaplectic representation) of the
metaplectic group. The smooth vectors of the unitary representation
$\mathfrak{m}:\mathsf{Mp}(2n)\to\mathsf{U}(L^{2}(\mathbb{R}^{n}))$ coincide
with the Schwartz space $\mathcal{S}(\mathbb{R}^{n})$ and are a model for the
symplectic spinors $\mathbb{S}^{\infty}$. Due to Stone-Von Neumann theorem the
representation is unique (up to unitary equivalence).

###  2.2. The symplectic Clifford algebra and the related Dirac operator

Let $(V,\omega)$ be a symplectic vector space. The symplectic Clifford algebra
$\mathsf{Cl}_{s}(V,\omega)$ is defined as the quotient algebra of the tensor
algebra $T(V)$ of $V$, by the two-sided ideal

$\mathcal{I}_{\omega}:=\\{v\otimes u-u\otimes v+\omega(v,u):u,v\in V\\}.$

In other words $\mathsf{Cl}_{s}(V,\omega):=T(V)/\mathcal{I}_{\omega}$ is the
algebra generated by $V$ in terms of the relation $[v,u]=-\omega(v,u)$, where
we have omitted the tensor product symbols. We refer to the symplectic
Clifford algebra on $\mathbb{R}^{2n}$ as the $n$th Weyl algebra
$\mathcal{W}_{n}$ with generators
$iq_{1},\dots,iq_{n},\partial_{q_{1}},\dots,\partial_{q_{n}}$ satisfying the
commutation relations $[q_{j},q_{k}]=0$ and
$[\partial_{q_{j}},q_{k}]=\delta_{jk}$.

Denote by $\mathcal{F}$ a suitable function space (e.g. the space of
polynomials, or smooth funtions).

######  Definition 2.1.

The symplectic Dirac operator on $(\mathbb{R}^{2n},\omega_{0})$ is the first-
order (in the base variables ${x}$ and ${y}$) differential operator acting on
a symplectic spinor-valued functions space
$\mathcal{F}\otimes\mathbb{S}^{\infty}$ given by
$D_{s}=\sum_{j=1}^{n}(iq_{j}\partial_{y_{j}}-\partial_{q_{j}}\partial_{x_{j}}).$
With respect to the symplectic Fischer inner product (see [5]), we obtain the
dual operator $X_{s}=\sum_{j=1}^{n}(iq_{j}{x_{j}}+\partial_{q_{j}}{y_{j}}).$

These operators satisfy the relations:

$\displaystyle[\mathbb{E}+n,X_{s}]=X_{s},$
$\displaystyle[\mathbb{E}+n,D_{s}]=-D_{s},$
$\displaystyle[D_{s},X_{s}]=-i(\mathbb{E}+n)$

where $\mathbb{E}=\sum_{j=1}^{n}(x_{j}\partial_{x_{j}}+y_{j}\partial_{y_{j}})$
is the Euler operator. In other words, the three operators give rise to a copy
of the Lie algebra $\mathfrak{sl}(2)$.

##  3\. The interaction with a complex structure

###  3.1. Definition of the twisted symplectic Dirac operators

We will now introduce a complex structure $\mathbb{J}$ on the symplectic
manifold $(\mathbb{R}^{2n},\omega_{0})$ which is compatible with the
symplectic form $\omega_{0}$. This means that $\omega_{0}(x,\mathbb{J}y)$
defines a Riemannian metric $g$. Otherwise said, we will be working with the
canonical Kähler manifold $(\mathbb{R}^{2n},\omega_{0},g,\mathbb{J})$. By
Darboux’s theorem, we obtain, with respect to the canonical symplectic basis
$\\{e_{j}\\}_{j=1}^{2n}$ the following complex structure
$\mathbb{J}=\left(\begin{smallmatrix}0&-I_{n}\\\
I_{n}&0\end{smallmatrix}\right).$ The action of the complex structure
$\mathbb{J}$ on $\mathbb{R}^{2n}$ is given by

$\displaystyle(x_{1},\dots,x_{n},y_{1},\dots,y_{n})\mapsto(y_{1},\dots,y_{n},-x_{1},\dots,-x_{n}).$

######  Definition 3.1.

The new differential operators acting on symplectic spinor-valued functions

$\displaystyle\tilde{D}_{s}$
$\displaystyle=\sum_{j=1}^{n}iq_{j}\partial_{x_{j}}+\partial_{y_{j}}\partial_{q_{j}}$
$\displaystyle\tilde{X}_{s}$
$\displaystyle=\sum_{j=1}^{n}x_{j}\partial_{q_{j}}-iy_{j}q_{j}$
$\displaystyle\mathbb{E}$
$\displaystyle=\sum_{j=1}^{n}x_{j}\partial_{x_{j}}+y_{j}\partial_{y_{j}}$

also give rise to a copy of the Lie algebra $\mathfrak{sl}(2)$. We call these
first two operators the twists of $D_{s}$ and $X_{s}$.

Both sets of operators, i.e. $(D_{s},X_{s})$ and
$(\tilde{D}_{s},\tilde{X}_{s})$, are symplectic invariant, albeit under the
following two different realisations of the symplectic Lie algebra given by:

$\displaystyle\begin{cases}X_{jk}=x_{j}\partial_{x_{k}}-y_{k}\partial_{y_{j}}-(q_{k}\partial_{q_{j}}+\frac{1}{2}\delta_{jk})&1\leq
j\leq k\leq n\\\
Y_{jk}=x_{j}\partial_{y_{k}}+x_{k}\partial_{y_{j}}+i\partial_{q_{j}}\partial_{q_{k}}&j<k=1,\dots,n\\\
Z_{jk}=y_{j}\partial_{x_{k}}+y_{k}\partial_{x_{j}}+iq_{j}q_{k}&j<k=1,\dots,n\\\
Y_{jj}=x_{j}\partial_{y_{j}}+\frac{i}{2}\partial_{q_{j}}^{2}&j=1,\dots,n\\\
Z_{jj}=y_{j}\partial_{x_{j}}+\frac{i}{2}q_{j}^{2}&j=1,\dots,n\end{cases}$

and

$\displaystyle\begin{cases}\tilde{X}_{jk}=x_{j}\partial_{x_{k}}-y_{k}\partial_{y_{j}}+q_{k}\partial_{q_{j}}+\frac{1}{2}\delta_{jk}&1\leq
j\leq k\leq n\\\
\tilde{Y}_{jk}=x_{j}\partial_{y_{k}}+x_{k}\partial_{y_{j}}-iq_{j}q_{k}&j<k=1,\dots,n\\\
\tilde{Z}_{jk}=y_{j}\partial_{x_{k}}+y_{k}\partial_{x_{j}}-i\partial_{q_{j}}\partial_{q_{k}}&j<k=1,\dots,n\\\
\tilde{Y}_{jj}=x_{j}\partial_{y_{j}}-\frac{i}{2}q_{j}^{2}&j=1,\dots,n\\\
\tilde{Z}_{jj}=y_{j}\partial_{x_{j}}-\frac{i}{2}\partial_{q_{j}}^{2}&j=1,\dots,n\end{cases}$

Of course, it is not very useful that $D_{s}$ and $\widetilde{D}_{s}$ are
invariant under different (yet isomorphic)
$\mathfrak{sp}(2n,\mathbb{R})$-realisations. Therefore, we will perform a
symmetry reduction so that both operators become invariant under one and the
same Lie algebra. To that end, we need to find the symplectic matrices which
commute with the complex structure.

######  Lemma 3.2.

We have that

$\mathsf{Sp}_{\mathbb{J}}(2n,\mathbb{R}):=\\{M\in\mathsf{Sp}(2n,\mathbb{R})\mid
M\mathbb{J}=\mathbb{J}M\\}$

defines a realisation for the unitary Lie group.

###### Proof.

In order to see this, assume that $M$ is of the block-form:
$\left(\begin{smallmatrix}A&B\\\ C&D\end{smallmatrix}\right),$ where $A,B,C$
and $D$ are $(n\times n)-$matrices. The condition that $M$ is symplectic is
equivalent to the one of the following conditions: the matrices $A^{T}C$ and
$B^{T}D$ are symmetric and $A^{T}D-C^{T}B=I$. So, in order to determine
$\operatorname{Sp}_{\mathbb{J}}(2n,\mathbb{R})$ we need to determine the
symplectic matrices $M$ which commute with the complex structure $\mathbb{J}$.
The latter conditions means that

$\displaystyle M\mathbb{J}=\mathbb{J}M$
$\displaystyle\iff\mathbb{J}^{-1}M\mathbb{J}=M$
$\displaystyle\iff\begin{pmatrix}0&I\\\ -I&0\end{pmatrix}\begin{pmatrix}A&B\\\
C&D\end{pmatrix}\begin{pmatrix}0&-I\\\ I&0\end{pmatrix}=\begin{pmatrix}A&B\\\
C&D\end{pmatrix}$ $\displaystyle\iff\begin{pmatrix}D&-C\\\
-B&A\end{pmatrix}=\begin{pmatrix}A&B\\\ C&D\end{pmatrix}$

This implies $A=D$ and $B=-C$. In other words the matrix $M$ is of the form:
$M=\left(\begin{smallmatrix}A&B\\\ -B&A\end{smallmatrix}\right).$ Next, we
still have the condition that $M$ is symplectic, i.e.

$\displaystyle M^{T}\Omega M=\Omega$
$\displaystyle\iff\begin{pmatrix}A^{T}&C^{T}\\\
B^{T}&D^{T}\end{pmatrix}\begin{pmatrix}0&I\\\
-I&0\end{pmatrix}\begin{pmatrix}A&B\\\ C&D\end{pmatrix}=\begin{pmatrix}0&I\\\
-I&0\end{pmatrix}$ $\displaystyle\iff\begin{pmatrix}-C^{T}&A^{T}\\\
-D^{T}&B^{T}\end{pmatrix}\begin{pmatrix}A&B\\\
C&D\end{pmatrix}=\begin{pmatrix}-C^{T}A+A^{T}C&-C^{T}B+A^{T}D\\\
-D^{T}A+B^{T}C&-D^{T}B+B^{T}D\end{pmatrix}=\begin{pmatrix}0&I\\\
-I&0\end{pmatrix}$ $\displaystyle\iff\begin{cases}A^{T}C=C^{T}A\\\
A^{T}D-C^{T}B=I\\\ B^{T}C-D^{T}A=-I\\\ B^{T}D=D^{T}B\end{cases}$

This means that $A^{T}C$ and $B^{T}D$ should be symmetric matrices and
$A^{T}D-C^{T}B=I$. But now, due to the first condition this reduces to
$B^{T}A=A^{T}B$ and $A^{T}A+B^{T}B=I$. In other words, the matrices we are
looking for must be of the form $M=\left(\begin{smallmatrix}A&B\\\
-B&A\end{smallmatrix}\right)$ with $B^{T}A=A^{T}B$ and $A^{T}A+B^{T}B=I$, i.e.

$\displaystyle\begin{pmatrix}A^{T}&-B^{T}\\\
B^{T}&A^{T}\end{pmatrix}\begin{pmatrix}A&B\\\
-B&A\end{pmatrix}=\begin{pmatrix}A^{T}A+B^{T}B&A^{T}B-B^{T}A\\\
B^{T}A-A^{T}B&B^{T}B+A^{T}A\end{pmatrix}=\begin{pmatrix}I&0\\\
0&I\end{pmatrix}$

Which is exactly the condition for a unitary matrix. The map
$\Phi:\mathsf{Sp}_{\mathbb{J}}(2n,\mathbb{R})\to\mathsf{U}(n):M\mapsto A+iB$
gives the wanted isomorphism. ∎

###  3.2. Unitary invariant symplectic Dirac operators

One can now check that the symplectic Dirac operator and its twist, are
unitary invariant differential operators. This can be done by verifying that
the operators commute with the following realisation of unitary Lie algebra
$\mathfrak{u}(n)$:

######  Lemma 3.3.

We have the following operator realisation of the Lie algebra
$\mathfrak{u}(n)$:

$\displaystyle\begin{cases}A_{jk}=y_{j}\partial_{x_{k}}+y_{k}\partial_{x_{j}}-x_{j}\partial_{y_{k}}-x_{k}\partial_{y_{j}}+i(q_{j}q_{k}-\partial_{q_{j}}\partial_{q_{k}})&\quad
1\leq j<k\leq n\\\
B_{jj}=y_{j}\partial_{x_{j}}-x_{j}\partial_{y_{j}}+\frac{i}{2}\left(q_{j}^{2}-\partial_{q_{j}}^{2}\right)&\quad
1\leq j\leq n\\\
C_{jk}=x_{j}\partial_{x_{k}}-x_{k}\partial_{x_{j}}+y_{j}\partial_{y_{k}}-y_{k}\partial_{y_{j}}+q_{j}\partial_{q_{k}}-q_{k}\partial_{q_{j}}&\quad
1\leq j<k\leq n\end{cases}$

This means that we can refine the $\mathfrak{sp}(2n)$-invariant PDE $D_{s}f=0$
into two $\mathfrak{u}(n)$-invariant PDEs given by $D_{s}f=0$ and
$\widetilde{D}_{s}f=0$, for a symplectic spinor valued polynomial
$f\in\mathcal{P}(\mathbb{R}^{2n},\mathbb{C})\otimes\mathcal{S}(\mathbb{R}^{n})$.
In analogy with the orthogonal case, we call the solutions hermitian
symplectic monogenics (or $h$-symplectic monogenics in short).

###  3.3. Symplectic Dolbeault operators

Moreover, there is a second way of introducing the twist of the symplectic
Dirac operator. Let us define the following operators which are known in the
literature as the symplectic Dolbeault operators [4] defined by means of

$D_{z}=\frac{D_{s}+i\widetilde{D}_{s}}{2}$

and

$D_{z}^{\dagger}:=\frac{D_{s}-i\widetilde{D}_{s}}{2}.$

One easily verifies that

$\displaystyle\frac{1}{2}(D_{s}+i\widetilde{D}_{s})$
$\displaystyle=-\sum_{j=1}^{n}\mathfrak{F}_{j}\partial_{z_{j}}\quad$
$\displaystyle\ \frac{1}{2}(D_{s}-i\widetilde{D}_{s})$
$\displaystyle=\sum_{j=1}^{n}\mathfrak{F}_{j}^{\dagger}\partial_{\overline{z}_{j}}.$

where we have introduced the symbols
$\mathfrak{F}_{j}=(q_{j}+\partial_{q_{j}})$ and
$\mathfrak{F}_{j}^{\dagger}=(q_{j}-\partial_{q_{j}})$. The structure of the
operators $D_{z}$ and $D_{z}^{\dagger}$ is similar to the orthogonal case.
However, the raising/lowering operators $\mathfrak{F}_{j}$ and
$\mathfrak{F}_{j}^{\dagger}$ are used instead of isotropic Witt vectors
$\mathfrak{f}_{j}$ and $\mathfrak{f}_{j}^{\dagger}$ (see for instance [6] and
the references therein).

###  3.4. Class of simultaneous solutions of $D_{s}$ and $\widetilde{D}_{s}$

We will now describe a wide class of examples of $h$-symplectic monogenics, by
making the link with holomorphic functions in several variables. Let
$f:\Omega\subset\mathbb{C}^{n}\to\mathbb{C}$ be a complex-valued function in
several complex variables which is of the class $\mathcal{C}^{1}(\Omega)$
(i.e. continuously differentiable). We say that $f$ is holomorphic (in several
variables) if $\partial_{\overline{z}_{j}}f(z_{1},\dots,z_{n})=0$ for all
$1\leq j\leq n$, where
$\partial_{\overline{z}_{j}}:=\frac{1}{2}(\partial_{x_{j}}+\partial_{y_{j}})$
is the Cauchy-Riemann operator in the relevant variable. Moreover, we denote
the set of holomorphic functions in $\Omega$ by $\operatorname{Hol}(\Omega)$.

In order not to overload notations, we use the summation convention. Suppose
that we have a function of the form
$F({x},{y},{q})=e^{-\frac{1}{2}|{q}|^{2}}H({x},{y})$. Letting the symplectic
Dirac operator act on $F$ gives:

$\displaystyle D_{s}\left(e^{-\frac{1}{2}|{q}|^{2}}H({x},{y})\right)$
$\displaystyle=(iq_{k}\partial_{y_{k}}-\partial_{x_{k}}\partial_{q_{k}})\left(e^{-\frac{1}{2}|{q}|^{2}}H({x},{y})\right)$
$\displaystyle=iq_{k}e^{-\frac{1}{2}|{q}|^{2}}\partial_{y_{k}}H({x},{y})+e^{-\frac{1}{2}|{q}|^{2}}q_{k}\partial_{x_{k}}H({x},{y})$
$\displaystyle=e^{-\frac{1}{2}|{q}|^{2}}q_{k}(\partial_{x_{k}}+i\partial_{y_{k}})H({x},{y})$

We note that this equals zero if
$(\partial_{x_{k}}+i\partial_{y_{k}})H({x},{y})=0$ for all $k=1,\dots,n$, i.e.
if $H({x},{y})$ is a holomorphic function in several variables. Completely
similar,

$\displaystyle\widetilde{D_{s}}\left(e^{-\frac{1}{2}|{q}|^{2}}H({x},{y})\right)$
$\displaystyle=(iq_{k}\partial_{x_{k}}+\partial_{y_{k}}\partial_{q_{k}})\left(e^{-\frac{1}{2}|{q}|^{2}}H({x},{y})\right)$
$\displaystyle=e^{-\frac{1}{2}|{q}|^{2}}q_{k}(i\partial_{x_{k}}-\partial_{y_{k}})H({x},{y})$
$\displaystyle=ie^{-\frac{1}{2}|{q}|^{2}}q_{k}(\partial_{x_{k}}+i\partial_{y_{k}})H({x},{y}),$

which is zero for holomorphic $H({x},{y})$.

This means that every function of the form
$e^{-\frac{1}{2}(q_{1}^{2}+\dots+q_{n}^{2})}H({x},{y})$ with $H$ an
holomorphic function is several variables is a solution of both $D_{s}$ and
$\widetilde{D}_{s}$. This observation generalises the class of solutions
obtained by Habermann in [2] for $n=1$. As a matter of fact, it turned out
that there are much more solutions than the ones of this form. In order to
describe these systematically, we will need the notion of Howe dualities and
corresponding Fischer decompositions. This will be done in full detail in our
upcoming paper [7]. In the following section we will reveal the algebraic
structures required for this approach.

##  4\. A unitary Howe duality associated with $D_{s}$ and $\widetilde{D}_{s}$

Recall that the Lie algebra $\mathfrak{su}(1,2)$ is a quasi-split real form of
the complex Lie algebra $\mathfrak{sl}(3)$ and is defined in terms of matrices
as

$\displaystyle\mathfrak{su}(1,2)=\left\\{\begin{pmatrix}\alpha&\beta&ic\\\
\gamma&\overline{\alpha}-\alpha&-\overline{\beta}\\\
id&-\overline{\gamma}&-\overline{\alpha}.\end{pmatrix}\mid c,d\in\mathbb{R}\
\&\ \alpha,\beta,\gamma\in\mathbb{C}\right\\}$

Some calculations now lead to observation that the Lie algebra generated by
the symplectic Dirac operators $D_{s},\widetilde{D}_{s}$ and their duals
$X_{s},\widetilde{X_{s}}$ gives rise to a copy of the Lie algebra
$\mathfrak{su}(1,2)$. In order to close the algebra, we introduce the
following differential operators:

$\displaystyle\mathcal{O}:=\sum_{j=1}^{n}i(x_{j}\partial_{y_{j}}-y_{j}\partial_{x_{j}})+\partial_{q_{j}}^{2}-q_{j}^{2},$
$\displaystyle\Delta:=\sum_{j=1}^{n}\partial_{x_{j}}^{2}+\partial_{y_{j}}^{2}$
and $\displaystyle r^{2}:=\sum_{j=1}^{n}{x_{j}}^{2}+{y_{j}}^{2}.$

Together with symplectic Dirac operators $D_{s},\widetilde{D}_{s}$ and their
duals $X_{s},\widetilde{X}_{s}$ these operators satisfy the commutator
relations of $\mathfrak{su}(1,2)$. Moreover, their commutators are given by:

$[\cdot,\cdot]$ | $D_{s}$ | $\widetilde{D}_{s}$ | $\Delta$ | $X_{s}$ | $\mathbb{E}$ | $\mathcal{O}$ | $\widetilde{X}_{s}$ | $r^{2}$
---|---|---|---|---|---|---|---|---
$D_{s}$ | 0 | $\Delta$ | 0 | -$\mathbb{E}$ | $D_{s}$ | $-3\widetilde{D}_{s}$ | $-\mathcal{O}$ | $-2\widetilde{X}_{s}$
$\widetilde{D}_{s}$ | $-\Delta$ | 0 | 0 | $\mathcal{O}$ | $\widetilde{D}_{s}$ | $3D_{s}$ | $-\mathbb{E}$ | $2X_{s}$
$\Delta$ | 0 | 0 | 0 | $\widetilde{D}_{s}$ | $\Delta$ | 0 | $-2{D}_{s}$ | $\mathbb{E}$
$X_{s}$ | $\mathbb{E}$ | $-\mathcal{O}$ | $-\widetilde{D}_{s}$ | 0 | $X_{s}$ | $-3\widetilde{X}_{s}$ | $-r^{2}$ | 0
$\mathbb{E}$ | $-D_{s}$ | $-\widetilde{D}_{s}$ | $-\Delta$ | $-X_{s}$ | 0 | 0 | $-\widetilde{X}_{s}$ | 0
$\mathcal{O}$ | $3\widetilde{D}_{s}$ | $-3D_{s}$ | 0 | $3\widetilde{X}_{s}$ | 0 | 0 | $-3X_{s}$ | 0
$\widetilde{X}_{s}$ | $\mathcal{O}$ | $\mathbb{E}$ | $2D_{s}$ | $r^{2}$ | $\widetilde{X}_{s}$ | $-3X_{s}$ | 0 | 0
$r^{2}$ | $2\widetilde{X}_{s}$ | $-2X_{s}$ | $-\mathbb{E}$ | 0 | $-2r^{2}$ | 0 | 0 | 0

  1. (1)

We have two copies of the Heisenberg algebra:

$\operatorname{Alg}\\{D_{s},\tilde{D}_{s},\Delta\\}\cong\operatorname{Alg}\\{X_{s},\widetilde{X}_{s},r^{2}\\}\cong\mathfrak{h}_{3}.$

  2. (2)

We have two copies of the Lie algebra $\mathfrak{sl}(2)$

$\operatorname{Alg}\\{D,D^{\dagger},\mathbb{E}\\}\cong\operatorname{Alg}\\{\tilde{D},\tilde{D}^{\dagger},\mathbb{E}\\}\cong\operatorname{Alg}\\{\Delta,N,\mathbb{E}\\}\cong\mathfrak{sl}(2).$

This means that there is a canonical $\mathfrak{su}(1,2)$-action on the space
of spinor valued polynomials
$\mathcal{P}(\mathbb{R}^{2n},\mathbb{C})\otimes\mathcal{S}(\mathbb{R}^{n})$
where restricting to the subalgebra $\operatorname{Alg}(D_{s},X_{s})$
corresponds to $\mathfrak{sl}(2)$-copy obtained in the Howe duality for
symplectic Clifford analysis (see [5] for more details). Now, taking into
account the symplectic Dirac operators and its twists, we obtain the dual pair
$\mathsf{U}(n)\times\mathfrak{su}(1,2)$ (i.e. the underlying group of
invariance, together with the algebra generated by the operators and their
duals).

We now focus on the reduction of the symplectic spinor space. In the
orthogonal case, the spinor space $\mathbb{S}$ decomposes under the action of
the unitary group $\mathsf{U}(n)$ as
$\mathbb{S}=\bigoplus_{r}\mathbb{S}_{(r)},$ with $\mathbb{S}_{(r)}$
inequivalent irreducible pieces, which are eigenspaces of the fermionic
quantum oscillator (also called spin-Euler operator, see for instance [6]). In
the symplectic case, the relevant operator for decomposing the infinite
dimensional spinor space is the bosonic quantum oscillator. The hamiltonian of
the quantum oscillator, the so-called Hermite operator, is given by

$\displaystyle\mathcal{H}:\mathcal{S}(\mathbb{R}^{n})\to\mathcal{S}(\mathbb{R}^{n}),\quad
f({q})\mapsto\frac{1}{2}\sum_{j=1}^{n}(\partial_{q_{j}}^{2}-q_{j}^{2})f({q}).$

Note that we can write

$\mathcal{O}=\sum_{j=1}^{n}i(x_{j}\partial_{y_{j}}-y_{j}\partial_{x_{j}})+2\mathcal{H},$

so that the Hermite operator is in fact the spinor-valued part of the operator
$\mathcal{O}$, i.e. the differential operator in $\partial_{q_{j}}$ and the
variables $q_{j}$. Moreover, the eigenspaces can be identified with the
irreducible decomposition of $\mathbb{S}^{\infty}$ into
$\mathfrak{u}(n)$-irreducible representations. This means that the symplectic
spinor space $\mathbb{S}^{\infty}$ decomposes into
$\mathfrak{u}(n)$-irreducible representations ${\widetilde{S}}^{\infty}_{(k)}$
of dimension ${n+k-1\choose k}$ which can be thought of as $k$-homogeneous
polynomials or the eigenspaces of the Hermite operator $\mathcal{H}$.

Moreover, the solutions of the corresponding Dirac operators, called
monogenics, can be introduced from a purely representation theoretical
viewpoint. In general, this boils down to determining the decomposition (this
is called a Fischer decomposition)
$\mathcal{P}_{k}(\mathbb{R}^{m},\mathbb{C})\otimes\mathbf{S}$ where
$\mathbf{S}$ is the spinor space, which is $\mathbb{S}$ in the orthogonal case
and $\mathbb{S}^{\infty}$ in the symplectic case, where we take $m=2n$ in
particular. Moreover, the space of $k$-homogeneous polynomials
$\mathcal{P}_{k}$ coincides with the $k$-symmetric power of the fundamental
representation of resp. the orthogonal or symplectic algebra. We denote by
$\mathcal{M}_{k}$ the $k$-homogeneous solutions of the Dirac operator
$\partial_{x}$, these are called monogenics. They can be defined as follows:

$\displaystyle\mathcal{M}_{k}\leftrightarrow(k,0,\dots,0)\boxtimes\mathbb{S}=(k)\boxtimes\left(\frac{1}{2},\dots,\frac{1}{2}\right)\cong\left(k+\frac{1}{2},\dots,\frac{1}{2}\right),$

where $\boxtimes$ denotes the Cartan product of the
$\mathfrak{so}(m)$-representations. In the symplectic case, we analoguously
obtain:

$\displaystyle\mathcal{M}_{k}^{s}\leftrightarrow(k,0,\dots,0)_{s}\boxtimes\mathbb{S}^{\infty}$
$\displaystyle=(k)_{s}\boxtimes\left(\left(-\frac{1}{2},\dots,-\frac{1}{2}\right)\oplus\left(-\frac{1}{2},\dots,-\frac{3}{2}\right)\right)$
$\displaystyle\cong\left(k-\frac{1}{2},\dots,-\frac{1}{2}\right)\oplus\left(k-\frac{1}{2},\dots,-\frac{3}{2}\right).$

In order to obtain an algebraic characterisation of the space of
$h$-symplectic monogenics, one proceeds as follows. First of all, we note that
we need to consider the symplectic spinors $\mathbb{S}^{\infty}$ from an
unitary viewpoint. We saw that $\mathbb{S}^{\infty}$ decomposes as an infinite
direct sum of finite dimensional $\mathfrak{u}(n)$-modules
${\widetilde{S}}^{\infty}_{(k)}$ which are in fact eigenspaces of the Hermite
operator. We denote the branched spinor space (which is in fact a direct sum
of $\mathfrak{u}(n)$-irreps) by $\widetilde{\mathbb{S}^{\infty}}$. However,
the space of $k$-homogeneous polynomials is not irreducible as a
$\mathfrak{u}(n)$-module and we denote the branched module by
$\widetilde{(k)}$. This means that we are left with the following Cartan
product
$\mathcal{M}_{k}^{hs}\leftrightarrow\widetilde{(k)}\boxtimes\widetilde{\mathbb{S}^{\infty}}$
as a representation theoretical definition of the $h$-symplectic monogenics.
Recall that this are the symplectic spinor-valued polynomial functions
$f\in\mathcal{P}(\mathbb{R}^{2n},\mathbb{C})\otimes\mathcal{S}(\mathbb{R}^{n})$
that satisfy the system of unitary unitary-invariant partial differential
equations

$\begin{cases}D_{s}f&=0\\\ \widetilde{D}_{s}f&=0\end{cases}$

The explicit calculation of the Cartan product (and more generally the tensor
product) will be done in [7]. Moreover, as an application we will prove a
Fischer decomposition for the Howe dual pair we obtained in this paper.

##  5\. Conclusion

In this paper we investigated a new Howe dual pair occurring in symplectic
Clifford analysis by allowing a compatible complex structure. This Howe
duality is of the form $(G,\mathfrak{g}^{\prime})$ where $G$ is the underlying
invariance group for which the relevant Dirac operators are invariant and
$\mathfrak{g}^{\prime}$ is the algebra generated by the Dirac operators and
their duals. Depending on the orthogonal or symplectic framework, we have the
following ‘types’ of Clifford analysis and refinements thereof:

  1. (1)

Orthogonal geometry (giving rise to a Clifford algebra)

    1. (a)

Clifford analysis: $\mathsf{SO}(m)\times\mathfrak{osp}(1|2)$

    2. (b)

Hermitian Clifford analysis $\mathsf{U}(m)\times\mathfrak{osp}(2|2)$

    3. (c)

Quaternionic Clifford analysis $\mathsf{USp}(m)\times\mathfrak{osp}(4|2)$

  2. (2)

Symplectic geometry (giving rise to a Weyl algebra)

    1. (a)

Symplectic Clifford analysis: $\mathsf{Sp}(2m)\times\mathfrak{sl}(2)$

    2. (b)

Hermitian symplectic Clifford analysis:
$\mathsf{U}(m)\times\mathfrak{su}(1,2)$

    3. (c)

Quaternionic symplectic Clifford analysis: $\mathsf{USp}(m)\times$?

Thus far, we extended the framework of hermitian Clifford analysis in the
presence of a symplectic structure in the case of the (flat) Kähler manifold
$\mathbb{R}^{2n}$. It is an interesting question to further reduce the
symmetry to the compact symplectic group $\mathsf{USp(n)}$ so that we have the
chain $\mathsf{Sp}(2n)\supset\mathsf{U}(n)\supset\mathsf{USp}(n)$. In our
furture work [7], we will describe the Fischer decomposition accompanying this
new Howe dual pair.

### Funding information

The author is supported by the FWO-EoS project ‘Symplectic Techniques in
Differential Geometry’ G0H4518N.

## References

  * [1] D. Eelbode, G. Muarem, The orthogonal branching problem for symplectic monogenics, Advances in Applied Clifford Algebras 32, (2022)
  * [2] K. Habermann, L. Habermann, Introduction to symplectic Dirac operators, In Lecture Notes in Mathematics. Springer Berlin Heidelberg (2006)
  * [3] F. Brackx, H. De Schepper, D. Eelbode, R. Lávička, V. Soucek, $\mathfrak{osp}(4|2)$–monogenicity in Clifford analysis, In 15th International Conference on Computational and Mathematical Methods in Science and Engineering (CMMSE 2015) (Vol. 1, pp. 240-243) (2015)
  * [4] E. Korman, Symplectic Dolbeault operators on Kähler manifolds, Annals of Global Analysis and Geometry 44(3), 339-358 (2013)
  * [5] H. De Bie, P. Somberg, V. Soucek, The metaplectic Howe duality and polynomial solutions for the symplectic Dirac operator, Journal of Geometry and Physics 75, 120-128 (2014)
  * [6] D. Eelbode, Stirling numbers and spin-Euler polynomials, Experimental Mathematics, 16(1), 55-66 (2007)
  * [7] Eelbode, D., Muarem G. (2022). A unitary refinement of the symplectic Fischer decomposition. In preparation.


# Capturing cross-session neural population variability through self-
supervised identification of consistent neuron ensembles

Justin Jude
University of Edinburgh
<EMAIL_ADDRESS>&Matthew G. Perich
Icahn School of Medicine at Mount Sinai
New York, NY 10029
<EMAIL_ADDRESS>
&Lee E. Miller
Feinberg School of Medicine
Northwestern
Chicago, IL 60611
<EMAIL_ADDRESS>
&Matthias H. Hennig
University of Edinburgh
<EMAIL_ADDRESS>
###### Abstract

Decoding stimuli or behaviour from recorded neural activity is a common
approach to interrogate brain function in research, and an essential part of
brain-computer and brain-machine interfaces. Reliable decoding even from small
neural populations is possible because high dimensional neural population
activity typically occupies low dimensional manifolds that are discoverable
with suitable latent variable models. Over time however, drifts in activity of
individual neurons and instabilities in neural recording devices can be
substantial, making stable decoding over days and weeks impractical. While
this drift cannot be predicted on an individual neuron level, population level
variations over consecutive recording sessions such as differing sets of
neurons and varying permutations of consistent neurons in recorded data may be
learnable when the underlying manifold is stable over time. Classification of
consistent versus unfamiliar neurons across sessions and accounting for
deviations in the order of consistent recording neurons in recording datasets
over sessions of recordings may then maintain decoding performance. In this
work we show that self-supervised training of a deep neural network can be
used to compensate for this inter-session variability. As a result, a
sequential autoencoding model can maintain state-of-the-art behaviour decoding
performance for completely unseen recording sessions several days into the
future. Our approach only requires a single recording session for training the
model, and is a step towards reliable, recalibration-free brain computer
interfaces.

##  1 Introduction

Neural decoders require stable neurons in a recorded population in order to
accurately predict behaviour such as movement or to allow decoding of stimuli.
However, over time instabilities in the recording equipment and drift in
neural activity lead to instabilities that prevent re-using a decoder trained
on one day for a session recorded on another day (Huber et al., 2012; Ziv et
al., 2013; Driscoll et al., 2017). At the same time, neural population
activity is highly structured and often confined to low-dimensional manifolds
(Cunningham and Byron, 2014) that can be recovered using latent variable
modelling approaches (Hurwitz et al., 2021). Importantly, recent work showed
that movement-related latent neural dynamics in population activity from the
primate motor cortex is stable and could be recovered over intervals as long
as two years (Gallego et al., 2020). This suggests that despite the
variability at the level of single neurons, in each session a subset of
neurons will remain informative about behaviour. A stable cross-session
decoder therefore has to be able to identify these neurons and utilise them
for decoding. Therefore, here we focus on identifying known recording neurons
in unseen sessions. In particular, we hypothesised that a latent encoding of
neural activity can be augmented by information about which neurons were seen
during training, and at which position in the input. We show that this is
sufficient to decode behaviour (in our case different cued arm movements by a
monkey with simultaneous motor cortex recordings) with high accuracy across
unseen sessions.

We achieve this with a self-supervised approach through training a recurrent
neural network (RNN) to predict original neuron positions following data
perturbation in a manner mirroring session to session variability. In essence,
the closer our perturbations mimic real inter-session variability (as shown in
Figure 1), the higher our behaviour prediction performance on an unseen
session. These perturbations include adding spikes to existing neurons from
randomly generated neurons, removing spikes from existing neurons, shifting
the entire neuron population by a constant amount, slightly shifting neurons
in time, replacing neurons with randomly generated neurons and eliminating
neurons entirely.

Figure 1: Inter-session ensemble variability possible when recording from
neural populations. Neurons from the original recording session can be lost to
the recording array, new neurons can become visible, neurons can move between
electrodes, original neurons can be replaced by unseen neurons and the entire
probe array can shift, causing a systematic change in neuron position. In
addition, spike sorting can induce variability as the signal to noise ratio of
individual neurons changes between sessions. The perturbations we apply to
each trial of recordings is in response to each of these sources of
variability. We model each unseen test trial as an instance of a perturbed
seen train trial and subsequently, our sequential autoencoder model attempts
to map each unseen trial to a known trial.

This neuron locator RNN is trained to predict original neuron position within
a single recording session from many perturbed variations of trials of this
training session. Once trained to predict original neuron positions, a
separate network, which in this case is a sequential autoencoder based on
Latent Factor Analysis via Dynamical Systems (LFADS) (Pandarinath et al.,
2017), is trained to predict original unperturbed neural recording trials from
perturbed variations of trials from the same session. The encoder of this
sequential autoencoder receives as additional input the embedding of the
neuron locator RNN activations, conditioning the encoder to produce latent
variables which are informative enough to accurately reconstruct the original
recording. The encoder produces latent variables which are separated by
behaviour (arm movement direction) in a self-supervised manner, from which
behaviour can be predicted without the model being explicitly trained on
behaviour.

Importantly, the joint neuron locator RNN and LFADS encoder ensemble can
predict behaviourally relevant latent variables for unseen recording sessions
that yield high decoding accuracy. Currently, there are no existing approaches
to accurately predict behaviour from an unseen recording session when training
on just one single session. We not only show this is possible with our method,
but that our approach is robust to inter-session variability for up to 8 days
when a sufficient number of neurons are persistent across sessions.

##  2 Related Work

There have been many recent approaches to creating robust behaviour decoders
of neural activity (Gallego et al., 2020; Farshchian et al., 2019; Sussillo et
al., 2016; Wen et al., 2021; Karpowicz et al., 2022; Wimalasena et al., 2021).
However these methods are not capable of decoding behaviour from a previously
unseen recording session if the recorded activity is subject to random
fluctuations.

Recent work in modelling neural activity shows the consequences of selectively
perturbing neural data in order to learn relevant latent variables in a self-
supervised way using an autoencoder (Liu et al., 2021; Azabou et al., 2021;
Zhu et al., 2021). These models take different views of the same neural data
and align the latent spaces of these views once passed through an encoder,
with the ultimate aim of reconstructing these views. We utilise a similar
technique to train our sequential autoencoder by aligning the latent variables
of perturbed versions of the same data and aim to generate the activity of the
original unperturbed trial. Importantly, Liu et al. (2021) propose a model
which is invariant to the specific neurons used to represent the neural state
within training data; in this work we look at unseen sessions and so do not
aim to produce a model invariant to new neurons, but one that is able to
identify and utilise seen neurons to reconstruct unperturbed trials.

Gonschorek et al. (2021) and Jude et al. (2022) use domain adaptation to align
data across recording sessions. In both studies the authors use an autoencoder
model and a domain classifier. However these models require training on many
days of recording sessions for good behaviour decoding accuracy. For instance,
Jude et al. (2022) requires as many as 12 training sessions and training on
behaviour explicitly in order to produce high behaviour decoding accuracy on
an unseen test session. In this work we achieve state-of-the-art behaviour
decoding performance on an unseen test recording session using just one
training recording session, and show that this decoding accuracy can be
maintained many days into the future without recalibration.

We train an RNN to predict original neuron position from perturbed trials and
utilise this network to inform the sequential autoencoder model. This is
considered self-supervised learning as we do not train our model on behaviour
explicitly but instead train on the subtasks of predicting original neuron
positions and reconstructing unperturbed trials from perturbed ones. This
approach is similar to that used in Noroozi and Favaro (2016), where authors
form 9 subsets of images and randomly permute these subsets, then task the
model with predicting the permutation.

##  3 M1 Recordings

Figure 2: Experimental setup: In each trial one randomly chosen target
direction (indicated by one of 8 coloured circles) appears on screen, and the
monkey is instructed to control the cursor (white circle) by moving the
manipulandum. The monkey moves the cursor to the target location after a go
cue. The collected data for each trial consists of the neural spikes and
monkey hand position across all timesteps. We predict hand position from
neural spikes at each timestep.

We apply our model to data from a previously published experiment (Gallego et
al., 2020). In this experiment, two monkeys were trained to perform a center-
out reach task towards eight outer targets. On a go cue, each monkey moves a
manipulandum along a 2D plane to guide a cursor on a screen to the target
location (Figure 2). On successful trials a liquid reward is given. Spiking
activity from the motor cortex (M1) along with the 2D hand position were
recorded during each trial. Spike trains were converted into spike counts in
10ms bins, and behaviour variables are used at the same resolution. In this
work, only successful trials are used, all trials are aligned to movement
onset and cut from movement onset to the shortest reach time across all
trials.

For our analysis, we train our model on one session of recorded data from a
single day which we denote day 0 (containing 173 trials for both monkeys) and
test on subsequent held out days of recordings for each monkey. A comparison
of the activity between sessions shows considerable variability, caused by
shifts in the order neurons appear as well as disappearance of neurons and the
appearance of new ones (see Appendix B, Figure 8). These changes are
particularly pronounced for longer time intervals, but are already significant
in recordings one day apart. In total we used 5 days of recordings for both
monkeys, with 55 recorded neurons across all sessions for Monkey C and 17 for
Monkey M. Each day for each monkey consists of one recording session.

##  4 Data Perturbations

Figure 3: Perturbations applied simultaneously to each trial of neural data,
demonstrated with a simple 5 neuron system. A) Replace entire spike train with
a randomly generated neuron of the same firing rate as the original neuron. B)
Spikes randomly added to spike train proportional to average firing rate of
all neurons in a given trial, to mirror influence of nearby newly added
unknown neurons. C) Spikes randomly removed to mimic removal or movement of
nearby known neurons. D) Deletion of entire neurons to simulate neuron loss
between sessions, with randomly generated neurons introduced as the first or
last neuron of the trial to keep neuron number consistent. E) Small random
time jitter of all neuron spike trains to simulate experimental variation
between sessions. F) Constant random shift of the order of all neurons to
mirror probe shift.

Fig. 3 outlines the perturbations forming each variation of a single trial
during the training of our model. Perturbations A) to D) in Fig. 3 are applied
with equal probability to a given neuron of a given trial. Perturbation E) is
applied to all neurons, time jitter is chosen randomly between -30ms and
+30ms. Perturbation F) is applied to all trials, the amount of this neuron
shift is chosen randomly between 0 and 25% of the total number of neurons. We
hypothesise that this combination of transformations sufficiently mirrors the
real day to day changes of recorded neuron ensembles.

##  5 Model

Figure 4: Our model consists of a neuron locator RNN (1) combined with a
sequential variational autoencoding approach (2). The neuron locator (1) is
trained first to identify original neuron position (or if the neuron is
randomly generated) in each trial after perturbations have been applied. Then
the neuron locator’s weights are frozen and its activations are given as
additional input to condition the encoder of the sequential autoencoder (2).
Notably, we perturb recording trials when training both the neuron locator and
sequential variational autoencoder. The sequential autoencoder is tasked with
reconstructing the original unperturbed recording trials. The encoder of the
sequential autoencoder maps perturbed versions of the same trial to similar
latent variables. This is accelerated by imposing an alignment loss across the
latent variables of variations of the same trial. The generator RNN of the
sequential autoencoder predicts original trials from latent variables produced
by the encoder RNN.

Our modelling approach is based on the hypothesis that the perturbations
mentioned above can capture the substantial variability between recording
sessions from the same animal. We also expect neural activity $x$ is related
to the latent variables $z$ through a simple function, however, this function
will differ between recording sessions as we expect to observe different
neurons in each session. The problem is thus to find the correct encoding
function $z=f(x)$ to transform perturbed neural activity into a consistent
latent space which then allows decoding of behaviour. In addition, for the
same behaviour we require $z_{i}$ for each trial $i$ to be similar despite
variations in the activity $x_{i}$.

We first train a fully connected layer and an RNN to predict original neuron
position in perturbed trials. We apply the perturbations from Figure 3 to each
trial, then task the network to predict the original position of each neuron
in the recording data or whether it was previously unseen. As shown in Figure
4, for each neuron in the recording data we project a softmax linear read-out
layer from the RNN which each form a probability distribution of predicted
original neuron position across all possible positions (plus an extra position
indicating that the neuron was randomly generated). Each of these is compared
against a one hot encoding of the original neuron position before any
perturbations have been applied. If the neuron is randomly generated then the
one-hot encoding is one at the dedicated extra position. Predictions of
original neuron position are made as follows:

$\displaystyle\bar{x}_{i,1:T}=\text{Perturb}(x_{i,1:T}),$ (1)
$\displaystyle\text{acts}_{i}=\text{GRU}_{\theta_{\textrm{pos}}}(f_{\textrm{pos}}(\bar{x}_{i,1:T})),$
(2)
$\displaystyle\text{pos}_{i,n}=\text{softmax}(W^{n}_{\textrm{neuron}}.\text{acts}_{i})$
(3) $\displaystyle\text{The predicted position for trial $i$ and neuron $n$ is
then:}\operatornamewithlimits{argmax}\text{pos}_{i,n}$

Perturb is the simultaneous application of all perturbations outlined in
Section 4 to a given trial. $f_{\textrm{pos}}$ is a fully connected layer and
$\theta_{\textrm{pos}}$ are the parameters of the locator network used to
predict original neuron position. $W^{n}_{\textrm{neuron}}$ is the set of
linear layers used to predict original neuron position, producing a
probability distribution when combined with a softmax layer for each neuron.

Once trained, the weights of this neuron locator network are frozen, and the
activations of the RNN are used as additional input to the encoder of an
LFADS-inspired sequential autoencoder. This input conditions the encoder in
predicting latent variables used to generate original trials from perturbed
trials. As proposed by Pandarinath et al. (2017) we assume that the latent
dynamics evolve autonomously provided a set of initial conditions $z_{i}$ that
are modelled as Gaussian random variables. These latent variables are produced
for each trial by an encoder network consisting of bidirectional Gated
Recurrent Units (Cho et al., 2014) (GRU). They are used to reconstruct the
original trial-specific neural activity from the perturbed trials. A further
bidirectional GRU is used as a generator for neural reconstruction of
unperturbed trials from latent variables $z_{i}$. Training is based on Poisson
likelihood for unperturbed neural activity reconstruction (as in (Pandarinath
et al., 2017)). The model is trained using real neural activity which
corresponds to consistent behaviours (movement directions in a centre-out
reach task). The generative process of our model is as follows:

$\displaystyle
z_{i}=f_{\textrm{enc}}(\text{GRU}_{\theta_{\textrm{enc}}}(\bar{x}_{i,1:T};\text{acts}_{i})),$
(4) $\displaystyle g_{1:T}=\text{GRU}_{\theta_{\textrm{gen}}}(z_{i}),$ (5)
$\displaystyle r_{t}=exp(W_{\textrm{rate}}.f_{\textrm{fac}}(g_{t})),$ (6)
$\displaystyle\hat{x}_{t}\sim\text{Poisson}(r_{t})$ (7)

where $\theta_{\textrm{enc}}$ and $\theta_{\textrm{gen}}$ are the parameters
of the GRUs used to encode perturbed spike trains into latent variables and
subsequently generate original unperturbed spike trains from the latent
variables. $f_{\textrm{enc}}$ and $f_{\textrm{fac}}$ are fully connected
layers which produce latent variables and neural activity factors
respectively. $W_{\textrm{rate}}$ is a linear transformation used to generate
firing rates at each time step per trial. At each training iteration the
following three losses are optimised with Adam (Kingma and Ba, 2015):

$\displaystyle
L_{{\textrm{rec}}}=-\sum_{t=1}^{t}\log(\text{Poisson}(x_{i,t}|r_{t}))$ (8)
$\displaystyle
L_{kl}=D_{KL}[\text{GRU}_{\theta_{\textrm{enc}}}(z_{i}|\bar{x}_{i};\text{acts}_{i})||\mathcal{N}(0,I)]=-\frac{1}{2}[\log(z_{i,\sigma}^{2})-z_{i,\mu}^{2}-z_{i,\sigma}^{2}+1]$
(9) $\displaystyle L_{\textrm{align}}=\frac{1}{P}\sum_{j=1}^{p}\sum_{k\neq
j}^{p}(z_{i,j}-z_{i,k})^{2}$ (10)

Together $L_{\textrm{rec}}$ and $L_{kl}$ are the usual evidence lower-bound of
the marginal log-likelihood in a VAE (Kingma and Welling, 2014).
$L_{\textrm{rec}}$ is minimised by the encoder network and the neural
generator network. As in Liu et al. (2021), we apply an alignment loss
($L_{\textrm{align}}$) across latent variables produced from perturbed trials
(where $P$ is the number of perturbations of a given trial) of the same
original trial $z_{i}$ which reduces training duration. We form 2 perturbed
variations of each trial in a given batch at each training iteration.
Kullback–Leibler ($L_{kl}$) divergence loss (between a multivariate standard
Gaussian distribution and the encoder-generated latent variables) and
$L_{align}$ are minimised by just the encoder network. We name our model
CAPTure and Identify Variability at Target Ensembles (CAPTIVATE). Further
implementation details can be found in Appendix A.

###  5.1 Comparison models

We compare the ability of CAPTIVATE to predict behaviour from sessions of
unseen spike data against existing methods and against a variation of our own
model where we do not use the locator network trained on original neuron
position to aid in aligning perturbed trials. We denote this model variation
CAPTIVATE-noLoc. In addition, we look at vanilla LFADS (Pandarinath et al.,
2017) in autoencoding trials without any perturbations. We also compare
against a baseline RNN (GRU) with a linear readout layer explicitly trained to
reconstruct movement behaviour from neural activity.

For all autoencoding models we use a separately trained GRU network to predict
behaviour from the day 0 training session latent space. We do not include ADAN
(Farshchian et al., 2019), NoMAD (Karpowicz et al., 2022) or the generative
model by Wen et al. (2021) as all require training data from a held out
session or subject to be effective. We also do not test against Gonschorek et
al. (2021) or (Jude et al., 2022) as these approaches require many training
sessions to be effective in predicting behaviour from an unseen session
whereas we aim to do this with just one training session.

##  6 Results

Figure 5 shows behaviour decoding performance of CAPTIVATE for an unseen
session that was recorded the day after the training session for different
total rates of perturbation. A total perturbation rate of 40% (i.e a rate of
10% for each perturbation A) - D) in section 4) for both monkeys appears to be
optimal. At perturbation rates above 40%, neural activity from perturbed day 0
train trials with a particular target movement direction begin to resemble
original trials of other movement directions, and thus hurt alignment.
Perturbation rates below 40%, particularly for Monkey C, are not sufficient to
simulate the inter-session variability between day 0 and day 1. Training the
neuron locator RNN on a total perturbation rate of 40% for both monkeys yields
85% and 93% accuracy on predicting original neuron position from day 0
perturbed trials from Monkey C and Monkey M respectively. Indeed, the neuron
locator network is 76% accurate at identifying original neuron position in a
simulated unseen session created with a total perturbation rate of 80% (see
Appendix D, Figure 10).

Figure 5: Behaviour decoding performance on an immediately subsequent unseen
session (day 1) of CAPTIVATE at different rates of total perturbation. Total
perturbation rate is the sum of the rates of perturbations A) - D) outlined in
section 4, each of which are applied at equal rates.

Using the optimal rate of 40% of perturbation to trials from both monkeys when
training CAPTIVATE leads to the results summarised in Figure 6. For both
monkeys we see high behaviour decoding performance on the unseen session from
day 1, surpassing previous methods. CAPTIVATE maintains high behaviour
decoding performance for Monkey C on an unseen session up to 8 days after the
day 0 training session was recorded. CAPTIVATE also accurately maps neurons
from trials across unseen sessions of Monkey C up to 8 days into the future to
known neurons from trials of the day 0 train session (see Appendix C, Figure
9). Notably, behaviour decoding for Monkey C is much more stable for future
unseen sessions than for Monkey M. This is likely due to sessions from Monkey
C containing more than 3 times as many neurons as Monkey M. However, we see in
Appendix E, Figure 11 that training CAPTIVATE with 20 neurons from the Monkey
C day 0 session is sufficient to achieve an $R^{2}$ of 0.68 when testing on 20
neurons of the day 8 session, indicating our model can be robust to a low
number of neurons.

Figure 6: Behaviour prediction performance when testing all models on 30% of
held-out trials from day 0 and subsequent days of completely unseen recording
sessions. We report the mean $R^{2}$ between the inferred and true x,y
positions. Each model is tested on held out trials from day 0 and trials from
unseen sessions recorded an increasing number of days into the future from the
original training session (day 0) for both monkeys. Each day 0 train session
is run 10 times with different random seeds, with error bars showing standard
deviation when applied to each unseen session.

Notably in the case of Monkey M, day 1 decoding performance is high at all
levels of perturbation from 0.1 to 0.4 (Figure 5), therefore it is likely that
the session to session variability between day 0 and day 1 is small. Thus, for
a subject with fewer neurons in recorded data, CAPTIVATE may only require a
low rate of total perturbation when aligning nearby unseen sessions.

CAPTIVATE-noLoc, Vanilla LFADS or an RNN model cannot capture session-to-
session variability even for the day 1 unseen session, as shown in Figure 7.
CAPTIVATE-noLoc cannot accurately reconstruct original trials from perturbed
variations of the day 0 train session, but has a similar day 0 and day 1
session behaviour decoding accuracy, implying our perturbations closely mirror
inter-session variability. This indicates that poor performance of CAPTIVATE-
noLoc on both monkeys is due to the inability of the encoder of this model to
recognise known neurons and thus, shows how crucial the neuron locator network
is in recognising known neuron ensembles in unseen recordings.

LFADS is trained solely on unperturbed trials and so cannot recognise the
shifts that occur between sessions as variations of the day 0 training
session, and thus cannot create an appropriate latent encoding. The RNN model
is also trained on unperturbed trials and is even less robust to later unseen
sessions than LFADS, however, this RNN baseline can recognise some behaviour
in both monkeys for the day 1 unseen session, indicating a relatively low
level of variability in adjacent day recordings. Similarly, we see a mean
$R^{2}$ of 0.37 ($\pm$ 0.02) when training an RNN on day 7 for monkey C and
testing on day 8 and an $R^{2}$ of 0.42 ($\pm$ 0.04) when training an RNN on
day 9 and testing on day 10 for monkey M. Importantly, none of these models
overfit as they yield high decoding accuracies for a held-out portion of day 0
trials for both monkeys and for all models, especially the RNN. Therefore the
performance drop of the RNN model when applied to unseen sessions is a clear
indication of substantial variations between sessions.

Figure 7: For each monkey, Top row: t-SNE embeddings of latent space for
CAPTIVATE when applied to each unseen session. In each embedding, points
denoted by a circle are trials from the day 0 training session. Points denoted
by a triangle are trials from the named unseen session. Each colour represents
a target direction for the centre-out reach task. Bottom row: Predicted 2D
monkey hand position of trials using a separately trained RNN decoder trained
only on the day 0 latent space of CAPTIVATE when applied to each unseen
session, with mean $R^{2}$ between all positions of each predicted and ground
truth trajectory shown across all trials in a given unseen session.

Figure 7 shows t-SNE visualisations of the latent space and behaviour
predictions made from the latent space of CAPTIVATE when trained on the day 0
session and applied to unseen sessions. For Monkey C, the majority of trials
from all unseen sessions are correctly aligned with the corresponding trials
in the training data set (Figure 7, compare dots and triangles in the t-SNE
plots where colour indicates movement direction; note that the latent space is
well partitioned by behaviour although the model is only trained on neural
activity). Occurrences where the unseen trials correctly overlap known train
trials, in turn, yields correctly decoded behaviour. The alignment becomes
progressively worse for later sessions, and as the alignment is less precise,
behaviour predictions also become worse. In contrast, for Monkey M the
alignment of trials beyond day 2 becomes increasingly worse, a consequence of
the smaller number of neurons in the recording.

Table 1: Mean decoding performance effects of ablating individual perturbations when training on day 0 session and tested on immediately subsequent (day 1) unseen session for both monkeys. For reference, the full CAPTIVATE model trained on day 0 achieves mean decoding $R^{2}$ performance of 0.84 ($\pm$ 0.02) on monkey C and 0.86 ($\pm$ 0.03) on monkey M when applied to the day 1 session. Ablation | No-Replace | No-Add | No-Remove | No-Delete | No-Jitter | No-Reorder
---|---|---|---|---|---|---
C Mean $R^{2}$ | 0.66 | 0.79 | 0.74 | 0.49 | 0.77 | 0.70
| ($\pm$ 0.03) | ($\pm$ 0.01) | ($\pm$ 0.02) | ($\pm$ 0.04) | ($\pm$ 0.01) | ($\pm$ 0.03)
M Mean $R^{2}$ | 0.71 | 0.75 | 0.81 | 0.63 | 0.81 | 0.77
| ($\pm$ 0.03) | ($\pm$ 0.02) | ($\pm$ 0.01) | ($\pm$ 0.03) | ($\pm$ 0.01) | ($\pm$ 0.02)

Ablations of individual perturbations (as outlined in Figure 3) applied when
training on the day 0 session reveal that perturbations which introduce
randomly generated neurons and alter the continuous ordering of neurons have
the highest impact on unseen session behaviour decoding performance. This
analysis is summarised in Table 1 and shows that neuron deletions,
replacements and probe shifts cause the majority of inter-session neuron
ensemble variability. Nonetheless, a combination of all perturbations are
necessary for the decoding performance achieved by CAPTIVATE in Figure 6. We
also train CAPTIVATE without the alignment loss in Eq. 10, which produces a
behaviour decoding mean $R^{2}$ of 0.82 on Monkey C and 0.85 on Monkey M on
trials from the day 1 unseen session. This minimal drop in decoding
performance when training without an explicit alignment loss is consistent
with results from (Liu et al., 2021). Additionally, decoding performance
across unseen sessions when training on day 0 and day 11 sessions separately
is almost symmetrical (as shown in Appendix F, Figure 12), indicating that our
model can effectively capture neural variability from unseen sessions both
forwards and backwards in time. We further assess robustness by testing
CAPTIVATE on a variable number of neurons across sessions (similar to a real
BCI setting) and show good generalisation, even surpassing performance of the
model trained with 55 neurons (as in Figure 6) for some unseen sessions (see
Appendix G, Figure 13).

##  7 Discussion

In this paper we use a self-supervised approach, CAPTIVATE, to train a model
to recognise and correct for session-to-session variability in neural
recordings. We then show that the combination of this approach with a latent
variable model that identifies low-dimensional dynamics in neural activity
yields a model that is now robust variability between recordings sessions. The
model is capable of successfully predicting behaviour with high accuracy from
unseen sessions, surpassing previous work by Jude et al. (2022) when comparing
against subsequent day decoding performance. Furthermore, our approach leads
to relatively high and stable behaviour decoding performance on unseen
sessions many days into the future when a sufficient number of neurons are
persistent across sessions. As a result, this method performs better for data
sets with more recorded neurons (Monkey C), while for fewer neurons the
performance degrades more quickly, only producing good results for sessions
close in time to the training session (Monkey M).

With CAPTIVATE we achieve stable behaviour decoding performance for up to 8
days, which is followed by a slow decline in performance. The decline is due
to an increase in variability that could no longer be compensated. This would
require a model to correct even stronger perturbations, but training a model
this way leads to an overall decrease in performance even for short time
intervals (Figure 5). Therefore long-term stable decoding currently still
requires re-training of the components of a latent variable encoder model such
that the altered neural dynamics are re-aligned with the latent dynamics (Wen
et al., 2021; Karpowicz et al., 2022; Farshchian et al., 2019). Equally, our
model fails to successfully decode behaviour from recordings from an unseen
animal (not illustrated) as this requires a more complex mapping function
between activity and latent space (Wen et al., 2021).

## References

  * Azabou et al. (2021) M. Azabou, M. G. Azar, R. Liu, C.-H. Lin, E. C. Johnson, K. Bhaskaran-Nair, M. Dabagia, B. Avila-Pires, L. Kitchell, K. B. Hengen, W. Gray-Roncal, M. Valko, and E. L. Dyer.  Mine your own view: Self-supervised learning through across-sample prediction, 2021.  URL https://arxiv.org/abs/2102.10106.
  * Cho et al. (2014) K. Cho, B. Van Merriënboer, C. Gulcehre, D. Bahdanau, F. Bougares, H. Schwenk, and Y. Bengio.  Learning phrase representations using RNN encoder-decoder for statistical machine translation.  In _EMNLP 2014 - 2014 Conference on Empirical Methods in Natural Language Processing, Proceedings of the Conference_ , 2014.  ISBN 9781937284961.  doi: 10.3115/v1/d14-1179.
  * Cunningham and Byron (2014) J. P. Cunningham and M. Y. Byron.  Dimensionality reduction for large-scale neural recordings.  _Nature Neuroscience_ , 17(11):1500–1509, 2014\.
  * Driscoll et al. (2017) L. N. Driscoll, N. L. Pettit, M. Minderer, S. N. Chettih, and C. D. Harvey.  Dynamic reorganization of neuronal activity patterns in parietal cortex.  _Cell_ , 170(5):986–999, 2017.
  * Farshchian et al. (2019) A. Farshchian, J. A. Gallego, L. E. Miller, S. A. Solla, J. P. Cohen, and Y. Bengio.  Adversarial domain adaptation for stable brain-machine interfaces.  In _7th International Conference on Learning Representations, ICLR 2019_ , 2019.
  * Gallego et al. (2020) J. A. Gallego, M. G. Perich, R. H. Chowdhury, S. A. Solla, and L. E. Miller.  Long-term stability of cortical population dynamics underlying consistent behavior.  _Nature Neuroscience_ , 23(2), 2020.  ISSN 15461726.  doi: 10.1038/s41593-019-0555-4.
  * Gonschorek et al. (2021) D. Gonschorek, L. Höfling, K. P. Szatko, K. Franke, T. Schubert, B. A. Dunn, P. Berens, D. A. Klindt, and T. Euler.  Removing inter-experimental variability from functional data in systems neuroscience.  In _Thirty-Fifth Conference on Neural Information Processing Systems_ , 2021.  URL https://openreview.net/forum?id=lVmIjQiJJSr.
  * Huber et al. (2012) D. Huber, D. A. Gutnisky, S. Peron, D. H. O’connor, J. S. Wiegert, L. Tian, T. G. Oertner, L. L. Looger, and K. Svoboda.  Multiple dynamic representations in the motor cortex during sensorimotor learning.  _Nature_ , 484(7395):473–478, 2012.
  * Hurwitz et al. (2021) C. Hurwitz, N. Kudryashova, A. Onken, and M. H. Hennig.  Building population models for large-scale neural recordings: Opportunities and pitfalls, 2021.  ISSN 18736882.
  * Jude et al. (2022) J. Jude, M. G. Perich, L. E. Miller, and M. H. Hennig.  Robust alignment of cross-session recordings of neural population activity by behaviour via unsupervised domain adaptation, 2022.  URL https://arxiv.org/abs/2202.06159.
  * Karpowicz et al. (2022) B. M. Karpowicz, Y. H. Ali, L. N. Wimalasena, A. R. Sedler, M. R. Keshtkaran, K. Bodkin, X. Ma, L. E. Miller, and C. Pandarinath.  Stabilizing brain-computer interfaces through alignment of latent dynamics.  _bioRxiv_ , 2022.  doi: 10.1101/2022.04.06.487388.  URL https://www.biorxiv.org/content/early/2022/04/08/2022.04.06.487388.
  * Kingma and Ba (2015) D. P. Kingma and J. L. Ba.  Adam: A method for stochastic optimization.  In _3rd International Conference on Learning Representations, ICLR 2015 - Conference Track Proceedings_ , 2015.
  * Kingma and Welling (2014) D. P. Kingma and M. Welling.  Auto-encoding variational bayes.  In _2nd International Conference on Learning Representations, ICLR 2014 - Conference Track Proceedings_ , 2014.
  * Liu et al. (2021) R. Liu, M. Azabou, M. Dabagia, C.-H. Lin, M. G. Azar, K. B. Hengen, M. Valko, and E. L. Dyer.  Drop, swap, and generate: A self-supervised approach for generating neural activity.  In _Thirty-Fifth Conference on Neural Information Processing Systems_ , 2021.
  * Noroozi and Favaro (2016) M. Noroozi and P. Favaro.  Unsupervised learning of visual representations by solving jigsaw puzzles.  In _Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics)_ , volume 9910 LNCS, 2016.  doi: 10.1007/978-3-319-46466-4{\\_}5.
  * Pandarinath et al. (2017) C. Pandarinath, D. O’Shea, J. Collins, R. Jozefowicz, S. Stavisky, J. Kao, E. Trautmann, M. Kaufman, S. Ryu, L. Hochberg, J. Henderson, K. Shenoy, and D. Sussillo.  Inferring single-trial neural population dynamics using sequential auto-encoders.  _Inferring single-trial neural population dynamics using sequential auto-encoders_ , 2017.  ISSN 1548-7091.  doi: 10.1101/152884.
  * Sussillo et al. (2016) D. Sussillo, S. D. Stavisky, J. C. Kao, S. I. Ryu, and K. V. Shenoy.  Making brain-machine interfaces robust to future neural variability.  _Nature Communications_ , 7, 2016.  ISSN 20411723.  doi: 10.1038/ncomms13749.
  * Wen et al. (2021) S. Wen, A. Yin, T. Furlanello, M. G. Perich, L. E. Miller, and L. Itti.  Rapid adaptation of brain–computer interfaces to new neuronal ensembles or participants via generative modelling.  _Nature Biomedical Engineering_ , 2021.  ISSN 2157846X.  doi: 10.1038/s41551-021-00811-z.
  * Wimalasena et al. (2021) L. N. Wimalasena, J. F. Braun, M. R. Keshtkaran, D. Hofmann, J. Á. Gallego, C. Alessandro, M. C. Tresch, L. E. Miller, and C. Pandarinath.  Estimating muscle activation from emg using deep learning-based dynamical systems models.  _bioRxiv_ , 2021.  doi: 10.1101/2021.12.01.470827.
  * Zhu et al. (2021) F. Zhu, A. Sedler, H. A. Grier, N. Ahad, M. Davenport, M. Kaufman, A. Giovannucci, and C. Pandarinath.  Deep inference of latent dynamics with spatio-temporal super-resolution using selective backpropagation through time.  In M. Ranzato, A. Beygelzimer, Y. Dauphin, P. Liang, and J. W. Vaughan, editors, _Advances in Neural Information Processing Systems_ , volume 34, pages 2331–2345. Curran Associates, Inc., 2021.  URL https://proceedings.neurips.cc/paper/2021/file/1325cdae3b6f0f91a1b629307bf2d498-Paper.pdf.
  * Ziv et al. (2013) Y. Ziv, L. D. Burns, E. D. Cocker, E. O. Hamel, K. K. Ghosh, L. J. Kitch, A. El Gamal, and M. J. Schnitzer.  Long-term dynamics of ca1 hippocampal place codes.  _Nature Neuroscience_ , 16(3):264, 2013.

## Checklist

  1. 1.

For all authors…

    1. (a)

Do the main claims made in the abstract and introduction accurately reflect
the paper’s contributions and scope? [Yes] See Results, Section 6.

    2. (b)

Did you describe the limitations of your work? [Yes] See Discussion, Section
7.

    3. (c)

Did you discuss any potential negative societal impacts of your work? [N/A]

    4. (d)

Have you read the ethics review guidelines and ensured that your paper
conforms to them? [Yes]

  2. 2.

If you are including theoretical results…

    1. (a)

Did you state the full set of assumptions of all theoretical results? [N/A]

    2. (b)

Did you include complete proofs of all theoretical results? [N/A]

  3. 3.

If you ran experiments…

    1. (a)

Did you include the code, data, and instructions needed to reproduce the main
experimental results (either in the supplemental material or as a URL)? [Yes]
Perturbation and model instructions in sections 4 and 5 respectively. Code
included in Supplement.

    2. (b)

Did you specify all the training details (e.g., data splits, hyperparameters,
how they were chosen)? [Yes] In Section 5 and further implementation details
in Appendix A.

    3. (c)

Did you report error bars (e.g., with respect to the random seed after running
experiments multiple times)? [Yes] See figures in Section 6.

    4. (d)

Did you include the total amount of compute and the type of resources used
(e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] Yes, see
details of compute used in Section 7.

  4. 4.

If you are using existing assets (e.g., code, data, models) or
curating/releasing new assets…

    1. (a)

If your work uses existing assets, did you cite the creators? [Yes] See
Section 3.

    2. (b)

Did you mention the license of the assets? [N/A]

    3. (c)

Did you include any new assets either in the supplemental material or as a
URL? [Yes] Code included in supplemental material.

    4. (d)

Did you discuss whether and how consent was obtained from people whose data
you’re using/curating? [N/A]

    5. (e)

Did you discuss whether the data you are using/curating contains personally
identifiable information or offensive content? [N/A]

  5. 5.

If you used crowdsourcing or conducted research with human subjects…

    1. (a)

Did you include the full text of instructions given to participants and
screenshots, if applicable? [N/A]

    2. (b)

Did you describe any potential participant risks, with links to Institutional
Review Board (IRB) approvals, if applicable? [N/A]

    3. (c)

Did you include the estimated hourly wage paid to participants and the total
amount spent on participant compensation? [N/A]

##  Appendix A Implementation and training details

Below are implementation details for the CAPTIVATE model.

CAPTIVATE
---
Parameter |  Value |  Notes
Neuron Locator Network |  |  Layer Normalisation on all layers
\- RNN Units |  784 X 3 |  Stacked Gated Recurrent Unit
\- $W_{pos}$ Units |  1024 X 3 |  Non-linear layer
\- $W_{pos}$ Dropout |  0.5 |
\- $W_{pos}$ L2 Regularisation |  100.0 |
Sequential Autoencoder Encoder |  |
\- RNN Units |  784 X 3 |  Stacked Gated Recurrent Unit
\- RNN L2 Kernel Regularisation |  0.1 |
\- RNN L2 Recurrent Regularisation |  0.1 |
\- $W_{enc}$ Units |  1024 X 3 |  Non-linear layer
\- $W_{enc}$ L2 Regularisation |  0.1 |
\- Latent space dimension |  64 |
Sequential Autoencoder Generator |  |
\- RNN Units |  512 X 3 |  Stacked Gated Recurrent Unit
\- RNN L2 Kernel Regularisation |  1.0 |
\- RNN L2 Recurrent Regularisation |  1.0 |
\- $W_{fac}$ Units |  512 |  Non-linear layer
Training |  |
\- KL divergence weighting ($\lambda_{kl}$) |  0.02 to 1.0 |  Rising exponentially
\- Batch size (Train Neuron Locator) |  16 |
\- Batch size (Train Seq. Autoencoder) |  4 |
\- Learning rate (Train Neuron Locator) |  0.0001 |  Adam Optimizer
\- Learning rate (Train Seq. Autoencoder) |  0.00001 |  Adam Optimizer

##  Appendix B Changes in recorded neural activity across sessions

Figure 8: A, Pairwise correlations of trial-averaged activity of single
neurons between two sessions. For each neuron, the average firing rate was
computed for each of the eight movement directions (see part B for examples)
at 10ms resolution. The activity for the eight movement directions was
concatenated and the Pearson correlation coefficients computed between all
neuron pairs. Each plot shows the correlation matrix for activity from session
from a different day and activity from the first day (day zero, the training
data set in Figure 6). This analysis shows that some neurons from the first
session can be matched to neurons recorded at subsequent days, but the
relative position of these matched neurons in the recording tends to shift
(see high off-diagonal correlations). As the average correlations do not
change systematically over this period of time (not illustrated), the gradual
changes in neuron identity is a main factor that prevents reliable decoding
from unseen sessions in previous models. B, Examples of trial-averaged firing
rates of three neurons that were tracked over all recording sessions. This
matching is based on the similarity of the firing rates, experimentally it is
hard to determine if these are indeed the same neurons. In all cases, the time
course of the activity is similar and shows consistent differences between
trial type (indicated by colour) across sessions. Also note that while these
neurons appear to reliably encode movement direction, the activity of a single
neuron alone is too noisy to allow for reliable direction decoding from single
trials, instead a population decoding approach is required. All data
illustrated here is from Monkey C.

##  Appendix C CAPTIVATE accurately maps perturbed neurons and neurons from
unseen sessions to known neurons from the Day 0 training session

CAPTIVATE is trained by mapping perturbed trials to known trials. If trials
from unseen sessions are similar to the perturbed trials then generalisation
to these sessions is possible. Therefore, we aim for the encoder network of
CAPTIVATE to map perturbed trials and trials from unseen sessions to day 0
trials. This entails that neurons across unseen sessions (even after neural
drift and ensemble change) are mapped directly to neuron positions of the day
0 session at the session. For trials of each movement direction from unseen
sessions, we expect that the trial average firing rates of these neurons will
map to the day 0 average firing rates for each neuron. As seen below for 4
neurons across 3 sessions (2 unseen), the CAPTIVATE generator network produces
trial average firing rates matching the day 0 train session firing rates.

Figure 9: CAPTIVATE is trained on the day 0 session of Monkey C. On the left
we show real trial averaged firing rates for each movement direction across 4
randomly selected neurons across the day 0 session and 2 unseen sessions. On
the right we show predictions from the generator network of CAPTIVATE. If
generalisation is achieved the generator should accurately map neurons across
unseen sessions to the neurons of day 0. We see that this is the case as the
predicted firing rates are closely matched in the unseen sessions to the day 0
firing rates across movement directions.

##  Appendix D Neuron Locator performance over simulated neural variation

As we do not have ground truth neuron identities from unseen sessions (with
respect to the day 0 train session), we simulate inter-session variability by
increasing perturbation rate and testing against CAPTIVATE trained on the day
0 session from monkey C with a total perturbation rate of 0.4 (as in the
results shown in Figure 6). We see that the neuron locator network of
CAPTIVATE can predict neuron identity with 68% accuracy even at a very high
total perturbation rate of 1.0.

Figure 10: Neuron locator network accuracy when predicting neuron identity
(with respect to unperturbed day 0 monkey C train session) as the total rate
of perturbation is increased. We are simulating neural drift and ensemble
shift across sessions. As we know the ground truth neuron identities, we can
assess how well the neuron locator can predict neuron identity.

##  Appendix E Training and testing CAPTIVATE with different numbers of
original neurons

Here we test the varying numbers of neurons across sessions of Monkey C when
using CAPTIVATE. We see that only 20 neurons are required across sessions for
good generalisation for up to 8 days.

Figure 11: Behaviour prediction performance when training CAPTIVATE on varying
numbers of neurons of the day 0 session recorded from Monkey C and testing on
all other unseen sessions of monkey C, using the same number of neurons as
used in the training session. We also test all neuron number variations of
CAPTIVATE on a held out portion of trials from day 0. We report the mean
$R^{2}$ between the inferred and true x,y positions for the entire movement
trajectory of each trial. Each day 0 train session is run 10 times with
different random seeds, with error bars showing standard deviation when
applied to each unseen session.

##  Appendix F Changing calibration session

Here we show that by training our model on perturbed trials we can generalise
to neural drift and recording array movement. CAPTIVATE accounts not only for
session variability in the future but also in the past.

Figure 12: Behaviour prediction performance when separately training CAPTIVATE
on the day 0 and day 11 sessions of Monkey C and testing on all other unseen
sessions. We see that performance across unseen sessions when training on
these sessions is almost symmetrical, indicating that our model can
effectively capture neural variability from sessions both backwards and
forwards in time.

##  Appendix G Variable neuron number per session

When using an implanted recording array we may lose electrodes or neurons due
to spike sorting error over a period of time. Here we show our model can
account for this variable neuron number.

Figure 13: We test CAPTIVATE with a variable number of neurons per session of
recording from Monkey C. In our original experiment we only utilise the first
55 neurons of each recording session as this is the minimum number across all
sessions. Here we use every neuron available per session (number of neurons
per session shown in Figure) and train CAPTIVATE on the day 0 session with 67
neurons. For all other sessions we add randomly generated neurons to
compensate. We see that CAPTIVATE is robust to the number of original neurons
being variable across sessions. Note the increase in generalisation
performance when the model is applied to the day 8 session. This is due to
this session having a relatively high number (60) of original neurons, and is
thus easier for the model to map trials from this session to known trials from
the day 0 training session than it is from other later unseen sessions.

##  Appendix H Testing trained model on known neural variability

We test the resilience of our whole model against an increasing total rate of
perturbation in order to ascertain how much variability the model can account
for. For the results below, CAPTIVATE is trained with a total perturbation
rate of 0.4 on the day 0 session of Monkey C. Note that the model is trained
to map perturbed trials to original unperturbed trials.

Figure 14: We train CAPTIVATE on trials from the day 0 session of Monkey C
with a 0.4 total rate of perturbation on each trial. We then test the trained
model on trials of the same session but with increasing rates of total
perturbation applied to trials. A) Mean r-squared error of movement predicted
from the latent space of the model vs. real movement trajectory of each trial.
B) Mean Poisson log-likelihood for neural activity reconstruction by the model
generator of original day 0 unperturbed trials. C) Mean squared error of model
predicted firing rates vs. real firing rates of original unperturbed day 0
trials.

##  Appendix I Testing trained model on known neural variability across held-
out trials

We test the resilience of our whole model against an increasing total rate of
perturbation in order to ascertain how much variability the model can account
for. For the results below, CAPTIVATE is trained with a total perturbation
rate of 0.4 on 70% of the trials of the day 0 session of Monkey C. We show
test performance on 30% of the trials of the day 0 session which are withheld
from training. Note that the model is trained to map perturbed trials to
original unperturbed trials.

Figure 15: We train CAPTIVATE on 70% of the trials from the day 0 session of
Monkey C with a 0.4 total rate of perturbation on each trial. We then test the
trained model on the remaining 30% of trials of the same session but with
increasing rates of total perturbation applied to these held-out trials. A)
Neuron locator network accuracy when predicting neuron identity (with respect
to unperturbed day 0 Monkey C train session) as the total rate of perturbation
is increased. We are simulating neural drift and ensemble shift across
sessions. As we know the ground truth neuron identities, we can assess how
well the neuron locator can predict neuron identity. B) Mean r-squared error
of movement predicted from the latent space of the model vs. real movement
trajectory of each trial. C) Mean Poisson log-likelihood for neural activity
reconstruction by the model generator of original day 0 unperturbed trials. D)
Mean squared error of model predicted firing rates vs. real firing rates of
original unperturbed day 0 trials.