# Sparse Linear Dynamical System with Its Application in Multivariate Clinical
Time Series

Zitao Liu
Department of Computer Science
University of Pittsburgh
Pittsburgh, PA 15213
<EMAIL_ADDRESS>
&Milos Hauskrecht
Department of Computer Science
University of Pittsburgh
Pittsburgh, PA 15213
<EMAIL_ADDRESS>

###### Abstract

Linear Dynamical System (LDS) is an elegant mathematical framework for
modeling and learning multivariate time series. However, in general, it is
difficult to set the dimension of its hidden state space. A small number of
hidden states may not be able to model the complexities of a time series,
while a large number of hidden states can lead to overfitting. In this paper,
we study methods that impose an $\ell_{1}$ regularization on the transition
matrix of an LDS model to alleviate the problem of choosing the optimal number
of hidden states. We incorporate a generalized gradient descent method into
the Maximum a Posteriori (MAP) framework and use Expectation Maximization (EM)
to iteratively achieve sparsity on the transition matrix of an LDS model. We
show that our Sparse Linear Dynamical System (SLDS) improves the predictive
performance when compared to ordinary LDS on a multivariate clinical time
series dataset.

##  1 Introduction

Developing accurate models of dynamical systems is critical for their
successful applications in outcome prediction, decision support, and optimal
control. A large spectrum of models have been developed and successfully
applied for these purposes in the literature [3, 11, 20]. In this paper we
focus on a popular model for time series analysis: the Linear Dynamical System
(LDS) [15] and its application to clinical time series [18, 19]. We aim to
develop a method to learn an LDS that performs better on future value
predictions when learned from a small amount of complex multivariate time
series dataset.

LDS is a widely used model for time series analysis of real-valued sequences.
The model is Markovian and assumes the dynamic behaviour of the system is
captured well using a small set of real-valued hidden-state variables and
linear-state transitions corrupted by a Gaussian noise. The observations in
LDS, similarly to hidden states, are real-valued. Briefly, the observations at
time $t$ are linear combinations of hidden state values for the same time.
While in some LDS applications the model parameters are known a priori, in the
majority of real-world applications the model parameters are unknown, and we
need learn them from data that consists of observation sequences we assume
were generated by the LDS model. While this can be done using standard LDS
learning approaches, the problem of learning an LDS model gives rise to
numerous important questions: Given the multivariate observation sequences,
how many hidden states are needed to represent the system dynamics well?
Moreover, since transition and observation matrices depend on the number of
hidden states, how do we prevent the overfit of the model parameters when the
number of examples is small?

In this work we address the above issues by presenting a method based on the
sparse representation of LDS (SLDS) that is able to adjust (depending on the
observation sequences in the data) the number of hidden states and at the same
time prevent the overfit of the model. Our approach builds upon the
probabilistic formulation of the LDS model, and casts the optimization of its
parameters as a maximum a posteriori (MAP) estimate, where the choice of the
parameter priors biases the model towards sparse solutions.

Our SLDS approach is distinctly different from previous work [4, 6, 9]. [4]
formulates the traditional Kalman filter as a one-step update optimization
procedure and incorporates sparsity constraints to achieve sparsity in the
hidden states. [9] trains an LDS for each training example and tries to find a
sparse linear combination of coefficients in order to combine the ensemble of
models. Neither [4] nor [9] directly achieve sparsity on the parameters of the
LDS, and furthermore, the performance of their resulting models still depends
on the optimal number of hidden states. [6] introduces a Bayesian
nonparametric approach to the identification of observation-only linear
systems, where no hidden states are involved. The underlining assumption is
that the observations are obtained from linear combinations of previous
observations and some system inputs, which may be too restrictive to model
complex multivariate time series and makes the model more sensitive to noisy
observations and outliers.

We test our sparse solution on the problem of modeling the dynamics of
sequences of laboratory test results. We show that it improves the learning of
the LDS model and leads to better accuracy in predicting future time-series
values.

Our paper is organized as follows. In Section 2 we review the basics of the
linear dynamical system. In Section 3 we describe SLDS – our method of
sparsifying the LDS parameters. Inference and learning details of SLDS are
explained in Section 3. Experimental results that compare SLDS method to
ordinary LDS are presented in Section 4. In Section 5, we summarize the work
and outline possible future extensions.

##  2 Linear Dynamical System (LDS)

The Linear Dynamical System (LDS) is a real-valued time series model that
represents observation sequences indirectly with the help of hidden states.
Let $\\{{z}_{t}\\}$, $\\{{y}_{t}\\}$ define sequences of hidden states and
observations respectively. The LDS models the dynamics of these sequences in
terms of the state transition probability $p({z}_{t}|{z}_{t-1})$, and state-
observation probability $p({y}_{t}|{z}_{t})$. These probabilities are modeled
using the following equations:

${z}_{t}=A{z}_{t-1}+{e}_{t};\hskip 14.22636pt{y}_{t}=C{z}_{t}+{v}_{t},$ (1)

where ${y}_{t}$ is a $d\times 1$ observation vector made at (current time)
$t$, and ${z}_{t}$ an $l\times 1$ hidden states vector. The transitions among
the current and previous hidden states are linear and captured in terms of an
$l\times l$ transition matrix $A$. The stochastic component of the transition,
${e}_{t}$, is modeled by a zero-mean Gaussian noise
${e}_{t}\sim\mathcal{N}(0,Q)$ with an $l\times 1$ zero mean and an $l\times l$
covariance matrix _Q_. The observations sequence is derived from the hidden
states sequence. The dependencies in between the two are linear and modeled
using a $d\times l$ emission matrix _C_. A zero mean Gaussian noise
${v}_{t}\sim\mathcal{N}(0,R)$ models the stochastic relation in between the
states and observation. In addition to $A,C,Q,R$, the LDS is defined by the
initial state distribution for ${z}_{1}$ with mean $\boldsymbol{\pi}_{1}$ and
covariance matrix $V_{1}$,
${z}_{1}\sim\mathcal{N}(\boldsymbol{\pi_{1}},V_{1})$. The complete set of the
LDS parameters is $\Omega=\\{A,C,Q,R,\boldsymbol{\pi_{1}},V_{1}\\}$. The
parameters of the LDS model can be learned using either the Expectation-
Maximization (EM) algorithm [8] or spectral learning algorithms [16, 21].

##  3 Sparse Linear Dynamical System (SLDS)

In this section, we propose a sparse representation of LDS that is able to
adjust the number of hidden states and at the same time prevents the overfit
of the model. More specifically, we impose $\ell_{1}$ regularizers on every
element of the transition matrix $A_{ij}$, which leads to zero entries in the
transition matrix _A_. The zero entries in the transition matrix of LDS indeed
reduce the actual number of parameters of LDS, sparsify the hidden states, and
avoid the overfitting problem from the real data, even if we set the number of
hidden states originally picked is too large.

To achieve sparsity on the transition matrix, we introduce a Laplacian prior
to each element of _A_ , $A_{ij}$, since Laplacian priors are equivalent to
$\ell_{1}$ regularizations [5, 10, 24]. In general, the Laplacian distribution
has the following form:
$p(x|\mu,\lambda)=\frac{1}{2\lambda}\exp(-\frac{|x-\mu|}{\lambda})$,
$\lambda\geq 0$ where $\mu$ is the location parameter and $\lambda$ is the
scale parameter. Here, we assume every element $A_{ij}$ is independent to each
other and has the following Laplacian density ($\mu=0$ and $\lambda=1/\beta$),
$p(A_{ij}|\beta)=\frac{\beta}{2}\exp(-\beta|A_{ij}|)$. Hence, the prior
probability for _A_ is
$p(A|\beta)=\prod_{i=1}^{l}\prod_{j=1}^{l}p(A_{ij}|\beta)$ and the log joint
distribution for SLDS is:

$\displaystyle\log p(\mathbf{z},\mathbf{y},A)=\log p(A)+\log
p(z_{1})+\sum_{t=2}^{T}\log p(z_{t}|z_{t-1},A)+\sum_{t=1}^{T}\log
p(y_{t}|z_{t})$ (2)

where _T_ is the observation sequence length.

###  3.1 Learning

In this section we develop an EM algorithm for the MAP estimation of the SLDS.
Let $\hat{z}_{t|T}\equiv\mathbb{E}[z_{t}|\mathbf{y}]$,
$M_{t|T}\equiv\mathbb{E}[z_{t}z_{t}^{{}^{\prime}}|\mathbf{y}]$,
$M_{t,t-1|T}\equiv\mathbb{E}[z_{t}z_{t-1}^{{}^{\prime}}|\mathbf{y}]$ and
define the $\mathcal{Q}$ function as
$\mathcal{Q}=\mathbb{E}_{\mathbf{z}}\Big{[}\log
p(\mathbf{z},\mathbf{y},A)|\Omega\Big{]}$, where

$\displaystyle\mathcal{Q}=\mathbb{E}_{\mathbf{z}}\Big{[}\log
p(z_{1})\Big{]}+\log p(A)+\mathbb{E}_{\mathbf{z}}\Big{[}\sum_{t=2}^{T}\log
p(z_{t}|z_{t-1},A)\Big{]}+\mathbb{E}_{\mathbf{z}}\Big{[}\sum_{t=1}^{T}\log
p(y_{t}|z_{t})\Big{]}$ (3)

In the E-step (Inference), we follow the backward algorithm in [8] to compute
$\mathbb{E}[z_{t}|\mathbf{y}]$,
$\mathbb{E}[z_{t}z_{t}^{{}^{\prime}}|\mathbf{y}]$ and
$\mathbb{E}[z_{t}z_{t-1}^{{}^{\prime}}|\mathbf{y}]$, which are sufficient
statistics of the expected log likelihood. In the M-step (Learning), we try to
find $\Omega$ that maximizes the likelihood lower bound $\mathcal{Q}$. In the
following, we derive the M-step for gradient based optimization of the
parameters $\Omega$. We omit the explicit conditioning on $\Omega$ for
notational brevity. Since the $\mathcal{Q}$ function is non-differentiable
with respect to _A_ , but differentiable with respect to all the other
variables ($C,R,Q,\pi_{1},V_{1}$ ), we separate the optimization into two
parts.

Optimization of $A$. In each iteration in the M-step, we need to maximize
$\mathbb{E}_{\mathbf{z}}\Big{[}\sum_{t=2}^{T}\log
p(z_{t}|z_{t-1},A)\Big{]}+\log p(A)$ with respect to _A_ , which is equivalent
to minimizing a function $f(A)$ that

$f(A)=\underbrace{\frac{1}{2}\sum_{t=2}^{T}\mathbb{E}_{\mathbf{z}}\Big{[}(z_{t}-Az_{t-1})^{\prime}Q^{-1}(z_{t}-Az_{t-1})\Big{]}}_{\text{g(A)}}+\underbrace{\beta||A||_{1}}_{\text{h(A)}}$
(4)

where $||A||_{1}$ is the $\ell_{1}$ norm on every element of matrix _A_ ,
$||A||_{1}=\sum_{i=1}^{l}\sum_{j=1}^{l}||A_{ij}||_{1}$.

As we can see $f(A)$ is convex but non-differentiable and we can easily
decompose $f(A)$ into two parts: $f(A)=g(A)+h(A)$, as shown in eq.(4). Since
$g(A)$ is differentiable, we can adopt the generalized gradient descent
algorithm to minimize $f(A)$. The update rule is:
$A^{(k+1)}=\mbox{prox}_{\alpha_{k}}(A^{(k)}-\alpha_{k}\bigtriangledown
g(A^{(k)}))$ where $\alpha_{k}$ is the step size at iteration _k_ and the
proximal function $\mbox{prox}_{\alpha_{k}}(A)$ is defined as the soft-
thresholding function $S_{\beta\alpha_{k}}(A)$

$[S_{\beta\alpha_{k}}(A)]_{ij}=\begin{cases}A_{ij}-\beta\alpha_{k}&\text{if
}A_{ij}>\beta\alpha_{k}\\\ 0&\text{if }-\beta\alpha_{k}\leq
A_{ij}\leq\beta\alpha_{k}\\\ A_{ij}+\beta\alpha_{k}&\text{if
}A_{ij}<-\beta\alpha_{k}\end{cases}$

######  Theorem 1.

Generalized gradient descent with a fixed step size $\alpha\leq
1/(||Q^{-1}||_{F}\cdot||\sum_{t=2}^{T}M_{t-1|T}||_{F})$ for minimizing eq.(4)
has convergence rate $O(1/k)$, where k is the number of iterations.

###### Proof.

$g(A)$ is differentiable with respect to _A_ , and its gradient is
$\bigtriangledown
g(A)=Q^{-1}(A\sum_{t=2}^{T}M_{t-1|T}-\sum_{t=2}^{T}M_{t,t-1|T})$. Using simple
algebraic manipulation we arrive at $||\bigtriangledown g(X)-\bigtriangledown
g(Y)||_{F}=||Q^{-1}(X-Y)\sum_{t=2}^{T}M_{t-1|T}||_{F}\leq||Q^{-1}||_{F}\cdot||\sum_{t=2}^{T}M_{t-1|T}||_{F}\cdot||X-Y||_{F}$
where $||\cdot||_{F}$ is the Frobenius norm and the inequality holds because
of the sub-multiplicative property of Frobenius norm. Since we know from
eq.(4), $f(A)=g(A)+\beta||A||_{1}$, and $g(A)$ has Lipschitz continuous
gradient with constant $||Q^{-1}||_{F}\cdot||\sum_{t=2}^{T}M_{t-1|T}||_{F}$,
according to [7, 22],
$f(A^{(k)})-f(A^{*})\leq||A^{(0)}-A^{*}||^{2}_{F}/2\alpha k$, where $A^{(0)}$
is the initial value and $A^{*}$ is the optimal value for $A$; _k_ is the
number of iterations. ∎

Theorem 1 gives us a simple way to set the step size during the generalized
gradient updates and also guarantees the fast convergence rate.

Optimization of $\Omega\backslash A=\\{C,R,Q,\pi_{1},V_{1}\\}$. Each of these
parameters is estimated similarly to the approach in [8] by taking the
corresponding derivative of the eq.(3), setting it to zero, and by solving it
analytically. Update rules for $\Omega\backslash A=\\{C,R,Q,\pi_{1},V_{1}\\}$
are shown in Algorithm 1.

The M-Step for optimizing $\Omega$ is summarized in Algorithm 1 in Appendix.

##  4 Experiments

We test our approach on time series data obtained from electronic health
records of 4,486 post-surgical cardiac patients stored in PCP database [2, 13,
12, 23]. To test the performance of our prediction model, we randomly select
600 patients that have at least 10 _Complete Blood Count_ (CBC) tests 111CBC
panel is used as a broad screening test to check for such disorders as anemia,
infection, and other diseases. ordered during their hospitalizations. The
three tests used in this experiment are Mean Corpuscular Hemoglobin
Concentration (MCHC), Mean Corpuscular Hemoglobin (MCH) and Mean Corpuscular
Volume (MCV). These time series data are noisy, their signals fluctuate in
time, and observations are obtained with varied time-interval period.

In order to get regularly sampled multivariate time series dataset, we apply
an 8-hour discretization on our original multivariate time series dataset and
use linear interpolation to fill the missing gaps from discretization. We
compare our sparse LDS (SLDS) with ordinary LDS (OLDS) on the above
multivariate time series dataset.

To evaluate the performance of our SLDS approach we split our time series for
600 patients into the training and testing sets, such that 50/100 times series
form the training data, and 500 are used for testing.

Evaluation Metric. Our objective is to test the predictive performance of our
approach by its ability to predict the future value of an observation for a
patient for some future time t given a sequence of patient’s past
observations. We judge the quality of the prediction using the Average Mean
Absolute Error (AMAE) on multiple test data predictions. More specifically,
the AMAE is defined as follows:

$AMAE=m^{-1}n^{-1}\sum_{i=1}^{m}\sum_{j=1}^{n}|y_{i,j}-\hat{y}_{i,j}|$ (5)

where $y_{i,j}$ is the _j_ th true observation from time series _i_ ,
$\hat{y}_{i,j}$ is the corresponding predicted value of $y_{i,j}$. _m_ is the
number of time series and _n_ is the length of each time series.

To conduct the evaluation, we use the test dataset to generate various
prediction tasks as follows. For each patient $p$ and complete time series _i_
for that patient, we calculate the number of observations $n_{i}^{p}$ in that
time series _i_. We use $n_{i}^{p}$ to generate different pairs of indices
$(\psi,\phi)$ for that patient, such that $1\leq\psi<\phi\leq n_{i}^{p}$,
where $\psi$ is the index of the last observation assumed to be seen, and
$\phi$ is the index of the observation we would like to predict. By adding
time stamp reading to each index, the two indices help us define all possible
prediction tasks that we can formulate on that time series. For each time
series _i_ from patient _p_ , we proceed by randomly picking 5 different pairs
of indices (or 5 different prediction tasks) for the total of 2500 predictions
tasks (500 x 5 = 2500). For each method, we repeat this random sampling
predictions 10 times and we use the Average Mean Absolute Error (AMAE) on
these tasks to judge the quality of test predictions. The prediction results
are shown in Table 1, Figure 2 and Figure 2.

From Figure 2 and Figure 2, we can see that OLDS achieves its lowest
prediction error _AMAE_ when the number of hidden states is 5. By varying the
number of states, the errors for OLDS first improve (till the optimal number
of states is reached) and then increase when the number of hidden states
exceeds the optimal point. This clearly shows the overfitting problem. The
SLDS performs similarly; its performance first impoves and after that it
deteriorates. However, its errors deteriorate at slower pace which shows it is
more robust to the overfitting problem. Comparing the two methods, the SLDS
always outperfroms the OLDS, indicating that the additional sparsity term
included in optimization helps it to better fit the underlying structure of
the transition matrix.

Table 1: Average mean absolute error for OLDS and SLDS with different hidden states sizes. # of states | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 12 | 15
---|---|---|---|---|---|---|---|---|---|---
OLDS(50) | 1.1653 | 1.0553 | 1.0561 | 0.6667 | 0.7268 | 0.9209 | 0.9502 | 0.9402 | 1.2291 | 1.2993
SLDS(50) | 0.6256 | 0.5858 | 0.5821 | 0.5684 | 0.6506 | 0.6200 | 0.8854 | 0.9236 | 1.0134 | 1.0430
OLDS(100) | 1.1527 | 1.0427 | 1.0039 | 0.6406 | 0.7153 | 0.8364 | 0.9210 | 0.9327 | 1.1427 | 1.2427
SLDS(100) | 0.5709 | 0.5429 | 0.5889 | 0.6379 | 0.6897 | 0.6949 | 0.7643 | 0.7811 | 0.7874 | 0.8309

Figure 1: AMAE on 50 training examples.

Figure 2: AMAE on 100 training examples.

##  5 Conclusion

In this paper, we have presented a sparse linear dynamical system (SLDS) for
multivariate time series predictions. Comparing with the traditional linear
state-space systems, SLDS model tries to (1) prevent the overfitting problem
and (2) represent additional structure in the transition matrix. Experimental
results on real world clinical data from electronic health records systems
demonstrated that this novel model achieves errors that is statistically
significantly lower than errors of ordinary linear dynamical system. We would
like to note that the results presented in this work are preliminary and
include only three time series. Further investigation of more complex higher
dimensional time-series data is needed and will be conducted in the future. In
addition, we would like to study group lasso regularization techniques which
we believe would be able to better control the dimensionality of the hidden
state space. Finally, we plan to study extensions of our model to switching-
state and controlled dynamical systems [14, 17].

Acknowledgement: This research work was supported by grants R01LM010019 and
R01GM088224 from the National Institutes of Health. Its content is solely the
responsibility of the authors and does not necessarily represent the official
views of the NIH. We would like to thank Eric Heim and Mahdi Pakdaman for
useful discussions and comments on this work.

## References

  * [1] Iyad Batal, Dmitriy Fradkin, James Harrison, Fabian Moerchen, and Milos Hauskrecht.  Mining recent temporal patterns for event detection in multivariate time series data.  In SIGKDD, pages 280–288, 2012.
  * [2] Iyad Batal, Hamed Valizadegan, Gregory F Cooper, and Milos Hauskrecht.  A pattern mining approach for classifying multivariate temporal data.  In BIBM, pages 358–365. IEEE, 2011.
  * [3] Iyad Batal, Hamed Valizadegan, Gregory F Cooper, and Milos Hauskrecht.  A temporal pattern mining approach for classifying electronic health record data.  TIST, Special Issue on Health Informatics, 2013.
  * [4] Adam Charles, Muhammad Salman Asif, Justin Romberg, and Christopher Rozell.  Sparsity penalties in dynamical system estimation.  In Information Sciences and Systems (CISS), 2011 45th Annual Conference on, pages 1–6. IEEE, 2011.
  * [5] Scott Shaobing Chen, David L Donoho, and Michael A Saunders.  Atomic decomposition by basis pursuit.  SIAM journal on scientific computing, 20(1):33–61, 1998.
  * [6] Alessandro Chiuso and Gianluigi Pillonetto.  Learning sparse dynamic linear systems using stable spline kernels and exponential hyperpriors.  In NIPS, pages 397–405, 2010.
  * [7] Massimo Fornasier and Holger Rauhut.  Iterative thresholding algorithms.  Applied and Computational Harmonic Analysis, 25(2):187–208, 2008\.
  * [8] Zoubin Ghahramani and Geoffrey E Hinton.  Parameter estimation for linear dynamical systems.  Technical report, Technical Report CRG-TR-96-2, University of Totronto, 1996.
  * [9] Bernard Ghanem and Narendra Ahuja.  Sparse coding of linear dynamical systems with an application to dynamic texture recognition.  In Pattern Recognition (ICPR), 2010 20th International Conference on, pages 987–990. IEEE, 2010.
  * [10] Yue Guan and Jennifer G Dy.  Sparse probabilistic principal component analysis.  In Proceedings of AISTATS, volume 5, pages 185–192, 2009.
  * [11] James Douglas Hamilton.  Time series analysis, volume 2.  Cambridge Univ Press, 1994.
  * [12] M. Hauskrecht, M. Valko, I. Batal, G. Clermont, S. Visweswaran, and G.F. Cooper.  Conditional outlier detection for clinical alerting.  In AMIA Annual Symposium Proceedings, volume 2010, page 286. American Medical Informatics Association, 2010.
  * [13] Milos Hauskrecht, Iyad Batal, Michal Valko, Shyam Visweswaran, Gregory F Cooper, and Gilles Clermont.  Outlier detection for patient monitoring and alerting.  Journal of Biomedical Informatics, 2012.
  * [14] Milos Hauskrecht and Hamish Fraser.  Modeling treatment of ischemic heart disease with partially observable markov decision processes.  Proceedings of the AMIA Symposium, page 538, 1998.
  * [15] Rudolph Emil Kalman et al.  A new approach to linear filtering and prediction problems.  Journal of basic Engineering, 82(1):35–45, 1960.
  * [16] T. Katayama.  Subspace methods for system identification.  Springer, 2005.
  * [17] Branislav Kveton and Milos Hauskrecht.  Solving factored mdps with exponential-family transition models.  In ICAPS, pages 114–120, 2006.
  * [18] Zitao Liu and Milos Hauskrecht.  Clinical time series prediction with a hierarchical dynamical system.  In Artificial Intelligence in Medicine, pages 227–237. Springer, 2013.
  * [19] Zitao Liu, Lei Wu, and Milos Hauskrecht.  Modeling clinical time series using gaussian process sequences.  In SIAM International Conference on Data Mining (SDM), pages 623–631, 2013.
  * [20] Lennart Ljung and Torkel Glad.  Modeling of dynamic systems.  1994\.
  * [21] P.V. Overschee, BLD Moor, D.A. Hensher, J.M. Rose, W.H. Greene, K. Train, W. Greene, E. Krause, J. Gere, and R. Hibbeler.  Subspace Identification for the Linear Systems: Theory–Implementation.  Boston: Kluwer AcademicPublishers, 1996.
  * [22] NZ Shor.  The rate of convergence of the generalized gradient descent method.  Cybernetics and Systems Analysis, 4(3):79–80, 1968.
  * [23] Michal Valko and Milos Hauskrecht.  Feature importance analysis for patient management decisions.  In MEDINFO, 2010.
  * [24] Peter M Williams.  Bayesian regularization and pruning using a laplace prior.  Neural computation, 7(1):117–143, 1995.

## Appendix

Algorithm 1 EM: M-step for the ($k+1$)th iteration. (We omit the explicit
superscript $(k+1)$ for notational brevity.)

INPUT:

  * •

Observation sequence $y_{t}$s, $t=1,\ldots,T$.

  * •

Sufficient statistics $\hat{z}_{t|T}$, $M_{t|T}$, $M_{t,t-1|T}$,
$i=1,\ldots,T$ from the ($k+1$)th iteration in E-Step.

PROCEDURE:

1: Update $\Omega\backslash A$:
$C=(\sum_{t=1}^{T}y_{t}\hat{z}_{t|T}^{{}^{\prime}})(\sum_{t=1}^{T}M_{t|T})^{-1}$,
$Q=\frac{1}{T-1}(\sum_{t=2}^{T}M_{t|T}-A\sum_{t=2}^{T}M_{t-1,t|T})$,
$R=\frac{1}{T}\sum_{t=1}^{T}(y_{t}y_{t}^{{}^{\prime}}-C\hat{z}_{t|T}y_{t}^{{}^{\prime}})$,
$\pi_{1}=\hat{z}_{1|T}$,
$V_{1}=M_{1|T}-\hat{z}_{1|T}\hat{z}_{1|T}^{{}^{\prime}}$.

2: Initialize $A$,
$A=(\sum_{t=2}^{T}M_{t,t-1|T})(\sum_{t=2}^{T}M_{t-1|T})^{-1}$.

3: repeat

4: Compute the fixed step size $\alpha$,
$\alpha=1/(||{Q}^{-1}||_{F}\cdot||\sum_{t=2}^{T}M_{t-1|T}||_{F})$.

5: Compute gradient of $g(A)$, $\bigtriangledown
g(A)={Q}^{-1}A\sum_{t=2}^{T}M_{t-1|T}-{Q}^{-1}\sum_{t=2}^{T}M_{t,t-1|T}$.

6: Update $A$, $A=S_{\beta\alpha}(A-t\bigtriangledown g(A))$.

7: until Convergence

OUTPUT: $\Omega^{(k+1)}=\\{A,C,Q,R,\boldsymbol{\pi}_{1},V_{1}\\}$.


# Interaction of oxygen with pristine and defective $\rm MoS_{2}$ monolayers

Murilo Kendjy Onita  Flavio Bento de Oliveira  Andréia Luisa da Rosa
<EMAIL_ADDRESS>Instituto de Física, Universidade Federal de Goiás, Campus
Samambaia, 74690-900, Goiânia, Goiás, Brazil.

###### Abstract

Atom controlled sub-nanometer MoS2 pores have been recently fabricated.
Oxidative environments are of particular interest for MoS2 applications in
electronics, sensing and energy storage. In this work we carried out first-
principles calculations of oxygen adsorption in plain and sub-nanometer MoS2
nanopores. The chemical stability of the layers and pores towards oxygen was
verified using density-functional theory. Dissociation and diffusion barriers
have been calculated in order to understand surface and pore oxidation and its
electronic properties at the atomic scale.

##  I Introduction

Owing to their fascinating properties, two dimensional transition metal
dichalcogenides (TMDs) have been explored for a variety of applications,
including electronics and optoelectronics, photonics, catalysis and energy
storage Feng et al. (2015, 2016a, 2016b); Karmodak et al. (2021); Chamlagain
and Khondaker (20x). In particular, molybdenum disulfide (MoS2), the most
promising TMDCs, is efficiently exfoliated in monolayer or multilayersSantosh
et al. (2015); Tsai et al. (2017). Recently the fabrication of MoS2 sub-
nanometer pores offer several opportunities being promising candidates for
several technological applications such as membranes for DNA
translocation,Feng et al. (2015); Sen and Gupta (2021); Graf et al. (2019a),
water filtration and desalination Cao et al. (2020); Kou et al. (2016);
Heiranian et al. (2015); Wang et al. (2021), energy harvesting Graf et al.
(2019b) and hydrogen evolution reaction Wu et al. (2019); Li et al. (2019); Xu
et al. (2022).

Control at atomic scale nanopores can be achieved by using electron or ion
beams to fabricate nanopores of specific sizes Chen et al. (2019); Shi et al.
(2017); Garoli et al. (2019); Lin et al. (2014); Vancsó et al. (2016). Sulfur
and molibdenum vacancies, the most abundant defect species in MoS2, may serve
as the nucleation sites for the nanopore formation Tsai et al. (2017); Trainer
et al. (2022). Pore sizes down to 0.5-1.2 nm have been created, corresponding
to a few atoms missing Zhao et al. (2018); Wang et al. (2017); Thiruraman et
al. (2020); Yin et al. (2018).

In particular, the interaction of small molecules with MoS2 is of particular
importance, since it plays a role in the properties and performance of two-
dimensional-based devices. Concurrently, structural defects such as vacancies
and edges are particularly susceptible to attacks in reactive environments
Gruenleitner et al. (2022). Adsorption energy and diffusion of oxygen on MoS2
is controversial and not completely understood. Moreover, oxygen on surfaces
either forms undesired alloys entering a sulphur site or possesses a low
binding energy. One important, not completely understood question, is the role
of point defects on the adsorption of oxygen on MoS2. Previous experimental
investigations suggested that the presence of defects significantly alters the
surface stability Santosh et al. (2015); Longo et al. (2017). Oxygen leads to
the formation of Mo oxide in the two-dimensional MoS2 lattice, yielding an
overall disordered and fragmented structure. On the other hand, other
investigations suggested that oxygen is incorporated in the MoS2 lattice at
substitutional sites Jeon et al. (2015). X-ray photoelectron spectroscopy
measurements have shown that the basal plane of MoS2 monolayers, when
subjected to long-term ambient exposure, spontaneously undergoes oxygen
substitution reactions, giving rise to a highly crystalline two-dimensional
molybdenum oxysulfide phase Petö et al. (2018); Chhowalla et al. (2013).

In this work we have investigated the interaction of oxygen atoms and
molecules with MoS2 monolayers and subnanometer pores. The pore size and
termination play a crucial role on the interaction with oxygen molecules.

##  II Methodology

All calculation have been performed using the density functional theory
Hohenberg and Kohn (1964); Kohn and Sham (1965) as implemented in the VASP
Kresse and Joubert (1999). For the exchange-correlation functional we use the
GGA approximation Perdew et al. (1996). Basis set consisting of an expansion
in plane waves and an energy cutoff of 300 eV were used. The Brillouin zone
was sampled according to the Monkhorst-Pack scheme Monkhorst and Pack (1976)
with $\Gamma$ point only in a $(6\times 6)$ unit cells. To avoid interaction
between neighbouring wires, a vacuum space of 15 Å was set in the direction
perpendicular to the two-dimensional MoS2 layers. The structural relaxation
was performed until the atomic forces were less than 10-3 eV/Å. In order to
determine the diffusion and reaction barriers we have performed CI-
NEBHenkelman (2000a, b) calculations with at least seven images to search the
minimum-energy reaction paths and saddle points between the initial state and
final state configurations.

##  III Results

Among the various arrangements for MoS2, the energetically most stable one at
room temperature is known as 2H which has a honeycomb structure as shown in
Fig. 1 (a) and (b). This arrangement has in its unit cell one molybdenum and
two sulfur atoms. The optimized Mo-S (S-S) distances are 2.42 (3.14) Å. The
lattice parameter $a$ = $b$ equals to 3.16 Å. The electronic band structure
and total density-of-states (DOS) of bare a MoS2 monolayer are shown in Fig.
1(b). The direct band gap at K–point is 1.76 eV, in agreement with previous
results Komsa and Krasheninnikov (2015).

The electronic band gap and formation energies of bare and defective
monolayered MoS2 are reported in our previous publication Kendjy et al.
(2022). The formation enthalpy of a bare MoS2 monolayer is -2.63 eV.

Figure 1:  a) Top view, b) side view and c) electronic band structure and DOS
of bare MoS2 monolayers.

Here the adsorption energy of oxygen is calculated at high oxygen chemical
potential according to: ${\rm
E_{b}=E_{MoS_{2}/O}-E_{MoS_{2}-bare}-\sum_{i}{\mu_{O}}}$, where ${\rm
E_{MoS_{2}/O}}$ is the total e nergy of the MoS2 the with oxygen adsorbed,
${\rm E_{MoS_{2}}^{bare}}$ is the total energy of the bare sheet, $\mu_{\rm
O}$ is the chemical potential of oxygen, which has ben chosen as the total
energy of the oxygen molecule.

Figure 2:  Relaxed sructure of oxygen substitutional in ${\rm MoS_{2}}$
nanopores. a) OS (11%), b) OS (22%), c) OS (25%), d) OS (11%) + VMo, e) OS
(11%) + VS (far) and f) OS (11%) + VS (close). Purple, yellow and red are Mo,
S and O atoms, respectively. The oxygen concentration is given within
brackets.

As possible configurations for defective MoS2, we have considered oxygen
substituting sulphur atoms, namely OS, with concentrations of 11%, 22% and 25%
as shown in Figs.2 (a), (b) and (c), respectively. We notice that substitution
at Mo sites is highly unfavorable. Furthermore, defects have been incorporated
considering Mo vacancies namely OS (11%) + VMo, Fig.2 (d) and S vacancies,
Figs. 2 (e) and (f) labeled as OS (11%) + VS (far) and OS (11%) + VS (close).
These concentrations were chosen based on experimental results reported in
Ref. Tang et al. (2020). We have included disorder effects by adding Mo
vacancies far and close to oxygen substitutional atoms, as mentioned above.

Figure 3:  DOS and PDOS of OS in $\rm MoS_{2}$ nanosheets at different oxygen
concentrations: a) OS at (11%), b) OS at (22%), c) OS (25%), d) O${}_{\rm
S}^{\rm close}$ \+ VMo (11%), e) O${}_{\rm S}^{\rm far}$ \+ VMo (11%), f)
O${}_{\rm S}^{\rm close}$ \+ VS (11%).

In Fig. 3 the DOS and PDOS of the oxygen substitutional in MoS2 is shown. a)
OS at (11%), b) OS at (22%), c) OS (25%), d) O${}_{\rm S}^{\rm close}$ \+ VMo
(11%), e) O${}_{\rm S}^{\rm far}$ \+ VMo (11%), f) O${}_{\rm S}^{\rm close}$
\+ VS (11%). The electronic structure of a) OS at (11%), b) OS at (22%) is
very similar. The later two systems exhibit semiconductor behavior.

At higher concentration as in OS (25%), Fig.3(c), additional state in the gap
appear, but the structure remain semiconducting. For d) O${}_{\rm S}^{\rm
close}$ \+ VMo (11%) and O${}_{\rm S}^{\rm far}$ \+ VMo (11%) states cross the
Fermi level, indicating a metalic behavior. Finally, for O${}_{\rm S}^{\rm
close}$ \+ VS (11%) shown in Fig. 3(f) the system is a semiconductor as well.

Table 1: Binding energies of oxygen molecules in ${\rm MoS_{2}}$ nanopores. Structure | Eb (eV)
---|---
pristine | -0.55
V${}_{\rm 1Mo-O_{2}}$ (non-diss.) | -0.37
V${}_{\rm 1S-O_{2}}$ (diss.) | -0.90
V${}_{\rm 2S-O_{2}}$ (diss.) | -3.44
V${}_{\rm 1Mo6S-3O_{2}}$(diss.) | -3.24
V${}_{\rm 1Mo6S-O_{2}}$ (non-diss.) | -0.45

Figure 4:  Top and rotated views of relaxed sructures of oxygen adsorption in
${\rm MoS_{2}}$ monolayers. a) and e) pristine, b) and f) V${}_{\rm 1Mo-
O_{2}}$, c) and g) V${}_{\rm 1Mo6S-3O_{2}}$ (dissociated) and d) and h)
V${}_{\rm 1Mo6S-O_{2}}$ (non-dissociated) defects.

Figure 5:  DOS and PDOS of oxygen adsorbed in ${\rm MoS_{2}}$ monolayers: a)
pristine and with b) ${\rm V_{1Mo-O_{2}}}$, c) ${\rm V_{1Mo6S-3O_{2}}}$
(dissociated), d) V${}_{\rm 1Mo6S-O_{2}}$ (non-dissociated) defects.

On the monolayer pristine the oxygen molecule spontaneously dissociates and
the oxygen atoms diffuses further apart to finally bind on the sulfur atoms.
This means that the barrier for O diffusion should be low, as suggested in
Ref. Rao et al. (2017). The binding energy of the final configuration is -1.13
eV/atom. In addition we have investigated the interaction of O2 with other
small defects which we call subnanometer pores: a Mo single vacancy, V${}_{\rm
Mo-O_{2}}$, shown in Fig. 4 (b), a Mo single vacancy plus S hexavacancy,
$V_{\rm 1Mo6S-O_{2}}$, shown in Fig. 4 (c) and d) V${}_{\rm 1Mo6S-O_{2}}$
(non-dissociated). Such structures can serve as prototypes for nanopores in
MoS2. Additionally, we have considred a Mo triple vacancy plus a S divacancy,
V${}_{\rm 3Mo2S-O_{2}}$, and a single Mo vacancy plus a S tetravacancy,
V${}_{\rm 1Mo4S-O_{2}}$, a Mo triple vacancy plus a S hexavacancy, V${}_{\rm
3Mo6S-O_{2}}$ (to be published).

Binding energies of oxygen for the relaxed structures are shown in Table 1. On
pristine MoS2, shown in Fig. 4 (a), adsorption of oxygen leads to a
dissociative configuration with binding energy is -0.55 eV/O atom. Oxygen
spontaneously dissociates with the oxygen atoms diffusing further apart with a
6.3 Å being the smallest O-O distance. The O sits right above the S atom with
S-O distance equal to 1.48 Å.

On V${}_{\rm 1Mo-O_{2}}$ oxygen chemisorbs in a non-dissociate manner with
binding energy of -0.37 eV/O atom. The O-O distance is 2.78 Å with each oxygen
bound to a sulphur atom in different layers. The relaxed geometry is shown in
Fig. 4 (b).

In a single sulfur vacancy defect, V${}_{\rm 1S-O_{2}}$, shown in Fig. 4 (c),
the oxygen molecule dissociates. Upon dissociation, the first atom goes to a
substitutional sulphur site, while the second oxygen atom moves to ontop of a
sulphur atom. This is in agreement with results reported previously Wang and
Paulus (2021). The binding energy is -0.90 eV/O atom.

On a double sulphur vacancy, V${}_{\rm 2S-O_{2}}$ both oxygen atoms go to
substitutional sulphur positions. The Mo-O bond length is 1.7 Å forming a
slightly distorted hexagon. One notes that typical Mo-O bond lengths range
from 1.69 to 1.73 Å Hardcastle and Wachs (1990). The O-O distance is our case
2.38 Å. The binding energy amounts -3.44 eV/atom, meaning that the reaction is
exothermic.

The V${}_{\rm 1Mo6S-O_{2}}$ interacts with the oxygen molecule via a
perpendicular orientation of the molecule with respect to the MoS2 surface in
its final configuration. The oxygen distance to the molybdenum atoms is around
3 Å.

For other large nanopores such as V${}_{\rm 3Mo6S-O_{2}}$, V${}_{\rm
1Mo4S-O_{2}}$ and V${}_{\rm 3Mo2S-O_{2}}$ (to be published)), the binding
energy indicates physisorption, with the O2 molecule lying in the middle of
the pore with O-O bond distance equals to 1.0 Å. We may argue here that there
is an interplay between the number of oxygen atoms and the number of dangling
bonds in the pore.

The density of states of the above discussed nanopores are show in Fig. 5.
V${}_{\rm 1S-O_{2}}$ shows states within the gap, as seen in Fig. 5(a). On the
other hand, V${}_{\rm 2S-O_{2}}$ does not have states in the gap as shown in
Fig. 5(b). Furthermore, V${}_{\rm 1Mo-O_{2}}$ and V${}_{\rm 1Mo6S-O_{2}}$ have
states within the gap due to sulphur dangling bonds.

Recently it has been suggested that the formation of oxygen at substitutional
sites with no band gap states can be formed in MoS2 layers Tang et al. (2020).
A possible explanation for the ordered oxygen incorporation is that due to the
higher strength of the Mo-O bonds compared to M-S, the substitutional
oxidation of the 2D MoS2 basal plane could be also thermodynamically
favourable Longo et al. (2017); Martincová et al. (2017). One can explain the
dissociation considering that SO2 has cohesive energy of -3.1 eV, which is
larger than the formation enthalpy of MoS2 which is -2.48 eV.

Figure 6:  Diffusion barrier (top panel) and corresponding paths (lower panel)
of an oxygen molecule on a molibdenum vacancy in MoS2 with perpendicular
orientation to the MoS2 basal plane.

Experiments reveal relatively low activation energy values for MoS2 oxidation,
ranging from 0.54 eV to 0.98 eV Rao et al. (2017). Recent theoretical
calculations Wang and Paulus (2021) report a barrier of 0.33 eV in the
presence of sulphur defects. Here we considered the presence of a molibdenum
vacancy. The reaction path is done by following the molecule along its minimum
energy path. We have used between five and eight intermediate states to
calculate the diffusion. While oxygen atom diffuses through the pore with a
barrier of 2.53 eV. On the other hand, an oxygen molecule diffuses barrierless
along the path shown in Fig. 6 respectively.

Figure 7:  Projected charge density on LUMO-1 of a) V${}_{\rm 1Mo-O_{2}}$, b)
V${}_{\rm 1Mo6S-O_{2}}$ and c) V${}_{\rm 3Mo6S-O_{2}}$. Isosurface value is
0.003 e/bohr3

In order to have further insight on the interaction between the oxygen
molecules and the nanopores we have calculated the projected partial charge
density on the first level above the Fermi level. These are shown in Fig. 7
(a) for V${}_{\rm 1Mo-O_{2}}$ which has the unnocupied orbitals of the pore
oriented towards the molecule. However, as the preferred orientation of the
molecule in the middle of the pore is perpendicular to the pore, this could
explain the small overlap between HOMO of the molecule and LUMO of the pore
and therefire the relative small berrier for V${}_{\rm 1Mo-O_{2}}$ where the
molecule easily diffuses through the pore.

For the V${}_{\rm 1Mo-O_{2}}$ the molecule diffuses without barrier.
Furthermore, the V${}_{\rm 3Mo6S-O_{2}}$ pore is large enough for the molecule
to diffuse either parallel or perpendicular to the pore. However, if the
molecule is close enough to the edges, it is likely to dissociate. On the
other hand, the V${}_{\rm 1Mo6S-O_{2}}$ pore have its Mo orbitals hybridizing
forming metallic-like bonds.

##  IV Conclusions

In this work we carried out first-principles calculations of oxygen on
monolayered MoS2 and sub-nanometer MoS2 nanopores. The dissociation and
diffusion of oxygen in MoS2 reveal two main features: the orientation of the
molecule with respect to the surface plays a role on the molecule diffusion
barrier; the reactivity of the pore plays a role on the molecule dissociation.

##  V Acknowledgements

We acknowledge financial support from Conselho Nacional de Pesquisa e
Desenvolvimento (CNPq) under grants 313081/2017-4 and 305335/2020-0. The
calculations have been performed using the computational facilities of Santos
Dumont Supercomputar at LNCC, CENAPAD at Unicamp and LaMCAD/UFG.

## References

  * Feng et al. (2015) J. Feng, K. Liu, R. D. Bulushev, S. Khlybov, D. Dumcenco, A. Kis, and A. Radenovic, Nat. Nanotechnol. 10, 1070 (2015).
  * Feng et al. (2016a) J. Feng, M. Graf, K. Liu, and D. Ovchinnikov, Nature 536, 7615 (2016a).
  * Feng et al. (2016b) J. Feng, M. Graf, and K. e. a. Liu, Nature 536, 197 (2016b).
  * Karmodak et al. (2021) N. Karmodak, L. Bursi, and O. Andreussi, The Journal of Physical Chemistry Letters 13, 58 (2021).
  * Chamlagain and Khondaker (20x) B. Chamlagain and S. I. Khondaker, ACS Omega 6, 24075–24081 (20x).
  * Santosh et al. (2015) K. C. Santosh, R. C. Longo, R. M. Wallace, and K. Cho, J. Appl. Phys. 117, 135301 (2015).
  * Tsai et al. (2017) C. Tsai, H. Li, S. Park, J. Park, H. S. Han, J. K. Norskov, X. Zheng, and F. A.-Pedersen, Nat. Commun. 8, 15113 (2017).
  * Sen and Gupta (2021) P. Sen and M. Gupta, RSC Advances 11, 6114 (2021).
  * Graf et al. (2019a) M. Graf, M. Lihter, D. Altus, S. Marion, and A. Radenovic, Nano Lett. p. 9075 (2019a).
  * Cao et al. (2020) Z. Cao, V. Liu, and A. B. Farimani, ACS Energy Lett. 5, 2217 (2020).
  * Kou et al. (2016) J. Kou, J. Yao, L. Wu, X. Zhou, H. Lu, F. Wu, and J. Fan, Phys. Chem. Chem. Phys. 18, 22210 (2016).
  * Heiranian et al. (2015) M. Heiranian, A. B. Farimani, and N. R. Aluru, Nat. Comm. 6, 8616 (2015).
  * Wang et al. (2021) L. Wang, Z. Wang, S. K. Patel, S. Lin, and M. Elimelech, ACS Nano 15, 4093 (2021).
  * Graf et al. (2019b) M. Graf, M. Lihter, D.Unuchek, A. Sarathy, J.-P. Leburton, A. Kis, and A. Radenovic, Joule 3, 1549 (2019b).
  * Wu et al. (2019) L. Wu, A. Longo, N. Y. Dzade, A. Sharma, M. M. R. M. Hendrix, A. A. Bol, N. H. de Leeuw, E. J. M. Hensen, and J. P. Hofmann, ChemSusChem 12, 4383 (2019).
  * Li et al. (2019) W. Li, G. Liu, J. Li, Y. Wang, L. Ricardez-Sandoval, Y. Zhang, and Z. Zhang, Appl. Surf. Sci. 498, 143869 (2019).
  * Xu et al. (2022) J. Xu, G. Shao, X. Tang, F. Lv, H. Xiang, C. Jing, S. Liu, S. Dai, Y. Li, J. Luo, et al., Nature Communications 13, 219 (2022).
  * Chen et al. (2019) C.-H. Chen, X. Chang, and C.-S. Wu, Sci. Rep. 9, 18663 (2019).
  * Shi et al. (2017) W. Shi, A. K. Friedman, and L. A. Baker, Anal. Chem 89, 157 (2017).
  * Garoli et al. (2019) D. Garoli, H. Yamazaki, N. Maccaferri, and M. Wanunu (2019).
  * Lin et al. (2014) J. Lin, O. Cretu, W. Zhou, K. Suenaga, D. Prasai, K. I. Bolotin, N. Y. Cuong, M. Otani, S. Okada, A. R. Lupini, et al., Nat. Nanotech. (2014).
  * Vancsó et al. (2016) P. Vancsó, G. Z. Magda, J. Petö, J.-Y. Noh, Y.-S. Kim, C. Hwang, and dL. Tapasztó, Sci. Rep. 6, 29726 (2016).
  * Trainer et al. (2022) D. J. Trainer, J. Nieminen, F. Bobba, B. Wang, X. Xi, A. Bansil, and M. Iavarone, 2D Materials and Applications 6, 13 (2022).
  * Zhao et al. (2018) X. Zhao, D. Fu, Z. Ding, Y. Zhang, D. Wan, S. J. R. Tan, Z. Chen, K. Leng, J. Dan, W. Fu, et al., Nano Lett. 18, 482 (2018).
  * Wang et al. (2017) S. Wang, H. Li, H. Sawada, C. S. Allen, A. I. Kirkland, J. C. Grossman, and J. H. Warner, Nanoscale 9, 6417 (2017).
  * Thiruraman et al. (2020) J. P. Thiruraman, P. M. Das, and M. Drndicć, ACS Nano 14, 11831 (2020).
  * Yin et al. (2018) K. Yin, S. Huang, X. Chen, X. Wang, J. Kong, Y. Chen, and J. Xue, ACS Appl. Mater. Interfaces p. 28909 (2018).
  * Gruenleitner et al. (2022) T. Gruenleitner, A. Henning, M. Bissolo, M. Zengerle, L. Gregoratti, M. Amati, P. Zeller, J. Eichhorn, A. V. Stier, A. W. Holleitner, et al., ACS Nano 16, 20364–20375 (2022).
  * Longo et al. (2017) R. C. Longo, R. Addou, K. C.Santosh, J.-Y. Noh, C. M. Smyth, D. Barrera, C. Zhang, J. W. P. Hsu, R. M. Wallace, and K. Cho, 2D Materials 4, 025050 (2017).
  * Jeon et al. (2015) M. H. Jeon, C. Ahn, H. Kim, K. N. Kim, T. Z. Li, H. Qin, Y. Kim, S. Lee, T. Kim, and G. Y. Yeo, Nanotechnology 26, 355706 (2015).
  * Petö et al. (2018) J. Petö, T. Ollár, P. Vancsó, Z. I. Popov, G. Z. Magda, G. Dobrik, C. Hwang, P. B. Sorokin, and L. Tapasztó, Nature Chem. 10, 1246 (2018).
  * Chhowalla et al. (2013) M. Chhowalla, H. S. Shin, G. Eda, L.-J. Li, K. P. Loh, and H. Zhang, Nature Chem. 5, 263 (2013).
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Perdew et al. (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Monkhorst and Pack (1976) H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
  * Henkelman (2000a) H. Henkelman, G.; Jónsson, J. Chem. Phys. 113 (2000a).
  * Henkelman (2000b) H. Henkelman, G.; Jónsson, J. Chem. Phys. 113, 9978–9985 (2000b).
  * Komsa and Krasheninnikov (2015) H.-P. Komsa and A. V. Krasheninnikov, Phys. Rev. B (2015).
  * Kendjy et al. (2022) M. Kendjy, A. L. da Rosa, and T. Frauenheim, 34, 044005 (2022).
  * Tang et al. (2020) J. Tang, Z. Wei, Q. Wang, Y. Wang, B. Han, X. Li, B. Huang, M. Liao, J. Liu, N. Li, et al., Small 16, 2004276 (2020).
  * Rao et al. (2017) R. Rao, A. E. Islam, P. M. Campbell, E. M. Vogel, and B. Maruyama, 2D materials 4, 025058 (2017).
  * Wang and Paulus (2021) K. Wang and V. Paulus, J. Phys. Chem. C 125, 19544–19550 (2021).
  * Hardcastle and Wachs (1990) F. D. Hardcastle and I. E. Wachs, J. Raman Spectroscopy 21, 683 (1990).
  * Martincová et al. (2017) J. Martincová, M. Otyepka, and P. Lazar, 23, 13233 (2017).


11institutetext: School of Biomedical Engineering & Imaging Sciences, KCL,
London, UK
11email<EMAIL_ADDRESS>22institutetext: Athinoula A. Martinos
Center for Biomedical Imaging,
MGH and HMS, Boston, USA
33institutetext: Wellcome Center for Human Neuroimaging, UCL, London, UK
44institutetext: Institute of Cognitive Neuroscience, UCL, London, UK
55institutetext: Institute of Neurology, UCL, London, UK

# Fitting Segmentation Networks on Varying Image Resolutions using Splatting

Mikael Brudfors 11 Yaël Balbastre 22 John Ashburner 33 Geraint Rees 44
Parashkev Nachev 55 Sébastien Ourselin 11 M. Jorge Cardoso 11

###### Abstract

Data used in image segmentation are not always defined on the same grid. This
is particularly true for medical images, where the resolution, field-of-view
and orientation can differ across channels and subjects. Images and labels are
therefore commonly resampled onto the same grid, as a pre-processing step.
However, the resampling operation introduces partial volume effects and
blurring, thereby changing the effective resolution and reducing the contrast
between structures. In this paper we propose a _splat layer_ , which
automatically handles resolution mismatches in the input data. This layer
pushes each image onto a _mean space_ where the forward pass is performed. As
the splat operator is the adjoint to the resampling operator, the mean-space
prediction can be pulled back to the native label space, where the loss
function is computed. Thus, the need for explicit resolution adjustment using
interpolation is removed. We show on two publicly available datasets, with
simulated and real multi-modal magnetic resonance images, that this model
improves segmentation results compared to resampling as a pre-processing step.

###### Keywords:

Image segmentation Splatting Resampling Pre-processing Image resolution.

##  1 Introduction

Automatic semantic segmentation of medical images is widely done using deep-
learning-based segmentation networks. To apply these networks, a pre-
processing step that resamples all images into the same space is currently
performed, as the images can have different orientation, field-of-view and
resolution. Choosing the common space can be done in many ways, _e.g._ , based
on the median voxel size of the training population [1]. This step is required
for stacking channel dimensions when working on multi-modal(channel) data, but
also if a batch size larger then one is required. This type of pre-processing
is performed in the majority of biomedical challenges, _e.g_ , BRATS [2],
Medical Segmentation Decathlon [3] and the WMH Segmentation Challenge [4].

Pre-processing images by resampling to a common space can be seen as a
normalisation step, intended to decrease data variance and facilitate both
model fitting and generalisability. However, resampling introduces values not
present in the original image through interpolation. Furthermore, it has a
smoothing effect that, unless coordinates fall exactly at voxel centres,
reduces the observed noise variance: let $y=ax_{1}+(1-a)x_{2}$ be the
interpolation of two values $x_{1}\sim\mathcal{N}(m_{1},~{}v)$ and
$x_{2}\sim\mathcal{N}(m_{2},~{}v)$, with $a\in[0,1]$ ; then,
$\text{Var}[y]=(1-2a(1-a))v\leq v$. In addition, interpolation algorithms do
not embed prior knowledge about the objects being interpolated, resulting in
overly smooth images that can bias analyses and cause false positives [5, 6].
These limitations can make it challenging to generalise segmentation networks
to a wide array of voxel sizes. Furthermore, the networks have no way to know
how confident they should be about a particular voxel value, _i.e._ , whether
it has been highly interpolated or preserves the raw value. This could be
particularly problematic when working with routine clinical MRIs, where thick-
sliced (high in-plane, low through-plane resolution), multi-modal acquisitions
are the default.

The simplest method for fusing modalities of different image resolution is
perhaps to fit separate networks to each modality and then combine their
individual predictions. This can be done by integrating multi-modal features
from the top layer of each network [7, 8], or by fusing their outputs via
averaging and majority voting [9]. However, such output-fusion strategies
learn only modality-specific features and ignore the complex interactions
between modalities. To better account for correlations across modalities,
several layer-level fusion methods have been proposed. For example, Hi-Net
[10], which learns individual modality features via a modality-specific
network and a layer-wise fusion strategy, or HyperDense-Net [11], which
employs dense connections across multi-modal streams. However, none of these
methods model the fact that the images are defined on different grids in their
native spaces. SynthSeg [12], on the other hand, introduced a convolutional
neural network (CNN) that learns a mapping between multi-modal scans, defined
on different grids, by simulating high-resolution scans from the training data
by interpolating the low-resolution scans to 1 mm isotropic. Since the
interpolation is simulated from the training data, the network becomes robust
to variable image resolutions. The method presented in this paper would avoid
interpolation, instead using the proposed splat layers. Finally, CNN-based
models exist that take irregularly sampled inputs [13], but they are currently
not easily extended to multi-modal data.

In this paper, we propose a method for directly fitting segmentation networks
to the raw image data. Our method is based on the splatting operation, which
pushes images, across subjects and channels, into a mean space. The network
produces its predictions in this space, which are then pulled back into the
native space of each input subject where the loss function is computed. If
multiple modalities are provided, the loss is computed on the image on which
the target segmentation was annotated. The splat layer avoids interpolating
the input data, allowing the network to instead infer on the raw voxel values.
We validate our proposed method on two semantic segmentation tasks, on
publicly available multi-modal brain MR images. Our validation shows that
extending a UNet with our proposed splat layers gives improved segmentations,
compared to fitting to data that have been pre-processed by resampling. Our
implementation uses a PyTorch backend, with custom splatting and resampling
layers written in C++/CUDA, publicly available at
https://github.com/balbasty/nitorch.

##  2 Methods

The idea of our method is quite simple; when fitting a segmentation network,
instead of resampling the input images as a pre-processing step, we instead
add two new layers to the network, one at its head and one at its tail. These
two layers are based on the splatting operation [14], which allow the network
to learn on the raw voxel data of the images. This avoids interpolation that
could introduce partial-volume effects, and for the loss to be computed on the
native space data. The idea is that the network implicitly interpolates the
data whilst training. To conceptualise the idea of splatting, we next show a
simple 1D toy example. The methodology is then extended to $D$-dimensional
input tensors in the subsequent section.

###  2.1 1D Toy Example

Let’s assume we have a training population of $M$ sets of native-space input
vectors ($D=1$), where each set of input vectors represents $C$ channels and
can be of different length (_i.e._ , it is not always possible to stack the
$C$ vectors of training data). For training a segmentation network, we want to
be able to concatenate all input vectors across $C$ and $M$ onto a common grid
(_i.e._ , having equal length). Let us define one of these vectors as
$\mathbf{f}_{mc}=[10,~{}11,~{}12,~{}13]^{\text{T}}$, with $N_{mc}=4$ elements
and the affine mapping111For medical images, the affine mapping can be read
from the image header. In general, it can be defined from knowledge of
orientation, pixel size and field-of-view.
$\mathbf{A}_{mc}=\begin{pmatrix}2.5&0\\\ 0&1\end{pmatrix}$. This vector’s
identity grid is given by $\mathbf{i}_{mc}=[0,~{}1,~{}2,~{}3]^{\text{T}}$. To
resize the input vector on a common training grid, we define its length
($N_{t}$) and affine mapping ($A_{t}$). This could be done in a number of
ways, in this paper we use a _mean space_. The mean space is defined later; in
this example, for simplicity, we assume $N_{t}=8$ and
$\mathbf{A}_{t}=\begin{pmatrix}1&0\\\ 0&1\end{pmatrix}$.

There are two ways of resizing the input vector onto the mean-space grid. The
standard method is resampling, in which we define an identity grid in the mean
space: $\mathbf{i}_{t}=[0,~{}1,~{}2,~{}3,~{}4,~{}5,~{}6,~{}7]^{\text{T}}$. We
then compose the affine mappings as
$\mathbf{A}=\mathbf{A}_{mc}^{-1}\mathbf{A}_{t}=\begin{pmatrix}0.4&0\\\
0&1\end{pmatrix}$ and transform the identity grid with $\mathbf{A}$ to get the
deformation
$\bm{\phi}_{t}=[0,~{}0.4,~{}0.8,~{}1.2,~{}1.6,~{}2,~{}2.4,~{}2.8]^{\text{T}}$.
Using $\bm{\phi}_{t}$ we can then pull values from $\mathbf{f}_{mc}$ onto a
grid $\mathbf{f}_{t}$ using some form of interpolation. With linear
interpolation we get
$\mathbf{f}_{t}=[10,~{}10.4,~{}10.8,~{}11.2,~{}11.6,~{}12,~{}12.4,~{}12.8]^{\text{T}}$.
This operation can be conceptualised as multiplying $\mathbf{f}_{mc}$ with a
sparse matrix $\mathbf{f}_{t}=\mathbf{\Psi}_{t}\mathbf{f}_{mc}$, where:

$\displaystyle\mathbf{\Psi}_{t}=\begin{bmatrix}1.0,&0.0,&0.0,&0.0\\\
0.6,&0.4,&0.0,&0.0\\\ 0.2,&0.8,&0.0,&0.0\\\ 0.0,&0.8,&0.2,&0.0\\\
0.0,&0.4,&0.6,&0.0\\\ 0.0,&0.0,&1.0,&0.0\\\ 0.0,&0.0,&0.6,&0.4\\\
0.0,&0.0,&0.2,&0.8\end{bmatrix}$

Figure 1: 1D examples of expansive and compressive splatting.

Splatting, on the other hand, does not interpolate the values in
$\mathbf{f}_{mc}$. Instead the affine mapping is defined as
$\mathbf{A}=\mathbf{A}_{t}^{-1}\mathbf{A}_{mc}=\begin{pmatrix}2.5&0\\\
0&1\end{pmatrix}$ and the identity grid of the input vector
($\mathbf{i}_{mc}$) transformed as
$\bm{\phi}_{mc}=[0,~{}2.5,~{}5,~{}7.5]^{\text{T}}$. This deformation is then
used to push each value in $\mathbf{f}_{mc}$ to a location in the mean-space
grid, giving us
$\mathbf{f}_{t}=[10,~{}0,~{}5.5,~{}5.5,~{}0,~{}12,~{}6.5,~{}6.5]^{\text{T}}$.
As multiple values can be pushed to the same location, a count image is also
computed
$\mathbf{c}_{t}=[1,~{}0,~{}0.5,~{}0.5,~{}0,~{}1,~{}0.5,~{}0.5]^{\text{T}}$,
which we additionally provide as input to the segmentation network. Note that
the splatting is expansive in this example, as there are fewer input voxels
than output voxels, so the count image has lots of zeros. Splatting can also
be compressive, meaning that multiple values can be pushed onto the same
output voxel. As resampling, splatting can be conceptualised as a matrix
multiplication $\mathbf{f}_{t}=\mathbf{\Lambda}_{t}\mathbf{f}_{mc}$, where:

$\displaystyle\mathbf{\Lambda}_{t}=\begin{bmatrix}1.0,&0.0,&0.0,&0.0,&0.0,&0.0,&0.0,&0.0\\\
0.0,&0.0,&0.5,&0.5,&0.0,&0.0,&0.0,&0.0\\\
0.0,&0.0,&0.0,&0.0,&0.0,&1.0,&0.0,&0.0\\\
0.0,&0.0,&0.0,&0.0,&0.0,&0.0,&0.5,&0.5\end{bmatrix}$

This matrix is in fact the conjugate transpose of the resampling from the mean
space to the input space, that is
$\mathbf{\Lambda}_{t}=\mathbf{\Psi}_{mc}^{\text{T}}$. The expansive and
compressive splatting operations are visualised in Fig. 1.

The splatting operation, and its conjugate transpose, are what we add as
layers to the head and tail of a segmentation network, respectively. When
splatting and its transpose are conceptualised as layers, we now have the
first and last layer of the network as adjoint operations, similarly to how
convolutions in the encoder and transposed convolutions in the decoder of a
UNet are adjoint operations. An additional parallel between splatting
(followed by a convolution) and transposed convolutions can be drawn, as they
both let the network learn how to optimally invent missing data.

Having the splat operation at the head of the network allows any sized input
images to be provided to a network, with multiple channels that may have
different sizes or orientations, because the image’s voxels will be pushed
onto the common training grid. Having the conjugate transpose of the splat
operator at the tail of the network allows the loss function (_e.g._ , Dice)
to be computed in the image’s native space, where the target segmentations
were originally drawn. Note that we apply the conjugate transpose to the
‘logit’ output of the network, and afterwards apply the softmax operation.
Furthermore, the splatting operations assume that voxels outside the field-of-
view are zero, and the conjugate transposes assume that data outside of the
field of view are part of the background class. In the next section we extend
the operations to an arbitrary number of dimensions.

###  2.2 Splatting and Resampling

Let us write as $\mathcal{N}:\mathbb{R}^{N\times
C_{i}}\rightarrow\mathbb{R}^{N\times C_{o}}$ a CNN that maps $N$ voxels and
$C_{i}$ input channels to $N$ voxels and $C_{o}$ output channels. Implicitly,
each voxel of the grid is associated with a spatial index $\mathbf{x}_{n}$,
which is identical in the input and output images, and in all channels within
them. Furthermore, voxel coordinates can be linked to a common coordinate
system through the mapping $\mathbf{A}:\mathbf{x}_{n}\mapsto\mathbf{y}$, which
can be encoded in an affine orientation matrix
$\mathbf{A}\in\mathbb{R}^{D+1,D+1}$ ($D$ is the dimensionality of the input).
These matrices are saved in the headers of most medical imaging file formats.
In practice, multiple MR contrasts
$\left\\{\mathbf{f}_{c}\in\mathbb{R}^{N_{c}}\right\\}_{c=1}^{C_{i}}$ may be
defined on different grids, with different coordinate systems
$\left\\{\mathbf{A}_{c}\in\mathbb{R}^{D+1,D+1}\right\\}_{c=1}^{C_{i}}$.
Because CNNs require all images to ‘live’ on the same grid, they are commonly
resampled to the same space. For segmentation, this is typically the space in
which the manual labels were drawn.

In general, resampling can be written as
$\hat{f}_{m}=\sum_{n=1}^{N}f_{n}w(\bm{\phi}(\hat{\mathbf{x}}_{m}),\mathbf{x}_{n})$,
where $\bm{\phi}$ is a mapping from the new grid ($\hat{\mathbf{x}}$) to the
old ($\mathbf{x}$) and $w$ is a weight that depends on the distance between a
voxel $\mathbf{x}_{m}$ and the sampled location $\bm{\phi}(\mathbf{x}_{n})$.
This operation can be conceptualised as a large linear operation
$\hat{\mathbf{f}}=\mathbf{\Phi}\mathbf{f}$, although in practice, the support
of $w$ is small and $\mathbf{\Phi}$ is sparse. In this paper, we use trilinear
interpolation weights.

Let us write the loss function as $\mathcal{L}$ and the labels as
$\mathbf{f}_{l}$, the forward pass of the CNN should therefore really be
written as:

$\left\\{\mathbf{f}_{c}\in\mathbb{R}^{N_{c}},\mathbf{A}_{c}\in\mathbb{R}^{D+1,D+1}\right\\}_{c=1}^{C_{i}}\mapsto\mathcal{L}\left(\mathcal{N}\left(\left[\mathbf{\Phi}_{c}\mathbf{f}_{c}\right]_{c=1}^{C_{i}}\right),\mathbf{\Phi}_{l}\mathbf{f}_{l}\right)~{},$
(1)

although in general the common space is chosen to be that of the labels so
that $\mathbf{\Phi}_{l}$ is the identity. When labels have a lower resolution
than some of the input images, a different formulation could be to re-slice
all images to the higher resolution space (_e.g._ , 1 mm isotropic), and
resample the output of the network to the label space:

$\left\\{\mathbf{f}_{c}\in\mathbb{R}^{N_{c}},\mathbf{A}_{c}\in\mathbb{R}^{D+1,D+1}\right\\}_{c=1}^{C_{i}}\mapsto\mathcal{L}\left(\mathbf{\Psi}_{l}\mathcal{N}\left(\left[\mathbf{\Phi}_{c}\mathbf{f}_{c}\right]_{c=1}^{C_{i}}\right),\mathbf{f}_{l}\right)~{},$
(2)

where $\mathbf{\Psi}_{l}$ maps from the common space to the native label
space, whereas $\mathbf{\Phi}_{l}$ was used to map from the native label space
to the common space (the underlying transformations $\bm{\psi}_{l}$ and
$\bm{\phi}_{l}$ are inverse of each other). However, this does not solve the
issues related to the resampling of the input images raised earlier.

In this paper, we propose to replace the initial resampling with the adjoint
operation of its inverse, as part of the forward pass. Since resampling is a
linear operation, its adjoint is simply its transpose
$\mathbf{\Psi}^{\text{T}}\mathbf{f}$. In practice, it means that native data
are _splatted_ onto the mean space:
$\hat{f}_{m}=\sum_{n=1}^{N}f_{n}w(\bm{\psi}(\mathbf{x}_{n}),\hat{\mathbf{x}}_{m})$.
Importantly, it means that if the resolution of the common space is higher
than that of the native space, the splatted image has many zeros (the data are
_not_ interpolated). The output of the network is then resampled to the native
label space, where the loss is computed:

$\left\\{\mathbf{f}_{c}\in\mathbb{R}^{N_{c}},\mathbf{A}_{c}\in\mathbb{R}^{D+1,D+1}\right\\}_{c=1}^{C_{i}}\mapsto\mathcal{L}\left(\mathbf{\Psi}_{l}\mathcal{N}\left(\left[\mathbf{\Psi}_{c}^{\text{T}}\mathbf{f}_{c},\mathbf{\Psi}_{c}^{\text{T}}\mathbf{1}\right]_{c=1}^{C_{i}}\right),\mathbf{f}_{l}\right)~{},$
(3)

where we have let the network know which zeros are missing and which are
native values, by concatenating splatted images of ones
($\mathbf{\Psi}_{c}^{\text{T}}\mathbf{1}$) to the input. We note that
$\mathbf{\Psi}_{c}^{\text{T}}\mathbf{f}_{c}$ can be seen as the gradient of
the resampling operation with respect to its input, while
$\mathbf{\Psi}_{c}^{\text{T}}\mathbf{1}$ can be seen as a diagonal
approximation of its Hessian.

###  2.3 The Mean Space

What is the best way of defining the common space, in which the training and
inference takes place? Using one of the input images to define this space, for
the complete dataset, is not optimal [5]. A more principled solution is to
compute a mean space from all input orientation matrices [15, 16]. Briefly,
this involves (1) extracting all linear components from the input orientation
matrices; (2) computing their barycentric mean in an iterative fashion by
alternately projecting them to the tangent space of GL(3) about the current
barycentre and updating the barycentre by zero-centering the tangent data; (3)
finding the closest matrix, in the least square sense, that can be encoded by
the product of a rotation matrix (the orientation) and an anisotropic scaling
matrix (the voxel size). In this work, we compute the mean space once, from
the entire training set, although one mean space per mini-batch could
alternatively be used. We constrain the mean-space dimensions to be a power of
two or three, to facilitate fitting encoding/decoding architectures. Finally,
we use a voxel size of 1 mm isotropic. This could be customised however,
_e.g._ , by using larger voxels for a more lightweight model.

##  3 Experiments and Results

This section investigates whether the splat layer can improve multi-modal MRI
brain segmentation in the scenario where, for each subject, we have multiple
MR contrasts of differing resolution and the target labels are defined on one
of the contrasts. We use a simple baseline network that we fit in two ways:
(1) to images that, for each subject, have been resampled to the grid of the
target labels; and (2), to native space images, by extending the baseline
network with our proposed splat layers. The number of learnable parameters in
both networks are the same.

###  3.1 The Baseline Network

We use a fairly light-weight UNet architecture [17] with
$(16,~{}32,~{}64,~{}128)$ channels in the encoder layer and
$(128,~{}64,~{}32,~{}32)$ in the decoder layer, where kernel size $3\times
3\times 3$ and stride two is used throughout. This is followed by a layer of
$3\times 3\times 3$ stacked convolutions with $(32,~{}16,~{}16)$ channels
each, and a stride of one. The last layer then outputs the $K$ segmentations
labels, which are passed through a softmax. All layers use ReLU activations
and batch normalisation. The final network has about 1 million parameters.
This is the baseline network, denoted _UNet_. The UNet is then extended with
our proposed splat layers as described in the Methods section. We denote this
network _MeanSpaceNet_. The mean-space has dimensions $(192,~{}192,~{}192)$
with 1 mm isotropic voxels. Note that the mean-space is defined on only the
training data. Both networks are optimised using the Dice loss and the ADAM
optimiser (lr=$10^{-3}$). During training, we augment with multiplicative
smooth intensity non-uniformities and random diffeomorphic deformations. For
the mean-space model, any spatial augmentation needs to be defined in the
mean-space and then composed to each image’s native space using the affine
matrices. We train for a fixed number of 100 epochs, with a batch size of one.

###  3.2 Simulated Data: Brain Tumour Segmentation

####  3.2.1 TCGA-GBM Dataset.

In this experiment, we use the pre-operative, multi-institutional scans of The
Cancer Genome Atlas (TCGA) Glioblastoma Multiforme (GBM) collection [18, 19],
publicly available in The Cancer Imaging Archive [20]. The dataset was
acquired from different MRI scanners. Each subject has skull-stripped and co-
registered multi-modal (T1w, T1-Gd, T2w, T2-FLAIR) MRIs and segmentation
labels of the enhancing part of the tumor core (ET), the non-enhancing part of
the tumor core (NET), and the peritumoral edema (ED). All MRIs have been
resampled to 1 mm isotropic voxel size and the same dimensions. In this
experiment, we use only the subjects with manually-corrected segmentation
labels, which gives in total $N_{\text{gbm}}=97$ subjects, each with four MR
modalities and three tumour labels.

####  3.2.2 Experiment.

We simulate two datasets from TCGA-GBM. The first dataset, denoted
$\mathcal{D}_{\text{nat}}^{\text{gbm}}$, is created by downsampling the
T2-FLAIR image and the segmentation by a factor of three, in a randomly
selected dimension. This emulates the situation where manual labels have been
drawn on one modality (here the T2-FLAIR), and the other modalities have
different voxel size (here the T1w, T1-Gd and T2w). The second dataset,
denoted $\mathcal{D}_{\text{res}}^{\text{gbm}}$, is created by trilinearly re-
slicing the T1w, T1-Gd and T2w images to the space of the downsampled
T2-FLAIR, so that all images have the same dimensions. This in turn emulates
the situation where all modalities are resampled to the space of the modality
on which the manual labels were drawn. We split the two datasets into equal
(train, validation, test) sets as $(40,~{}17,~{}40)$. We then fit the UNet to
$\mathcal{D}_{\text{res}}^{\text{gbm}}$ and the MeanSpaceNet to
$\mathcal{D}_{\text{nat}}^{\text{gbm}}$. Note that it would not be possible to
train the UNet model on the $\mathcal{D}_{\text{nat}}^{\text{gbm}}$ dataset,
as the subjects’ input modalities have different dimensions. After training we
apply the two trained networks to their test sets and compute pairwise Dice
scores between predicted and target segmentations.

Figure 2: Pairwise Dice scores computed on the TCGA-GBM test set
($N_{\text{gbm}}^{\text{test}}=40$), for three tumour labels (ED, ET, NET),
and two CNN models (MeanSpaceNet, UNet). On each box, the central mark
indicates the median, and the bottom and top edges of the box indicate the
25th and 75th percentiles, respectively. The whiskers extend to the most
extreme data points not considered outliers. Asterisks indicate statistical
significance of paired Wilcox tests after Holm–Bonferroni correction: 0.05
($\ast$), 0.01 ($\ast\ast$), 0.001 ($\ast\ast\ast$) & 0.0001
($\ast\ast\ast\ast$).

####  3.2.3 Results.

The experimental results are shown in the boxplot in Figure 2. The
MeanSpaceNet model achieves the best median Dice score over all classes 0.677
vs 0.528; as well as for all individual classes: 0.705 vs 0.639 for ED, 0.731
vs 0.626 for ET and 0.513 vs 0.340 for NET. Paired Wilcoxon tests with Holm-
Bonferroni correction shows that the segmentation results are all significant
($p<0.0001$).

###  3.3 Real Data: Brain Tissue Segmentation

####  3.3.1 MRBrains18 Dataset.

In this experiment, we use the original scans (before any pre-processing) from
the MICCAI MRBrainS18 grand segmentation challenge
(https://mrbrains18.isi.uu.nl). The dataset was acquired from the same MRI
scanner. Each subject has multi-modal (T1w, T1-IR, T2-FLAIR) MRIs and
segmentation labels of ten brain structures: Cortical gray matter (CGM),
subcortical gray matter (SGM), white matter (WM), white matter lesions (WML),
cerebrospinal fluid in the extracerebral space (ECSF), ventricles (VEN),
cerebellum (CBM), brain stem (BS), infarction (INF) and other (OTH). In this
experiment, we do not use the INF and OTH labels, and we combine WM and WML
into a single label. The images’ voxel sizes (mm) are: T1w
$(1.0,~{}1.0,~{}1.0)$, T1-IR $(0.96,~{}0.96,~{}3.0)$ and T2-FLAIR
$(0.96,~{}0.96,~{}3.0)$. The segmentations were drawn in resampled
$(0.96,~{}0.96,~{}3.0)$ space. In total, we have $N_{\text{brain}}=7$
subjects, each with three MR modalities and seven brain labels.

Figure 3: Pairwise Dice scores computed on the MRBrainS18 datasets from leave-
one-out cross-validation ($N=7$), for seven brain labels (BS, CBM, CGM, ECSF,
SGM, VEN, WM), and two CNN models (MeanSpaceNet, UNet).

####  3.3.2 Experiment.

As in the TCGA-GBM experiment, we construct two datasets. The first dataset,
denoted $\mathcal{D}_{\text{nat}}^{\text{brain}}$, is simply the original
scans from the MRBrainS18 dataset. The second dataset, denoted
$\mathcal{D}_{\text{res}}^{\text{brain}}$, is created by, for each subject in
the MRBrainS18 dataset, trilinearly re-slice the T1w image to the space of the
segmentation, so that all images have the same dimensions. As the MRBrainS18
dataset has a small number of subjects we perform leave-one-out cross-
validation, fitting the UNet on $\mathcal{D}_{\text{res}}^{\text{brain}}$ and
the MeanSpaceNet on $\mathcal{D}_{\text{nat}}^{\text{brain}}$, seven times in
total, and for each fold computing pairwise Dice scores on the subject held
out from training.

####  3.3.3 Results.

The experimental results are shown in Figure 3. The MeanSpaceNet model
achieves the best median Dice score over all classes 0.864 vs 0.813; as well
as for all individual classes: 0.878 vs 0.827 for BS, 0.908 vs 0.836 for CBM,
0.816 vs 0.798 for CGM, 0.783 vs 0.770 for ECSF, 0.824 vs 0.749 for SGM, 0.943
vs 0.910 for VEN and 0.871 vs 0.836 for WM. Paired Wilcoxon tests with Holm-
Bonferroni correction shows that the segmentation results are significant
($p<0.05$) for four out of seven segmentation classes.

##  4 Conclusions

In this paper, we described a splat layer that allows a segmentation network
to automatically handle resolution mismatches in input data. The idea is to
splat each input channel into a mean-space, without interpolating. The forward
pass is then computed in this mean-space and its output prediction pulled back
into the original resolution of the target labels. The splat layer therefore
removes the need for explicit resolution adjustment. We showed on two multi-
modal MRI segmentation tasks that splatting was preferred over resampling.
Besides allowing segmentation networks to work on the raw image voxels, and
computing the loss function in the space of the original target labels, the
splat model could also streamline model deployment as a user does not need to
ensure that input images have a specific size. The dimension of the mean-space
can additionally be defined to allow optimal application of operations, such
as strided convolutions, and/or be made small for faster inference.

Splatting and resampling have the same complexity, which is linear in the
number of native voxels. In practice, the loop over native voxels is
parallelised (whether running on a CPU or GPU), which makes splatting slightly
slower than resampling because multiple native voxels may be pushed to the
same output voxel, necessitating the use of atomic assignment (this also makes
splatting non-deterministic, as the order in which values are summed-up in an
output voxel is architecture-dependent). The cost is therefore somewhat
equivalent, compared to resampling all images to the same grid. However, we
introduce an additional input channel with the count image, which increases
the number of convolution filters in the first layer. However, if the mean
space is known a priori, input images can be splatted offline, as a
preprocessing step. In this case, only resampling of the prediction to the
loss space needs to happen online.

## References

  * [1] F. Isensee, P. F. Jaeger, S. A. Kohl, J. Petersen, and K. H. Maier-Hein, “nnU-Net: a self-configuring method for deep learning-based biomedical image segmentation,” Nature methods, vol. 18, no. 2, pp. 203–211, 2021\.
  * [2] B. H. Menze, A. Jakab, S. Bauer, J. Kalpathy-Cramer, K. Farahani, J. Kirby, Y. Burren, N. Porz, J. Slotboom, R. Wiest, et al., “The multimodal brain tumor image segmentation benchmark (BRATS),” IEEE transactions on medical imaging, vol. 34, no. 10, pp. 1993–2024, 2014.
  * [3] M. Antonelli, A. Reinke, S. Bakas, K. Farahani, B. A. Landman, G. Litjens, B. Menze, O. Ronneberger, R. M. Summers, B. van Ginneken, et al., “The medical segmentation decathlon,” arXiv preprint arXiv:2106.05735, 2021\.
  * [4] H. J. Kuijf, J. M. Biesbroek, J. De Bresser, R. Heinen, S. Andermatt, M. Bento, M. Berseth, M. Belyaev, M. J. Cardoso, A. Casamitjana, et al., “Standardized assessment of automatic segmentation of white matter hyperintensities and results of the WMH segmentation challenge,” IEEE transactions on medical imaging, vol. 38, no. 11, pp. 2556–2568, 2019.
  * [5] P. A. Yushkevich, B. B. Avants, S. R. Das, J. Pluta, M. Altinay, C. Craige, A. D. N. Initiative, et al., “Bias in estimation of hippocampal atrophy using deformation-based morphometry arises from asymmetric global normalization: an illustration in ADNI 3 T MRI data,” Neuroimage, vol. 50, no. 2, pp. 434–445, 2010.
  * [6] W. K. Thompson, D. Holland, A. D. N. Initiative, et al., “Bias in tensor based morphometry Stat-ROI measures may result in unrealistic power estimates,” NeuroImage, vol. 57, no. 1, pp. 1–4, 2011.
  * [7] H.-I. Suk, S.-W. Lee, D. Shen, A. D. N. Initiative, et al., “Hierarchical feature representation and multimodal fusion with deep learning for AD/MCI diagnosis,” NeuroImage, vol. 101, pp. 569–582, 2014\.
  * [8] D. Nie, L. Wang, Y. Gao, and D. Shen, “Fully convolutional networks for multi-modality isointense infant brain image segmentation,” in 2016 IEEE 13Th international symposium on biomedical imaging (ISBI), pp. 1342–1345, IEEE, 2016.
  * [9] K. Kamnitsas, W. Bai, E. Ferrante, S. McDonagh, M. Sinclair, N. Pawlowski, M. Rajchl, M. Lee, B. Kainz, D. Rueckert, et al., “Ensembles of multiple models and architectures for robust brain tumour segmentation,” in International MICCAI brainlesion workshop, pp. 450–462, Springer, 2017\.
  * [10] T. Zhou, H. Fu, G. Chen, J. Shen, and L. Shao, “Hi-net: hybrid-fusion network for multi-modal MR image synthesis,” IEEE transactions on medical imaging, vol. 39, no. 9, pp. 2772–2781, 2020.
  * [11] J. Dolz, K. Gopinath, J. Yuan, H. Lombaert, C. Desrosiers, and I. B. Ayed, “HyperDense-Net: a hyper-densely connected CNN for multi-modal image segmentation,” IEEE transactions on medical imaging, vol. 38, no. 5, pp. 1116–1126, 2018.
  * [12] B. Billot, D. N. Greve, O. Puonti, A. Thielscher, K. Van Leemput, B. Fischl, A. V. Dalca, and J. E. Iglesias, “SynthSeg: Domain randomisation for segmentation of brain MRI scans of any contrast and resolution,” arXiv preprint arXiv:2107.09559, 2021.
  * [13] A. B. Szczotka, D. I. Shakir, D. Ravì, M. J. Clarkson, S. P. Pereira, and T. Vercauteren, “Learning from irregularly sampled data for endomicroscopy super-resolution: a comparative study of sparse and dense approaches,” International journal of computer assisted radiology and surgery, vol. 15, pp. 1167–1175, 2020.
  * [14] L. Westover, “Interactive volume rendering,” in Proceedings of the 1989 Chapel Hill workshop on Volume visualization, pp. 9–16, 1989.
  * [15] J. Ashburner and G. R. Ridgway, “Symmetric diffeomorphic modeling of longitudinal structural MRI,” Frontiers in neuroscience, vol. 6, p. 197, 2013.
  * [16] X. Pennec and V. Arsigny, “Exponential barycenters of the canonical Cartan connection and invariant means on Lie groups,” in Matrix Information Geometry, pp. 123–166, Springer, 2013.
  * [17] O. Ronneberger, P. Fischer, and T. Brox, “U-Net: Convolutional networks for biomedical image segmentation,” in International Conference on Medical image computing and computer-assisted intervention, pp. 234–241, Springer, 2015\.
  * [18] S. Bakas, H. Akbari, A. Sotiras, M. Bilello, M. Rozycki, J. S. Kirby, J. B. Freymann, K. Farahani, and C. Davatzikos, “Advancing the cancer genome atlas glioma MRI collections with expert segmentation labels and radiomic features,” Scientific data, vol. 4, no. 1, pp. 1–13, 2017.
  * [19] S. Bakas, H. Akbari, A. Sotiras, M. Bilello, M. Rozycki, J. S. Kirby, J. B. Freymann, K. Farahani, and C. Davatzikos, “Segmentation labels for the pre-operative scans of the TCGA-GBM collection,” 2017.  Data retrieved from the Cancer Imaging Archive, https://doi.org/10.7937/K9/TCIA.2017.KLXWJJ1Q.
  * [20] K. Clark, B. Vendt, K. Smith, J. Freymann, J. Kirby, P. Koppel, S. Moore, S. Phillips, D. Maffitt, M. Pringle, et al., “The Cancer Imaging Archive (TCIA): maintaining and operating a public information repository,” Journal of digital imaging, vol. 26, no. 6, pp. 1045–1057, 2013.