# Immunophenotypes of Acute Myeloid Leukemia
From Flow Cytometry Data Using Templates

Ariful Azad1 Bartek Rajwa 2 and Alex Pothen 3111to whom correspondence should
be addressed<EMAIL_ADDRESS>Computational Research Division, Lawrence Berkeley
National Laboratory, California, USA
<EMAIL_ADDRESS>Bindley Bioscience Center, Purdue University, West
Lafayette, Indiana, USA
<EMAIL_ADDRESS>Department of Computer Science, Purdue University, West
Lafayette, Indiana, USA

###### Abstract

##  1 Motivation:

We investigate whether a template-based classification pipeline could be used
to identify immunophenotypes in (and thereby classify) a heterogeneous disease
with many subtypes. The disease we consider here is Acute Myeloid Leukemia,
which is heterogeneous at the morphologic, cytogenetic and molecular levels,
with several known subtypes. The prognosis and treatment for AML depends on
the subtype.

##  2 Results:

We apply flowMatch, an algorithmic pipeline for flow cytometry data created in
earlier work, to compute templates succinctly summarizing classes of AML and
healthy samples. We develop a scoring function that accounts for features of
the AML data such as heterogeneity to identify immunophenotypes corresponding
to various AML subtypes, including APL. All of the AML samples in the test set
are classified correctly with high confidence.

##  3 Availability:

flowMatch is available at www.bioconductor.org/
packages/devel/bioc/html/flowMatch.html; programs
specific to immunophenotyping AML are at www.cs.purdue.edu/
homes/aazad/software.html.

##  4 Contact:

apothen<EMAIL_ADDRESS>
##  5 Introduction

Can Acute Myeloid Leukemia (AML) samples be distinguished from healthy ones
using flow cytometry data from blood or bone marrow with a template-based
classification method? This method builds a template for each class to
summarize the samples belonging to the class, and uses them to classify new
samples. This question is interesting because AML is a heterogeneous disease
with several subtypes and hence it is not clear that a template can succinctly
describe all types of AML. Furthermore, we wish to identify immunophenotypes
(cell types in the bone marrow and blood) that are known to be characteristic
of subtypes of AML. Pathologists use these immunophenotypes to visualize AML
and its subtypes, and a computational procedure that can provide this
information would be more helpful in clinical practice than a classification
score that indicates if an individual is healthy or has AML.

In earlier work, we have developed a template-based classification method for
analyzing flow cytometry (FC) data, which consists of measurements of
morphology (from scattering) and the expression of multiple biomarkers (from
fluorescence) at the single-cell level. Each FC sample consists of hundreds of
thousands or more of such single-cell measurements, and a study could consist
of thousands of samples from different individuals at different time points
under different experimental conditions (Aghaeepour et al., 2013; Shapiro,
2005). We have developed an algorithmic pipeline for various steps in
processing this data (Azad et al., 2013, 2010, 2012). We summarize each sample
by means of the cell populations that it contains. (These terms are defined in
Table 5 and illustrated in Fig. 1.) Similar samples belonging to the same
class are described by a template for the class. A template consists of meta-
clusters that characterize the cell populations present in the samples that
constitute the class. We compute templates from the samples, and organize the
templates into a template tree. Given a sample to classify, we compare it with
the nodes in the template tree, and classify it to the template that it is
closest to. A combinatorial measure for the dissimilarity of two samples or
two templates, computed by means of a mixed edge cover in a graph model
(described in the next section), is at the heart of this approach.

We have applied our algorithmic pipeline for template-based classification to
various problems: to distinguish the phosphorylation state of T cells; to
study the biological, temporal, and technical variability of cell types in the
blood of healthy individuals; to characterize changes in the immune cells of
Multiple Sclerosis patients undergoing drug treatments; and to predict the
vaccination status of HIV patients. However, it is not clear if the AML data
set can be successfully analyzed with this scheme, since AML is a
hetereogeneous disease at the morphologic, cytogenetic and molecular levels,
and a few templates may not describe all of its subtypes.

AML is a disease of myeloid stem cells that differentiate to form several
types of cells in the blood and marrow. It is characterized by the profusion
of immature myeloid cells, which are usually prevented from maturing due to
the disease. The myeloid stem cell differentiates in several steps to form
myeloblasts and other cell types in a hierarchical process. This hierarchical
differentiation process could be blocked at different cell types, leading to
the multiple subtypes of AML. Eight different subtypes of AML based on cell
lineage are included in the French-American-British Cooperative Group (FAB)
classification scheme (Bennett et al., 1985). (A different World Health
Organization (WHO) classification scheme has also been published.) Since the
prognosis and treatment varies greatly among the subtypes of AML, accurate
diagnosis is critical.

We extend our earlier work on template-based classification here by developing
a scoring function that accounts for the subtleties of FC data of AML samples.
Only a small number of the myeloid cell populations in AML samples are
specific to AML, and there are a larger number of cell populations that these
samples share with healthy samples. Furthermore, the scoring function needs to
account for the diversity of the myeloid cell populations in the various
subtypes of AML.

Our work has the advantage of identifying immunophenotypes of clinical
interest in AML from the templates. Earlier work on the AML dataset we work
with has classified AML samples using methods such as nearest neighbor
classification, logistic regression, matrix relevance learning vector
quantization, etc., but they have not identified these immunophenotypes; e.g.,
(Biehl et al., 2013; Manninen et al., 2013; Qiu, 2012).

Template-based classification has the advantage of being more robust than
nearest neighbor classification since a template summarizes the characteristic
properties of a class while ignoring small sample-to-sample variations. It is
also scalable to large numbers of samples, since we compare a sample to be
classified only against a small number of templates rather than the much
larger number of samples. The comparisons with the templates can be performed
efficiently using the structure of the template tree. It also reduces the data
size by clustering the data to identify cell populations and then working with
the statistical distributions characterizing the cell populations, in contrast
to some of the earlier approaches that work with data sets even larger than
the FC data by creating multiple variables from a marker (reciprocal, powers,
products and quotients of subsets of the markers, etc.).

Summary of terminology  Terms Meaning Sample From FC data. Characterized by
the collection of cell populations included within it. Cell population
(cluster) A group of cells with identical morphology and expressing similar
biomarkers, e.g., helper T cells, B cells. Computed from a sample. Meta-
cluster A set of similar cell clusters from different samples. Computed from
similar clusters in samples. Template A collection of meta-clusters from
samples of the same class.

Template-based classification has been employed in other areas such as
character, face, and image recognition, but its application to FC is
relatively recent. In addition to our work, templates have been used for
detecting the effects of phosphorylation (Pyne et al., 2009), evaluating the
efficiency of data transformations (Finak et al., 2010), and labeling clusters
across samples (Spidlen et al., 2013).

##  6 Methods

###  6.1 The AML Dataset

We have used an FC dataset on AML that was included in the DREAM6/FlowCAP2
challenge of 2011. The dataset consists of FC measurements of peripheral blood
or bone marrow aspirate collected from $43$ AML positive patients and $316$
healthy donors over a one year period. Each patient sample was subdivided into
eight aliquots (“tubes”) and analyzed with different biomarker combinations,
five markers per tube (most markers are proteins). In addition to the markers,
the forward scatter (FS) and side scatter (SS) of each sample was also
measured in each tube. Hence, we have $359\times 8=2,872$ samples and each
sample is seven-dimensional (five markers and the two scatters). Tube 1 is an
isotype control used to detect non-specific antibody binding and Tube 8 is an
unstained control for identifying background or autofluorescence of the
system. Since the data has been compensated for autofluorescence and spectral
overlap by experts, we omit these tubes from our analysis. The disease status
(AML/healthy) of $23$ AML patients and $156$ healthy donors are provided as
training set, and the challenge is to determine the disease status of the rest
of the samples, $20$ AML and $157$ healthy, based only on the information in
the training set. The complete dataset is available at
http://flowrepository.org/.

The side scatter (SS) and all of the fluorescence channels are transformed
logarithmically, but the forward scatter (FS) is linearly transformed to the
interval [0,1] so that all channels have values in the same range. This
removes any bias towards FS channel in the multi-dimensional clustering phase.
After preprocessing, an FC sample is stored as an $n\times p$ matrix $A$,
where the element $A(i,j)$ quantifies the $j^{th}$ feature in the $i^{th}$
cell, and $p$ is the number of features measured in each of $n$ cells. In this
dataset, $p=7$ for each tube and $n$ varies among the samples.

###  6.2 Identifying cell populations in each sample

We employ a two-stage clustering approach for identifying phenotypically
similar cell populations (homogeneous clusters of cells) in each sample. At
first, we apply the $k$-means clustering algorithm for a wide range of values
for $k$, and select the optimum number of clusters $k^{*}$ by simultaneously
optimizing the Calinski-Harabasz and S_Dbw cluster validation methods (Halkidi
et al., 2001). Next, we model the clusters identified by the $k$-means
algorithm with a finite mixture model of multivariate normal distributions. In
the mixture model, the $i^{th}$ cluster is represented by two distribution
parameters $\boldsymbol{\mu_{i}}$, the $p$-dimensional mean vector, and
$\Sigma$, the $p\times p$ covariance matrix. The distribution parameters for
each cluster are then estimated using the Expectation-Maximization (EM)
algorithm. The statistical parameters of a cluster are used to describe the
corresponding cell population in the rest of the analysis.

###  6.3 Dissimilarity between samples

We calculate the dissimilarity between a pair of cell populations by the
Mahalanobis distance between their distributions. Let
$c_{1}(\boldsymbol{\mu_{1}},\Sigma_{1})$ and
$c_{2}(\boldsymbol{\mu_{2}},\Sigma_{2})$ be two normally distributed clusters
and $\Sigma_{p}$ be the pooled covariance of $\Sigma_{1}$ and $\Sigma_{2}$.
The Mahalanobis distance $d(c_{1},c_{2})$ between the clusters is computed as
follows:

$\begin{array}[]{lcl}d(c_{1},c_{2})=\frac{1}{2}(\boldsymbol{\mu_{1}}-\boldsymbol{\mu_{2}})^{\top}\Sigma_{p}^{-1}(\boldsymbol{\mu_{1}}-\boldsymbol{\mu_{2}}),\
\ \text{where}\\\ \\\
\Sigma_{p}={((n_{1}-1)*\Sigma_{1}+(n_{2}-1)*\Sigma_{2}))}/{(n_{1}+n_{2}-2)}.\end{array}$
(1)

We compute the dissimilarity between a pair of samples by optimally matching
(in a graph-theoretic model) similar cell clusters and summing up the
dissimilarities of the matched clusters. In earlier work, we have developed a
robust variant of a graph matching algorithm called the Mixed Edge Cover (MEC)
algorithm that allows a cluster in one sample to be matched with zero, one, or
more clusters in the second sample (Azad et al., 2010). The cell population in
the first sample could be either absent, or present, or split into two or more
cell populations in the second sample. These can happen due to changes in
biological conditions or due to artifactual errors in clustering.

Consider two FC samples $A$ and $B$ consisting of $k_{a}$ and $k_{b}$ cell
populations such that $A=\\{a_{1},a_{2},...,a_{k_{a}}\\}$, and
$B=\\{b_{1},b_{2},...,b_{k_{b}}\\}$ where $a_{i}$ is the $i^{th}$ cluster from
sample $A$ and $b_{j}$ is the $j^{th}$ cluster from $B$. The mixed edge cover
computes a mapping mec, of clusters across $A$ and $B$ such that ${\tt
mec}(a_{i})\in\mathcal{P}(B)$ and ${\tt mec}(b_{j})\in\mathcal{P}(A)$, where
$\mathcal{P}(A)$ ($\mathcal{P}(B)$) is the power set of $A$ ($B$). When a
cluster $a_{i}$ (or $b_{j}$) remains unmatched under mec , i.e., ${\tt
mec}(a_{i})=\emptyset$, we set $d(a_{i},-)=\lambda$, where the fixed cost
$\lambda$ is a penalty for leaving a vertex unmatched. We set $\lambda$ to
$\sqrt{p}$ so that a pair of clusters get matched only if the average squared
deviation across all dimensions is less than one. The cost of a mixed edge
cover mec is the sum of the dissimilarities of all pairs of matched clusters
and the penalties due to the unmatched clusters. A minimum cost mixed edge
cover is a mixed edge cover with the minimum cost. We use this minimum cost as
the dissimilarity $D(A,B)$ between a pair of samples $A$ and $B$:

$\min_{\begin{subarray}{c}\text{mixed}\ \text{edge}\\\ \text{covers},\ {\tt
mec}\end{subarray}}\ (\sum_{\begin{subarray}{c}1\leq i\leq k_{a}\\\
b_{j}\in{\tt mec}(a_{i})\end{subarray}}d(a_{i},b_{j})\
+\sum_{\begin{subarray}{c}1\leq i\leq k_{b}\\\ a_{j}\in{\tt
mec}(b_{i})\end{subarray}}d(b_{i},a_{j})),$ (2)

where $d(a_{i},b_{j})$ is computed from Equation (1). A minimum cost mixed
edge cover can be computed by a modified minimum weight perfect matching
algorithm in $O(k^{3}\log k)$ time where $k$ is the maximum number of clusters
in a sample (Azad et al., 2010). The number of cell clusters $k$ is typically
small (fewer than fifty for the AML data), and the dissimilarity between a
pair of samples can be computed in less than a second on a desktop computer.

###  6.4 Creating templates from a collection of samples

We have designed a hierarchical matching-and-merging (HM&M) algorithm that
arranges a set of similar samples into a binary _template tree_ data structure
(Azad et al., 2012). A node in the tree represents either a sample (leaf node)
or a template (internal node). In both cases, a node is characterized by a
finite mixture of multivariate normal distributions each component of which is
a cluster or meta-cluster. Fig. 1 shows an example of a template-tree created
from four hypothetical samples, $S_{1},S_{2},S_{3}$, and $S_{4}$.

Figure 1: (a) A hierarchical template tree created by the HM&M algorithm from
four hypothetical samples $S_{1},S_{2},S_{3}$ and $S_{4}$. Cells are denoted
with dots, and clusters with solid ellipses in the samples are at the leaves
of the tree. An internal node represents a template created from its children,
and the root represents the template of these four samples. A meta-cluster is
a homogeneous collection of clusters and is denoted by a dashed ellipse inside
the template. (b) One phase of the HM&M algorithm creating a sub-template
$T(S_{3},S_{4})$ from samples $S_{3}$ and $S_{4}$. At first, corresponding
clusters across $S_{3}$ and $S_{4}$ are matched by the MEC algorithm, and then
the matched clusters are merged to construct new meta-clusters.

Let a node $v_{i}$ (representing either a sample or a template) in the
template tree consist of $k_{i}$ clusters or meta-clusters $c^{i}_{1}$,
$c^{i}_{2}$, $\ldots$, $c^{i}_{k_{i}}$. A node $v_{i}$ is called an “orphan”
if it does not have a parent in the template-tree. Consider $N$ flow cytometry
samples $S_{1},S_{2},\ldots,S_{N}$ belonging to a class. Then the HM&M
algorithm for creating a template tree from these samples can be described by
the following three steps.

1\. _Initialization_ : Create a node $v_{i}$ for each of the $N$ samples
$S_{i}$. Initialize all these nodes to the set of orphan nodes. Repeat the
matching and merging steps until a single orphan node remains.

2\. _Matching_ : Compute the dissimilarity $D(v_{i},v_{j})$ between every pair
of nodes $v_{i}$ and $v_{j}$ in the current ${\tt Orphan}$ set with the mixed
edge cover algorithm. (using Equation (2))

3\. _Merging:_ Find a pair of orphan nodes $(v_{i},v_{j})$ with minimum
dissimilarity $D(v_{i},v_{j})$ and merge them to create a new node $v_{l}$.
Let ${\tt mec}$ be a function denoting the mapping of clusters from $v_{i}$ to
$v_{j}$. That is, if $c_{x}^{i}\in v_{i}$ is matched to $c_{y}^{j}\in v_{j}$,
then $c_{y}^{j}\in{\tt mec}(c_{x}^{i})$, where $1\leq x\leq k_{i}$ and $1\leq
y\leq k_{j}$. Create a new meta-cluster $c^{l}_{z}$ from each set of matched
clusters, $c^{l}_{z}=\\{c_{x}^{i}\cup{\tt mec}(c^{i}_{x})\\}$. Let $k_{l}$ be
the number of the new meta-clusters created above. Then the new node $v_{l}$
is created as a collection of these newly created meta-clusters, i.e.,
$v_{l}=\\{c^{l}_{1},c^{l}_{2},...,c^{l}_{k_{l}}\\}$. The distribution
parameters, $(\mu^{l}_{z},\Sigma^{l}_{z})$, of each of the newly formed meta-
clusters $c^{l}_{z}$ are estimated by the EM algorithm. The height of $v_{l}$
is set to $D(v_{i},v_{j})$. The node $v_{l}$ becomes the parent of $v_{i}$ and
$v_{j}$, and the set of orphan nodes is updated by including $v_{l}$ and
deleting $v_{i}$ and $v_{j}$ from it. If there are orphan nodes remaining, we
return to the matching step, and otherwise, we terminate.

When the class labels of samples are not known a priori, the roots of well-
separated branches of tree give different class templates. However, if samples
belong to the same class – as is the case for the AML dataset studied in this
paper, the root of the template-tree gives the class-template. The HM&M
algorithm requires $O(N^{2})$ dissimilarity computations and $O(N)$ merge
operations for creating a template from a collection of $N$ samples. Let $k$
be the maximum number of clusters or meta-clusters in any of the nodes of the
template-tree. Then a dissimilarity computation takes $O(k^{3}\log k)$ time
whereas the merge operation takes $O(k)$ time when distribution parameters of
the meta-clusters are computed by maximum likelihood estimation. Hence, the
time complexity of the algorithm is $O(N^{2}k^{3}\log k)$, which is $O(N^{2})$
for bounded $k$. The complexity of the algorithm can be reduced to $O(N\log
N)$ by avoiding the computation of all pairwise dissimilarities between the
samples, for larger numbers of samples $N$, but we did not need to do this
here.

###  6.5 Classification score of a sample in AML dataset

Consider a sample $X$ consisting of $k$ cell populations
$S=\\{c_{1},c_{2},...,c_{k}\\}$, with the $i^{th}$ cluster $c_{i}$ containing
$|c_{i}|$ cells. Let $T^{-}$ and $T^{+}$ be the templates created from AML-
negative (healthy) and AML-positive training samples, respectively. We now
describe how to compute a score $f(X)$ in order to classify the sample $X$ to
either the healthy class or the AML class.

The intuition behind the score is as follows. An AML sample contains two kinds
of cell populations: (1) AML-specific myeloblasts and myeloid cells, and (2)
AML-unrelated cell populations, such as lymphocytes. The former cell
populations correspond to the immunophenotypes of AML-specific metaclusters in
the AML template, and hence when we compute a mixed edge cover between the AML
template and an AML sample, these clusters get matched to each other. (Such
clusters in the sample do not match to any metacluster in the healthy
template.) Hence we assign a positive score to a cluster in sample when it
satisfies this condition, signifying that it is indicative of AML. AML-
unrelated cell populations in a sample could match to meta-clusters in the
healthy template, and also to AML-unrelated meta-clusters in the AML template.
When either of these conditions is satisfied, a cluster gets a negative score,
signifying that it is not indicative of AML. Since AML affects only the
myeloid cell line and its progenitors, it affects only a small number of AML-
specific cell populations in an AML sample. Furthermore, different subtypes of
AML affect different cell types in the myeloid cell line. Hence there are many
more clusters common to healthy samples than there are AML-specific clusters
common to AML samples. (This is illustrated later in Fig. 3 (c) and (d).) Thus
we make the range of positive scores relatively higher than the range of
negative scores. This scoring system is designed to reduce the possibility of
a false negative (an undetected AML-positive patient), since this is more
serious in the diagnosis of AML. Additional data such as chromosomal
translocations and images of bone marrow from microscopy could confirm an
initial diagnosis of AML from flow cytometry.

In the light of the discussion above, we need to identify AML-specific
metaclusters initially. Given the templates $T^{+}$ and $T^{-}$, we create a
complete bipartite graph with the meta-clusters in each template as vertices,
and with each edge weighted by the Mahalanobis distance between its endpoints.
When we compute a minimum cost mixed edge cover in this graph, we will match
meta-clusters common to both templates, and such meta-clusters represent non-
myeloid cell populations that are not AML-specific. On the other hand, meta-
clusters in the AML template $T^{+}$ that are not matched to a meta-cluster in
the healthy template $T^{-}$ correspond to AML-specific metaclusters. We
denote such meta-clusters in the AML template $T^{+}$ by the set $M^{+}$.

Now we can proceed to compare a sample against the template for healthy
samples and the template for AML. We compute a minimum cost mixed edge cover
between a sample $X$ and the healthy template $T^{-}$, and let ${\tt
mec^{-}}(c_{i})$ denote the set of meta-clusters in $T^{-}$ mapped to a
cluster $c_{i}$ in the sample $X$. Similarly, compute a minimum cost mixed
edge cover between $X$ and the AML template $T^{+}$, and let ${\tt
mec^{+}}(c_{i})$ denote the set of meta-clusters in $T^{+}$ mapped to a
cluster $c_{i}$. These sets could be empty if $c_{i}$ is unmatched in the
mixed edge cover. We compute the average Mahalanobis distance between $c_{i}$
and the meta-clusters matched to it in the template $T^{-}$, and define this
as the dissimilarity $d(c_{i},{\tt mec^{-}}(c_{i}))$. From the formulation of
the mixed edge cover in (Azad et al., 2010), we have $d(c_{i},{\tt
mec^{-}}(c_{i}))\leq 2\lambda$. Hence we define the similarity between $c_{i}$
and ${\tt mec^{-}}(c_{i})$ as $s(c_{i},{\tt
mec^{-}}(c_{i}))=2\lambda-d(c_{i},{\tt mec^{-}}(c_{i}))$. By analogous
reasoning, the similarity between $c_{i}$ and ${\tt mec^{+}}(c_{i})$ is
defined as $s(c_{i},{\tt mec^{+}}(c_{i}))=2\lambda-d(c_{i},{\tt
mec^{+}}(c_{i}))$.

The score of a sample is the sum of the scores of its clusters. We define the
score of a cluster $c_{i}$, $f(c_{i})$, as the sum of two functions
$f^{+}(c_{i})$ and $f^{-}(c_{i})$ multiplied with suitable weights. A positive
score indicates that the sample belongs to AML, and a negative score indicates
that it is healthy.

The function $f^{+}(c_{i})$ contributes a positive score to the sum if $c_{i}$
is matched to an AML-specific meta-cluster in the mixed edge cover between the
sample $X$ and the AML template $T^{+}$, and a non-positive score otherwise.
For the latter case, there are two subcases: If $c_{i}$ is unmatched in the
mixed edge cover, it corresponds to none of the meta-clusters in the template
$T^{+}$, and we assign it a zero score. If $c_{i}$ is matched only to non-AML
specific meta-clusters in the AML template $T^{+}$, then we assign it a small
negative score to indicate that it likely belongs to the healthy class (recall
that $k$ is the number of clusters in sample $X$). Hence

$\displaystyle f^{+}(c_{i})$ $\displaystyle=$
$\displaystyle\begin{cases}s\left(c_{i},{\tt mec^{+}}(c_{i})\right),&\text{if\
\ }{\tt mec^{+}}(c_{i})\cap M^{+}\neq\emptyset,\\\
-\frac{1}{k}\left[s(c_{i},{\tt mec^{+}}(c_{i}))\right],&\text{if\ \ }{\tt
mec^{+}}(c_{i})\cap M^{+}=\emptyset,\\\ &{\text{and\ \ }\tt
mec^{+}}(c_{i})\neq\emptyset,\\\ 0,&\text{if\ \ }{\tt
mec^{+}}(c_{i})=\emptyset.\\\ \end{cases}$

The function $f^{-}(c_{i})$ contributes a negative score to a cluster $c_{i}$
in the sample $X$ if it is matched with some meta-cluster in the healthy
template $T^{-}$, indicating that it likely belongs to the healthy class. If
it is not matched to any meta-cluster in $T^{-}$, then we assign it a positive
score $\lambda$. This latter subcase accounts for AML-specific clusters in the
sample, or a cluster that is in neither template. In this last case, we
acknowledge the diversity of cell populations in AML samples. Hence we have

$\displaystyle f^{-}(c_{i})$ $\displaystyle=$
$\displaystyle\begin{cases}-\frac{1}{k}\left[s(c_{i},{\tt
mec^{-}}(c_{i}))\right],&\text{if\ \ }{\tt mec^{-}}(c_{i})\neq\emptyset,\\\
\lambda,&\text{if\ \ }{\tt mec^{-}}(c_{i})=\emptyset.\\\ \end{cases}$

Finally, we define

$f(X)=\sum_{c_{i}\in
X}\frac{|c_{i}|}{|X|}\frac{1}{2}(f^{+}(c_{i})+f^{-}(c_{i})).$ (3)

Here $|X|$ is the number of cells in the sample $X$. The score of a cluster
$c_{i}$ is weighted by the fractional abundance of cells in it.

Figure 2: Cell types identified on the side scatter (SS) and CD45 channels for
a healthy and an AML positive sample. Cell populations are discovered in the
seven-dimensional samples with the clustering algorithm and then projected on
these channels for visualization. A pair of clusters denoting the same cell
type is marked with the same color. The proportion of myeloid blast cells
(shown in red) increases significantly in the AML sample.

Figure 3: The healthy and AML templates created from Tube 6. (a) The template-
tree created from 156 healthy samples in the training set. (b) The template-
tree created from 23 AML samples in the training set. Samples in the red
subtree exhibit the characteristics of Acute Promyelocytic Leukemia (APL) as
shown in Subfigure (f). (c) Fraction of 156 healthy samples present in each of
the 22 meta-clusters in the healthy template. Nine meta-clusters, each of them
shared by at least 60% of the healthy samples, form the core of the healthy
template. (d) Fraction of 23 AML samples present in each of the 40 meta-
clusters in the AML template. The AML samples, unlike the healthy ones, are
heterogeneously distributed over the meta-clusters. (e) The expression levels
of markers in the meta-cluster shown with blue bar in Subfigure (d). (Each
horizontal bar in Subfigures (e) and (f) represents the average expression of
a marker and the error bar shows its standard deviation.) This meta-cluster
represents lymphocytes denoted by medium SS and high CD45 expression and
therefore does not express the AML-related markers measured in Tube 6. (f)
Expression of markers in a meta-cluster shown with red bar in Subfigure (d).
This meta-cluster denotes myeloblast cells as defined by the SS and CD45
levels. This meta-cluster expresses HLA-DR-CD117+CD34-CD38+, a characteristic
immunophenotype of APL. Five AML samples sharing this meta-cluster are similar
to each other as shown in the red subtree in Subfigure (b).

##  7 Results

###  7.1 Cell populations in healthy and AML samples

In each tube, we identify cell populations in the samples using the clustering
algorithm described in Section 6.2. Each sample contains five major cell types
that can be seen when cell clusters are projected on the side scatter (SS) and
CD45 channels, as depicted in Fig. 2. (Blast cells are immature progenitors of
myeloid cells or lymphocytes.) The side scatter measures the granularity of
cells, whereas CD45 is variably expressed by different white blood cells
(leukocytes). AML is initially diagnosed by rapid growth of immature myeloid
blast cells with medium SS and CD45 expressions (Lacombe et al., 1997) marked
in red in Fig. 2. According to the WHO guidelines, AML is initially confirmed
when the sample contains more than 20% blasts. This is the case for all,
except one of the AML samples in the DREAM6/FlowCAP2 training set, and the
latter will be discussed later.

###  7.2 Healthy and AML templates

From each tube of the AML dataset, using the training samples, we build two
templates: one for healthy samples, and one for AML. As described in Section
6.4, the HM&M algorithm organizes samples of the same class into a binary
template tree whose root represents the class template. The template trees
created from the healthy and AML training samples in Tube 6 are shown in
Subfigures 3(a) and 3(b) respectively. The height of an internal node in the
template tree measures the dissimilarity between its left and right children,
whereas the horizontal placement of a sample is arbitrary. In these trees, we
observe twice as much heterogeneity in the AML samples than among the healthy
samples (in the dissimilarity measure), despite the number of healthy samples
being five times as numerous as the AML samples. The larger heterogeneity
among AML samples is observed in other tubes as well. The template-tree for
AML partitions these samples into different subtrees that possibly denote
different subtypes of AML. For example, the subtree in Fig. 3(b) that is
colored red includes samples (with subject ids 37, 58, 67, 89, and 117) with
immunophenotypes of Acute Promyelocytic Leukemia (APL) (discussed later in
this section).

Together, the meta-clusters in a healthy template represent a healthy immune
profile in the feature space of a tube from which the template is created. We
obtained $22$ meta-clusters in the healthy template created from Tube 6. The
percentage of samples from the training set participating in each of these
meta-clusters is shown in Fig. 3(c). Observe that $60\%$ or more of the
healthy samples participate in the nine most common meta-clusters (these
constitute the core of the healthy template). The remaining thirteen meta-
clusters include populations from a small fraction of samples. These
populations could correspond to biological variability in the healthy samples,
variations in the FC experimental protocols, and possibly also from the
splitting of populations that could be an artifact of the clustering
algorithm.

The AML template created from Tube 6 consists of forty meta-clusters (almost
twice the number in the more numerous healthy samples). Fig. 3(d) shows that,
unlike the healthy samples, the AML samples are heterogeneous with respect to
the meta-clusters they participate in: There are $21$ meta-clusters that
include cell populations from at least $20\%$ of the AML samples. Some of the
meta-clusters common to a large number of AML samples represent non-AML
specific cell populations. For example, Fig. 3(e) shows the average marker
expressions of the meta-cluster shown in the blue bar in Fig. 3(d). This meta-
cluster has low to medium side scatter and high CD45 expression, and therefore
represents lymphocytes (Fig. 2). Since lymphocytes are not affected by AML,
this meta-cluster does not express any AML-related markers, and hence can be
described as HLA-DR-CD117-CD34-CD38-, as expected. Fig. 3(f) shows the
expression profile of another meta-cluster shown in the red bar in Fig. 3(d).
This meta-cluster consists of five cell populations from five AML samples
(with subject ids 37, 58, 67, 89, and 117) and exhibits medium side scatter
and CD45 expression and therefore, represents myeloid blast cells.
Furthermore, this meta-cluster is HLA-DR-CD117+CD34-CD38+, and represents a
profile known to be that of Acute Promyelocytic Leukemia (APL) (Paietta,
2003). APL is subtype M3 in the FAB classification of AML (Bennett et al.,
1985)) and is characterized by chromosomal translocation of retinoic acid
receptor-alpha (RAR$\alpha$) gene on chromosome 17 with the promyelocytic
leukemia gene (PML) on chromosome 15, a translocation denoted as t(15;17). In
the feature space of Tube 6, these APL samples are similar to each other while
significantly different from the other AML samples. Our template-based
classification algorithm groups these samples together in the subtree colored
red in the AML template tree shown in Fig. 3(b).

###  7.3 Identifying meta-clusters symptomatic of AML

In each tube, we register meta-clusters across the AML and healthy templates
using the mixed edge cover (MEC) algorithm. Meta-clusters in the AML template
that are not matched to any meta-clusters in the healthy template represent
the abnormal, AML-specific immunophenotypes while the matched meta-clusters
represent healthy or non-AML-relevant cell populations. Table 7.3 lists
several unmatched meta-clusters indicative of AML from different tubes. As
expected, every unmatched meta-cluster displays medium side scatter and CD45
expression characteristic of myeloid blast cells, and therefore we omit FS,
SS, and CD45 values in Table 7.3. We briefly discuss the immunophenotypes
represented by each AML-specific meta-cluster in each tube, omitting the
isotype control Tube 1 and unstained Tube 8.

Some of the meta-clusters characteristic of AML for the 23 AML samples in the
training set. In the second column, ‘$-$’, ‘low’, and ‘$+$’ denote very low,
low and high, abundance of a marker, respectively, and $\pm$ denotes a marker
that is positively expressed by some samples and negatively expressed by
others. The number of samples participating in a meta-cluster is shown in the
third column. The average fraction of cells in a sample participating in a
meta-cluster, and the standard deviation, are shown in the fourth column.
Tube Marker expression #Samples Fraction of cells 2
Kappa${}^{\text{low}}$Lambda${}^{\text{low}}$CD19+CD20- 5 $6.3\%(\pm 6.8)$ 3
CD7+CD4-CD8-CD2- 4 $18.0\%(\pm 4.8)$ 4 CD15-CD13+CD16-CD56- 17 $16.6\%(\pm
6.9)$ 4 CD15-CD13+CD16-CD56+ 8 $11.1\%(\pm 5.7)$ 5 CD14-CD11c-CD64-CD33+ 10
$13.5\%(\pm 5.2)$ 5 CD14-CD11c+CD64-CD33+ 18 $10.8\%(\pm 3.8)$ 5
CD14${}^{\text{low}}$CD11c+CD64${}^{\text{low}}$CD33+ 6 $13.8\%(\pm 4.3)$ 6
HLA-DR+CD117+CD34+CD38+ 11 $13.3\%(\pm 2.6)$ 6 HLA-DR+CD117±CD34+CD38+ 13
$17.3\%(\pm 6.6)$ 6 HLA-DR-CD117±CD34-CD38+ 5 $12.9\%(\pm 4.7)$ 7
CD5-CD19+CD3-CD10- 3 $12.3\%(\pm 2.4)$ 7 CD5+CD19-CD3-CD10- 3 $10.0\%(\pm
8.5)$ 7 CD5-CD19-CD3-CD10+ 1 $9.9\%$

Tube 6 is the most important panel for diagnosing AML since it includes
several markers expressed by AML blasts. HLA-DR is an MHC class II cell
surface receptor complex that is expressed on antigen-presenting cells, e.g.,
B cells, dendritic cells, macrophages, and activated T cells. It is expressed
by myeloblasts in most subtypes of AML except M3 and M7 (Campana and Behm,
2000). CD117 is a tyrosine kinase receptor (c-KIT) expressed in blasts of some
cases ($30-100\%$) of AML (Campana and Behm, 2000). CD34 is a cell adhesion
molecule expressed on different stem cells and on the blast cells of many
cases of AML ($40\%$) (Mason et al., 2006). CD38 is a glycoprotein found on
the surface of blasts of several subtypes of AML but usually not expressed in
the M3 subtypes of AML (Keyhani et al., 2000). In Tube 6, we have identified
two meta-clusters with high expressions of HLA-DR and CD34. One of them also
expresses CD117 and CD34, and Fig. 4(c) shows the bivariate contour plots of
the cell populations contained in this meta-cluster. The second meta-cluster
expresses positive but low levels of CD117 and CD34. These two HLA-DR+CD34+
meta-clusters together are present in $18$ out of the $23$ training AML
samples. The remaining five samples (subject id: 5, 7, 103, 165, 174) express
HLA-DR-CD117±CD34-CD38+ myeloblasts, which is an immunophenotype of APL
(Paietta, 2003) as was discussed earlier. Fig. 4(d) shows the bivariate
contour plots of this APL-specific meta-cluster.

Tube 5 contains several antigens typically expressed by AML blasts, of which
CD33 is the most important. CD33 is a transmembrane receptor protein usually
expressed on immature myeloid cells of the majority of cases of AML ($91\%$
reported in (Legrand et al., 2000)). The AML specific meta-clusters identified
from markers in Tube 5 (see Table 7.3) include CD33+ myeloblasts from every
sample in the training set. Several of the CD33+ populations also express
CD11c, a type I transmembrane protein found on monocytes, macrophages and
neutrophils. CD11c is usually expressed by blast cells in acute myelomonocytic
leukemia (M4 subclass of AML), and acute monocytic leukemia (M5 subclass of
AML) (Campana and Behm, 2000). Therefore CD14-CD11c+CD64-CD33+ meta-cluster
could represent patients with M4 and M5 subclasses of AML. We show the
bivariate contour plots of this meta-cluster in Fig. 4(b) .

Tube 4 includes several markers usually expressed by AML blasts, of which CD13
is the most important. CD13 is a zinc-metalloproteinase enzyme that binds to
the cell membrane and degrades regulatory peptides (Mason et al., 2006). CD13
is expressed on the blast cells of the majority of cases of AML (95% as
reported in (Legrand et al., 2000)). Table 7.3 shows two AML-specific meta-
clusters detected from the blast cells in Tube 4. In addition to CD13, eight
AML samples express CD56 glycoprotein that is naturally expressed on NK cells,
a subset of CD4+ T cells and a subset of CD8+ T cells. Raspadori et al.
(Raspadori et al., 2001) reported that CD56 was more often expressed by
myeloblasts in FAB subclasses M2 and M5, which covers about 42% of AML cases
in a study by Legrand et al. (Legrand et al., 2000). In this dataset, we
observe more AML samples expressing CD13+CD56- blasts than expressing
CD13+CD56+ blasts, which conforms to the findings of Raspadori et al.
(Raspadori et al., 2001). Fig. 4(a) shows the bivariate contour plots of the
CD13+CD56- meta-cluster.

Tube 2 is a B cell panel measuring B cell markers CD19 and CD20, and Kappa
($\kappa$) and Lambda ($\lambda$), immunoglobulin light chains present on the
surface of antibodies produced by B lymphocytes. B-cell specific markers are
occasionally co-expressed with myeloid antigens especially in FAB M2 subtype
of AML (with chromosomal translocation t(8;21)) (Campana and Behm, 2000;
Walter et al., 2010). In Tube 2, we have identified a meta-cluster in the
myeloblasts that expresses high levels of CD19 and low levels of Kappa and
Lambda. The five samples with subject ids 5, 7, 103, 165, and 174
participating in this meta-cluster possibly belong to the FAB-M2 subtype of
AML. Tube 3 is a T cell panel measuring T cell specific markers CD4, CD8, CD2,
and CD7. Tube 7 is a lymphocyte panel with several markers expressed on T and
B lymphocytes and is less important in detecting AML since they are
infrequently expressed by AML blasts.

Figure 4: Bivariate contour plots (side scatter vs. individual marker) for two
meta-clusters (one in each row) indicative of AML. The ellipses in a subplot
denote the 95th quantile contour lines of cell populations included in the
corresponding meta-cluster. Myeloblast cells have medium side scatter (SS) and
CD45 expressions. The red lines indicate approximate myeloblast boundaries
(located on the left-most subfigures in each row and extended horizontally to
the subfigures on the right) and confirm that these meta-clusters represent
immunophenotypes of myeloblast cells. Blue vertical lines denote the +/-
boundaries of a marker. Gray subplots show contour plots of dominant markers
defining the meta-cluster in the same row. (a) HLA-DR+CD117+CD34+CD38+ meta-
cluster shared by 11 AML samples in Tube 6. (b) HLA-DR-CD117±CD34-CD38+ meta-
cluster shared by 5 AML samples in Tube 6. This meta-cluster is indicative of
acute promyelocytic leukemia (APL). These bivariate plots are shown for
illustration only, since the populations of specific cell types are identified
from seven-dimensional data.

Four statistical measures evaluating the performance of the template-based
classification in the training set and test set of the AML data. The
statistical measures are computed for each tube separately and two
combinations of tubes.  Tubes Training set Test set Precision Recall
Specificity F-value Precision Recall Specificity F-value 4 0.94 0.74 0.99 0.83
1.00 0.75 1.00 0.86 5 0.75 0.91 0.96 0.82 0.65 0.85 0.94 0.74 6 1.00 0.70 1.00
0.82 1.00 0.80 1.00 0.89 All (2-7) 1.00 0.74 1.00 0.85 1.00 0.85 1.00 0.92
4,5,6 1.00 0.96 1.00 0.98 1.00 1.00 1.00 1.00

Figure 5: Average classification score from Tubes 4,5,6 for each sample in the
(a) training set and (b) test set. Samples with scores above the horizontal
line are classified as AML, and as healthy otherwise. The actual class of each
sample is also shown. An AML sample (subject id 116) is always misclassified
in the training set, and this is discussed in the text.

###  7.4 Impact of each tube in the classification

As discussed in the methods section, we build six independent classifiers
based on the healthy and AML templates created from Tubes 2-7 of the AML
dataset. A sample is classified as an AML sample if the classification score
is positive, and as a healthy sample otherwise. Let true positives (TP) be the
number of AML samples correctly classified, true negatives (TN) be the number
of healthy samples correctly classified, false positives (FP) be the number of
healthy samples incorrectly classified as AML, and false negatives (FN) be the
number of AML samples incorrectly classified as healthy. Then, we evaluate the
performance of each template-based classifier with the well-known four
statistical measures: Precision, Recall(Sensitivity), Specificity, and
F-value, defined as $\text{Precision}=\frac{\text{TP}}{\text{TP}+\text{FP}}$,
$\text{Recall(Sensitivity)}=\frac{\text{TP}}{\text{TP}+\text{FN}}$,
$\text{Specificity}=\frac{\text{TN}}{\text{FP}+\text{TN}}$, and
$\text{F-value}=\frac{2(\text{Precision}\times\text{Recall})}{\text{Precision}+\text{Recall}}$.
These four measures take values in the interval [0,1], and the higher the
values the better the classifier.

First, we evaluate the impact of each tube in the classification of the
training samples. For a training sample $X$, the classification score is
computed by comparing it with the healthy and AML templates created from the
training set after removing $X$. The predicted status of $X$ is then compared
against true status to evaluate the classification accuracy. Table 7.3 (left
panel) shows various statistical measures for the classifiers defined in Tubes
2-7 of the training set. The classifiers based on Tubes 4, 5, and 6 have the
highest sensitivity because these tubes include several markers relevant to
AML diagnosis (Campana and Behm, 2000; Paietta, 2003). The number of true
negatives TN is high in every tube since the identification of healthy samples
does not depend on the detection of AML-specific markers. Hence specificity is
close to one for all tubes. Analogously, FP is low for most tubes, and we
observe high precision for most tubes. The F-value is a harmonic mean of
precision and recall, and denotes the superior classification ability of
markers in Tubes 4-6. Averaging scores from all tubes does not improve the
sensitivity and F-value dramatically. However, combining Tubes 4-6 gives
almost perfect classification with one misclassification for the training set.
We plot the average classification scores from Tubes 4-6 for the training
samples in Fig. 5(a). The class labels of samples are also shown (blue circles
for healthy and red triangles for AML samples).

In Fig. 5(a), we observe an AML sample (subject id 116) with score below the
classification boundary. In this subject, the proportion of myeloid blasts is
4.4%, which is lower than the minimum 20% AML blasts necessary to recognize a
patient to be AML-positive according to the WHO guidelines (Estey and Döhner,
2006) (the FAB threshold is even higher, at $30\%$). Hence this is either a
rare case of AML, or one with minimal residual disease after therapy, or
perhaps it was incorrectly labeled as AML in the training set. Subject 116 was
classified with the healthy samples by methods in other published work (Biehl
et al., 2013).

###  7.5 Classifying test samples

Now we turn to the test samples. For each tube, we compute the classification
score for each sample in the test set using templates created from the
training set and applying Eq. 3. Since the average classification score from
Tubes 4-6 performs best for the training set, we use it as a classifier for
the test set as well. Since the status of test samples was released after the
DREAM6/FlowCAP2 challenge, we can determine the classification accuracy of the
test samples. Fig. 5(b) shows the classification scores of the test samples,
where samples are placed in ascending order of classification scores. In Fig.
5(b), we observe perfect classification in the test set. Similar to the
training set, we tabulate statistical measures for the classifiers in Table
7.3.

When classifying a sample $X$, we assume the null hypothesis: $X$ is healthy
(non-leukemic). The sample $X$ receives a positive score if it contains AML-
specific immunophenotypes, and the higher the score, the stronger the evidence
against the null hypothesis. Since Tube 1 (isotype control) does not include
any AML-specific markers, it can provide a background distribution for the
classification scores. In Tube 1, 174 out of 179 training samples have
negative classification scores, but five samples have positive scores, with
values less than $0.2$. In the best classifier designed from Tubes 4, 5, 6, we
observe that two AML-positive samples in the training set and three AML-
positive samples in the test set have scores between 0 and 0.2. The classifier
is relatively less confident about these samples; nevertheless, the p-values
of these five samples (computed from the distribution in Tube 1) are still
small ($<0.05$), so that they can be classified as AML-positive. The rest of
the AML samples in the training and test sets have scores greater than $0.2$
and the classifier is quite confident about their status (p-value zero).

Four AML samples in the test set (ids 239, 262, 285, and 326) were
subclassified as APL by comparing against distinct template trees for APL and
the other AML samples in the training set (cf. Fig. 3 (b)).

Finally, we state the computational times required on an iMac with four 2.7
GHz cores and 8 GB memory. Our code is in R. Consider a single tube with $359$
samples in it. The $k$-means clustering of all samples took one hour,
primarily because we need to run the algorithm multiple times (about ten on
the average) to find the optimal value of the number of clusters. Creating the
healthy template from $156$ samples in the training set required $10$ seconds
(s) on one core, and the AML template for $23$ AML samples took $0.5$s on one
core. Cross validation (leave one out) of the training set took $30$ minutes,
and computing the classification score for the $180$ test samples took $15$s,
both on four cores. We could have reduced the running time by executing the
code in parallel on more cores. We have made the dominant step, the $k$-means
clustering of all the samples with an optimal number of clusters, faster using
a GPU, reducing the total time to a few minutes.

##  8 Conclusions

We have demonstrated that an algorithmic pipeline for template-based
classification can successfully identify immunophenotypes of clinical interest
in AML. These could be used to differentiate the subtypes of AML, which is
advantageous since prognosis and treatment depends on the subtype. The
templates enable us to classify AML samples in spite of their heterogeneity.
This was accomplished by creating a scoring function that accounts for the
subtleties in cell populations within AML samples. We are currently applying
this approach to a larger AML data set, and intend to analyze other
heterogeneous data sets.

## Acknowledgments

This research was supported by NIH grant IR21EB015707-01, NSF grant
CCF-1218916, and DOE grant 13SC-003242.

## References

  * Aghaeepour et al. (2013) Aghaeepour, N., Finak, G., Hoos, H., Mosmann, T. R., Brinkman, R., Gottardo, R., Scheuermann, R. H., et al. (2013).  Critical assessment of automated flow cytometry data analysis techniques.  Nature Methods, 10(3), 228–238.
  * Azad et al. (2010) Azad, A., Langguth, J., Fang, Y., Qi, A., and Pothen, A. (2010).  Identifying rare cell populations in comparative flow cytometry.  Lecture Notes in Computer Science, 6293, 162–175.
  * Azad et al. (2012) Azad, A., Pyne, S., and Pothen, A. (2012).  Matching phosphorylation response patterns of antigen-receptor-stimulated T cells via flow cytometry.  BMC Bioinformatics, 13(Suppl 2), S10.
  * Azad et al. (2013) Azad, A., Khan, A., Rajwa, B., Pyne, S., and Pothen, A. (2013).  Classifying immunophenotypes with templates from flow cytometry.  In Proceedings of the International Conference on Bioinformatics, Computational Biology and Biomedical Informatics (ACM BCB), page 256. ACM.
  * Bennett et al. (1985) Bennett, J. M., Catovsky, D., Daniel, M. T., Flandrin, G., Galton, D. A., Gralnick, H. R., and Sultan, C. (1985).  Proposed revised criteria for the classification of acute myeloid leukemia: A report of the French-American-British Cooperative Group.  Annals of Internal Medicine, 103(4), 620–625.
  * Biehl et al. (2013) Biehl, M., Bunte, K., and Schneider, P. (2013).  Analysis of flow cytometry data by matrix relevance learning vector quantization.  PLoS One, 8(3), e59401.
  * Campana and Behm (2000) Campana, D. and Behm, F. G. (2000).  Immunophenotyping of leukemia.  Journal of Immunological Methods, 243(1), 59–75.
  * Estey and Döhner (2006) Estey, E. and Döhner, H. (2006).  Acute myeloid leukaemia.  The Lancet, 368(9550), 1894–1907.
  * Finak et al. (2010) Finak, G., Perez, J., Weng, A., and Gottardo, R. (2010).  Optimizing transformations for automated, high throughput analysis of flow cytometry data.  BMC Bioinformatics, 11(1), 546.
  * Halkidi et al. (2001) Halkidi, M., Batistakis, Y., and Vazirgiannis, M. (2001).  On clustering validation techniques.  Journal of Intelligent Information Systems, 17(2-3), 107–145.
  * Keyhani et al. (2000) Keyhani, A., Huh, Y. O., Jendiroba, D., Pagliaro, L., Cortez, J., Pierce, S., Pearlman, M., Estey, E., Kantarjian, H., and Freireich, E. J. (2000).  Increased CD38 expression is associated with favorable prognosis in adult acute leukemia.  Leukemia Research, 24(2), 153–159.
  * Lacombe et al. (1997) Lacombe, F., Durrieu, F., Briais, A., Dumain, P., Belloc, F., Bascans, E., Reiffers, J., Boisseau, M., and Bernard, P. (1997).  Flow cytometry CD45 gating for immunophenotyping of acute myeloid leukemia.  Leukemia, 11(11), 1878–1886.
  * Legrand et al. (2000) Legrand, O., Perrot, J.-Y., Baudard, M., Cordier, A., Lautier, R., Simonin, G., Zittoun, R., Casadevall, N., and Marie, J.-P. (2000).  The immunophenotype of 177 adults with acute myeloid leukemia: proposal of a prognostic score.  Blood, 96(3), 870–877.
  * Manninen et al. (2013) Manninen, T., Huttunen, H., Ruusuvuori, P., and Nykter, M. (2013).  Leukemia prediction using sparse logistic regression.  PLoS One, 8(8), e72932.
  * Mason et al. (2006) Mason, K. D., Juneja, S. K., and Szer, J. (2006).  The immunophenotype of acute myeloid leukemia: is there a relationship with prognosis?  Blood Reviews, 20(2), 71–82.
  * Paietta (2003) Paietta, E. (2003).  Expression of cell-surface antigens in acute promyelocytic leukaemia.  Best Practice & Research Clinical Haematology, 16(3), 369–385.
  * Pyne et al. (2009) Pyne, S., Hu, X., Wang, K., Rossin, E., Lin, T., Maier, L., Baecher-Allan, C., McLachlan, G., Tamayo, P., Hafler, D., et al. (2009).  Automated high-dimensional flow cytometric data analysis.  Proceedings of the National Academy of Sciences, 106(21), 8519–8524.
  * Qiu (2012) Qiu, P. (2012).  Inferring phenotypic properties from single-cell characteristics.  PLoS One, 7(5), e37038.
  * Raspadori et al. (2001) Raspadori, D., Damiani, D., Lenoci, M., Rondelli, D., Testoni, N., Nardi, G., Sestigiani, C., Mariotti, C., Birtolo, S., Tozzi, M., et al. (2001).  CD56 antigenic expression in acute myeloid leukemia identifies patients with poor clinical prognosis.  Leukemia, 15(8), 1161–1164.
  * Shapiro (2005) Shapiro, H. M. (2005).  Practical Flow Cytometry.  Wiley-Liss.
  * Spidlen et al. (2013) Spidlen, J., Barsky, A., Breuer, K., Carr, P., Nazaire, M.-D., Hill, B. A., Qian, Y., Liefeld, T., Reich, M., Mesirov, J. P., et al. (2013).  Genepattern flow cytometry suite.  Source Code for Biology and Medicine, 8(1), 1–8.
  * Walter et al. (2010) Walter, K., Cockerill, P., Barlow, R., Clarke, D., Hoogenkamp, M., Follows, G., Richards, S., Cullen, M., Bonifer, C., and Tagoh, H. (2010).  Aberrant expression of CD19 in AML with t (8; 21) involves a poised chromatin structure and pax5.  Oncogene, 29(20), 2927–2937.


# Switching of ferrotoroidal domains via an intermediate mixed state in the
multiferroic Y-type hexaferrite Ba0.5Sr1.5Mg2Fe12O22

Jiahao Chen  Department of Physics, University of Oxford, Oxford OX1 3PU,
England, United Kingdom  Francis Chmiel  Department of Physics, University of
Oxford, Oxford OX1 3PU, England, United Kingdom  Roger D. Johnson  Department
of Physics and Astronomy, University College London, Gower Street, London WC1E
6BT, United Kingdom  Department of Physics, University of Oxford, Oxford OX1
3PU, England, United Kingdom  Jieyi Liu  Department of Physics, University of
Oxford, Oxford OX1 3PU, England, United Kingdom  Dharmalingam Prabhakaran
Department of Physics, University of Oxford, Oxford OX1 3PU, England, United
Kingdom  Paolo G. Radaelli<EMAIL_ADDRESS>Department of
Physics, University of Oxford, Oxford OX1 3PU, England, United Kingdom

###### Abstract

We report a detailed study of the magnetic field switching of
ferrotoroidal/multiferroic domains in the Y-type hexaferrite compound
Ba_0.5Sr_1.5Mg_2Fe_12O_22. By combining data from SQUID magnetometry, magneto-
current measurements, and resonant X-ray scattering experiments, we arrive at
a complete description of the deterministic switching, which involves the
formation of a temperature-dependent mixed state in low magnetic fields. This
mechanism is likely to be shared by other members of the hexaferrite family,
and presents a challenge for the development of high-speed read-write memory
devices based on these materials.

##  I Introduction

Harnessing the intrinsic functionality of magneto-electric (ME) and
multiferroic (MF) materials in which (anti)ferromagnetic order can be
controlled by electric fields, and conversely ferroelectric order controlled
by magnetic fields [1, 2, 3, 4], may enable the development of high-speed, low
energy magnetic storage devices. A particularly promising group of materials
is the hexaferrite family, which includes several “type-II” multiferroics with
significant electric polarisation induced via the inverse-Dzyaloshinskii-Moria
interaction (i-DMI) upon non-collinear magnetic ordering [5, 6]. Remarkably
low fields are required for switching, and this class of materials could well
be integrated into in MF storage devices provided ordering temperatures are
raised above room temperature [7, 8].

Within the wider hexaferrite family, Y-type hexaferrites were the first found
to show significant ME effects, and have remained of great interest given that
their low-temperature ME coefficients are among the largest observed in any
material. The Y-type compound Ba_0.5Sr_1.5Mg_2Fe_12O_22 (BSMFO) discussed in
this paper displays switching behaviour typical of most Y-type hexaferrites:
at low temperatures and under moderate magnetic fields, BSMFO adopts a complex
magnetic structure, known as ‘two-fan transverse conical’ (TFTC - see below),
which has a sizeable ferrimagnetic moment M and a magnetically-induced
polarisation P that is orthogonal to M and to the 6-fold axis of the crystal
structure. The development and switching of the macroscopic polarisation under
remarkably low magnetic fields are most readily described in terms of a
ferrotoroidal order parameter $\mathbf{T}\propto\mathbf{P}\times\mathbf{M}$,
which has become known as the fourth primary ferroic order parameter [9, 10,
11]. In essence, cooling in a cross-field configuration tends to stabilise a
single or at least a dominant TFTC toroidal domain (which we will refer to as
the ‘mono-domain state’), as demonstrated in our previous work (Figure 7 in
reference 12). The ferrotoroidal domain population is unaffected by subsequent
reversals of the magnetic field at low temperatures, except very near $B=0$
(see below). Indeed, reversing the magnetic or electric fields alone will not
change the sign of the ferrotoroidicity, which can only be reversed by the
application of cross fields. This behaviour leads to the deterministic
switching of the electric polarisation upon reversal of the sample’s
magnetisation for a given, preconditioned, ferrotoroidal state. At the
microscopic level, within a single ferrotoroidal domain, it was proposed that
reversing the magnetic field led to a coherent rotation of the TFTC structure
about an in-plane axis, such that the deterministic switching is mediated via
an intermediate longitudinal-conical (LC) phase, defined below [13]. We later
showed that the rotation actually occurs about the $c$ axis, and therefore not
via an LC phase [12]. Despite this progress, anomalous temperature-dependent
features in the magnetisation, electric polarisation, and resonant X-ray
scattering (RXMS) intensities observed upon domain switching at low fields had
not been explained [14, 15, 12], but were thought to indicate a more complex
switching mechanism; perhaps a two-step process rather than a simple rotation.

In this paper, we present a combined analysis of magnetometry, ME
measurements, and linearly polarised RXMS to show that, while an intermediate
longitudinal-conical phase does not mediate ferrotoroidal domain switching,
its presence in a mixed state at low fields can account for all anomalous
features observed. Hence, we are now able to provide a fully comprehensive
description of deterministic multiferroic domain switching in BSMFO based on a
simple rotation of the TFTC structure within ferrotoridal domains of a unique
and controllable sign.

##  II Crystal and magnetic structures

Figure 1: The crystal structure (a) and magnetic structures (b)-(d) of Y-type
hexaferrites. The magnetic structures are internally complex but are depicted
using the so-called ‘block model’ approximation (see text), in which the
hexagonal unit cell is divided into quasi-collinear ferrimagnetic blocks
stacked along the $c$ axis. (b) At room temperature the magnetic structure is
described by a planar helix in the $ab$ plane. (c) On cooling, the spins cant
slightly towards the $c$ axis and the longitudinal conical (LC) structure
emerges with a small net ferrimagnetic moment parallel to the $c$ axis (less
than 1 $\mu_{B}$/f.u.)[15, 16, 17]. (d) Upon application of a magnetic field
in the $ab$ plane, the cone ‘flops’ so that its axis, and ferrimagnetic
magnetisation, is directed along the field, leading to the two-fan transverse
cone (TFTC) magnetic structure.

Ba_0.5Sr_1.5Mg_2Fe_12O_22 (BSMFO) crystallises in the $R\bar{3}m$ space group,
with lattice parameters $a=b\approx 5.8$ Å and $c\approx 43$ Å. The primitive
unit cell contains six symmetry-inequivalent Fe ions, four in octahedral
coordination and two in tetrahedral coordination. A simple model of the
magnetic structure can be obtained by grouping the Fe sites into quasi-
collinear ferrimagnetic blocks, whereby the intra-block exchange interactions
are largely unfrustrated and predominantly antiferromagnetic. By contrast, the
inter-block interactions are severely frustrated by geometry, which is the
main cause of the hexaferrites’ complex magnetic phase diagram [18]. This
simplification leads to the so called “block model”, in which all the spins in
each block are added together into a “super-spin” proportional to the block’s
ferrimagnetic moment. This approximation can be used to construct
phenomenological models, and provide a reasonable approximation for the
neutron scattering structure factor used to determine the magnetic structure
of each phase [19, 20, 21]. The simplest version of the block model is the
two-block model, in which block L (large “super-spin”) and block S (small
“super-spin”) are taken as the basic units of the magnetic structure.

Figure 1 shows the crystal structure (Figure 1(a)) and three two-block
magnetic structures (Figure 1(b)(c)(d)), which are most relevant for Y-type
hexaferrite BSMFO. At high temperatures, BSMFO is ferrimagnetic, while at room
temperature and zero magnetic field it exhibits a planar helical structure
with proper-screw axis along the c axis and magnetic propagation vector
$k_{0}$ = (0,0,0.9), as shown in Figure 1(b). On cooling below 50 K, a
longitudinal cone (LC) structure with an out-of-plane component of the moment
is formed [16, 22], as shown in Figure 1(c), giving rise to a bulk
ferrimagnetic moment parallel to the $c$ axis. Upon the application of a small
magnetic field within the $ab$ plane, the ferrimagnetic magnetisation, and
hence axes of the cone, rotates to coincide with the field direction, forming
the so called transverse cone (TC) magnetic structure [23, 13]. Although
different variants of the TC structures are reported throughout the Y-type
hexaferrite family, the one relevant to BSMFO is the two-fan TC phase (TFTC),
described by two magnetic propagation vectors $k_{1}$ = (0,0,0) and $k_{2}$ =
(0,0,1.5) [2], whereby $k_{1}$ describes the $\Gamma$ point component of the
two block structure giving rise to the ferrimagnetic magnetisation, while
$k_{2}$ describes the oscillating ‘fan’ component. The TFTC phase is
ferroelectric due to the i-DMI effect, the polarisation being within the $ab$
plane and perpendicular to the ferrimagnetic magnetisation (i.e. to the
applied magnetic field). By contrast, the in-plane incommensurate helical
phase and the LC phase are non-polar and therefore have zero polarisation.

##  III Experimental techniques

The in-plane electric polarisation was measured as a function of perpendicular
in-plane magnetic field at fixed temperature using the magneto-current
technique. A single crystal sample was prepared in a parallel-plate capacitor
geometry, and mounted on a custom probe within a Quantum Design PPMS. Magnetic
fields were swept at a fixed rate of 5 mT/s while logging the perpendicular
magneto-current, which was integrated to obtain the electric polarisation, P.
The in-plane magnetisation was measured using a Quantum Design MPMS-3
magnetometer, with the magnetic field oriented along the $a^{*}$ axis of the
single crystal (the [2,1,0] direction in real space). Field sweeps at fixed
temperature (MvH) were performed at 10, 25, 50, and 125 K. We found that, in
general, the transitions observed in the magnetisation were significantly
_sharper_ than those observed in the magneto-current, likely due to the fact
that the polarisation switching curve is broadened both by the sweep rate of
the magnetic field and by the impedance of the capacitive/inductive component
in the effective circuit. We also note that the accuracy of the above magneto-
current measurements depends both on the conductivity of the sample (which is
very low at low temperature but increases on warming) and on the background
noise.

RXMS measurements at the Fe $L_{3}$ edge (707 eV) were performed on the RASOR
diffractometer at beamline I10, Diamond Light Source (UK). A cleaved single
crystal sample was mounted with the $(00l)$ zone axis surface normal, on the
end of a 4He cryostat equipped with a small resistive electromagnet capable of
applying $\pm 250$ mT within the plane of the sample (full details are given
in reference 12). $\theta/2\theta$ scans that spanned the $(0,0,3)$ Bragg peak
and its magnetic satellites were performed at different magnetic fields and
temperatures, with the helical and TFTC phases identified by their
characteristic propagation vectors (satellite positions) of $(0,0,\pm 0.9)$
and $(0,0,\pm 1.5)$, respectively. The presence of circular dichroism
(normalized difference in diffraction intensity for right- and left-handed
circularly polarised light) provided quick confirmation of the magnetic origin
of the satellite intensities. No circular dichroism was observed at $(0,0,3)$
in the helical phase, while for the TFTC phase this peak displayed weak
circular dichroism since it includes a significant contribution from the
$\Gamma$ point ferrimagnetic order. Magnetic moment vectors of a given Fourier
component of each magnetic structure could be decomposed using linear
polarisation analysis of the diffraction intensity, which was measured in four
channels; $\sigma-\sigma^{\prime}$,$\sigma-\pi^{\prime}$,
$\pi-\sigma^{\prime}$, and $\pi-\pi^{\prime}$, where $\sigma$ denotes the
polarisation vector perpendicular to the scattering plane, $\pi$ denotes the
polarisation vector within the scattering plane, and the prime indicates the
scattered radiation. Detailed calculations of the linear polarisation
dependence on the diffraction intensities can be found in reference [12].

Figure 2: Two-stage magnetic-field-driven electrical polarisation switching of
BSFMO at 10 K. (a)-(d) Current-field (I-H) and polarisation-field (P-H)
hysteresis loops, based on magneto-current measurement collected after cooling
with poling electric and magnetic fields of different orientations (see
boxes): (a) (-E,-H) field; (b) (-E,+H) field; (c) (+E,+H) field; (d) (+E,-H)
field. The poling fields were removed after cooling and prior to measuring.
The starting point of each magnetic field sweep is indicated by red crosses.
The second, smaller peak associated with the second stage of switching is
emphasised by a dark (green) filling. The insets show the directions of the
initial polarisation P and magnetisation M after poling, and the resulting
direction of the toroidal moment along the $c$ axis, which is never changed
during the magnetoelectric switching process.

##  IV Results and discussion

###  IV.1 Ferrotoroidal domains

The measurements reported in the following sections were performed having
preconditioned the sample into a monodomain state by cross-field cooling. The
field poling process began at 200 K when the sample is sufficiently insulating
to maintain the bias electric field. A magnetic field H of $\pm$250 mT was
applied along the $a$ axis (deep into the TFTC phase at low temperature), with
a $\pm$100 V bias voltage simultaneously applied along $b^{*}$ (i.e.
orthogonal to $a$ within $ab$ plane). The sample was cooled to 10 K and then
the electric field bias was removed. Depending on the sign of the two fields
one of four combinations of P and M could be established;
$\\{-\mathbf{P}|-\mathbf{M}\\}$, $\\{-\mathbf{P}|+\mathbf{M}\\}$,
$\\{+\mathbf{P}|+\mathbf{M}\\}$, or $\\{+\mathbf{P}|-\mathbf{M}\\}$ (see
Figure 2 a to d).

The conservation of ferrotoroidicity through magnetic field reversals is
clearly demonstrated in Figure 2. Upon reversing the magnetic field at low
temperatures, a current in the orthogonal direction (black line) is induced
due to the reversal of polarisation. Furthermore, the full magnitude of the
polarisation (blue lines) is completely recovered at high fields after
reversal, regardless of the initial configuration, and is retained even after
multiple cycles. These observations are consistent with the ferrotoroidal
domain population always being restored at high applied magnetic fields, while
both P and M simultaneously switch upon magnetic field reversal.

The magneto-current data and resultant P-H loops indicate that there is
anomalous behaviour in the low-field region where switching occurs. For
example, starting from the (-E, -H) state (in panel a), a sharp positive peak
emerges when the field approaches zero, followed by a second, smaller positive
peak at around +75 mT, giving rise to a low-field plateau in the electric
polarisation within an apparent 2-step switching process. By reversing the
direction of the sweep, the pattern is very similar but all signs are
reversed, with a first large peak at B = 0 and smaller peak near -75 mT.
Different bias field configurations in Figure 2 (b-d) yield the same
qualitative picture, though the magnitude of the saturated polarisation varies
slightly, implying a small bias towards one ferrotoroidal domain over the
other. The reason for this bias is still unknown but similar effects in
electric field poling were identified in Ba_0.4Sr_1.6Mg_2Fe_12O_22 and
attributed to clamping between ferroelectric/two-fan cone domains [24, 23].

###  IV.2 Low-field intermediate mixed state

Figure 3: TFTC/LC phase coexistence through the magnetoelectric switching of
BSFMO, as determined from RXMS data collected at low temperature along the
$(0,0,l)$ zone axis. (a). Peaks at (0,0,3) (left), (0,0,3.9) (center) and
(0,0,4.5) (right) correspond, respectively, to the $\Gamma$ point (both
phases, Thomson and magnetic scattering), the incommensurate magnetic peak of
the LC phase and the commensurate $Z$-point magnetic peak of the TFTC phase.
The red line is a guide to the eye. (b). Intensity and FWHM (full width at
half maximum) as a function of magnetic field of the commensurate (0,0,4.5)
and incommensurate (0,0,3.9) peaks, extracted from Lorentzian fits on peaks by
subtracting the background. The dashed line indicates the field region where
the intensity of the (0,0,3.9) ICM peak is zero.

RXMS was used to investigate the microscopic origin of the anomalous behaviour
observed at low fields during switching at low temperature. We recall that the
RXMS signatures of the TFTC phase are the purely magnetic commensurate
satellites at (0,0,1.5) and (0,0,4.5) (here only the latter was measured),
while the (0,0,3) Bragg peak acquires a magnetic contribution. Figure 3 shows
the evolution of the diffraction pattern as a function of magnetic field. As
the magnetic field is swept from -250 mT towards zero, we observe a decrease
in the intensity of the (0,0,4.5) satellite starting from around -10 mT, while
a broad incommensurate peak at (0,0,3.9) emerges, roughly corresponding to the
propagation vector of the LC phase. We note that the large width of the
(0,0,3.9) peak indicates that this incommensurate magnetic order has a
relatively short correlation length. As highlighted by the dotted guide line,
the (0,0,4.5) peak intensity remains low in the field range between -6 mT and
100 mT, throughout which the ICM peak persists. Interestingly, in this field
regime, the width of the (0,0,4.5) peak also increases, indicating a shorter
correlation length, most likely associated to the formation of ICM domains
that break up larger CM domains. Above 100 mT the ICM peak disappears and the
(0,0,4.5) peak intensity quickly recovers. The (0,0,3) peak (not shown in
full) also exhibits a behaviour indicative of phase coexistence in the low-
field range, which will be discussed in more detail in section IV.4. Similar
measurements performed on sweeping the magnetic field from +250 mT to -250 mT
(not shown here) showed the same qualitative behaviour, with the ICM peak
appearing in the field range +6 mT and -100 mT. If we assume that the
diffraction intensities characteristic of the TFTC and LC structures are
proportional to the respective phase fractions, our RXMS results are
consistent with an intermediate mixed TFTC/LC state occurring in magnetic
fields that exactly coincide with the anomalous features observed in the
electric polarisation switching. Considering these data together, we can
propose the following scenario: as the field passes through zero, the TFTC
phase rotates rapidly, but is also partially converted to the non-polar ICM LC
phase (N.B. this is the only phase observed upon cooling in zero magnetic
field). The mixed state persists up to 100 mT, when the non-polar LC phase
disappears and a pure TFTC state is fully recovered in the same ferrotoroidal
domain as at -250 mT. This results in a sudden increase of the polarisation at
$\sim$ 100 mT, consistent with the second ’peak’ observed in the magneto-
current measurements.

###  IV.3 Field-temperature phase diagram

The field dependent magnetisation (MvH) was measured at 10 K and 50 K, and is
shown in Figure 4 (panels b and e) alongside the respective PvH response
(panels a and d). Intermediate plateaus between two steps in the PvH curves,
which were found to originate in the intermediate mixed state, are mirrored in
the MvH data. This observation is consistent with the expectation that only
the TFTC phase has a finite magnetisation within the hexagonal basal plane,
which is reduced (in the bulk) throughout the mixed state. In this case, the
magnetisation may serve as an indirect measure of the TFTC domain population.
The remaining panels in Figure 4 (c and f) show the MvH data for an
intermediate temperature (25 K) and a higher temperature (125 K). At all
temperatures a zero-field step is present regardless of the sweeping
directions, followed by a second hysteretic step at higher fields. However,
the relative height of the zero-field step decreases upon warming, while the
relative height of the hysteretic second step increases and is pushed to
higher fields. The hysteresis of the second steps also narrows at high
temperatures, indicating a reduced coercivity. These observations imply that
on warming, deterministic switching of the ferroelectric polarisation due to
the rotation of the TFTC phase is maintained at zero-field, but in a reduced
ferrotoridal domain fraction, while the phase fraction of the intermediate
mixed state increases correspondingly.

Figure 4: (Colour online) Magnetoelectric switching of BSFMO hexaferrite at
different temperatures. (b,c,e, and f) MvH loops and (a and d) PvH loops,
measured at different temperatures with the magnetic field applied along the
along the $a$ axis and the polarisation measured along the $b^{*}$ direction.
Open/closed symbols correspond to virgin/subsequent cycles, respectively (see
arrows for the field directions). A field-derivative curve is provided at the
bottom for convenience. In panels (a) and (d), the bottom curve (in yellow)
actually corresponds to the absolute value of the magneto-current data, which
are integrated to yield the polarisation (green).

The field-temperature phase diagram shown in Figure 5 was derived from a
series of MvH measurements at different temperatures up to 120 K (the maximum
temperature was limited to ensure accurate electric field poling). The red and
blue dots connected by dashed lines of the same colors correspond,
respectively, to the positions of the low-field and high-field steps (peaks in
the first derivative of the magnetisation – see Fig. 4), with the area in
between (shaded in light blue) representing the hysteretic region. The TFTC
structure is always fully restored at high field, as shown in gray filling,
and the region where the LC and TFTC structures always coexist in a mixed
state regardless of history is shown in red.

Figure 5: Field-temperature phase diagram of BSFMO. Red and blue symbols
connected by dashed lines are extracted from magnetisation data (see Fig. 4)
and indicate the fields at which the secondary hysteretic steps occur (see
text). The central red shaded area correspond to the region where TFTC/LC
phase coexistence always occurs, while the light blue shaded area is the
hysteretic region. In the solid grey area only the TFTC phase is observed.

###  IV.4 Modelling RXMS data through the field switching transition

Figure 6: (a) Temperature- and field- dependent phase fraction of the TFTC
structure, calculated using the $\delta(H,T)$ function in Eq. IV.4 and
parameters extracted from the fitting to the magnetisation data (Table 1).
Solid (dot) curves depict the fraction values change by increasing
(decreasing) the field. (b) The models of magnetisation change by considering
the $\delta(H,T)$ function at different temperature. Inset: Temperature
dependence of the TFTC phase fraction $\delta_{min}(T)$ that ‘survives’ during
the switching process in the mixed phase at $H=0$.

To corroborate the above switching behaviour, we demonstrate here that a
‘minimal’ model directly based on the magnetisation data is capable of fitting
some of the fine details of our RXMS data in all polarisation channels. In
order to achieve this, our model must include the relative fractions of the
TFTC and LC phases. In the following we assume that no other phases are
present, and that only the TFTC phase contributes to the in-plane
magnetisation, which for a given magnetic domain we take to be independent of
field (i.e. no change in the TFTC cone angle).

We introduce the phase fraction parameter $\delta=V_{TFTC}/(V_{TFTC}+V_{LC})$
where $V_{TFTC}$ and $V_{LC}$ are the volumes of the two phases such that
$\delta=1(0)$ for pure TFTC (LC) phases, respectively. The parameter $\delta$
is both field- and temperature dependent, and it is strongly hysteretic as a
function of field – in other words, the up-sweep (represented by $+$) and
down-sweep (represented by $-$) functions $\delta^{+}(H,T)$ and
$\delta^{-}(H,T)$ are different, being mirror image of each other through
$H=0$ (see below).

The phase diagram in Figure 5 implies that the form of $\delta^{-}(H,T)$ vs
$H$ is that of an inverted top hat, dropping from 1 at high fields to a
temperature-dependent minimum value $\delta_{min}(T)$ upon decreasing the
field, and rising again to 1 when the opposite field is increased. We have
modelled the $\delta(H,T)$ function by a linear combination of back-to-back
error functions:

$\displaystyle\delta(H,T)$
$\displaystyle=1+\frac{(1-\delta_{min}(T))}{2}\left(\mathrm{erfc}\left[\frac{\mu_{0}H-\mu_{0}H_{2}(T)}{w_{c}(T)}\right]\right.$
$\displaystyle\quad\left.-\mathrm{erfc}\left[\frac{\mu_{0}H-\mu_{0}H_{1}(T)}{w_{c}(T)}\right]\right)$
(1)

where $w_{c}(T)$ represents the temperature dependent ‘speed of inter-
conversion’ of the phases (which may also depend on the field sweep rate), and
$H_{1}(T)$ and $H_{2}(T)$ are the characteristic fields defining the sides of
the top hat, which are also temperature dependent. We note that the symmetry
of the hysteresis loop is maintained when $H_{1}^{-}\simeq-H_{2}^{+}\ll
H_{1}^{+}\simeq-H_{2}^{-}$, consistent with Figure 5. The $\delta_{min}(T)$
and $w_{c}(T)$ functions can be estimated from the magnetometry measurements,
while the characteristic fields can be read directly from the phase diagram in
Figure 5.

The function $\delta(H,T)$ thus modelled is plotted in Figure 6 (a) for
different temperatures and sweep directions, with the corresponding parameters
listed in Table 1 – only the values for field up-sweep (solid curves) are
listed.

Temperature |  $\delta_{min}$ |  $\mu_{0}H_{1}^{+}$(mT) |  $\mu_{0}H_{2}^{+}$(mT) |  $w_{c}$(mT)
---|---|---|---|---
10 K |  0.75 |  70 |  10 |  50
25 K |  0.64 |  80 |  5 |  50
50 K |  0.18 |  100 |  -25 |  20
125 K |  0.18 |  205 |  -175 |  18
Table 1: Parameters for the TFTC fraction function $\delta(H,T)$.

Having modelled the field-dependent fractions of the TFTC and LC phases, we
now consider the magnetisation of the TFTC phase alone, which switches upon
rapid rotation about the $c$ axis. Following reference 12, we describe the
TFTC rotation angle about $c$ through a simple non-hysteretic function of
field. Since the magnetometer measures the magnetization _parallel_ to the
magnetic field, the signal from the TFTC phase can be expressed as:

$m_{\parallel}(H,T)=\lvert m(H,T)\rvert\sin\zeta$ (2)

where

$\zeta=\frac{\pi}{2}\tanh\left(\frac{\mu_{0}H}{v_{c}(T)}\right)$ (3)

is the cone rotation angle, whose low-field switching behaviour has been
modelled with a hyperbolic cotangent function as we had done previously in
ref. 12.

In Eq. 3, $m(H,T)$ is the magnetisation of the TFTC phase, while the parameter
$v_{c}(T)$ described the steepness of the magnetisation vs field curves (the
‘velocity of switching’) and is itself temperature dependent. Note that, in
fitting the magnetisation data, we have allowed $\lvert m(H,T)\rvert$ to be
field dependent through the linear relation $m(H,T)=m(0,T)+\chi\ H$, which is
necessary for an accurate fit of the high-field behaviour in fig. 6. At the
microscopic level, this would correspond to a field-induced reduction of the
TFTC cone angle. By combining the magnetisation rotation and phase fraction
information, we obtain a single equation to fit the temperature- and field-
dependent measured magnetisation data $m_{\mathrm{meas.}}(H,T)$:

$m_{\mathrm{meas.}}(H,T)=\delta(H,T)\;m_{\parallel}(H,T)$ (4)

where $m_{\parallel}(H,T)$ is from Eq. 2.

This model provides a satisfactory fit to the magentometry data, as shown in
Figure 6 (b), with the fitting parameters of equation 4 reported in Table 2.

Temperature |  10 K |  25 K |  50 K |  125 K
---|---|---|---|---
$v_{c}(T)$ (mT) |  25 |  25 |  50 |  55
$\lvert m(0,T)\rvert$ ($\mu_{B}$/f.u) |  0.78 |  0.77 |  0.75 |  0.66
$\chi$ (10${}^{-4}\mu_{B}/mT$) |  3.53 |  3.63 |  3.87 |  4.60
Table 2: Temperature dependent parameters for SQUID switching model in Eq. 2.

The minimum TFTC phase fraction, $\delta_{min}$, is plotted as a function of
temperature in the inset of panel 6 (b) (see also Table 1). $\delta_{min}$
drops suddenly around 50 K, indicating that the TFTC becomes unstable (rather
than metastable) in zero fields above this temperature. The susceptibility
$\chi$ can be taken as temperature-independent, while $\lvert m(0,T)\rvert$
changes little in this temperature range, since the magnetic ordering
temperature is above room temperature.

Figure 7: Calculated (dots) and measured (lines) scattering intensity for the
(0,0,3) (panels (a)-(c))and (0,0,4.5) (panels (d)-(f)) reflections, in three
linear polarization channels during the magnetic polarity switching process at
34 K. Yellow/blue lines and symbols are for field up-sweeps/down-sweeps,
respectively. Dashed lines are calculated using the model in reference 12,
while solid lines are calculated using the model described in the present
paper, which takes into account the TFTC phase fraction.

The parameters obtained from fitting the magnetometry data can now be used as
a starting point to refine the modelling of our RXMS data through the magnetic
field switching transition. Very similar data were previously presented in
ref. 12, which also contains a detailed analysis of the magnetic scattering.
Very briefly, RXMS data were collected at the Fe $L_{3}$ edge (E=707 eV) on
the $(0,0,3)$ and $(0,0,4.5)$ reflections, with linear incident X-ray
polarisation, $\sigma$ or $\pi$, and with an analyser selecting the final
linear polarisation to be detected, $\sigma^{\prime}$ or $\pi^{\prime}$. We
measured field-dependent reflection intensities at T = 40 K in the three
channels $\sigma-\pi^{\prime}$, $\pi-\sigma^{\prime}$, and $\pi-\pi^{\prime}$,
which are most sensitive to the components of the magnetic structure (see ref.
12 for details). The $(0,0,4.5)$ reflection is purely magnetic, while the
$(0,0,3)$ reflection has a Thomson (structural) component in the
$\pi-\pi^{\prime}$ channel but not in the other two channels which are purely
magnetic. In our geometry, the $\sigma-\pi^{\prime}$, $\pi-\sigma^{\prime}$
channels are sensitive to the angular momentum components in the scattering
plane, i.e., $a^{*}$ (along the magnetic field) and $c$, whereas the
$\pi-\pi^{\prime}$ channel is sensitive to the $b$ component (orthogonal to
the field), and its magnetic contribution is only activated during the field
rotation. In ref. 12, the magnetic scattering was modelled using a modified
version of the ‘standard’ TFTC phase, which was found to be necessary to
explain the RXMS. This m-TFTC model is also used as a basis to fit the data in
the present paper. In ref. 12, RXMS data were analysed under the assumption of
a rigid rotation of the TFTC phase and no change in phase fraction, which
yielded a qualitatively acceptable agreement with the experimental
observations. In particular, the rigid TFTC rotation model successfully
describes the fast change in the RXMS diffraction intensity in the different
$I_{channels}$ near $H=0$, as the magnetic cone structure rotates rapidly and
non-hysterically through 180∘ about the $c$ axis.

In spite of the qualitative similarities, there are obvious differences
between the experimental and calculated RXMS profiles in ref. 12. In
particular, the data display clear anomalies at intermediate magnetic fields
that are not reproduced by the model. The low-temperature, intermediate-field
anomalies observed in both magnetisation and polarisation data appear at very
similar magnetic fields to those found in the RXMS data, strongly suggesting
that they have the same origin. In principle, the LC phase is weakly
ferromagnetic and therefore contributes to the $\Gamma$-point magnetic signal,
but in practice the magnetic moment from magnetometry is extremely small[15,
16, 17] ($\approx 0.2\mu_{B}$/f.u. from page 44 in [25]) and no such
scattering, with either linear or circular polarisations, was previously
detected. Hence, to good approximation, the TFTC phase is the only one
contributing to magnetic scattering at $(0,0,3)$, and uniquely contributes to
the magnetic scattering at $(0,0,4.5)$. With this assumption, a minimal
revised model is obtained simply by multiplying the calculated intensities in
each channel from reference 12 by the TFTC phase fraction $\delta(H,T)$ as
deduced above (N.B. we use $\delta(H,T)$ determined at 25 K). The intensities
in each channel, for $\Gamma$ or $Z$ point diffraction peaks, are therefore:

$\begin{split}I_{\Gamma/Z}^{\sigma-\pi}(H,T)&=\delta(H,T)\left|F^{(1)}[(M_{\Gamma/Z}^{x}\cos\zeta-
M_{\Gamma/Z}^{y}\sin\zeta)\cos\theta+M_{\Gamma/Z}^{z}\sin\theta]\right|^{2}\\\
I_{\Gamma/Z}^{\pi-\sigma}(H,T)&=\delta(H,T)\left|F^{(1)}[(M_{\Gamma/Z}^{x}\cos\zeta-
M_{\Gamma/Z}^{y}\sin\zeta)\cos\theta-M_{\Gamma/Z}^{z}\sin\theta]\right|^{2}\\\
I_{\Gamma/Z}^{\pi-\pi}(H,T)&=\delta(H,T)\left|F^{(0)}\cos 2\theta-
iF^{(1)}[(M_{\Gamma/Z}^{x}\sin\zeta+M_{\Gamma/Z}^{y}\cos\zeta)\sin
2\theta]\right|^{2}\\\ \end{split}$ (5)

where $\zeta$ is the rotation angle of the m-TFTC moment around the $c$ axis
and is strongly field-dependent, $M_{\Gamma/Z}^{x}$, $M_{\Gamma/Z}^{y}$ and
$M_{\Gamma/Z}^{z}$ are the components of the magnetic interaction vectors of
the m-TFTC phase for the $(0,0,3)$ ($\Gamma$ point) and $(0,0,4.5)$
(commensurate $Z$ point) – see [12], eq. B10, and $2\theta$ is the scattering
angle. The internal parameters of the m-TFTC phase and the magnetic field
dependence of $\zeta$ were kept as in [12], while the magnetic field
dependence of $\delta$ was taken from the magnetometry data at 25 K, since no
corresponding data at 40 K was available.

Figure 7 displays a comparison between the experimental RXMS data (filled
circles), the original rigid-rotation model from reference 12 (short dash) and
our new model including LC/TFTC intermediate mixed state (solid line). The top
and bottom rows are the RXMS results obtained at the $(0,0,3)$ ($\Gamma$
point) and $(0,0,4.5)$ (commensurate $Z$ point) positions, respectively in the
different cross-polarisation channels $\sigma-\pi^{\prime}$,
$\pi-\sigma^{\prime}$ and $\pi-\pi^{\prime}$. The new model is much closer to
the experimental data, with the intermediate-field anomalies being well
reproduced.

As explained above, this ‘minimal’ model has no adjustable parameters other
than an overall scale factor, since all its parameters were taken either from
our previous work [12] or from our magnetometry data. An obvious refinement of
the model would entail introducing the field dependence of the main fan angle
(referred to as $\phi$ in [12], Eq. B10), consistent with the non-zero
magnetic susceptibility we observed in the magnetometry. Although a refinement
of these and other parameters would undoubtedly produce a better fit to the
data, including the high-field slopes of panels (a) and (b) in 7, we believe
that our minimal model provides sufficient confirmation of the switching model
without introducing additional complexity.

##  V Conclusions

In conclusion, we have reported a detailed study of the field switching
behaviour in the multiferroic Y-type hexaferrite Ba_0.5Sr_1.5Mg_2Fe_12O_22
(BSMFO), using SQUID magnetometry, magneto-current measurements and resonant
X-ray magnetic scattering (RXMS) on single crystals that were poled using
crossed magnetic and electric fields. Both magnetisation and electrical
polarisation display an apparent multi-step hysteretic behaviour upon sweeping
the external magnetic field at low temperatures. Similar behaviour was
previously reported, but not fully understood, in related materials. By
combining these data sets, we are able to construct a model that closely
reproduces all our observations. In particular, we demonstrated that the
multiferroic two-fan transverse cone (TFTC) phase, which is stable at low
temperatures and at high magnetic fields, partially transforms into the
longitudinal-cone (LC) phase upon field reversal, giving rise to an
intermediate mixed state with reduced magnetisation and polarisation. At
higher magnetic fields, the LC phase transforms back into TFTC yielding once
again a pure state with maximum magnetisation and polarisation – remarkably,
without any need to re-apply an external bias. We have also shown that this
phenomenology is highly temperature-dependent, with the field of stability and
phase fraction of the LC phase notably broadening/increasing upon warming.
Combining this new insight with our previous model of the TFTC cone rotation
through $B=0$ produces a much better agreement with the RXMS data, yielding an
essentially complete picture of the field-driven polarity switching in Y-type
BSMFO, which can be likely extended to other multiferroic hexaferrites.

## Acknowledgments

We acknowledge Beamline I10, Diamond Light Source (DLS) under Proposal No.
SI17388 and No. SI20608, and thank P. Steadman and R. Fan for assisting the
RXMS experiment. We thank J-C Lin and H. Jani for useful discussion. Jiahao
Chen acknowledges the financial support from China Scholarship Council (CSC).

## References

  * Eerenstein  _et al._ [2006] W. Eerenstein, N. Mathur, and J. F. Scott, Multiferroic and magnetoelectric materials, Nature 442, 759 (2006).
  * Kimura [2012] T. Kimura, Magnetoelectric hexaferrites, Annual Review of Condensed Matter Physics 3, 93 (2012).
  * Nan  _et al._ [2008] C.-W. Nan, M. Bichurin, S. Dong, D. Viehland, and G. Srinivasan, Multiferroic magnetoelectric composites: Historical perspective, status, and future directions, Journal of Applied Physics 103, 1 (2008).
  * Hill [2000] N. A. Hill, Why are there so few magnetic ferroelectrics, The Journal of Physical Chemistry B 104, 6694 (2000).
  * Katsura  _et al._ [2005] H. Katsura, N. Nagaosa, and A. V. Balatsky, Spin current and magnetoelectric effect in noncollinear magnets, Phys. Rev. Lett. 95, 057205 (2005).
  * Mostovoy [2006] M. Mostovoy, Ferroelectricity in spiral magnets, Phys. Rev. Lett. 96, 067601 (2006).
  * Kocsis  _et al._ [2019] V. Kocsis, T. Nakajima, M. Matsuda, A. Kikkawa, Y. Kaneko, J. Takashima, K. Kakurai, T. Arima, F. Kagawa, Y. Tokunaga, _et al._ , Magnetization-polarization cross-control near room temperature in hexaferrite single crystals, Nature Communications 10, 1 (2019).
  * Kitagawa  _et al._ [2010] Y. Kitagawa, Y. Hiraoka, T. Honda, T. Ishikura, H. Nakamura, and T. Kimura, Low-field magnetoelectric effect at room temperature, Nature Materials 9, 797 (2010).
  * Van Aken  _et al._ [2007] B. B. Van Aken, J.-P. Rivera, H. Schmid, and M. Fiebig, Observation of ferrotoroidic domains, Nature 449, 702 (2007).
  * Sawada and Nagaosa [2005] K. Sawada and N. Nagaosa, Optical magnetoelectric effect in multiferroic materials: evidence for a lorentz force acting on a ray of light, Phys. Rev. Lett. 95, 237402 (2005).
  * Fiebig and Spaldin [2009] M. Fiebig and N. A. Spaldin, Current trends of the magnetoelectric effect, The European Physical Journal B 71, 293 (2009).
  * Chmiel  _et al._ [2019] F. P. Chmiel, D. Prabhakaran, P. Steadman, J. Chen, R. Fan, R. D. Johnson, and P. G. Radaelli, Magnetoelectric domains and their switching mechanism in a $\mathrm{Y}$-type hexaferrite, Phys. Rev. B 100, 104411 (2019).
  * Ishiwata  _et al._ [2008] S. Ishiwata, Y. Taguchi, H. Murakawa, Y. Onose, and Y. Tokura, Low-magnetic-field control of electric polarization vector in a helimagnet, Science 319, 1643 (2008).
  * Taniguchi  _et al._ [2008] K. Taniguchi, N. Abe, S. Ohtani, H. Umetsu, and T. H. Arima, Ferroelectric polarization reversal by a magnetic field in multiferroic $\mathrm{Y}$-type hexaferrite $\mathrm{Ba}_{2}\mathrm{Mg}_{2}\mathrm{Fe}_{12}\mathrm{O}_{22}$, Applied Physics Express 1, 031301 (2008).
  * Sagayama  _et al._ [2009] H. Sagayama, K. Taniguchi, N. Abe, T.-h. Arima, Y. Nishikawa, S.-i. Yano, Y. Kousaka, J. Akimitsu, M. Matsuura, and K. Hirota, Two distinct ferroelectric phases in the multiferroic $\mathrm{Y}$-type hexaferrite $\mathrm{Ba}_{2}\mathrm{Mg}_{2}\mathrm{Fe}_{12}\mathrm{O}_{22}$, Phys. Rev. B 80, 180419 (2009).
  * Ishiwata  _et al._ [2010] S. Ishiwata, D. Okuyama, K. Kakurai, M. Nishi, Y. Taguchi, and Y. Tokura, Neutron diffraction studies on the multiferroic conical magnet $\mathrm{Ba}_{2}\mathrm{Mg}_{2}\mathrm{Fe}_{12}\mathrm{O}_{22}$, Phys. Rev. B 81, 174418 (2010).
  * Kocsis  _et al._ [2020] V. Kocsis, T. Nakajima, M. Matsuda, A. Kikkawa, Y. Kaneko, J. Takashima, K. Kakurai, T. Arima, Y. Tokunaga, Y. Tokura, and Y. Taguchi, Stability of multiferroic phase and magnetization-polarization coupling in $\mathrm{Y}$-type hexaferrite crystals, Phys. Rev. B 101, 075136 (2020).
  * Utsumi  _et al._ [2007] S. Utsumi, D. Yoshiba, and N. Momozawa, Superexchange interactions of ($\mathrm{Ba}_{1-x}\mathrm{Sr}_{x})_{2}\mathrm{Zn}_{2}\mathrm{Fe}_{12}\mathrm{O}_{22}$ system studied by neutron diffraction, Journal of the Physical Society of Japan 76, 034704 (2007).
  * Momozawa and Yamaguchi [1993] N. Momozawa and Y. Yamaguchi, Field-induced commensurate intermediate phases in helimagnet ($\mathrm{Ba}_{1-x}\mathrm{Sr}_{x})_{2}\mathrm{Zn}_{2}\mathrm{Fe}_{12}\mathrm{O}_{22}$(x=0.748), Journal of the Physical Society of Japan 62, 1292 (1993).
  * Lee  _et al._ [2012] H. B. Lee, S. H. Chun, K. W. Shin, B.-G. Jeon, Y. S. Chai, K. H. Kim, J. Schefer, H. Chang, S.-N. Yun, T.-Y. Joung, and J.-H. Chung, Heliconical magnetic order and field-induced multiferroicity of the $\mathrm{Co}_{2}$ $\mathrm{Y}$-type hexaferrite $\mathrm{Ba}_{0.3}\mathrm{Sr}_{1.7}\mathrm{Co}_{2}\mathrm{Fe}_{12}\mathrm{O}_{22}$, Phys. Rev. B 86, 094435 (2012).
  * Momozawa  _et al._ [2001] N. Momozawa, Y. Nagao, S. Utsumi, M. Abe, and Y. Yamaguchi, Cation distribution and helimagnetic structure of the $\mathrm{Ba}_{2}(\mathrm{Zn}_{1-x}\mathrm{Mg}_{x})_{2}\mathrm{Fe}_{12}\mathrm{O}_{22}$ system as revealed by magnetization measurements and neutron diffraction, Journal of the Physical Society of Japan 70, 2724 (2001).
  * Nakamura  _et al._ [2012] S. Nakamura, Y. Tsunoda, and A. Fuwa, Mössbauer study on $\mathrm{Y}$-type hexaferrite $\mathrm{Ba}_{2}\mathrm{Mg}_{2}\mathrm{Fe}_{12}\mathrm{O}_{22}$, Hyperfine Interactions 208, 49 (2012).
  * Zhai  _et al._ [2017] K. Zhai, Y. Wu, S. Shen, W. Tian, H. Cao, Y. Chai, B. C. Chakoumakos, D. Shang, L. Yan, F. Wang, _et al._ , Giant magnetoelectric effects achieved by tuning spin cone symmetry in $\mathrm{Y}$-type hexaferrites, Nature Communications 8, 1 (2017).
  * Zhai  _et al._ [2018] K. Zhai, Y. Chai, J. Cong, D. Shang, and Y. Sun, Electric control of exchange bias in multiferroic hexaferrite $\mathrm{Ba}_{0.4}\mathrm{Sr}_{1.6}\mathrm{Mg}_{2}\mathrm{Fe}_{12}\mathrm{O}_{22}$, Phys. Rev. B 98, 144405 (2018).
  * Chmiel [2018] F. Chmiel,  _Magnetic order in functional iron oxides_ , Ph.D. thesis, University of Oxford (2018).