\section{Introduction}
\label{sec:intro}
\input{sec-intro}

\subsection{Related Work}
\label{sec:relatedwork}
\input{sec-relatedwork}
\section{Methodology}
\label{sec:method}
\input{sec-method}
\section{Experimental Results}
\label{sec:experiment}
\input{sec-experiment}
\section{Conclusion}
This paper presents a novel topological biomarker, \emph{TopoTxR}, that leverages the rich geometric information embedded in structural MRI and enables improvement in downstream CNN processing. In particular, we compute 1D cycles and 2D bubbles from breast DCE-MRIs with the theory of persistent homology; these structures are then used to direct the attention of neural networks. We further demonstrate that \emph{TopoTxR} on treatment-naive imaging is predictive of pCR. 




\bibliographystyle{abbrv}


\subsection{Discussion and TopoTxR Feature Interpretation}
The topological structures extracted by TopoTxR capture the breast tissue structures.
Learning directly on these structures and their vicinity provides the opportunity for interpreting the learning outcomes and drawing novel biological insights. 
Here we provide some visual analysis as a proof of concept.

Fig.~\ref{fig:qualitative} shows the TopoTxR topographical structures from different views for a representative DCE-MRI scan from each group. We observe that the structures (1D and 2D) are sparse for the case exhibiting pCR, and are relatively dense for the non-pCR case. In the corresponding MRI images, we note that the pCR breast has scattered fibroglandular breast density with minimal background parenchymal enhancement. The non-pCR breast has a more heterogenous fibroglandular breast density with moderate background parenchymal enhancement. This possibly suggests that the \textit{TopoTxR} features capture the complex fibrogladular structure which can be a potential indicator of treatment response. 

We also compare the topological behavior of the two populations. Recall the birth time of a topological structure is the threshold at which a cycle appears. In our experiments, since we use the inverse image $f=-I$, the birth time essentially captures -1 times the brightness of a structure. In Fig.~\ref{fig:qualitative} (right), we plot the cumulative density function (CDF) of the birth time of topological structures for pCR (red) and non-pCR (blue) patients. The CDFs suggest that pCR patients' tissue structures are generally less bright (or less visible) compared with that of non-pCR patients. This is consistent with our observation on qualitative examples. 
A Kolmogorov-Smirnov test \cite{massey1951kolmogorov} is performed to compare these CDFs. The computed \textit{p-value} is 0.0002, indicating a significant difference between the distributions of birth times of the pCR and non-pCR patient groups.




\subsection{Background: Persistent Homology}
We review basics of persistent homology in this section. Interested readers may refer to \cite{edelsbrunner2010computational} for more details.
Persistent homology extracts the topological information of data observed via a scalar function. Given an image domain, $X$, and a real-valued function $f:X\rightarrow \mathbb{R}$, we can construct a sublevel set $X_{t} =
\{x \in X: f(x) \leq t\}$ where $t$ is a threshold controlling the ``progress"
of sublevel sets. The family of sublevel sets $\mathcal{X} = \{X_{t}\}_{t \in \mathbb{R}}$ defines a filtration, i.e., a family of subsets of $X$ nested with respect to the inclusion: $X_{\alpha} \subseteq X_{\beta}$ if $\alpha \leq \beta$. As the threshold $t$ increases from $-\infty$ to $+\infty$, topological structures such as connected components, handles, and voids appear and disappear. The birth time of a topological structure is the threshold $t$ at which the structure appears in the filtration. Similarly, the death time is the threshold at which the structure disappears. Persistent homology tracks the topological changes of sublevel sets $X_{t}$ and encodes them in a \emph{persistence diagram}, i.e., a point set in which each point $(b,d)$ represents a topological structure born at $b$ and killed at $d$.

\begin{figure}[btph]
\begin{center}
\includegraphics[width=.96\textwidth]{figs/ipmi21-phf.pdf}
\caption[ph]{From left to right: a synthetic image $f$, sublevel sets at thresholds $b_{1} < b_{2} < d_{2} < d_{1}$ (in black), and the persistence diagram of dimension 1. The red loop represents a 1D structure born at $b_1$ and killed at $d_1$. The green loop represents a 1D structure born at $b_{2}$ and killed at $d_2$. They correspond to the red and green dots in the diagram. }
\label{fig:ph}
\end{center}
\end{figure}

See Fig.~\ref{fig:ph} for an example function $f$ and its sublevel sets at different thresholds. At time $b_1$, a new handle (delineated by the red cycle $c_{1}$) is created. This handle is later destroyed at time $d_{1}$. Another handle delineated by the green cycle $c_{2}$ is created and killed at $b_{2}$ and $d_{2}$ respectively. The topological changes are summarized in a persistence diagram on the right. Each handle corresponds to a 2D dot in $\mathbb{R}^{2}$, whose $x$ and $y$ coordinates are birth and death times. The difference between a dot's birth and death times is called its \emph{persistence}. 


\subsection{Persistence Cycles and Their Computation}
In this section, we introduce cycles that represent topological structures discovered by persistent  homology. We also present an algorithm to compute these cycles.
Intuitively, a topological cycle of dimension $p$ is a $p$-manifold without boundary. A 1-dimensional (1D) cycle is a loop (or a union of a set of loops). A 2-dimensional (2D) cycle is a bubble (or a union of a set of bubbles). 
A cycle $z$ represents a persistent homology structure if it delineates the structure at its birth. For example, in Fig.~\ref{fig:ph}, the red and the green loops represent the handles born at time $b_1$ and $b_2$, respectively.

We assume a discretization of the image domain into distinct elements, i.e., vertices (corresponding to voxels), edges connecting adjacent vertices, squares, and cubes.
These elements are 0-, 1-, 2-, and 3-dimensional cells.
Any set of $p$-cells is called a \emph{$p$-chain}.
We define the \emph{boundary operator} of a $p$-cell, $\sigma$, as the set of its $(p-1)$ faces.
The boundary of an edge is its two vertices. The boundary of a square consists of the 4 edges enclosing it. The boundary of a cube consists of the 6 squares enclosing it. 
For any $p$-chain, $c$, its boundary is  
$\partial(c) = \sum_{\sigma\in c} \partial(\sigma)$, under mod-2 sum. 
For a set of edges forming a path, its boundary are the two end vertices. For any set of squares forming a patch, its boundary is the loop enclosing the patch. The boundary of a set of cubes is the bubble enclosing them.

\begin{figure}[btph]
\begin{center}
\includegraphics[width=.88\textwidth]{figs/ipmi21-cycle.pdf}
\caption[cycle]{(a) Example of a cubical complex whose cells are sorted according to the function values. (b) 2D Boundary matrix $\partial$. (c) Reduced boundary matrix. (d) Persistence diagram and resulting cycles corresponding to $\partial$. (e) 1D boundary matrix.}
\label{fig:cycle}
\end{center}
\end{figure}

A $p$-chain is a \emph{$p$-cycle} if its boundary is empty.
All $p$-cycles form the null space of the boundary operator, i.e. $\{c : \partial(c)=\emptyset\}$. Any topological structure, formally defined as a \emph{homology class}, can be delineated by different cycles, which are equivalent to each other in terms of topology. We can choose any of them to represent this class. In persistent homology, for each dot in the persistence diagram, we can represent it with one representative cycle at its birth. In Fig.~\ref{fig:ph}, the red and green cycles represent the two corresponding handles.
Note that the choice of representative cycle is not unique. A relevant question is to choose the shortest representative cycle (i.e., one with the least number of edges) for each dot in the diagram \cite{wu2017optimal,zhang2019heuristic}. In this paper, we focus on choosing a standard representative cycle, leaving the optimal cycle for future work.

\myparagraph{Computation of persistent homology and representative cycles.}
We assume a filtration function on a discretization of the image domain.
An example discretization of a 2D image is given in Fig.~\ref{fig:cycle}(a). We first sort all cells in increasing order according to their function values. The computation of persistence diagrams is then performed by encoding the $p$-dimensional boundary operator in binary matrices named \emph{boundary matrices}, $\partial_p$. $\partial_p$ maps $p$-cells to their boundaries. Fig.~\ref{fig:cycle} shows the 1D and 2D boundary matrices of the given complex and its filtration. The 1D boundary matrix is essentially the incidence matrix of the underlying graph (Fig.~\ref{fig:cycle}(e)). 
High dimensional boundary matrices are defined similarly (e.g., a 2D boundary matrix in Fig.~\ref{fig:cycle}(b)).

Persistence diagram is computed by reducing the boundary matrix similar to Gaussian elimination, but without row or column perturbation. We reduce through column operations performed on $\partial$ from left to right. 
Fig.~\ref{fig:cycle}(c) shows the reduced 2D boundary matrix. 
Once the boundary matrices are reduced. Each non-zero column corresponds to a persistent dot in the diagram. The reduced column itself is the cycle representing the corresponding topological structure. In this paper, we pay attention to both 1D and 2D cycles, corresponding to loops and bubbles. The extracted cycles will be used to explicitly guide 3D CNNs for analysis.
The computation of topological cycles is of the same complexity as the computation of persistent homology. In theory, it takes $O(n^\omega)$ time ($\omega \approx 2.37$ is the exponent in the matrix multiplication time, i.e., time to multiply two $n\times n$ matrices) \cite{milosavljevic2011zigzag}. Here $n$ is the number of voxels in an image. In practice, to compute all cycles of an input image ($256^3$), it takes approximately 5 minutes.
\subsection{Topological-Cycle-Driven 3D CNN}
An overview of our topological-cycle-driven 3D CNN has been provided in the beginning of Section \ref{sec:method}.
Here we describe the technical details.

To compute persistence and topological cycles, we invert the MRI image $f = -I$ so that the tissue structures correspond to low intensity.
After the computation, we select topological cycles representing dots of the diagram with high persistence. 
The general belief is that low-persistence dots tend to be due to noise. Thus we only select high-persistence cycles, which are considered more salient structures, and more likely to represent true tissue structures. The threshold is a hyperparamter tuned in practice. 

Next, we create two binary 3D masks representing 1D and 2D topological cycles, respectively. Both masks are dilated slightly in order to cover both the structures and their vicinity. 
Instead of directly using these binary masks for pCR prediction, we fill the foreground voxels with their original image intensity values. In other words, we mask the input image with the complement of the cycle mask. See\emph{Topological Structure Masking} step in Fig.~\ref{fig:overview} for the masked MRI image. We generate masked images for both 1D and 2D, and provide them to two CNNs.
All masked MRIs are padded to the same size of 256 $\times$ 256 $\times$ 256.

We use separate networks with the same architecture for cycles and bubbles. The CNN consists of 5 3D convolution layers, each followed by a batch normalization layer and a LeakyReLU. The output feature maps from these two 3D CNNs are reshaped and concatenated into a feature vector. This feature vector is sent into a fully connected (FC) network with three FC layers for final pCR prediction. Besides ReLU and batch normalization, a dropout layer is added to the second FC layer. The final output is a vector of size 2. All three networks are trained together in an end-to-end fashion with stochastic gradient descent (SGD) as the optimizer and cross-entropy as loss.
