# Reproducibility in Cytometry:
Signals Analysis and its Connection to Uncertainty Quantification

Paul N. Patrone<EMAIL_ADDRESS>Matthew DiSalvo  Anthony J. Kearsley
Geoffrey B. McFadden  Gregory A. Cooksey  National Institute of Standards and
Technology
100 Bureau Drive, Gaithersburg, MD 20899, USA

###### Abstract

Signals analysis for cytometry remains a challenging task that has a
significant impact on uncertainty. Conventional cytometers assume that
individual measurements are well characterized by simple properties such as
the signal area, width, and height. However, these approaches have difficulty
distinguishing inherent biological variability from instrument artifacts and
operating conditions. As a result, it is challenging to quantify uncertainty
in the properties of individual cells and perform tasks such as doublet
deconvolution. We address these problems via signals analysis techniques that
use scale transformations to: (I) separate variation in biomarker expression
from effects due to flow conditions and particle size; (II) quantify
reproducibility associated with a given laser interrogation region; (III)
estimate uncertainty in measurement values on a per-event basis; and (IV)
extract the singlets that make up a multiplet. The key idea behind this
approach is to model how variable operating conditions deform the signal shape
and then use constrained optimization to “undo” these deformations for
measured signals; residuals to this process characterize reproducibility.
Using a recently developed microfluidic cytometer, we demonstrate that these
techniques can account for instrument and measurand induced variability with a
residual uncertainty of less than 2.5% in the signal shape and less than 1% in
integrated area.

##  I Introduction

In the last 30 years, cytometry has evolved as a powerful technique for
clinical diagnostics, drug development, and biotechnology [1, 2]. This success
has recently motivated fundamental questions and studies designed to better
understand the ultimate capabilities of cytometers, as well as their potential
for quantitative, reproduceable, and even traceable measurements [3, 4, 5, 6].
However, in-depth uncertainty quantification (UQ) has yet to be fully
realized, limiting efforts to refine the metrology aspects of cytometry.

In this context, signals analysis is a challenging and often overlooked task
that has a significant impact on uncertainty. For example, conventional flow
cytometers assume that individual events are well characterized by simple
properties such as the signal area, height, and width [7, 8, 9, 10]. However,
such approaches discard the vast majority of information in the measurement.
Moreover, cytometers only interrogate a particle once per laser region, so
that population variability is inherently convolved with other sources of
uncertainty. As a result, it is difficult to determine whether an event
corresponds to a valid measurand, characterize reproducibility on a per-event
basis, and ultimately, quantify confidence in the measurement process.

The goal of this manuscript is to address such problems by incorporating UQ
directly into the signals analysis. We propose a collection of techniques that
use scale transformations to identify sources of variation in events arising
from changes in particle speed, size, and brightness. The main idea behind
these analyses is to represent signals in terms of low-order mathematical
interpolations and define the relevant transformations in terms of generic
physical models. We then use constrained optimization to map all events onto
one another, which also yields the physical parameters (e.g. size) associated
with the particles. Using a recently developed microfluidic cytometer (see
Ref. [11]), we illustrate how this analysis can: (i) quantify phenomena such
as noise in flow conditions and sample variability; and (ii) estimate the per-
event reproducibility associated with all remaining sources of uncertainty. We
also demonstrate that this analysis provides a foundation for more advanced
signal processing techniques by using it to extract the individual singlets
comprising a multiplet signal, i.e. an event composed of multiple overlapping
singlets.

The need to incorporate UQ directly into signals analysis arises from several
issues unique to cytometers. For example, cells are only measured once per
laser region; moreover, it is difficult to ensure identical flow and optical
conditions for all measurands. Thus, it is impossible to directly characterize
repeatability and reproducibility of any measurement. Signals analysis and
mathematical modeling therefore take on new roles, since they allow us to
answer the questions, “how does variation in measurement conditions deform the
signals, and how can these deformations be undone?” In this way, the theory
approximately reconstructs an imaginary scenario in which one could repeat
measurements on identical particles, thereby overcoming experimental
limitations.

Multiplets provide a similar motivation for our analysis techniques. Because
engineering solutions alone cannot prevent multiplets [12, 13, 14, 15], common
practice discards their information. However, this introduces uncertainty into
population variability estimates [12, 14]. Perhaps worse, common data analysis
strategies that rely on integrated areas cannot distinguish dim doublets from
bright singlets, for example. This leads to an uncomfortable situation in
which even the identity and number of measurands have potentially
unquantifiable uncertainties. Again, we demonstrate that signals analysis is
an appropriate avenue to address such problems, since it allows us to extract
individual singlets by understanding how they combine to make the multiplet.
More generally, these examples illustrate that signals analysis is a useful
tool for untangling sources of uncertainty that become intertwined by the
measurement process.

This example of multiplet deconvolution highlights another theme of our work:
there is a natural feedback between UQ and signals analysis that must be
exploited to realize the full potential of cytometry. While it is obvious that
uncertainty estimates arise from data analysis, we also illustrate the
reverse: new data analyses are enabled by UQ. For example, we recast multiplet
deconvolution as the task of finding the most probable singlets that
reconstruct the original multiplet. Thus, estimating the uncertainty in the
shape of singlets is a prerequisite step. Similar arguments apply to the
related but distinct problem of multiplet detection. More generally, UQ can
inform data analysis by ensuring that results are physically
meaningful.111These observations clarify our stance that UQ is the broad set
of tasks that increase confidence in a measurement.

Because the present work develops new analysis tools, it is useful to perform
validation measurements on particles with well characterized properties. For
this reason, the examples below are restricted to fluorescent microspheres.
However, extensions to cell-based measurements are straightforward, if not
trivial, and we point to relevant issues throughout. We note also that the
current manuscript does not consider issues associated with gating strategies
and quantification of population variability per se. Such topics are left for
future work.

The reader should also note that while our analysis strives to be as general
as possible – for example, we need not specify a form of the laser profile –
we make certain assumptions that may not apply to all cytometers. These relate
to uniformity of the laser profile perpendicular to the flow direction and
light-collection geometric factors. Importantly, these assumptions allow us to
formulate the scale transformations as relatively simple, closed-form
expressions, thereby facilitating downstream numerical optimization. An
alternative would be to tailor our approach to a specific cytometer through
more detailed modeling, but this may introduce significant computational
overhead. Such issues are discussed in Sec. V.

The rest of the manuscript is organized as follows. Section II presents the
main assumptions of our analysis in the context of a generic model of a
cytometry measurement (II.1); derives the associated scale transformations
(II.2); and formulates the optimization problem needed to undo deformations
associated with variable measurement conditions (II.3). Section III validates
these methods using experimental data. Section IV highlights the usefulness of
this analysis for UQ (IV.1) and multiplet deconvolution (IV.2). Section V
considers this work in the greater context of cytometry and signals analysis
thereof. An appendix presents technical aspects of the multiplet
deconvolution.

##  II Motivation and Main Ideas

A key observation motivates our work: under generic conditions and for a given
cytometer, the shapes of all signals are identical up to a set of linear
transformations that depend only on the particle properties. As will become
clear, these transformations can be interpreted as changes of units that bring
the numerical values of different measurements into agreement. This provides a
quantitative framework for comparing events, e.g. to determine if they are
valid measurands, characterize their properties, and estimate variation about
mean behavior. The goal of this section is to specify the conditions under
which this observation holds and develop the tools needed to compare signals.

A key challenge in formulating this analysis is that the measurements are
given in a form that does not permit direct comparison. Part of our task is
therefore to convert the signals we can acquire into ones we can use. The
multitude and complexity of steps involved motivates a separation of the
analysis into three subsections. In Sec. II.1, we propose a model that
characterizes deformation in signals arising from differences in particle
speed and size. Critically, the particles are assumed to have the same
trajectory. In Sec. II.2, we develop the machinery to convert signals acquired
over the same time interval (which is experimentally feasible) to ones having
the same trajectory length. Section II.3 combines these results to solve the
backwards problem: determine the particle properties by undoing the signal
deformations.

###  II.1 Physical Effects Causing Signal Variation

Consider the interplay between physical, chemical, and biological processes
during a cytometry measurement. We take an event to be the fluorescence signal
$f(t)$ collected as a function of time $t$ while a bead or cell passes a laser
excitation region; see Fig. 1. For particles moving with a constant velocity
and parallel to the fluid flow, a model of this process is

$\displaystyle f(t)=\\!\int_{D_{L}}\\!\\!\\!\\!{\rm d}x\,{\rm d}y\,{\rm
d}z\,C(x,y,z-vt)\Psi(x,y,z)\Phi(x,y,z)$ (1)

where $C(x,y,z-vt)$ is the concentration of fluorophores on the particle,
$\Psi(x,y,z)$ is the laser light intensity, $\Phi(x,y,z)$ is the amount of
emitted fluorescent light (per unit laser light and fluorophore number) coming
from $(x,y,z)$ and collected by the photodetector, $v$ is the constant
advection velocity of the particle, $D_{L}$ is the spatial domain over which
the signal is generated, and $z$ is parallel to the direction of flow, i.e.
the axial direction. The goal of our analysis is to quantify how $f(t)$
changes as a function of the bead radius $R$, velocity $v$, and concentration
$C$, so that we can then solve the reverse problem of determining these
parameters from $f(t)$.

Figure 1: Schematic of the signal generation process associated with a
cytometry event. As a bead or cell traverses the laser profile, the
concentration of fluorophores and object shape convolve with the laser profile
to create fluorescent light. This light is converted to a voltage by a
photodetector. See Eq. (1) for a mathematical description of this process.
Conceptually we can imagine the signal as beginning and ending when the
particle reaches an initial and final position as it moves with the fluid. We
refer to this as a constant-trajectory signal. In practice, it is easier to
acquire a signal over a constant time interval, i.e. by approximately
centering the peak in a time window of fixed duration, without regard to the
distance traveled.

We first simplify Eq. (1). Considering, for example, the microfluidic
cytometer in Ref. [11], we assume that Eq. (1) can be reduced to

$\displaystyle f(t)=\int_{D_{L}}{\rm d}x\,{\rm d}y\,{\rm
d}z\,C(x,y,z-vt)\psi(z)$ (2)

where $\psi(z)$ depends only on $z$ and accounts for both the laser and
collection optics. In an abuse of terminology, we henceforth refer to
$\psi(z)$ as the laser profile.

Equation (2) is valid for a wide range of operating conditions. For example,
letting $\boldsymbol{\rho}=(x,y)$ denote the radial direction, Eq. (2) applies
when: (i) differently sized, radially symmetric (in $\boldsymbol{\rho}$)
particles on a single streamline encounter a laser profile and geometric
factors whose product is linear in $\boldsymbol{\rho}$; and (ii) arbitrary
shaped particles on any streamline encounter a laser profile and geometric
factors that are independent of $\boldsymbol{\rho}$. Either case may happen
when the particles are small relative to the laser interrogation
region.222Case (i) arises from symmetry arguments applied to Eq. (1). See Ref.
[11] for relevant details associated with the device studied in this
manuscript, including in particular the focusing strategy used to ensure that
particles do not cross streamlines.

We model $C(x,y,z)$ as a step function of the form

$\displaystyle C(x,y,z)=c\Theta\left(R^{2}-x^{2}-y^{2}-z^{2}\right)$ (3)

where $c$ is a constant fluorophore concentration and the Heaviside function
is defined such that $\Theta(r)=1$ if $r>0$ and $\Theta(r)=0$ otherwise.
Equation (3) corresponds to a bead or cell with a uniform concentration of
fluorophores throughout its volume. Through appropriate modification of the
Heaviside function this expression can describe, for example, surface
concentrations. Moreover, Eq. (3) can be parameterized to model non-radial
deformations associated with deformable cells, although such tasks are beyond
the scope of this manuscript.

In light of these simplifications, Eq. (1) becomes

$\displaystyle f(t)$
$\displaystyle=c\pi\\!\\!\int_{-L_{0}}^{L_{0}}\\!\\!\\!\\!\\!{\rm
d}z\,\\!\left[R^{2}\\!-\\!(z\\!-\\!vt)^{2}\right]\\!\Theta\\!\left[R^{2}\\!-\\!(z\\!-\\!vt)^{2}\right]\\!\psi(z)$
(4)

where where we assume that the laser is fully contained inside the domain
$[-L_{0},L_{0}]$ for some length $L_{0}$. For later convenience, we assume
(without loss of generality) that $L_{0}$ satisfies the inequality

$\displaystyle 0<a<L_{0}-R,$ (5)

where the domain $[-a,a]$ is the smallest set for which $\psi(z)>0$.
Physically, inequality (5) implies an event in which we track a particle
starting before it enters the laser and ending after it fully exits; the
corresponding function $f(t)$ is padded on the left and right by zeros.333We
assume that any constant offset or background has been subtracted from the
fluorescence signal. Because $L_{0}$ is constant, the locations at which we
start and stop tracking the particle are fixed. Taking $t=0$ as the time when
the particle is at the (arbitrary) “center” $z=0$ of the laser, $f(t)$ is
defined on the interval $[-\tau,\tau]$, where $\tau=L_{0}/v$.

Equation (4) implies that increasing or decreasing $v$ “compresses” or
“stretches” the signal in the time domain. Thus, we should be able to solve
the reverse problem: estimate the relative velocities by undoing the
deformations in such a way that the signals coincide. To achieve this, it is
convenient to consider Fourier representations of the signals. We denote the
corresponding transforms of $\psi(z)$ and $[R^{2}-z^{2}]\Theta[R^{2}-z^{2}]$
by $\hat{\psi}(k)$ and $\hat{g}(k;R)$, where a simple computation yields444We
use the normalization
$\mathfrak{F}[\cdot]=\int_{-L_{0}}^{L_{0}}\cdot\exp(-ikz/L_{0}){\rm d}z$,
where $\mathfrak{F}[\cdot]$ is the Fourier transform acting on $\cdot$.

$\displaystyle\hat{g}(k;R)\\!$
$\displaystyle=\\!\begin{cases}\frac{-4RL_{0}^{2}\cos(kR/L_{0})}{k^{2}}+\frac{4L_{0}^{3}\sin(kR/L_{0})}{k^{3}}&\\!k\neq
0\\\ \frac{4}{3}R^{3}&\\!k=0\end{cases}$ (6)

and $k=\pi n$ for any integer $n$. In practice, we assume that $-M\leq n\leq
M$, where $M$ is a mode-cutoff associated with the noise floor of the signal;
see also Sec. III. Unless otherwise stated, all sums over $k$ range from
$-M\pi$ to $M\pi$. Note that $\hat{\psi}(k)$ is unknown.

In light of these assumptions, Eq. (4) becomes

$\displaystyle f(t)$
$\displaystyle=\frac{c\pi}{2L_{0}}\sum_{k}\hat{\psi}(k)\int_{-L_{0}}^{L_{0}}\\!\\!\\!\\!{\rm
d}z\,\left[R^{2}\\!-\\!(z\\!-\\!vt)^{2}\right]$
$\displaystyle\qquad\qquad\qquad\qquad\Theta\left[R^{2}\\!-\\!(z\\!-\\!vt)^{2}\right]e^{ikz/L_{0}},$
(7)

where we have used the fact that

$\displaystyle\psi(z)=\frac{1}{2L_{0}}\sum_{k}\hat{\psi}(k)e^{ikz/L_{0}}.$ (8)

To further simplify this, use inequality (5) to impose the periodicity
relationship $C(x,y,z-vt)=C(x,y,z-vt+2nL_{0})$, where $n\in\mathbb{Z}$ is an
integer. This assumption does not change the signal on the domain
$-L_{0}/v\leq t\leq L_{0}/v$. However, it does allow us to invoke the
convolution theorem [16], which yields

$\displaystyle f(t)$
$\displaystyle=c\pi\sum_{k}\hat{\psi}(k)\hat{g}(k;R)e^{ikt/\tau},$ (9)

Equations (6) and (9) illustrate how the measurand and measurement conditions
affect the signal. Varying the concentration $c$ and velocity $v$ (which
controls $\tau$) alters the height and width of the signal. Changes in
particle size $R$ have more complicated effects as characterized by Eq. (6),
since a particle is convolved with a larger fraction of the laser profile as
$R$ increases. Taking the Fourier transform of Eq. (9) with respect to time
(over the interval $-\tau$ to $\tau$) yields

$\displaystyle\hat{f}(k;\tau)=2c\pi\tau\hat{\psi}(k)\hat{g}(k;R).$ (10)

The inclusion of $\tau$ in the argument of $\hat{f}(k;\tau)$ is to emphasize
the time-domain over which the transform is taken, which is important in the
following sections. Letting $f(t)$ and $f_{0}(t)$ denote test and reference
signals, the laser profile $\hat{\psi}(k)$ is eliminated by taking the ratio

$\displaystyle\frac{v\hat{f}(k;\tau)}{v_{0}\hat{f}_{0}(k;\tau_{0})}=\frac{c\hat{g}(k;R)}{c_{0}\hat{g}(k;R_{0})},$
(11)

where quantities with $0$ subscripts are associated with $f_{0}(t)$ and the
trajectory length $L_{0}$ is assumed to be the same for both particles. In the
remainder of the manuscript, we assume that $\tau_{0}=L_{0}/v_{0}$, i.e. the
reference particle defines the length and time scales against which we compare
all events.

Equation (11) is the key result of this section: it defines the shape of one
signal in terms of another and as a function of the particle properties,
provided the trajectories are identical.

###  II.2 Converting between fixed-time and fixed-trajectory representations

In practice, it is difficult to measure an event over a fixed trajectory
length and variable time domain. It is more common for acquisition to occur in
fixed time windows, which we assume to be $[-\tau_{0},\tau_{0}]$. However, if
$\tau_{0}$ remains fixed, the corresponding trajectory lengths of each
particles must vary, which violates the key assumption of the previous
section, so that we cannot use Eq. (11).555Allowing trajectory length to vary
per-particle would yield varying spectral representations of $\psi(z)$, making
it difficult to eliminate the laser profile via a relationship akin to Eq.
(11).

The resolution to this problem is to identify a mapping that converts a fixed-
time Fourier representation to one on a fixed spatial domain. Observe first
that in Eq. (10), the dimensionless frequency $k$ is common to all signals;
moreover, $\hat{f}(k;\tau)$ is linear in the particle-dependent time-scale
$\tau=L_{0}/v$. Instead of $\hat{f}(k;\tau)$ we are given a measured signal
$\hat{f}_{m}(k;\tau_{0})$, where $\tau_{0}=L/v=L_{0}/v_{0}$ and $L\neq L_{0}$
if $v\neq v_{0}$. The linearity of $\hat{f}(k;\tau)$ in $\tau$ suggests that
there is a transformation matrix $\boldsymbol{\rm X}$ such that

$\displaystyle\hat{f}(k;\tau)=\sum_{k^{\prime}}X_{k,k^{\prime}}(\tau,\tau_{0})\hat{f}_{m}(k^{\prime};\tau_{0}),$
(12)

where $\hat{f}(k;\tau)$ and $\hat{f}_{m}(k^{\prime};\tau_{0})$ are the
constant-trajectory and constant-time representations of the event. The
subscript $m$ emphasizes that the latter is the measured signal. The time
domains associated with the right and left sides of Eq. (12) are
$D_{v}=[-\tau,\tau]$ and $D_{0}=[-\tau_{0},\tau_{0}]$. (The subscript $v$ on
$D_{v}$ is to emphasize that $v$, not $L_{0}$, controls $\tau$.) We
temporarily assume that $\tau$, $\tau_{0}$, and $L_{0}$ are known; in the next
section we show how to find these quantities.

Equation (12) does not itself specify the operator $\boldsymbol{\rm X}$. Its
definition depends on how we wish to extend or truncate a signal when mapping
it to a constant-trajectory domain. There are two equivalent perspectives one
may take: active or passive. In the former, we deform (i.e. expand or
contract) the time-series so that its constant-trajectory domain $D_{v}$
becomes $D_{0}$. This collapses all time-series onto one another. In the
latter perspective, we deform $D_{0}$ to become $D_{v}$, keeping the time-
series unchanged. While both approaches are equivalent, the corresponding
derivations have subtle differences.

To better understand these issues, consider the case in which $v>1$. Clearly
$D_{v}$ is a subset of $D_{0}$. In the passive approach, we are given
$\hat{f}_{m}(k;\tau_{0})$, i.e. the Fourier modes on $D_{0}$, and wish to find
the corresponding modes on the smaller domain $D_{v}$. In this sense,
$\boldsymbol{\rm X}$ is a “domain-restriction” operator. Assume that
inequality (5) holds, so that $f_{m}(t)\to 0$ as $t\to\pm L_{0}/v$. Defining
the restriction $f_{r}(t)=f_{m}(t)$ for $t\in[-\tau,\tau]$, we take the
Fourier transform of $f_{r}(t)$ on $D_{v}$ to find

$\displaystyle\hat{f}(k;\tau)$ $\displaystyle:=\hat{f}_{r}(k;\tau)$
$\displaystyle=\frac{v_{0}}{v}\sum_{k^{\prime}}{\rm
sinc}\\!\left(k\\!-\\!\frac{k^{\prime}v_{0}}{v}\right)\hat{f}_{m}(k^{\prime};\tau_{0}),$
(13)

so that $X_{k,k^{\prime}}:=(v_{0}/v){\rm sinc}(k-k^{\prime}v_{0}/v)$ and the
$:=$ symbol means that we define the left-hand side to be equal to the right-
hand side. We adopt the convention that

$\displaystyle{\rm sinc}(x)=\sin(x)/x.$ (14)

To arrive at this result via the active perspective, we scale the argument of
$f_{m}(t)$ by a factor $v_{0}/v$, so that $f_{m}(t)\to f_{m}(v_{0}t/v)$. Then,
restricting $f_{m}(v_{0}t/v)$ to the domain $D_{0}$ and taking the Fourier
transform yields an estimate that we denote $\hat{f}_{r}(k;\tau_{0})$ in a
slight abuse of notation.666The restriction $f_{r}(t)$ is only defined on
$D_{v}$, so that by the notation $\hat{f}(k;\tau_{0})$ we mean the Fourier
transform of the expanded version of $f_{r}(t)$, which has the constant-
trajectory domain $D_{0}$. To convert this to $\hat{f}(k;\tau)$ (i.e. the
mode-weight on the domain $D_{v}$), we invoke Eq. (11) with $R=R_{0}$ and
$c=c_{0}$, which yields
$\hat{f}_{r}(k;\tau)=(v_{0}/v)\hat{f}_{r}(k;\tau_{0})$. This derivation
suggests the interpretation that $\boldsymbol{\rm X}$ is also a “signal-
expansion” operator, since $v>v_{0}$. We leave it as an exercise to the reader
to show the the corresponding matrix $\boldsymbol{\rm X}$ is identical to that
given by Eq. (13).

When $v<1$, the domain $D_{v}=[-\tau,\tau]$ contains the interval
$[-\tau_{0},\tau_{0}]$, so that in the passive-interpretation,
$\boldsymbol{\rm X}$ is a “domain-extension” operator that extends $f_{m}(t)$
from $D_{0}$ onto $D_{v}$. While there is no unique extension $f_{e}(t)$, we
assume without loss of generality that $f_{m}(t)\to 0$ as
$t\to\pm\tau_{0}$.777Signals violating this condition can be treated as
multiplets, which are considered in Sec. IV. This together with inequality (5)
suggests the definition

$\displaystyle f_{e}(t):=\begin{cases}f_{m}(t)&t\in[-\tau_{0},\tau_{0}]\\\
0&{\rm otherwise}.\end{cases}$ (15)

Again, taking the Fourier transform on $D_{v}$, one finds

$\displaystyle\hat{f}(k;\tau):=\sum_{k^{\prime}}{\rm
sinc}\left(k\frac{v}{v_{0}}-k^{\prime}\right)\hat{f}_{m}(k^{\prime};\tau_{0}),$
(16)

so that $X_{k,k^{\prime}}:={\rm sinc}[k(v/v_{0})-k^{\prime}]$. In the active
interpretation, $\boldsymbol{\rm X}$ is a “signal-contraction” operator; we
leave derivation of the equivalence with Eq. (16) as an exercise for the
reader.

Combining these results, we define $f(t)$ on $D_{v}$ to be

$\displaystyle f(t):=\begin{cases}f_{r}(t)&v\geq v_{0}\\\ f_{e}(t)&v\leq
v_{0},\end{cases}$ (17)

with $\hat{f}(k;\tau)$ being the corresponding Fourier transform. In both Eqs.
(13) and (16), setting $v\to v_{0}$ yields the identity transformation, as
expected.

In practical settings, distinct signals $f_{m}(t)$ may not be identically
centered on the domain $D_{0}$ (i.e. before applying $\boldsymbol{\rm X}$).
However, Eq. (4) assumes that the center of any given particle is at $z=0$
when $t=0$, and moreover $\boldsymbol{\rm X}$ deforms the signal about $t=0$.
Thus, it is necessary to consider transformations of the form $f_{m}(t)\to
f_{m}(t+\Delta t)$. As before, assume that $f_{m}(t)$ is padded by zeros as
$t\to\pm\tau_{0}$, so that cyclic permutations by sufficiently small $\Delta
t$ only wrap zeros around the left and right sides of the peak. Assuming that
$f_{m}(t)$ is periodic with a period of $2\tau_{0}$, one finds that

$\displaystyle f_{m}(t)$
$\displaystyle=\sum_{k}\hat{f}_{m}(k;\tau_{0})e^{ikt/\tau_{0}}$
$\displaystyle\qquad\to\sum_{k}\hat{f}_{m}(k;\tau_{0})e^{ikt/\tau_{0}+ik\Delta
t/\tau_{0}}$ (18)

so that $\hat{f}_{m}(k;\tau_{0})\to\hat{f}_{m}(k;\tau_{0})\exp\left(ik\Delta
t/\tau_{0}\right)$ under time translation. While the padding assumption is
technically not necessary for the validity of Eq. (18), it ensures that
subsequent transformations by $\boldsymbol{\rm X}$ do not truncate non-zero
parts of the signal.

###  II.3 Signal Matching

Assume that we have reference and test signals $f_{0}(t)$ and $f_{m}(t)$, both
mapped to the interval $-\tau_{0}\leq t\leq\tau_{0}$ and sufficiently padded
by zeros on left and right. Both Fourier expansions use the same set of
frequencies. Let $\boldsymbol{\rm k}$ denote the vector of frequencies, and
$\hat{\boldsymbol{\rm f}}_{0}$ and $\hat{\boldsymbol{\rm f}}_{m}$ denote the
corresponding vectors of Fourier coefficients. Also define the matrix
operators $\boldsymbol{\rm X}$, $\boldsymbol{\rm T}$, and $\boldsymbol{\rm S}$
having elements

$\displaystyle X_{k,k^{\prime}}(v,v_{0})$
$\displaystyle:=\begin{cases}(v_{0}/v){\rm sinc}(k-k^{\prime}v_{0}/v)&v\geq
v_{0}\\\ \,{\rm sinc}(kv/v_{0}-k^{\prime})&v\leq v_{0}\end{cases}$ (19)
$\displaystyle T_{k,k^{\prime}}(\Delta t)$
$\displaystyle:=\delta_{k,k^{\prime}}e^{ik\Delta t/\tau_{0}}$ (20)
$\displaystyle S_{k,k^{\prime}}(R;R_{0})$
$\displaystyle:=\delta_{k,k^{\prime}}\hat{g}(k;R_{0})/\hat{g}(k;R)$ (21)

where $\delta_{k,k^{\prime}}=1$ if $k=k^{\prime}$ and
$\delta_{k,k^{\prime}}=0$ if $k\neq k^{\prime}$. Equations (19) and (20)
derive from matrix analogues of Eqs. (13), (16), and (18). The quantity
$\hat{g}(k;R)$ is given by Eq. (6) and defines the relative transformation
associated with changing particle radius. The vector analogue of
$\hat{f}(k;\tau)$ is denoted $\hat{\boldsymbol{\rm f}}(\boldsymbol{\rm
k};\tau)$.

For noise-less signals and a known set of transformation parameters $\Delta
t$, $R$, $R_{0}$, $v$, and $v_{0}$, Eq. (11) implies that

$\displaystyle\boldsymbol{\Lambda}\hat{\boldsymbol{\rm f}}_{m}(\boldsymbol{\rm
k};\tau_{0})-\hat{\boldsymbol{\rm f}}_{0}(\boldsymbol{\rm k};\tau_{0})=0$ (22)

where

$\displaystyle\boldsymbol{\Lambda}=\frac{(v/v_{0})}{(c/c_{0})}\boldsymbol{\rm
S}(R,R_{0})\boldsymbol{\rm X}(v,v_{0})\boldsymbol{\rm T}(\Delta t).$ (23)

The interpretation of Eq. (22) is straightforward. The product
$\boldsymbol{\Lambda}\hat{\boldsymbol{\rm f}}_{m}(\boldsymbol{\rm
k};\tau_{0})$: (i) centers $f_{m}$; (ii) scales it to the domain $D_{v}$;
(iii) matches the radius scaling to $R_{0}$; and (iv) normalizes the amplitude
to the reference signal. The order of operations of the matrices mirrors the
assumptions of the analysis above and cannot be changed.888The transformations
matrices do not all commute with one another. For example, centering a peak
and then truncating its edges can yield a different signal than truncating
followed by centering.

In practice, the transformation parameters must be determined from noisy
signals. This motivates the objective

$\displaystyle\mathcal{L}=\Big{|}\Big{|}\boldsymbol{\Lambda}\hat{\boldsymbol{\rm
f}}_{m}-\hat{\boldsymbol{\rm f}}_{0}\Big{|}\Big{|}^{2}$ (24)

where the square is interpreted as the sum over magnitude-squared of the
elements of the complex vector argument. Ostensibly minimizing $\mathcal{L}$
as a function of the $\Delta t$, $R$, $R_{0}$, $v$, $v_{0}$, $c$, and $c_{0}$
should yield their numerical values. However, many of these quantities only
appear in $\boldsymbol{\Lambda}$ via the ratios $v/v_{0}$, $R_{0}/L_{0}$, and
$c/c_{0}$. As a result, we can only determine their relative values. This
reflects a deeper freedom that we have not yet exploited: the ability to
arbitrarily pick the numerical scale of the coordinate system. For
convenience, we pick $c_{0}=1$, $L_{0}=1$ (i.e. the reference particle
traverses a distance of $2$), and $v_{0}=1$, which implies that $\tau_{0}=1$,
with all units now being dimensionless. The latter two choices fix the time
and length scales of the system, so that $R$, $R_{0}$, and $v$ are defined
relative to the reference trajectory and velocity. We use this convention
throughout the remainder of the manuscript.

With these choices, we define the true $c$, $v$, $R$, $R_{0}$, and $\Delta t$
to be solutions to the optimization problem

$\displaystyle\\{c^{*}\\!,v^{*}\\!,R^{*}\\!,R_{0}^{*},\Delta
t^{*}\\}=\mathop{\rm argmin}_{\\{c,v,R,R_{0},\Delta
t\\}}\\!\\!\\!\\!\\!\mathcal{L}(c,v,R,R_{0},\Delta t;\hat{\boldsymbol{\rm
f}}_{m}).$

The objective $\mathcal{L}$ is highly non-linear and may have local minima.
Thus, it is important to identify reasonable initial points as inputs to the
optimization. While these issues are discussed at greater length in Sec. III,
we note that simple properties of the signal such as its height, width at half
max, etc. yield initial conditions that should be sufficiently close to
optimal solutions to ensure convergence of the optimization.

##  III Validation with Experimental Data

To validate the analysis presented in Sec. II, we consider a collection of
1302 events measured in the optofluidic device described in Ref. [11]. The
measurands are $15.3\,\mu$m diameter polystyrene microspheres with dispersed
Dragon Green fluorophore.999Certain commercial equipment, instruments,
software, or materials are identified in this paper in order to specify the
experimental procedure adequately. Such identification is not intended to
imply recommendation or endorsement by the National Institute of Standards and
Technology, nor is it intended to imply that the materials or equipment
identified are necessarily the best available for the purpose. We use a high-
speed camera system to manually verify that all events correspond to singlets.
The fluorescence signals, denoted $\boldsymbol{\rm f}_{m}$, are discretized on
a 16-bit data acquisition card that records samples at 2 MHz. Thus, the
$\boldsymbol{\rm f}_{m}$ are vectors whose entries correspond to the discrete
time interval over which the data is collected. The digitizer outputs
fluorescence values in units of volts; we divide these measurements by $1$ V
to non-dimensionlize.

As a preprocessing step, we use a moving average filter to smooth peaks
followed by a polynomial fit peak finder to estimate the height $h_{m}$, width
$w_{m}$, and the time interval of the location of the peak maxima for each
event; see Ref. [11, 17, 18]. The first two quantities are used later in the
optimization. Using the third, we perform a cyclic permutation of the data
(corresponding to temporal shifts) to approximately align the peak maximum
with the midpoint of the time interval, which we rescale to be $[-1,1]$.

Having estimated these quantities, we return to the original (non-smoothed)
signal for subsequent analysis. Because the signals contain noise (likely from
the photodetector and/or counting statistics of photons), we perform a
discrete Fourier transform (DFT) to identify a noise-floor, which begins
around the tenth mode, corresponding to $M=9$ according to our indexing
convention; see also Ref. [19]. We take this value to be a cutoff for a low
pass filter and also keep the complex conjugate weights, which correspond to
the last 9 modes (the case $k=0$ is its own conjugate mode).101010The exact
cutoff will depend on application and the noise-floor of the signal.  Given
that the DFT yields a spectral representation that is a continuous
interpolation of the data, we can directly use the corresponding mode weights
and frequencies in calculations involving Eq. (24).

To test Eq. (24), we pick a reference event $\hat{\boldsymbol{\rm f}}_{0}$ at
random from among those having a height and width roughly equal to the median
value for this population. The specific choice of this event is unimportant,
since the goal is to demonstrate that all signals can be collapsed to any
chosen reference. Next, we choose $N_{s}$ additional samples
$\hat{\boldsymbol{\rm f}}_{m}$ and define the modified objective

$\displaystyle\mathcal{L}_{N_{s}}$
$\displaystyle=\sum_{m=1}^{N_{s}}\mathcal{L}(c_{m},v_{m},R_{m},R_{0},\Delta
t_{m};\hat{\boldsymbol{\rm f}}_{m})$ (25)

where $R_{0}$ is the unknown radius of the reference particle and
$(c_{m},v_{m},R_{m},\Delta t_{m})$ are the unknown transformation parameters
associated with $\hat{\boldsymbol{\rm f}}_{m}$. In Eq. (25), $m$ has been
upgraded from a subscript to an index. Because the beads have on the order of
$10^{9}$ fluorophores per particle, we anticipate that variations in $c_{m}$
(e.g. due to shot noise) are negligible. Thus, we set $c_{m}=1$ for all values
of $m$, so that it does not play a further role in the analysis.

In principle, minimizing $\mathcal{L}_{N_{s}}$ with respect to the unknown
parameters yields estimates of their values, as well as a “consensus” estimate
of $R_{0}$. However, several computational considerations limit our ability to
use Eq. (25) directly. First, the objective function is nonlinear and contains
on the order of $20N_{s}$ terms and $3N_{s}$ parameters. Second, the objective
has rapid oscillations due to the forms of $\boldsymbol{\rm X}(v)$ and
$g(k;R)$. Third, in the limit that $R\to 0$ for $R/R_{0}$ constant, the ratio
$g(k;R_{0})/g(k;R)\to(R_{0}/R)^{3}$, so that $\mathcal{L}_{N_{s}}$ develops a
connected set of solutions. That is, only the ratio of $R_{0}$ to $R$ can be
estimated from $\mathcal{L}_{N_{s}}$.

A resolution to these problems is to consider a regularized version of the
objective

$\displaystyle\mathcal{L}_{r}=\mathcal{L}_{N_{s}}+\epsilon(R_{0}-\bar{R}_{0})^{2},$
(26)

where $\epsilon\ll 1$ is a regularization parameter and $\bar{R}_{0}$ is an
estimate of the particle radius given by outside sources of information. For
the beads under consideration, manufacturer specifications indicate that the
average radius is approximately $7.625$ $\mu$m, which we use as the value for
$\bar{R}_{0}$ in dimensional units. To convert to dimensionless units, we note
that: (i) the digitizer samples at 2 MHz; (ii) each peak is 1200 samples long;
and (iii) the average velocity as estimated from time-of-flight between two
interrogation regions is 337 mm/s. Thus, the characteristic distance $d_{c}$
traversed by a particle during an event is

$\displaystyle d_{c}=\frac{1200\,{\rm samples}}{2{\rm MHz}}337{\rm
mm/s}\approx 202\,\mu{\rm m}$ (27)

The characteristic radius 7.625 $\mu$m normalized by $d_{c}$ and doubled (to
account for mapping to $-1\leq z\leq 1$) yields $\bar{R}_{0}=0.0755$. We also
set $\epsilon=10^{-2}$ and fix $N_{s}=5$.

Figure 2: Example of data collapse for 1302 time-series. Top: The original
curves to which we apply the scale transformations described in the main text.
The bold-yellow time-trace is used as the reference. The time and fluorescence
values are given in their original units, although subsequent analysis is done
after non-dimensionalizing. Bottom: All 1302 curves after data collapse and
conversion back to the original units. To achieve, collapse, we use Eq. (22)
to turn each measured curve into a realization of $f_{0}(t)$. For the purposes
of defining units in this and other figures, we adopt this active
interpretation of the transformations. The inset shows the point-wise
residuals between the transformed curves and reference. The difference is
normalized by the maximum value of the reference. Note that the residuals have
the characteristic frequency of the mode cutoff. Figure 3: Absolute values of
the mode-weights associated with the 1302 curves in Fig. 2 after
transformation and collapse via Eq. (22). The inset shows the signal-to-noise
ratio, which was computed by: (i) using the first 10 modes to determine the
transformation parameters; (ii) applying the transformation matrices to the
first 100 modes of each signal; and (iii) computing the ratio of the absolute
value of the mean to standard deviation of the resulting mode-weights. The
dotted-red vertical line shows the mode cutoff, which occurs when the mode-
weights approach a signal-to-noise ratio of unity.

To estimate transformation parameters associated with the remaining curves, we
minimize $\mathcal{L}$ for each of the remaining 1296 curves separately, using
the consensus value of $R_{0}$. The results of this exercise are shown in Fig.
2. The top subplot shows the original time-traces, while the bottom subplot
shows the data collapse; the inset shows the relative errors. See also Fig. 3.
We also compute the sample mean $\bar{R}$ and sample variance $\sigma_{R}^{2}$
estimates of the particle radius to estimate the coefficient of variation (CV)
given by $CV_{R}=\sigma_{R}/\bar{R}$ [20]. We find $CV_{R}=3.18$ %, which is
consistent with the manufacturer specified range of particle sizes ($4.4$%)
and the number of measurements considered.

Figure 4: Histograms of relative particle sizes (top) and velocities (bottom).
The coefficient of variation (CV) in particle radii is 3.18%, which matches
well with the manufacturer specification of 4.4% (Bangs Laboratories,
FSDG009).

##  IV Applications

We now consider two applications of the analysis discussed in Sec. III:
estimating per-event measurement reproducibility and doublet deconvolution. A
key theme of these examples is that UQ enables one to extract otherwise
inaccessible information from cytometry measurements. In particular, UQ allows
us to resolve a uniqueness problem of identifying the most probable singlets
that generate a multiplet. Our main task is to recast the previous results in
a probabilistic framework and demonstrate how this leads to new kinds of
signals analyses.

###  IV.1 Uncertainty Quantification

Equation (22) implies that the optimal transformation parameters undo signal
deformation due to effects such as particle speed, size, brightness, etc.
Thus, each transformed signal $\boldsymbol{\Lambda}\hat{\boldsymbol{\rm
f}}_{m}$ can be viewed as a realization of a measurement of identical
particles, i.e. a version of $\hat{\boldsymbol{\rm f}}_{0}$. Treating these as
independent and identically distributed random variables, we can estimate the
per-event reproducibility in terms of corresponding statistical estimators.
For example, if $A_{j}$ denotes the integrated area of the $j$th transformed
signal, then reproducibility in area measurements is given in terms of a
standard deviation $\sigma_{A}$ computed via

$\displaystyle\sigma_{A}^{2}$
$\displaystyle=\frac{1}{N_{e}-1}\sum_{j=1}^{N_{e}}(A_{j}-\bar{A})^{2}$ (28)
$\displaystyle\bar{A}$ $\displaystyle=\frac{1}{N_{e}}\sum_{j=1}^{N_{e}}A_{j}$
(29)

where $N_{e}$ is the number of events.111111In Sec. II we use the subscript
$m$ to denote a measured time-series, which we implicitly index from $1\leq
m\leq N_{e}$. Unless otherwise specified, in this section we use the subscript
$j$ when referring to transformed signals. Note that the number of transformed
signals is still $N_{e}$. Figure 5 shows the results of this analysis for the
1302 curves considered in Fig. 2. The inset to the bottom plot of Fig. 2
indicates that the point-wise (in time and relative to the reference
amplitude) reproducibility in shape of the filtered signals is on the order of
2 % or less.

We make this last observation more precise by considering the Fourier modes of
each realization of $\hat{\boldsymbol{\rm f}}_{0}$, which we denote by
$\hat{\boldsymbol{\rm f}}_{0,j}$. As a practical matter, it is easier to
consider the real and imaginary parts of $\hat{f}_{0,j}(k;\tau_{0})$
separately, since their fluctuations should be independent. Dropping explicit
reference to $\tau_{0}$, we decompose $\hat{f}_{0,j}(k)=\hat{f}_{j,{\rm
Re}}(k)+i\hat{f}_{j,{\rm Im}}(k)$, where

$\displaystyle{\rm Re}[\hat{f}_{0,j}(k)]=\hat{f}_{j,{\rm Re}}(k)=\bar{f}_{\rm
Re}(k)+\Delta\hat{f}_{j,{\rm Re}}(k)$ (30a) $\displaystyle{\rm
Im}[\hat{f}_{0,j}(k)]=\hat{f}_{j,{\rm Im}}(k)=\bar{f}_{\rm
Im}(k)+\Delta\hat{f}_{j,{\rm Im}}(k),$ (30b)

$\bar{f}_{{\rm Re}}(k)$ [$\bar{f}_{{\rm Im}}(k)$] are average real [imaginary]
Fourier modes, and $\Delta\hat{f}_{j,{\rm Re}}(k)$ $[\Delta\hat{f}_{j,{\rm
Im}}(k)]$ are random variables associated with uncertainty in each mode. The
average mode weights are constructed by analogy to Eq. (29). As a bookkeeping
step, it is useful to define a vector $\hat{\mathfrak{f}}_{j}$ whose first
$M+1$ elements (where $M$ is the mode cutoff) are the real parts of modes
$k=0,\pi,2\pi,...,M\pi$ of $\hat{\boldsymbol{\rm f}}_{0,j}$ and whose next $M$
elements are the corresponding imaginary parts for $k\neq 0$; see Appendix A
for more detailed motivation of this decomposition.121212Note that
$\hat{\mathfrak{f}}$ contains all of the same information as
$\hat{\boldsymbol{\rm f}}$ because $f(t)$ is real. There are $N_{e}$ such
realizations of $\hat{\mathfrak{f}}_{j}$, where $1\leq j\leq N_{e}$. The
corresponding average (sample mean) weights and perturbations are
$\bar{\mathfrak{f}}$ and $\Delta\hat{\mathfrak{f}}_{j}$. We also use the
$\hat{\mathfrak{f}}_{j}$ to construct a covariance matrix $\Xi$, as well as
the probability of a deviation $\Delta\hat{\mathfrak{f}}$

$\displaystyle
P_{1}(\Delta\hat{\mathfrak{f}})\propto\exp\left[-\frac{1}{2}\Delta\hat{\mathfrak{f}}^{\rm\,T}\,\Xi^{-1}\Delta\hat{\mathfrak{f}}\right],$
(31)

where the proportionality factor depends on the dimensionality and determinant
of $\Xi$. Equation (31) characterizes uncertainty in the shape of a
signal.131313The real basis underlying $\hat{\mathfrak{f}}$ in Eq. (31) allows
us to express correlations between modes in terms of the usual covariance
matrix. Likewise, we construct a probability density

$\displaystyle
Q_{1}(\Delta\chi)\propto\exp\left[-\frac{1}{2}\Delta\chi^{\rm\,T}\Upsilon^{-1}\Delta\chi\right]$
(32)

where $\Delta\chi^{\rm\,T}=(R-\bar{R},v-\bar{v})$ characterizes deviation from
the average radius and velocity, and $\Upsilon$ is the corresponding
covariance matrix. In the next section, we show how these probability
densities play a fundamental role in the tasks of doublet deconvolution.

Figure 5: Histogram of integrated areas of the transformed signals in Fig. 2.

###  IV.2 Doublet Deconvolution

As an illustration of how UQ can inform downstream analysis, consider the task
of doublet deconvolution. In typical cytometry protocols, such data is
identified indirectly via gating procedures and subsequently rejected. We
propose an alternate strategy based on constrained optimization.

First assume that a doublet $d(t)$ is a linear combination singlet signals. In
Fourier space this implies

$\displaystyle\hat{\mathfrak{d}}$
$\displaystyle=\boldsymbol{\Theta}^{-1}(c_{1},v_{1},R_{1},R_{0},\Delta
t_{1})[\bar{\mathfrak{f}}+\Delta\hat{\mathfrak{f}}_{1}]$
$\displaystyle\qquad+\boldsymbol{\Theta}^{-1}(c_{2},v_{2},R_{2},R_{0},\Delta
t_{2})[\bar{\mathfrak{f}}+\Delta\hat{\mathfrak{f}}_{2}],$ (33)

where the parameters $c_{j}$, $v_{j}$, $R_{j}$, and $\Delta t_{j}$ transform
the reference signal into the corresponding singlets that comprise the
doublet, and $\hat{\mathfrak{d}}$ and $\boldsymbol{\Theta}$ are the
representations of $\hat{\boldsymbol{\rm d}}$ and $\boldsymbol{\Lambda}$ in
the same basis as $\Delta\hat{\mathfrak{f}}$ (cf. the Appendix). We use the
inverse $\boldsymbol{\Theta}^{-1}$, since we are transforming from the
reference signal to the measured signal. Generalizing Eq. (33) to an arbitrary
multiplet $m(t)$, yields

$\displaystyle\hat{\mathfrak{m}}=\sum_{j}\boldsymbol{\Theta}^{-1}(c_{j},v_{j},R_{j},R_{0},\Delta
t_{j})[\bar{\mathfrak{f}}+\Delta\hat{\mathfrak{f}}_{j}].$ (34)

Assuming that a sample is known to be a multiplet comprised of $\mathcal{M}$
singlets, the probability that the singlets have a given set of deformations
$\Delta\mathfrak{f}_{1},\Delta\mathfrak{f}_{2},...,\Delta\mathfrak{f}_{\mathcal{M}}$
before transformation to $\hat{\mathfrak{m}}$ can be expressed as

$\displaystyle
P_{\mathcal{M}}(\Delta\hat{\mathfrak{f}}_{1},\Delta\hat{\mathfrak{f}}_{2},...,\Delta\hat{\mathfrak{f}}_{\mathcal{M}})=\prod_{j=1}^{\mathcal{M}}P_{1}(\Delta\hat{\mathfrak{f}}_{j}).$
(35)

Likewise, the probability of a set of transformation parameters is given by

$\displaystyle
Q_{\mathcal{M}}(\Delta\chi_{1},\Delta\chi_{2},...,\Delta\chi_{\mathcal{M}})=\prod_{j=1}^{\mathcal{M}}Q_{1}(\Delta\chi_{j})$
(36)

We assume that the true values of $\Delta\hat{\mathfrak{f}}_{j}$ and
$\Delta\chi_{j}$ are those that maximize the joint probability
$P_{\mathcal{M}}Q_{\mathcal{M}}$. However, the singlets must reconstruct the
multiplet exactly. Together this suggests the objective

$\displaystyle\mathcal{L}_{\mathcal{M}}=\sum_{j}\Delta\hat{\mathfrak{f}}_{j}^{\rm\,T}{\boldsymbol{\Xi}}^{-1}\Delta\hat{\mathfrak{f}}_{j}+\Delta\chi_{j}^{\rm\boldsymbol{\rm
T}}\Upsilon^{-1}\Delta\chi_{j}$ (37)

which we minimize with respect to the $\Delta\hat{\mathfrak{f}}_{j}$ and
$\Delta\chi_{j}$, subject to the constraint given by Eq. (34).

This optimization problem is performed over $\mathcal{M}(2M+1)+4\mathcal{M}$
parameters, where $\mathcal{M}$ is the number of singlets comprising the
signal, and $M$ is the mode cutoff. In practice, however, this can be reduced
to a $4\mathcal{M}$ dimensional problem, since the constraint is linear in the
$\Delta\hat{\mathfrak{f}}_{j}$. We leave details of the optimization to the
Appendix. We also set $c_{j}=1$, consistent with Sec. III.

Figure 6: Example of doublet deconvolution. Signals associated with two known
singlets were added together to make a synthetic doublet. Original signals are
shown in dotted lines. The solid lines show the recovered signals after
deconvolution (purple and yellow) as well as the filtered doublet (red). The
inset shows the error in the recovered singlets relative to the original
filtered singlets, normalized by the maximum values of the latter. Figure 7:
Impact of peak separation on doublet deconvolution. The inset shows synthetic
doublets constructed by varying the peak separation distance of the singlets
in Fig. 6. The main figure shows the time-dependent relative difference
between one of the recovered singlets and its true (filtered) signal. Colors
have the same interpretation in both the main plot and inset.  Figure 8:
Deconvolution of two visually verified doublets. The inset shows false color
images of the particles traversing the laser interrogation region. In the
inset, the horizontal axis is parallel to the flow direction, and the vertical
direction is parallel to the laser. The horizontal elongation is due to
blurring. In both plots, two lobes are visible. The main figure shows the
original signal (blue), as well as the reconstructed singlets (orange and
yellow). The sum of singlets is superimposed in purple over the blue curve but
is visually indistinguishable. In the bottom plot, one particle is likely
rotated out of the focal plane, which would account for the differences in
visual intensities relative to signal peaks. This assumes the geometric
factors associated with illumination and light collection are constant for
both particles; see Ref. [11] for the validity of this assumption.

Figure 6 illustrates the results of this analysis applied to a synthetic
doublet constructed from two singlets taken from Fig. 2. The analysis recovers
the original signals to within roughly 5% accuracy or better point-wise in
time. The actual transformation parameters associated with each particle are
recovered to within approximately 2%. Figure 7 shows the relative errors (with
doublets inset) associated with different peak spacings between the singlets.
All reconstructions have point-wise errors of less than 10%, even for closely
spaced singlets.

Figure 8 shows the results of this analysis applied to two measured signals
that were visually verified to be doublets. In the second event (bottom plot),
one of the beads was likely out-of-focus of the camera, leading to a reduced
visual signal. See Ref. [11] for more details. In both cases the analysis
recovers singlets that are consistent with the synthetic examples shown in
Fig. 6. While more work is needed to make this analysis fully robust for
commercial settings, the examples provided herein indicate that the
constrained optimization formulation is a useful tool for doublet
deconvolution.

##  V Discussion

###  V.1 The Role of Modeling in Cytometry

From a metrology perspective, the operation of a cytometer is at odds with the
fundamental assumptions used to characterize reproducibility and uncertainty.
The inherent separation of scales – thousands of distinct, micron-sized cells
traveling meters per second in a complicated fluid-dynamic system – coupled
with the native variability of biological systems means that it is challenging
to repeat an independent measurement on the same particle in the same optical
region. Moreover, as Fig. 2 demonstrates, even small variations in reference
materials can lead to dramatic changes in raw measurement signals. Thus, a key
challenge in characterizing the accuracy of a cytometer arises from an
inability to design operating conditions that isolate individual sources of
uncertainty. In other words, experimental analysis of uncertainty is
difficult.

The modeling herein provides a distinct approach to this problem through
synthesis, i.e. building up a prediction of the experimental result, taking
into account the cumulative effects of multiple sources of variation. Under
ideal circumstances, comparing this synthetic result with reality leads to a
unique quantification of the physical phenomena (e.g. particle size, speed,
etc.) that generate each signal. Any remaining variation that we cannot
account for is then treated as a reasonable proxy for reproducibility.

This distinction between analysis and synthesis highlights both the potential
roles and challenges of using mathematical modeling for device
characterization. In particular, the model permits us to infer the relative
magnitudes of coupled physical effects, thereby compensating for experimental
limitations. But critically, the accuracy of these estimates relies on the
validity of the underlying theoretical assumptions. For example, the
assumption that the particles are spherical may not be appropriate for
deformable cells, requiring revision of Eq. (7). Likewise, the properties of
laser profiles considered herein may not apply to all cytometers.

These observations suggest a need for deeper coordination between UQ, signals
analysis, and design of cytometers. While models can often be revised to
account for increasingly complex phenomena (e.g. deformable cells, non-uniform
lasers), computations invariably become too expensive to be useful. In such
cases, theory can instead inform design changes that may be experimentally
achievable and lead to improved accuracy through consistency with modeling
assumptions. See also Ref. [11] for related ideas.

###  V.2 Further Applications and Open Directions

The analysis presented herein offers routes to solving several outstanding
problems in signals analysis for cytometry. In particular, the ability to
quantify reproducibility of a measurement suggests a task wherein one seeks to
minimize this uncertainty as a function of operating conditions, e.g. flow
velocity and degree of focusing, particle density, etc. Moreover,
reproducibility estimates suggest the possibility for propagating uncertainty
into populations studies so as to inform best gating and classification
strategies. In this spirit, the work presented in Ref. [21] may be relevant.

Recent publications have also suggested the importance of cell deformability
during the measurement process. Equation (4) makes a key simplifying
assumption that requires modification in this case. To the extent that it is
possible to parameterize deformation modes (e.g. in terms of spherical
harmonics or empirically defined shapes), the transformations associated with
Eq. (6) can be generalized to account for non-spherical particles. This
ability to extract shape information from a relatively simple time-series
could yield significant improvements in throughput relative to imaging
cytometers while still characterizing more nuanced information about the cell
status.

The problem of doublet identification is also largely unresolved. Traditional
strategies address this task through subjective gating of populations.
However, the distribution of Fourier spectra may provide more objective,
probability-based methods. The left plot of Fig. 9 shows the envelope of
spectra (red region) associated with the distribution of transformed singlets
in Fig. 2. The right plot shows a collection of synthetic doublets; their
corresponding spectra (normalized to the same area as the mean singlet) are
shown relative to the singlet envelope on the left. Note that for almost any
peak separation at least one mode associated with the doublet falls outside
the admissible window for singlets. This suggests that signal shape may be an
exquisitely sensitive tool for identifying events that are candidates for our
doublet deconvolution algorithm. Such questions, however, are left for future
work.

Figure 9: Left: Fourier spectra of synthetic doublets (dotted) compared with
the uncertainty window for singlets (solid red). For a fixed mode number, the
uncertainty window for singlets is the set of all mode weights between the
minimum and maximum values computed in Fig. 2. Synthetic doublets are composed
of two copies of a singlet (taken at random from Fig. 2) whose peaks are
separated by 0 (green), 10 (blue), 20 (cyan), 30 (yellow), 40 (magenta), and
50 (black) microseconds. The doublet mode weights are divided by 2 to
normalize the them to scale of the singlets. Note that after a separation of
only 10 microseconds the doublet spectra fall outside the uncertainty window
of the singlets. Right: The doublets whose spectra are shown on the left plot.
Solid lines are the doublets, whereas dotted lines of the same color are the
corresponding singlets. The colors have the same interpretation as on the left
plot. Note that all doublets appear as a single peak and are visually
difficult to distinguish from singlets.

###  V.3 Metrics in the Context of Past Work

Signals analysis has been integral to cytometry since its inception [7, 8, 9,
10, 19]. However, virtually all techniques characterize cells in terms of
properties such as the signal area, height, width, etc., and the
justifications for such analyses are not universally valid. For example,
common practice dictates that forward-scatter (FSC) vs area or FSC-height
measurements are appropriate for detecting doublets, as cells should nominally
follow one another single file [12]. But as Fig. 8 illustrates, cells may pass
the laser interrogation region side-by-side or even at a diagonal to the flow
direction. A recent study of inertial effects also suggests that typical flow
focusing strategies are ineffective [11]. Furthermore, biological processes of
interest may interfere with measurements typically used to distinguish
doublets [22].

Many of these problems arise from the fact that quantities such as signal
area, width, and height are functionals of high-dimensional data (i.e. the
full event time-trace) that produce a scalar. Thus, much of the underlying
information that could be used to characterize measurands is lost before
signals are actually compared. This suggests a need to revisit the order of
operations, e.g. by directly comparing full signals before quantifying their
differences in terms of simple descriptors.

To better understand this point, it is useful to consider the concept of a
metric, which addresses the question: how “far apart” are the generic objects
$h$ and $g$? In other words, the metric, often denoted $d(h,g)$, defines a
notion of distance appropriate to the structure of $h$ and $g$. While a
complete treatment is beyond the scope of this work (see Ref. [23]), we note
two fundamental properties: (i) $d(h,g)\geq 0$ for any objects $h$ and $g$,
i.e. all distances are non-negative; and (ii) $d(h,g)=0$ implies $h=g$.

While seemingly abstract, these observations have important ramifications for
UQ of cytometry. Conventional analyses implicitly use the absolute value
metric

$\displaystyle
d(\mathcal{P}[f_{j}(t)],\mathcal{P}[f_{k}(t)])=|\mathcal{P}[f_{j}(t)]-\mathcal{P}[f_{k}(t)]|$
(38)

and related notions of Euclidean distance as the basis for comparing
measurements, where $\mathcal{P}[f_{j}(t)]$ is a functional that returns a
scalar value (e.g. area, height) associated with the time-series $f_{j}(t)$.
Since the dimensionality of $\mathcal{P}[f_{j}(t)]$ is low, the statement
$d(h,g)=0$ applied to Eq. (38) implies that two distinct time-series may still
be treated as equivalent. For example, a large, dim cell may yield the same
integrated area as a small bright one, despite the objects being very
different. Yet according to Eq. (38) both objects would be the same based on
intensity alone. Moreover, such approaches make it impossible to separate
effects (e.g. flow rate) that may change the signal shape while keeping scalar
properties such as the height constant. Such shortcomings have limited UQ
studies to those effects due to photodetectors and total uncertainties as
characterized by population histograms [24, 25, 26, 27].

The objectives considered herein address such problems by defining the metric
in terms of the time-series explicitly; that is, Eq. (24) considers the
situation $d(\hat{\boldsymbol{\rm f}}_{j},\hat{\boldsymbol{\rm
f}}_{k})=||\hat{\boldsymbol{\rm f}}_{j}-\hat{\boldsymbol{\rm f}}_{k}||$ in
which the notion of distance $||\cdot||$ is applied directly to the full time-
series (expressed as its Fourier transform). Equation (24) in particular is a
variation on the commonly used $\ell_{2}$ or “sum-of-squares” metric [23]. A
benefit of this approach is that it leverages the full information content of
each signal to characterize multiple sources of variation. However, the
resulting analysis is more computationally expensive and does not have a
simple interpretation. Nonetheless, the example of doublet deconvolution
highlights the usefulness of such techniques, and we speculate that more
advanced signals analyses in cytometry will require further development of
appropriate metrics.

###  V.4 Limitations

The modeling framework described Sec. II.1 sets forth the minimum assumptions
required for the validity of our analysis. A key goal of the theory is to
avoid the need for detailed device characterization, which can be costly.
However, our main assumptions on the laser profile and geometric factors may
not be applicable to all cytometers. In such cases, using our analysis would
decrease confidence in results by introducing additional model-form
uncertainty.

In principle, such limitations can be overcome by more detailed modeling of
the experimental system. The general structure of Eq. (24) can be maintained,
although the transformation matrix $\boldsymbol{\Lambda}$ may not exist in
closed form. Rather, it may be necessary to evaluate it on-the-fly in terms of
a simulation or other computational model of the measurement process. In this
case, a key challenge will be to formulate an optimization routine that can
minimize Eq. (24) in a reasonable amount of time. Reduced-order modeling or
computationally inexpensive approximations are possible routes for addressing
such problems.

##  Appendix A Optimization for Multiplet Deconvolution

In this appendix we describe the mathematical formulation and solution of
optimization for multiplet deconvolution. For convenience, we restate the key
equations underlying this problem. Specifically, the objective is given by

$\displaystyle\mathcal{L}_{d}$
$\displaystyle=\sum_{j=1}^{\mathcal{M}}\Delta\hat{\mathfrak{f}}_{j}^{\rm
T}\,\,\Xi^{-1}\,\Delta\hat{\mathfrak{f}}_{j}+\Delta\chi_{j}^{\rm\boldsymbol{\rm
T}}\Upsilon^{-1}\Delta\chi_{j}$ (39)

whereas the constraint in the complex basis is

$\displaystyle\hat{\boldsymbol{\rm m}}$
$\displaystyle=\sum_{j}\boldsymbol{\Lambda}^{-1}(c_{j},v_{j},R_{j},\Delta
t_{j})[\bar{\boldsymbol{\rm f}}+\Delta\hat{\boldsymbol{\rm f}}_{j}],$ (40)

where $\bar{\boldsymbol{\rm f}}$ is the complex representation of
$\bar{\mathfrak{f}}$. Our goal is to convert Eq. (40) to the same basis as
$\hat{\mathfrak{f}}$ and then simplify the optimization problem implied by Eq.
(39).

Define
$\boldsymbol{\Lambda}_{j}^{-1}=\boldsymbol{\Lambda}^{-1}(c_{j},v_{j},R_{j},\Delta
t_{j})$. Next, note that a vector $\hat{\boldsymbol{\rm f}}$ has $2M+1$
elements. Because we compute them via a DFT, we take the convention that the
first $M+1$ modes correspond to $k=0,\pi,2\pi,...,M\pi$, whereas the last $M$
modes correspond to $k=-M\pi,-(M-1)\pi,...,-\pi$. Moreover, because the
signals are real in the time-domain, we know that
$\hat{f}(k)=\hat{f}(-k)^{\star}$, where $\star$ denotes the complex conjugate.
This implies that all of the relevant information about $\hat{\boldsymbol{\rm
f}}$ is contained in $\hat{\mathfrak{f}}$, since only $M$ complex Fourier
modes are needed to describe the signal. To see this explicitly, express an
arbitrary $\hat{\boldsymbol{\rm f}}$ as

$\displaystyle\hat{\boldsymbol{\rm f}}$
$\displaystyle=\begin{pmatrix}{\mathfrak{a}}_{1}\\\
{\mathfrak{a}}_{2}+i{\mathfrak{a}}_{M+2}\\\
{\mathfrak{a}}_{3}+i{\mathfrak{a}}_{M+3}\\\ \vdots\\\
{\mathfrak{a}}_{M+1}+i{\mathfrak{a}}_{2M+1}\\\
{\mathfrak{a}}_{M+1}-i{\mathfrak{a}}_{2M+1}\\\
{\mathfrak{a}}_{M}-i{\mathfrak{a}}_{2M}\\\ \vdots\\\
{\mathfrak{a}}_{2}-i{\mathfrak{a}}_{M+2}\end{pmatrix}$ (41)

where ${\mathfrak{a}}_{j}$ is the $j$th element of $\hat{\mathfrak{f}}$.

To derive the transformed version of $\boldsymbol{\Lambda}$ in the basis of
$\hat{\mathfrak{f}}$, first decompose the matrix into blocks via

$\displaystyle\boldsymbol{\Lambda}$ $\displaystyle=\begin{bmatrix}A&B\\\
C&D\end{bmatrix}$ (42)

where $A$ is $(M+1)\times(M+1)$, $B$ is $(M+1)\times M$, $C$ is
$M\times(M+1)$, and $D$ is $M\times M$. Clearly $A$ couples the first $M+1$
modes of $\hat{\boldsymbol{\rm f}}$ into one another, $B$ couples the
remaining $M$ modes into the first $M+1$, and so forth. However, because the
signal remains real in the time-domain after transformation by
$\boldsymbol{\Lambda}$, knowledge of $A$ and $B$ is sufficient to determine
the transformation matrix in the basis of $\hat{\mathfrak{f}}$. Let
${\mathcal{F}}$ denote that operator that reverse the order of columns in a
matrix, $\mathcal{T}_{c}$ ($\mathcal{T}_{r}$) denote the operator that removes
the first column (row) of a matrix, and $\boldsymbol{0}_{M}$ denote a column
vector with $M$ zeros. Then it is straightforward to show that the operator
$\boldsymbol{\Lambda}$ transforms to

$\displaystyle\boldsymbol{\Theta}\\!$
$\displaystyle=\\!\begin{bmatrix}\tilde{A}&\tilde{B}\\\
\tilde{C}&\tilde{D}\end{bmatrix}$ (43)

where

$\displaystyle\tilde{A}$
$\displaystyle=\Re(A)+\left[\boldsymbol{0}_{M+1},{\mathcal{F}}(\Re(B))\right]$
$\displaystyle\tilde{B}$
$\displaystyle=-\mathcal{T}_{c}(\Im(A))+{\mathcal{F}}(\Im(B))$
$\displaystyle\tilde{C}$
$\displaystyle=\mathcal{T}_{r}(\Im(A))+\left[\boldsymbol{0}_{M},{\mathcal{F}}(\Im(B))\right]$
$\displaystyle\tilde{D}$
$\displaystyle=\mathcal{T}_{r}(\Re(A))-{\mathcal{F}}(\Re(B))$ (44)

and $\Re$ and $\Im$ denote the real and imaginary components. Thus we arrive
at an expression for $\hat{\boldsymbol{\rm m}}$ expressed in the basis of
$\hat{\mathfrak{f}}$, viz,

$\displaystyle\hat{\mathfrak{m}}=\sum_{j}\boldsymbol{\Theta}^{-1}_{j}[\bar{\mathfrak{f}}+\Delta\hat{\mathfrak{f}}_{j}].$
(45)

Minimizing Eq. (39) subject to Eq. (45) may entail optimizing over
$\mathcal{O}(100)$ or more variables corresponding to: (i) the scale
parameters $c$, $v$, $R$, and $\Delta t$; and (ii) the real and imaginary
parts of the mode-weights. The latter comprise the majority of variables,
although they only appear up to second order. In contrast, the transformation
variables, while few in number, appear in highly non-linear function
associated with the matrix $\boldsymbol{\Theta}^{-1}$ in Eq. (45). Further
compounding these issues is the fact that both $\boldsymbol{\Theta}^{-1}$ and
$\Xi$ are dense matrices, the latter possibly having eigenvalues close to
zero. This may yield a relatively large numerical problem that is poorly
scaled, and thus challenging to solve.

Fortunately, the constraint given by Eq. (40) is linear in the mode-weights,
which yields a key simplification. Without loss of generality, one finds

$\displaystyle\Delta\mathfrak{f}_{1}=\boldsymbol{\Theta}_{1}\left[\hat{\mathfrak{m}}-\boldsymbol{\Theta}_{1}^{-1}\bar{\mathfrak{f}}-\sum_{j=2}^{\mathcal{M}}\boldsymbol{\Theta}_{j}^{-1}[\bar{\mathfrak{f}}+\Delta\hat{\mathfrak{f}}_{j}]\right].$
(46)

Equation (46) can be substituted into Eq. (37) and minimization performed over
the remaining modes $\Delta\hat{\mathfrak{f}}_{j}$ for $j\geq 2$ for fixed
transformation parameters associated with the $\boldsymbol{\Theta}_{j}$. We
leave this exercise for the reader. For the case of doublets, one find that

$\displaystyle G$
$\displaystyle=\boldsymbol{\Theta}_{1}\left[\hat{\mathfrak{m}}-(\boldsymbol{\Theta}_{1}^{-1}+\boldsymbol{\Theta}_{2}^{-1})\bar{\mathfrak{f}}\right]$
$\displaystyle\Delta\hat{\mathfrak{f}}_{2}^{*}$
$\displaystyle=\Big{[}(\boldsymbol{\Theta}_{1}\boldsymbol{\Theta}_{2}^{-1})^{\rm
T}\Xi^{-1}\boldsymbol{\Theta}_{1}\boldsymbol{\Theta}_{2}^{-1}+\Xi^{-1}\Big{]}^{-1}\boldsymbol{\Theta}_{1}\boldsymbol{\Theta}_{2}^{-2}\Xi^{-1}G$
$\displaystyle\Delta\hat{\mathfrak{f}}_{1}^{*}$
$\displaystyle=\boldsymbol{\Theta}_{1}[G-\boldsymbol{\Theta}_{2}^{-1}\Delta\hat{\mathfrak{f}}_{2}],$

where $\Delta\hat{\mathfrak{f}}_{2}^{*}$ and
$\Delta\hat{\mathfrak{f}}_{1}^{*}$ are the optimal mode perturbations (the $*$
is distinct from the complex conjugate $\star$). Having the
$\Delta\hat{\mathfrak{f}}_{j}^{*}$ in terms of the transformation parameters
(via the $\boldsymbol{\Theta}_{j}$), we may then express the objective as

$\displaystyle\mathcal{L}_{d}=\sum_{j}\left[\Delta\hat{\mathfrak{f}}_{j}^{*}\right]^{\rm
T}\Xi^{-1}\Delta\hat{\mathfrak{f}}_{j}^{*}+\Delta\chi_{j}^{\rm\boldsymbol{\rm
T}}\Upsilon^{-1}\Delta\chi_{j}.$ (47)

This $\mathcal{L}_{d}$ can then be optimized as a function of the scale
transformations parameters.

## References

  * Adan  _et al._ [2017] A. Adan, G. Alizada, Y. Kiraz, Y. Baran, and A. Nalbant, Flow cytometry: basic principles and applications, Critical Reviews in Biotechnology 37, 163 (2017).
  * McKinnon [2018] K. M. McKinnon, Flow cytometry: An overview, Current Protocols in Immunology 120, 5.1.1 (2018).
  * Wang  _et al._ [2021] L. Wang, R. Bhardwaj, H. Mostowski, P. N. Patrone, A. J. Kearsley, J. Watson, L. Lim, J. Pichaandi, O. Ornatsky, D. Majonis, S. R. Bauer, and H. A. Degheidy, Establishing cd19 b-cell reference control materials for comparable and quantitative cytometric expression analysis, PLOS ONE 16, 1 (2021).
  * Laskowski  _et al._ [2020] T. J. Laskowski, A. L. Hazen, R. S. Collazo, and D. Haviland, Rigor and reproducibility of cytometry practices for immuno-oncology: A multifaceted challenge, Cytometry Part A 97, 116 (2020).
  * Tighe  _et al._ [2019] R. M. Tighe, E. F. Redente, Y.-R. Yu, S. Herold, A. I. Sperling, J. L. Curtis, R. Duggan, S. Swaminathan, H. Nakano, W. J. Zacharias, W. J. Janssen, C. M. Freeman, R. R. Brinkman, B. D. Singer, C. V. Jakubzick, and A. V. Misharin, Improving the quality and reproducibility of flow cytometry in the lung. an official american thoracic society workshop report, American Journal of Respiratory Cell and Molecular Biology 61, 150 (2019).
  * Kalina [2020] T. Kalina, Reproducibility of flow cytometry through standardization: Opportunities and challenges, Cytometry Part A 97, 137 (2020).
  * Xiong and Tanik [2015] F. Xiong and M. M. Tanik, Deciphering flow cytometry electronics, in  _SoutheastCon 2015_ (2015) pp. 1–5.
  * Zilmer  _et al._ [1995] N. A. Zilmer, M. Godavarti, J. J. Rodriguez, T. A. Yopp, G. M. Lambert, and D. W. Galbraith, Flow cytometric analysis using digital signal processing, Cytometry 20, 102 (1995).
  * Grimaldi  _et al._ [2004] D. Grimaldi, A. Palumbo, and S. Rapuano, Hardware and software improvements in flow-cytometry measurements, Measurement 36, 111 (2004).
  * Godavarti  _et al._ [1996] M. Godavarti, J. J. Rodriguez, T. A. Yopp, G. M. Lambert, and D. W. Galbraith, Automated particle classification based on digital acquisition and analysis of flow cytometric pulse waveforms, Cytometry 24, 330 (1996).
  * [11] M. DiSalvo, P. N. Patrone, and G. A. Cooksey, Serial flow cytometry in an optofluidic microchip for direct assessment of measurement variations, submitted .
  * Stadinski and Huseby [2020] B. D. Stadinski and E. S. Huseby, How to prevent yourself from seeing double, Cytometry Part A 97, 1102 (2020).
  * Donnenberg and Donnenberg [2015] V. S. Donnenberg and A. D. Donnenberg, Coping with artifact in the analysis of flow cytometric data, Methods 82, 3 (2015), flow Cytometry Methods and Approaches.
  * Ahmed  _et al._ [2019] R. Ahmed, Z. Omidian, A. Giwa, B. Cornwell, N. Majety, D. R. Bell, S. Lee, H. Zhang, A. Michels, S. Desiderio, S. Sadegh-Nasseri, H. Rabb, S. Gritsch, M. L. Suva, P. Cahan, R. Zhou, C. Jie, T. Donner, and A. R. A. Hamad, A public bcr present in a unique dual-receptor-expressing lymphocyte from type 1 diabetes patients encodes a potent t cell autoantigen, Cell 177, 1583 (2019).
  * Burel  _et al._ [2020] J. G. Burel, M. Pomaznoy, C. S. Lindestam Arlehamn, G. Seumois, P. Vijayanand, A. Sette, and B. Peters, The challenge of distinguishing cell–cell complexes from singlet cells in non-imaging flow cytometry and single-cell sorting, Cytometry Part A 97, 1127 (2020).
  * Oppenheim [1999] A. Oppenheim,  _Discrete-Time Signal Processing_ , Pearson education signal processing series (Pearson Education, 1999).
  * Wallace  _et al._ [2004] W. E. Wallace, A. J. Kearsley, and C. M. Guttman, An operator-independent approach to mass spectral peak identification and integration, Analytical Chemistry 76, 2446 (2004).
  * Kearsley  _et al._ [2005] A. J. Kearsley, W. E. Wallace, J. Bernal, and C. M. Guttman, A numerical method for mass spectral data analysis, Applied Mathematics Letters 18, 1412 (2005).
  * Sloot  _et al._ [1987] P. M. A. Sloot, P. Tensen, and C. G. Figdor, Spectral analysis of flow cytometric data: Design of a special-purpose low-pass digital filter, Cytometry 8, 545 (1987).
  * GUM [2008] Jcgm. jcgm 100: Evaluation of measurement data - guide to the expression of uncertainty in measurement, Joint Committee for Guides in Metrology (2008).
  * Patrone and Kearsley [2021] P. N. Patrone and A. J. Kearsley, Classification under uncertainty: Data analysis for diagnostic antibody testing, (2021).
  * Bohmer  _et al._ [2011] R. M. Bohmer, E. Bandala-Sanchez, and L. C. Harrison, Forward light scatter is a simple measure of t-cell activation and proliferation but is not universally suited for doublet discrimination, Cytometry Part A 79A, 646 (2011).
  * Kreyszig [1989] E. Kreyszig,  _Introductory Functional Analysis with Applications_ , Wiley Classics Library (Wiley, 1989).
  * Hoffman and Wood [2007] R. A. Hoffman and J. C. Wood, Characterization of flow cytometer instrument sensitivity, Current Protocols in Cytometry 40, 1.20.1 (2007).
  * Ubezio and Andreoni [1985] P. Ubezio and A. Andreoni, Linearity and noise sources in flow cytometry, Cytometry 6, 109 (1985).
  * Steen [1992] H. Steen, Noise, sensitivity, and resolution of flow cytometers, Cytometry 13, 822—830 (1992).
  * Chase and Hoffman [1998] E. S. Chase and R. A. Hoffman, Resolution of dimly fluorescent particles: A practical measure of fluorescence sensitivity, Cytometry 33, 267 (1998).