1

Accelerated Time-of-Flight Mass Spectrometry

arXiv:1212.4269v2 [math.OC] 28 Jul 2013

Morteza Ibrahimi, Student Member, IEEE, Andrea Montanari, Senior Member, IEEE,
and George Moore, Senior Member, IEEE

Abstract—We study a simple modification to the conventional
time of flight mass spectrometry (TOFMS) where a variable and
(pseudo)-random pulsing rate is used which allows for traces
from different pulses to overlap. This modification requires
little alteration to the currently employed hardware. However,
it requires a reconstruction method to recover the spectrum
from highly aliased traces. We propose and demonstrate an
efficient algorithm that can process massive TOFMS data using
computational resources that can be considered modest with
today’s standards. This approach can be used to improve duty
cycle, speed, and mass resolving power of TOFMS at the same
time. We expect this to extend the applicability of TOFMS to
new domains.
Index Terms—Time of flight mass spectrometry.

I. I NTRODUCTION
Mass spectrometry (MS) refers to a family of techniques
used to analyze the constituent chemical species in a sample.
The applications abound in science and technology and new
fields of scientific investigations have evolved around these
techniques. An example is proteomics which refers to the
science of analyzing peptides and proteins. Proteins are the
workhorse of many biological mechanism. Of great interest to
biological sciences, medical research and drug discovery and
developments is identifying and analyzing the composition and
structure of proteins and other large chemical compounds. It
has become possible only recently to analyze the composition
of proteins with high throughput and accuracy through mass
spectrometry techniques [1] [2]. Other applications include
measuring isotopic ratio, space exploration, testing for illegal
substances etc. Mass spectrometers are usually accompanied
with gas or liquid chromatography techniques and are used in
different configurations in tandem with other mass spectrometer of the same or different types. These configurations provide
a wide range of utility and performance criteria making mass
spectrometry relevant for many different applications.
A typical mass spectrometer consists of three main modules:
an ionizer, a mass analyzer, and a detector. The ionizer converts the species of interest, and possibly other compounds in
the sample to ions in gas phase. Recent advances in ionization
techniques, namely matrix assisted laser desorption ionization
[3], and electrospray ionization [2] has made it possible to
ionize and transform into gas phase large intact molecules
like proteins. These techniques provided new applications for
mass spectrometry and open new avenues for analyzing the
composition and structure of proteins [4].
M. Ibrahimi is with the Department of Electrical Engineering, Stanford
University, Stanford, CA, 94305 USA email: ibrahimi@stanford.edu.
A. Montanari is with the Department of Electrical Engineering and Department of Statistics, Stanford University, Stanford, CA, 94305 USA email:
montanari@stanford.edu.

Fig. 1: Different parts of a TOFMS.

The purpose of the mass analyzer module is to separate the
ions according to their mass to charge ratio (MCR). Today’s
common mass analyzers separate the ions by subjecting them
to electromagnetic fields. These fields exert different forces
to different ions. One class of instruments, broadly referred
to as sector instruments cause the ions with different MCR to
take different trajectories, effectively beamforming a stream of
flying ions of particular MCR toward a detector [5] Another
technique is to have all the ions travel a common trajectory
but with different velocities. This is the basis for time of flight
(TOF) mass spectrometers which we shall describe throughly
in the sequel. A third approach is to guide only a particular
MCR to have a stable trajectory. This is the basis for the
Quadrupole and ion trap mass analyzers [6]. These instruments
can act as a MCR filter or scan a wider range by sweeping
the filter pass band.
The detector module senses the ions by detecting the impact
of charged compounds with the detector surface or the charge
or current they induce by their particular motion.
In this paper we are concerned with time of flight mass
spectrometry (TOFMS) which is a simple yet powerful MS
technique. TOFMS was introduced in the 1940s by Stephens
[7]. TOFMS offers two major benefits over alternative techniques. It has essentially unlimited mass range and high
repetition rate. These properties along with the recent advances
in available hardware and ionization techniques have made
TOFMS an appealing choice for the analysis of samples with
wide mass range [8], biological macromolecules [9] and in
combination with other mass spectrometers [10].
A basic TOF mass spectrometer consists of four parts: an
ionization chamber, an acceleration chamber, a drift region,
and a detector (c.f., Fig 1). The sample is ionized in the
ionization chamber. These ions are then subjected to a very
strong electrical field in the acceleration chamber, effectively
firing them into the drift region. Ideally, the ions entering the
drift region have kinetic energy, K, proportional to their charge
z, i.e., if the potential difference in the acceleration chamber is
U then the following holds K =q
U z. This means that an ion
2U z
with mass m has velocity v =
m . Therefore, assuming

2

that the length of the drift region to be L, the time to reach
the detector is
r
m L
.
(1)
t=
2z 2U
In other words, the time p
that takes for an ion to reach the
detector is proportional to m/z where m is the mass of the
ion and z is its net charge.
TOFMS is a pulsed technique, i.e., ions are formed in an
ionization stage and subsequently accelerated as a packet into
the drift region with ions with different MCR traveling at
different speeds. As the ions impact the detector, they generate
a continuous electrical signal which is then sampled resulting
in a discrete signal. The result of thispprocess, which we
m/z spectrum. A
call a scan, is a noisy sample of the
single scan is often too noisy and this process is repeated
from hundreds to a few thousands
p times and averaged to
m/z spectrum. We will
obtain an accurate estimate
of
the
p
call an estimate of m/z, which is the result of processing
many scans, one acquisition. In many applications multiple
consecutive acquisitions are obtained: to construct a movie
of an evolving sample like an ongoing chemical reaction;
to analyze the output of a preceding chromatography stage;
or to mass analyze samples through an automated system
where TOFMS instrument in being fed automatically, e.g., in
pharmaceutical applications [11].
There are several metrics that describe the performance
of a mass spectrometer. Some of the widely used metrics
are: mass resolving power, mass accuracy, dynamic range,
sensitivity and speed. Mass resolving power is the minimum
difference in mass to charge ratio for two present species to
be distinguishable by the instrument. Mass accuracy is the
normalized precision with which the instrument can report the
MCR of present species, measured in MCR error divided by
MCR. Mass range refers to the range of MCRs the instrument
can detect. Sensitivity is the minimum concentration of an
specie to be detectable by the instruments. Finally, speed is
the number of acquisition the instrument can acquire per unit
time.
These metrics are not independent. Several trade offs exist
among these metrics based on how a TOFMS instrument is
designed and operated. For example, speed can be increased,
at the cost of mass resolving power, mass accuracy and
sensitivity, by decreasing the number of scans collected for
each acquisition. Another trade off exists between speed in
one hand and mass resolving power and accuracy on the other
hand, which is the focus of this paper and is described in
details below.
In conventional TOFMS, the time between consecutive
pulses is set to be long enough to avoid overlap between
different scans, i.e., for the slowest ion in an scan to arrive
at the detector before the fastest ion of the next scan. Hence,
acquisition time is lower bounded by,
∗
Tacquisition ≥ N × Tscan
p
L p
≥N×√
( (m/z)max − (m/z)min ), (2)
2U

where N is the number of scans collected for each acquisition.

Furthermore, the difference in time of arrival for two ions is
p
L p
(3)
t2 − t1 = √
( m2 /z2 − m1 /z1 ).
2U
Hence, increasing the length of the drift region L (i)
increases the acquisition time and therefore decreases the
speed (ii) spreads the ions further apart and therefore increases
the mass accuracy and resoling power. This issue is of fundamental importance because of the following.
Other factors that can improve the mass accuracy and
resolving power of a TOFMS, e.g., detector characteristics and
the speed of the electronics, have reached their limits while
new applications demand even better performance in terms of
mass accuracy and resolving power. One option that remains
available for improving the mass accuracy and resolving power
is to increase the length of the drift region.
However, first, there is the obvious desire for higher speed
and accuracy at the same time. Second, some applications have
stringent requirements in terms of speed, mass accuracy, resolving power and sensitivity. This could be due to exogenous
time restrictions, e.g., when monitoring a chemical reaction
or experimentation choice, e.g., when TOFMS is preceded
by a chromatography stage or used in tandem, with another
mass spectrometry stage. There is also significant economical
implications, a high end instrument costs at the order of
hundreds of thousands of dollars and improving the speed
and throughput while keeping or improving the accuracy can
result in significant savings. This is most clear in the case of
large scale automated experiments used in drug discovery and
development activities.
Therefore, simultaneous improvement of the speed and mass
accuracy and resolving power is of fundamental interest in
TOFMS [12].
Conventional TOFMS works by repeating the same experiment multiple times and averaging the results. The choice of
averaging was mainly due to its simplicity.
In particular, the volume and rate of data generated by
TOFMS instruments prohibited the use of more sophisticated
techniques. Our ability to commit more computational resources has increased significantly since the introduction of
TOFMS. At the same time, the data rate of these instruments
generate has also dramatically increased. In this paper, we
present an efficient, highly parallelizable algorithm that in
conjunction with a simple modification to the conventional
TOFMS can improve mass accuracy, mass resolving power
and speed at the same time.
There has been previous work trying to alleviate this
problem. One approach, called Fourier transform TOF, is
to modulate a continuous ion beam at the source using a
periodic waveform and subsequently accelerate it into the drift
region [13]. The detected signal is then demodulated to obtain
the spectrum. Another approach, called Hadamard transform
TOF (HT-TOF) [14], is based on modulation (gating) of a
continuous ion source by a 0/1 pulse. In this approach, an
ion beam is deflected according to a pseudorandom sequence
of pulses. If the pulse is 1, the beam is undeflected and will
reach the detector. In contrast, if the pulse is 0 the beam is
deflected away from the detector. The pulse sequence has the

3

same frequency as the detector. In an ideal case where there
is neither shot noise nor additive noise, the output can be
described as
y = Hx,
(4)
where y ∈ Rn is the observed signal at the detector and
x ∈ Rn is the TOF spectrum. H ∈ Rn×n is a 0/1 matrix
where each column is the pseudorandom sequence shifted
by the index of the column. As long as H is full rank and
the model is accurate the TOF spectrum can be recovered by
applying the inverse of H to y. The spectrum is obtained by
a deconvolution that can be implemented efficiently using the
fast Hadamard transform.
One drawback of these methods is that they treat the
reconstruction process as a deterministic inversion problem
and ignore the noisy nature of the observations. Furthermore,
they require substantial modification to the hardware of a
conventional TOFMS. In this paper, we describe a different
method, called accelerated time of flight mass spectrometry
(ATOF), which simultaneously achieves mass resolving power,
duty cycle, and speed improvement using essentially the
same hardware as a conventional TOFMS. Our reconstruction
scheme acknowledges the stochastic nature of the observation.
Simulation results using real data confirm the performance
improvement of this scheme.
Notations and Terminology: Let [n] = {1, 2, . . . , n} and
{x[i]}i=[n] be the output of the detector for a single scan.
With a slight abuse of notation we also refer to x as a
scan. Typically a TOFMS experiment consists of many scans
which are later processed (simply averaged) to obtain a more
accurate estimate of the spectrum. Let x(l) [i] be the lth scan.
Define the true spectrum, x̄[i], as the
PNaverage of infinitely
many scans, i.e., x̄[i] = limN →∞ N1 l=1 x(l) [i]. Each scan
x(l) [i] is therefore a noisy version of x̄[i]. We define the trace,
y[t], t = 1, . . . T , to be the observed detector response for
multiple, possibly overlapping scan. Given an observed trace
b of x̄.
y, the goal is to find a good estimate x
In what follows we treat the discrete signals as column
vectors. For u and v two vector of the same dimension, let v ∗
denote the transpose of v and hu, vi the scaler product of u
and v.
As a matter of convention we refer to each element of
the vectors that represent the spectrum (x(l) and x̄) as a bin
and to that of the trace as a sample, e.g., y[1] represents the
first sample of the trace. When an ion impacts the detector
it generates a bell-shaped pulse in the output of the detector.
We refer to an observed pulse in the trace as an impact event,
or event for short. Usually the sampling rate of the detector
response is such that an event spans multiple samples.
For pedagogical reasons, we first describe the algorithm as
if each event could occupy only one sample and there was no
time jitter, i.e., all the ions of the same species are associated
with the same bin. We then describe the algorithm without this
assumptions with some minor modifications. All the results
presented in this paper are obtained using real data from a
conventional TOFMS instrument which is used to simulate
the output of an ATOFMS. The algorithm used to obtain these
results is the generalized version of the algorithm.

Fig. 2: Difference between TOFMS and ATOFMS. In
ATOFMS different scans can overlap resulting in shorter
acquisition time for the same number of scans but also a
convoluted observed trace.

II. M EASUREMENT S CHEME AND THE DATA M ODEL
A TOF measurement from a single scan is commonly very
sparse (after removing the additive electrical noise through
preprocessing, c.f. Section IV). Furthermore, a single measurement of the whole spectrum is not expensive and it can
be viewed as being performed in parallel as all ions are flying
in the drift region at the same time. However, the observed
signal from a single scan is too noisy and many repetitions
of the same measurement are necessary to obtain an accurate
estimate of the spectrum. In a conventional P
TOFMS setting,
N
the observed trace can be expressed as y[t] = l=1 x(l) [t−ln]
(l)
th
where x is the detector response to the l scan and x[i] is
understood to be zero for i ≤ 0 or i > n.
We incorporate a simple, yet powerful, modification to this
scheme [15] (c.f. Fig. 2). Conventional TOF (TOF) requires
collection of many scans, each scan collected independently
with no overlap. ATOF idea is to increase the repetition rate
and allow the subsequent scans to overlap.
Define τl , the firing time, to be the starting time of the
lth scan, i.e., the time when the lth ion packet is accelerated
into the drift region. Define ∆τl ≡ τl+1 − τl . In TOF ∆τl =
∆τ ≥ n to avoid overlapping between consecutive scans. We
relax this condition and let ∆τl be a random variable with
E[∆τl ] = αn, for some α < 1. As is the case with the HTTOF we assume that the detector response to overlapping scans
is the superposition of the individual responses,
y[t] =

N
X
l=1

x(l) [t − τl ].

(5)

In this case at each time t, y[t] is the superposition of multiple
overlapping scans. Assume there are a total of N scans and
let 0 = τ1 < τ2 < · · · < τN be the firing times. For a given
τ = (τ1 , τ2 , . . . , τN ), define the matrix A ∈ RT ×n as

1 if ∃ l ∈ [N ] s.t. i = t − τl
A(t, i) =
(6)
0 Otherwise.
The matrix A can be considered as the adjacency matrix of a
bipartite graph (c.f., Fig 3), with rows of A corresponding to
the samples in the trace y and columns of A corresponding
to the bins on the spectrum x. Sample t on the spectrum is
connected to bin i on the spectrum, i.e., Ati = 1, if and only
if for some scan l ∈ [T ] the ions from bin i of the scan x(l)
arrive at time t in the trace y. In what follows, we will refer

4

35

0

i1

i2

n

i3

30
25
ADC count

t

0

T

Fig. 3: Adjacency matrix A and the corresponding bipartite
graph. The signal on the top represents the spectrum and the
bottom signal represents the trace. The trace is the overlapped
concatenation of noisy copies of the spectrum. Nodes are color
coded where blue correspond to the first scan, purple to the
second, and green to the third. Neighbors of sample t on the
trace are those bins on the spectrum who could potentially
contribute to an event at sample t (again color coded).
to the neighbors of sample t as ∂t = {i ∈ [n] | A(t, i) > 0},
and similarly to the neighbors of bin i as ∂i.
y[t] can be considered as a noisy version of linear measurements of x̄, hAt , x̄i, with At the tth row of A as a column
vector. In this notation, the TOF is a special case where each
row of A has only one nonzero element, i.e., measurement
y[t] corresponds to a noisy observation of x̄[i] for some bin
i. The structure of matrix A reveals the difference between
ATOF and TOF.




1 0 0 0
1 0 0 0
0 1 0 0
0 1 0 0




0 0 1 0
0 0 1 0




 0 0 0 1
1 0 0 1




1 0 0 0 
0 1 0 0




0 1 0 0



ATOF = 
 0 0 1 0 , AAT OF = 1 0 1 0 
 0 1 0 1




 0 0 0 1

1 0 1 1



1 0 0 0 
0 1 0 0




0 1 0 0


0 0 1 0
0 0 1 0
0 0 0 1
0 0 0 1
Given the trace y and adjacency matrix A one can attempt
to solve for x̄ using an ordinary least squares, x̂LS =
argminkAx − yk2 or `1 -regularized least squares x̂LASSO =
argminkAx − yk2 + λkxk1 [16]. However, simulation results
demonstrate poor performance for both these methods. The
reason lies in the choice of the objective function. Sum
of square residuals approximates the negative log likelihood
when the measurement noise is additive Gaussian. However,
TOFMS is dominated by shot-noise which is signal dependent
and non-additive. Similar issues arises in applications like
photon-limited imaging where the observations are again shotnoise limited. Regularized maximum likelihood approaches
proved effective in these settings [17]. Here we propose
a stochastic model for the observation y and present an
algorithm that optimizes the `1 -regularized log likelihood.

20
15
10
5
0
1.2907 1.2908 1.2908 1.2909 1.2909 1.291

1.291 1.2911

×105

Fig. 4: A sample event from the trace

Figure 4 shows a sample observed event plotted as the
output of ADC vs bin number. As in this figure, each such
event can span multiple samples. However, to simplify the
presentation of the algorithm we first assume that each ion
impact can occupy only one sample on the trace. The algorithm is extended to the realistic case of events spanning
multiple samples in the next section. All the experimental
results presented in this paper also corresponds to this general
case. However, we have not extended the theoretical results of
our paper to the general case.
Define w ∈ Rn such that w[i] is the average number of
ions that impact the detector at bin i for a single scan. In
each scan a large number of molecules of each species enter
the instrument. However, each molecule has a slight chance
of passing all the stages of the instrument and reaching the
detector. Therefore, it is a natural choice to assume that the
number of ions that impact the detector at time i follows
a Poisson distribution with mean w[i]. Figure 5 shows the
estimated empirical probability mass function (EPMF) for the
number of ions that impact the detector. Note that here we
mentioned estimated EPMF since we do not directly observe
the number of ion impacts. What we observe instead is the
current at the output of the detector which has an arbitrary
scaling. We estimate the number of ion impacts as follows.
We start with ten thousands acquisition of the same sample
and identify a set of rare ions as ions that are observed in
more than 0.1% but less than 1% of acquisitions. These ions
correspond to chemical species with low concentrations and
have small probability of having multiple ion impacts in any
acquisition. We take the median area under the pulse (weight)
for these ions as the estimate for the weight of a single ion
impact. We then normalize the weight of all events by this
estimate and round it to the closest integer. Figure 5 shows
the estimated empirical probability mass function (EPMF) for
the number of ions that impact the detector and its maximum
likelihood Poisson fit.
This result indicates that a Poisson model for y is inadequate. In particular, a Poisson random variable, having its

5

scribing the detector response. Given w and A the probability
density function of y[t] is

102
Estimated EPMF
Empirical PMF and the best (ML) Poisson fit

100

0

P(y[t]|w, A) = e−hAt ,wi−w δ(y[t])
∞
X
+
Ek,µ (y[t])PhAt ,wi+w0 (k),

Best Poisson fit

10−2
10−4

(7)
(8)

k=1

where δ( · ) is the Dirac delta representing a probability mass at
zero and Ek,µ ( · ) and Pλ ( · ) are the Erlang PDF and Poisson
PMF defined as

10−6
10−8
10−10
10−12
10−14
−2

0

2

4

6

10

8

12

Number of impacts

Fig. 5: Estimated empirical probability mass function (EPMF)
for the number of ions that impact the detector and the best
(ML) Poisson fit for this data. See Section V for details of
estimating EPMF. This figure shows that a Poisson random
variable is not adequate for modeling the number of impacts.

−1

10

y k−1 exp(− µy )

,
µk (k − 1)!
exp(−λ)λk
Pλ (k) =
.
k!

Ek,µ (y) =

(9)
(10)

w0 > 0 is a constant accounting for the spurious ion impacts
observed at the detector.
Assuming that each ion impact is observed in only one
sample of the trace, different trace samples are the result of
different ion impacts. Hence, given w and A the observed
responses in different samples can be considered independent,
namely
P(y|w, A) =

Data
Pois+Erlang
Exponential

T
Y

P(y[t]|w, A)

t=1

Density

Let I( · ) denote the indicator function, 1 ∈ RT the vector of
all ones, and define log(0)I(FALSE) = 0. Hence the negative
log-likelihood function takes the form
`(w; y, A) = h1, Awi + N w0
−

10−2

100

200

300

400 500 600
Event weight

700

800

900

Fig. 6: Empirical pdf of the normalized observed event weight
(area under pulse) and a Poisson+Erlang model fitted to this
data.

mean and variance tied together, cannot explain the observed
variation in the tail. We assume that the detector response
is additive for multiple concurrent impacts [14]. Furthermore,
as suggested in [18] we assume each ion impact generates
a cumulative ADC count which itself is an exponentially
distributed random variable. Therefore, conditioning on the
event that k ions impact the detector at a certain time the
cumulative ADC count has an Erlang distribution with the
shape parameter k. Figure 6 shows the empirical probability
density function of normalized event weight and its maximum
likelihood fit of the Poisson+Erlang model. Comparing Figs 6
and 5 reveals that the Poisson+Erlang is a much more satisfactory model for this data.
Let µ be the mean of the exponential random variable de-

T
X
t=1

log

∞
X
y[t]k (hAt , wi + w0 )k

k=1

(k − 1)! k! µk

!
I(y[t] > 0)

(11)

for w ≥ 0, and infinity otherwise. Define F ⊆ [T ] as
F = {t ∈ [T ]| y[t] > 0}. Alternatively, for w ≥ 0 the loglikelihood function can be written as
`(w; y, A) = N kwk1


X
p
2 p
−
log
hAt , wi + w0 I1 ( √
y[t]hAt , wi + w0 ) ,
µ
t∈F

(12)
where I1 is the modified Bessel Function of the first kind, and
we dropped the terms 21 (log(y[t]) − log(µ)) and N w0 which
do not depend on w.
We note that the negative log-likelihood function (12) is
strictly convex in w which is remarkable given the existence
of the hidden variable k.
Lemma II.1. The function `y,A (w) is strictly convex for w ≥
0.
Proof: The proof is immediate using the following theorem due to Findling [19].
Theorem 1. (Findling 95) The function xI(x) is strictly logconcave on {x ∈ R | x > 0}

6

A simple transformation of the negative log-likelihood funce = w/µ, then for
tion (12) is insightful. Let λ = 2N µ, w
e ≥ 0,
w
e w;
e y, A) + λkwk
e 1,
`(w; y, A) = `(

(13)

where
e w;
e y, A) = −
`(

X
t∈F




s
s
X
X
ei I1  y[t]
ei .
w
w
log  y[t]
i∈∂t

i∈∂t

A few remarks are in order. First, note that the scaling of
e is irrelevant for our purpose and only the
the variable w
relative values are important. Hence, this is a single-parameter
representation of the NLL function. This is of great importance
for practical systems where optimal tuning of multiple parameters in different operational scenarios can be complicated
and require additional expertise. Second, the single tuning
parameter appears as the multiplicative factor in front of the
`1 regularization term. Given our intuition about the effect of
`1 regularization [16], [20] the parameter governs the sparsity
e i.e., the number of species that appear in
of the estimate w,
the output. In what follows we free the parameter λ from
our original interpretation of it as the product 2N µ and refer
to it as the regularization parameter. Further, we define the
regularized negative log-likelihood cost function Cλ (w; y, A)
as
e
Cλ (w; y, A) = `(w;
y, A) + λkwk1
(14)
III. A LGORITHM
Given the regularized negative log-likelihood cost function
Cλ (w; y, A) our algorithm attempts to solve the following
optimization problem.
e
minimize `(w;
y, A) + λkwk1 ,
w

s.t. w ≥ 0,

(15)

We use the now standard method of iterative soft thresholding
to solve this convex but non-differentiable optimization problem. For a doubly differentiable function f : D ⊂ Rn → R let
∇f and ∇2 f be the gradient and Hessian of f respectively.
e
Let γ > 0 be such that k∇2 `(w;
y, A)k2 < γ −1 for w ≥ 0.
It is easy to see that γ exists because of the presence of the
chemical noise term w0 . We call the parameter γ the step size
as we use it to scale the steps the algorithm takes in each
iteration. Let w(k) be our estimate of w at step k. Then we
can compute an upper bound for the cost function Cy,A,λ (w)
as follows.
e (k) ; y, A) + w∗ ∇`(w
e (k) ; y, A)
Cλ (w; y, A) ≤ `(w
+ γ −1 kw − w(k) k22 + kwk1 .

(16)

Equation (16) provides an approximation for Cy,A,λ (w) when
w is in a small neighborhood of w(k) . Minimizing the righthand-side of Eq. (16) with respect to w as a surrogate for the
actual cost function results in


e (k) ; y, A) .
w(k+1) = ηθ w(k) − γ∇`(w
(17)
where ηθ ( · ) is the soft thresholding function, ηθ (x) = (|x| −
θ)+ with ( · )+ being the positive part and θ ∝ λ−1 . Note that
this is the one-sided soft thresholding function which differs

from the two-sided soft thresholding function by mapping all
negative values to zero. From Eq. (12), ∇`eλ (w; y, A) can be
calculated as
X
1
∇`e∗λ (w; y, A) = −
2hAt , wi
t∈F

 p

 p
I0 2 y[t]hAt , wi + I2 2 y[t]hAt , wi 
 p

+
At . (18)
p
2 y[t]hAt , wiI1 2 y[t]hAt , wi
Given the gradient of the log-likelihood the algorithm is as
follows.
Algorithm
Input: trace y, firing times τ , and constants (θ0 , θ1 , µ)
b
Output: estimated spectrum x
1: Calculate the adjacency matrix A as in Eq. (6).
2: Set w(0) = 0, θ = θ0 + θ1 .
2: Repeat until stopping criterion is met
θ ← θ0 + k12 θ
1

w(k+1) ← ηθ w(k) − γ∇`eλ (w(k) ; y, A) .
3:
4:

b=0
Set x
For t ∈ F
i∗ = arg max w[i]
i∈∂t

5:

b[i∗ ] = x
b[i∗ ] + y[t]
x
b
Return x.

Step 4 in the algorithm is worth noting. It was mentioned
that each observed event t has a set of possible bins on the
spectrum it can be caused by, ∂t (c.f., Fig. 3). The problem
of estimating the spectrum from the observed ATOF trace can
be thought of as assigning each observed event to one of its
neighbors which is the true cause of the event. This framing
of the problem enables us to terminate the slow first order
optimization method as soon as we are confident about the
likely assignment of an event. In the generalized algorithm
which is concerned with the case of multi-sample events this
technique proves instrumental in decreasing the bias in the
estimated spectrum.
IV. E XPERIMENTAL E VALUATION
In this section we present performance evaluation results for
the ATOF algorithm. We use a commercial TOFMS instrument for data collection and obtain 10, 000 scans for a high
concentration multimode chemical sample using conventional
TOFMS technique. Each scan is 100 µs in length sampled at
25 ps intervals. Therefore, in our notation n = 4 × 105 . The
average of these ten thousand scans is considered the ground
truth.
For evaluation, we use a random subset of 1, 000 scans
and simulate ATOF using these scans as follows. First, we
construct a vector of firing times, τ = (τ0 , τ1 , . . . , τN −1 ) by
setting τ0 = 0 and choosing ∆τi ≡ τi − τi−1 uniformly at
random from the interval [∆τmin , ∆τmax ]. Given τ , the trace
is constructed using unaliased scans as prescribed by Eq. (5).
In addition to the aliasing effect, the trace is corrupted by
noise. Henceforth, we preprocess the trace before applying
the reconstruction algorithm by setting the trace to zero

7

0.1

False negative rate

0.08
0.06
d d

dd

→
←1→↔2
←3→4
0.04
↓
↓
↓
↓
↓↓↓↓
hw
0.02
↓↓ ↓
↓↓↓↓
h0
0
4300 4350 4400 4450 4500 4550 4600 4650 4700
Fig. 7: Preprocessing the data. There are four pulses that
exceed the threshold hw from which three pass the minimum
width condition di ≥ dmin (1, 2, and 4). The markers at the
h0 level indicate the start and end of each marked event.

1
False negative rate

0.9
0.8
0.7
0.6
0.5
0.4
0

5

10

15
20
Iteration

×10−1
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
25
30

False discovery rate

Fig. 8: Illustration of true positive (TP), false positive (FP),
and false negative (FN).

Fig. 9: False negative rate (dashed) and false discovery rate
(solid) vs. iteration.

0.1
0.09
0.08
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0
0.4

0.5

0.7
0.6
0.8
False discovery rate

0.9

1

Fig. 10: False negative rate vs. false discovery rate.

unless it is likely to be the result of an ion impact. As
mentioned before, the detector response to each ion impact
is a bell shaped pulse which spreads across multiple samples.
These pulses are corrupted by electrical noise that can be
modeled as additive noise. However, the electrical noise level
is significantly smaller in magnitude compared to the detector
response to an ion impact and henceforth ion impacts can be
marked with high confidence.
We select the potential ion impacts through the following
procedure. Define three constants h0 , hw , and dmin and label a
pulse as a potential ion impact event if the width of the pulse at
hight hw is greater than dmin . If an event satisfies this criterion
the support interval of the event is determined by thresholding
the signal at level h0 . Figure 7 demonstrates this procedure
through an example. In this figure there are four pulses that
exceed the threshold hw in peak magnitude. From these pulses
d1 , d2 , and d4 satisfy the minimum width condition of di ≥
dmin at height hw . Henceforth, after the preprocessing there
are three valid events with the start and end times marked
at level h0 . We set the trace equal to zero wherever it does
not support a valid event. After the preprocessing step, the
trace can be represented as a list of events whereby each event
describes a single pulse in the trace. Note that observed traces
are outputs of an ADC and have an arbitrary scaling. We keep
this scaling but note that the absolute value of the amplitude
of the trace and the corresponding parameters like hw and h0
are irrelevant for our purposes.
Our procedure to identify valid events also enables us to
define metrics for quantitative evaluation of different techniques. We take the true spectrum x̄ to be the average of
all 10, 000 scans. Let x
b be an estimate of x̄. For some
constants h0 , hw , and dmin define Ē = {ē1 , ē2 , . . . , ēm̄ }, and
Eb = {b
e1 , eb2 , . . . , ebm
b
b } to be the set of events in x̄ and x
respectively, obtained through the procedure described above.
For two events ēi and ebj we say ebj matches ēi if ēi overlap
with at least 50% of the width of ebj . For ēi ∈ Ē we say ēi
is a false negative if none of the events in Eb matches ēi . For
ebj ∈ Eb we say ebj is a true positive if there exist ēi ∈ Ē such

8

35
TOF-10K
ATOF-10K

30

Intensity

25

20

15

10

5

0
1000

500

2000
1500
Mass (Dalton)

3000

2500

Fig. 11: Sample reconstructed spectrum for the acceleration
factor of 10.

1
TOF-10K
0.9

ATOF-10K

0.8
0.7

Intensity

0.6
0.5
0.4
0.3
0.2
0.1
0
820

821

822

823

824
826
825
Mass (Dalton)

827

828

829

830

Fig. 12: Sample reconstructed spectrum for the acceleration
factor of 10. (different scale of Fig 11)

that ebj matches ēi and we say ebj is a false positive if it does
not match any event in Ē. See Fig 8 for an illustration of these
concepts.
Let TP, FP, and FN be the number of true positives, false
positives and false negatives respectively. We consider false
negative rate (FNR = FN/(TP+FN)), true positive rate (TPR =
TP/(TP+FN)), and false discovery rate (FDR = FP/(FP+TP)) as
the metrics of interest. Note that the notion of a true negative
is ill-defined in this problem and hence we cannot use the
false positive rate metric. In particular, observed pulses are of
different width and shapes and they can overlap. Therefore,
given an estimated spectrum the question how many pulses
are not observed? is not well-posed.
Unless otherwise stated, the parameters used to obtain the
results of this section are as follows. µ = 225, γ = 2.5×10−3 ,
(trace)
(spectrum)
θ1 = 2 × 10−2 , θ0 = 5 × 10−4 , hw
= 2, hw
= 0.2,
5
3
n = 4 × 10 , N = 10 , ∆τmin = 0, ∆τmax = 2 × 105 .
n
We also define the acceleration factor as the ratio ≡ E[∆τ
].

For example, E∆τ = 14 n results in an acceleration factor of
4, which means that ATOF is four times faster compared to
conventional TOFMS in terms of the time it takes to collect
the same number of scans. At the same time, each event has on
average 4 different positions on the spectrum it can be assigned
to and the algorithm should be able to infer the correct position
with satisfactory accuracy.
Figure 9 shows the false negative rate and false discovery
rate as a function of the iteration for the ATOF with acceleration factor ten. The algorithm starts with the all-zero spectrum.
Hence, the false negative rate is equal to 1 at iteration 0 and
as the algorithm proceeds the false negative rate decreases,
converging to a final value of 0.47. The false discovery rate on
the other hand increases as the algorithm proceeds settling at
a final value of about 0.085. Inspecting Fig. 9 suggest that the
algorithm converges, in the sense of establishing the existence
of ions, in about 15 iterations.
Note that the large value of FDR is by design. Firstly, by
setting the hw threshold (c.f. Fig 7) very low we are requiring
the algorithm to discover ions with diminutive abundance in
the solutions that are observable in the ten thousand scans
ground truth but very rarely appear in a random one thousand
sample. Secondly, in a typical application of TOFMS declaring
the presence of an ion that does not exist is considered a more
costly mistake than missing an ion that is present. Hence, in
line with these type of applications of TOFMS instruments
we choose to operate in a high FNR and low FDR regime.
Figure 10 shows the same plot on the FNR vs. FDR plane.
This figure shows how the algorithm converges as the number
of iteration increases. Figures 11 and 12 show the ground truth
(solid) and reconstructed (dashed) spectrums at two different
scale. Visual inspection of the graphs indicates substantial
match between the reconstructed spectrum and the ground
truth.
The output of a TOFMS is usually used to generate a list
of peaks, i.e., the list of MCRs that are deemed present by the
instrument. Peak picking is an important task for a TOFMS
instrument that can have significant effect on the instrument
performance and the practice involves both publicly available
methods as well as patents and trade secretes. These estimated
peaks can be a few order of magnitudes more precise than the
width of the pulses at the estimated spectrum.
Another more practically important but less transparent
method for defining TP, FP, and FN is to use the list of
peaks generated from the estimated spectrum. We use the
peak picking software that ships with Agilent Technologies
TOFMS. Since the estimated peaks are real-valued variables
we also need to consider a precision. With a slight abuse of
notation we define two peaks to match if they are p
within ∆m
m1 /z1 −
distance
of
each
other
on
the
MCR
scale,
i.e.,
p
m2 /z2 ≤ ∆m.
Figure 13 shows the TPR vs FDR for different values of ∆m
and number of peaks in the ground truth. Each plot contains
four curves corresponding to TOF with one hundred scans
(TOF-100), ATOF with one thousand scans and acceleration
factor 10 (ATOF-1K), TOF with one thousand scans (TOF1000), and ATOF with ten thousand scans and acceleration
factor 10 (ATOF-10K). Note that (TOF-100) and (ATOF-

9

1.0

0.6

0.4

0.2

●

TOF−1K
TOF−100
ATOF−10K
ATOF−1K

0.8

True Positive rate

True Positive rate

0.8

1.0
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

0.6

0.4

0.2

0.0
0.0

0.2

0.4

0.6

0.8

1.0

●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

0.0
0.0

●

0.2

0.4

False Discovery rate

(a) ∆m = 0.1 Dalton, top 400 peaks.

0.8

1.0

(b) ∆m = 0.01 Dalton, top 400 peaks.
1.0

●

0.8

●

● ●

●●
●
●●

●●

●●
●●●●●
●●

●

●

TOF−1K
TOF−100
ATOF−10K
ATOF−1K

●

0.8

●
●
●
●
●

0.6

●
●
●
●

●●

●●●●●●●
● ● ●●●●

●
●
●
●
●
●
●
●

●

●

●

●

●

●

0.2

●

●

●

●

●

●

●

0.0 ●
0.0

●

●
●●

●

0.4

●

●

●

TOF−1K
TOF−100
ATOF−10K
ATOF−1K

●

0.6

●

0.4

●

True Positive rate

True Positive rate

0.6

False Discovery rate

1.0

0.2

TOF−1K
TOF−100
ATOF−10K
ATOF−1K

●

0.2

0.4

0.6

0.8

False Discovery rate

(c) ∆m = 0.1 Dalton, top 1000 peaks.

1.0

0.0 ●
0.0

0.2

0.4

0.6

0.8

1.0

False Discovery rate

(d) ∆m = 0.01 Dalton, top 1000 peaks.

Fig. 13: TPR vs. FDR for different values of ∆m and number of peaks. The effect of hw can be seen by comparing top and
bottom plots. Lowering hw admits more smaller peaks as valid peaks in the ground truth making it more difficult to detect
them. Comparing left and right plots shows the uncertainity in determinig MCR.

10

1

1.0

0.8

TOF−10K
TOF−1K
ATOF−10K
ATOF−1K

TOF-1K
TOF-250
ATOF-1K
Naı̈ve ATOF-1K

0.9
0.8

ECDF

True positive rate

0.6

0.4

0.2

0.7
0.6
0.5
0.4

0.0
1e−04

width (FWHM) to intensity ratio

1e−02

0.3

0

Fig. 14: Empirical CDF of the width (full width at half
maximum) to intensity ratio for the output of TOF (dashed)
and ATOF (solid). ATOF results in little pulse broadening
compared with TOF with the same number of scans but 10
times longer acquisition time.

1K) have the same acquisition time. Similarly, (TOF-1K)
and (ATOF-10K) have the same acquisition time. For these
plots we run the peack picking algorithm on the ground
truth spectrum and obtain a list of picks. These picks are
sorted based on their amplitude and we choose the k most
significant picks where k ∈ {400, 1000}. We then run the
peak picking algorithm on the reconstructed spectrums and
obtain the list of peaks. Similar to the process for the ground
truth we keep the k most significant peaks for each of these
spectrums. A peak in an estimated spectrum is considered
to match a peak in the ground truth if and only if their
estimated MCRs are within ∆m of each other. For the top
400 peaks and ∆m = 0.1 (Fig. 13a, ATOF-10K achieves a
nearly perfect reconstruction. TOF with 1K scans (TOF-1K)
demonstrates acceptable but significantly inferior performance
compared to ATOF-10K This is while ATOF-10K and TOF-1K
have the same acquisition time. ATOF-1K and TOF-1k have
comparable performance for low TPR but TOF-1k outperforms
ATOF-1K for high TPR. Note that ATOF-1K is ten times faster
than TOF-1K. TOF-100 scans do not have enough information
to achieve any significant accuracy in this regime. Figure 13b
shows the same curves when we decrease ∆m to 0.01, i.e.,
when we adopt a more restrictive definition of two peaks
matching. The overall trend is similar to that of Fig. 13b.
Figs 13c and 13d are similar but for the top one thousand
significant peaks. The TPR and FDR degrade for all the cases
since now we are expecting many more smaller peaks to
be detected. However, the relative performance of different
methods and configurations remain unchanged.
The goal of a TOFMS instrument is to precisely measure
the MCR of the present ions. Hence, for the instrument to
be able to perform accurate peak picking it is important that
ATOF does not significantly distort the shape of the pulses at
the estimated spectrum. Figure 14 shows, for ATOF and TOF
and for different number of scans and acceleration factor ten,

0.1

0.2

0.3
0.4
0.5
False discovery rate

0.6

0.7

0.8

Fig. 15: True positive rate vs false discovery rate

the empirical CDF of the width to intensity (hight) ratio of the
peaks (a measure of peak spikiness). As can be seen ATOF
introduces little broadening in the peaks.
To better understand the achieved accuracy, we compare
ATOF and three other cases in Fig. 15. We run multiple
randomized sample reconstruction and calculate bootstrap
confidence intervals. The figure shows the true positive rate
vs. false discovery rate, as parameterized by θ0 for the
bATOF obtained by the ATOF with one
estimated spectrum x
thousand scans and acceleration factor 4. Each data point is
obtained dividing the ten thousand scans into ten buckets of
one thousand scans each. The trace is constructed using one
thousand scans from one bucket. The average of the other nine
buckets is considered as the ground truth. Error bars indicate
2σ̂ confidence intervals where σ̂ is the standard error estimate.
We Also plot the corresponding curves for three other cases.
The red curve corresponds to the naı̈ve ATOF, obtained by
mapping each event to all possible positions on the spectrum,
X 1
bN [i] =
x
Atj y[t],
(19)
degt
t∈F

where A
Ptjn is the adjacency matrix defined in Eq. (6) and
degt = j=1 Atj is the number of positions on the spectrum
that event t can be mapped to. The naı̈ve ATOFuses the
simplest way of processing the ATOFtrace where one assumes
that each event is equally likely to be caused by an ion from
any of its potential locations on the spectrum.
The green (TOF-1K) and blue (TOF-250) curves correspond
to conventional TOFMS, i.e., the spectrum obtained by averaging the scans when there is no overlapping,
bave =
x

N
1 X (l)
x .
N

(20)

l=1

In TOF-250, the number of scans is chosen such that the time
which takes to perform the TOFMS and ATOF are the same.
TOF-1K corresponds to TOFMS with the same number of
scans as ATOF, which has an acquisition time 4 times longer

11

compared to the ATOF. Equivalently, one can think of TOF1K corresponding to the case where an oracle is available for
ATOF experiment that to each event assigns its true position
on the spectrum.
The naı̈ve ATOFand conventional TOFMS curves are parameterized by the thresholding parameter hw which is used
to identify the events in the estimated spectrum. Similar to θ0 ,
hw enables us to obtain a trade off between true positive rate
and false discovery rate.
This comparison shows that ATOF significantly outperforms
the conventional TOFMS. One the one hand, it allows for
a four-fold speed-up with essentially unchanged accuracy
(comparison with TOF-1K). On the other, it allows a two-fold
increase in true positive rate for FDR = 0.2 if the experiment
duration is unchanged (comparison with TOF-250).
V. E VENTS THAT SPAN MORE THAN ONE SAMPLE
In order to simplify the presentation of the main ideas in the
previous sections we assumed that all ion impacts are confined
to one sample in the observed trace. However, as was evident
from Fig. 4 this assumption is far from being true. In this
section we show how the algorithm presented in Section III
can be extended to the case where each ion impact event can
spread across multiple samples in the trace.
Consider the conditional probability of an observation y[t]
from Eq. (7) which is repeated below for convenience.
0

P(y[a]|w, A) = e−hAa ,wi−w δ(y[a])
∞
X
+
Ek,µ (y[a])PhAa ,wi+w0 (k).

(21)

k=1

ta
X

y[t].

(22)

t=ta

Recall that w[i] + w0 represents the arrival rate of ions at bin
i for a single scan, i.e., w[i] + w0 is the expected number of
ions that arrive at bin i in a single scan. Henceforth, assuming
that ions arrive independent of each other hAt , w + w0 1i is
the cumulative mean of the number of ions that arrive at time
t where 1 is the vector of all ones. Given the assumption that
an ion impact being confined to only one sample we used
hAt , w + w0 1i as the expected number of ions that caused
the observation y[t]. Here, we make the alternative assumption
that for an ion impact observation spanning the interval [ta , ta ]
the expected number of ions responsible, K, is given by
E[K] =

ta
X
t=ta

hAt , w + w0 i

(23)

Let the generalized neighborhood of an event be
∂a = {i ∈ [n] | ∃ t ∈ [ta , ta ] s.t. Ati > 0}.
Then, Eq. 23 can be written as
X
E[K] =
(w[i] + w0 ).
i∈∂a

i∈∂a

)

(

za
+ exp − −
µ
∞
X
k=1

X

0

(w[i] + w )

i∈∂a

X
zak−1
(
w[i] + w0 )k .
µk (k − 1)!k!

(26)

i∈∂a

For a given trace y, let F be the set of observed events.
Similar to section II assume that given w and A distinct
events (ion impacts) are independent. Then the probability of
observing the set of non-zero-weight events F can be written
as
(
)
Y
X
za
P(F |w, A) =
exp − −
(w[i] + w0 )
µ
a∈F
i∈∂a

∞
X
X
zak−1
0 k
(
w[i] + w ) .
µk (k − 1)!k!
k=1

i∈∂a

(27)

Assume instead of being confined to one sample the ion
impact event a is spread from time ta to ta . Define the weight
of ion impact a, za , as
za =

In other words, we are assuming that an observed event can be
caused by one or multiple ions impacting at any time during
the span of the event.
Similar to Section II we assume that z, the weight of an ion
impact, is distributed as an Erlang random variable with shape
parameter K where K is a Poisson random variable with mean
given in Eq. (23). Let a be an ion
P impact event with weight
za and expected number of ions i∈∂a w[i]. Then, similar to
Eq. (7) the probability of observing a given w and A can be
written as
(
)
X
0
P(a|w, A) = exp −
(w[i] + w ) δ(za )

(24)

(25)

Note that here we dropped the first term in Eq (26) since for
an observed ion impact event the weight za is always positive
and the first term vanishes.
The next step is to write the negative log-likelihood function. There is a subtle point to be noted here. F is the set of ion
impactP
events with non-zero weights. However, za can be zero
while i∈∂a w[i] > 0. We refer to such observations as zeroweight observations. Zero-weight observations are informative
and should be included in the log-likelihood function. In our
model the conditional probability of observing a zero-weight
event a in the interval [ta , ta ] is given by
(
)
X
P(a|w, A) = exp −
(w[i] + w0 )
 i∈∂a

ta
 X

= exp −
hAt , (w + w0 )i .



(28)

t=ta

The difficulty that seems to exist here is how to identify the
zero-weight ion impact events when they can span more than
one sample. However, as we shall see shortly, the particular
form of the log-likelihood function for the zero-weight events
eliminates the need to distinguish between adjacent zeroweight events for calculating the log-likelihood function.
Let F0 be the set of zero-weight observations and define
U0 = ∪a∈F0 [ta , ta ], i.e., the union of all time intervals we
observed zero-weight events at. Again, making the assumption

12

that distinct zero-weight ion impacts are independent events
given w and A we can write the joint probability of observing
F0 as
(
)
Y
X
P(F0 |w, A) =
exp −
w[i]
a∈F0

Y

=

i∈∂a

exp

a∈F0

= exp(−





X
t∈U0

−

ta
X
t=ta



hAt , wi


i∈∂a

hAt , wi).

(29)
e
And the gradient of the function `(w|y,
A) is

Note that P(F0 |w, A) does not depend on the number of zeroweight events or the beginning or end time of a particular
event. We need only to identify all the samples that are part
of a zero-weight event, i.e., not part of any observed ion impact
event.
Putting Equations (27) and (29) together and assuming that
all the events are independent we have the probability of
observing a trace y as
P(y|w, A) = P(F |w, A)P(F0 |w, A)
)
(
X X
0
= exp −
(w[i] + w )
exp −

XX
a∈F i∈∂a

∞
Y X
a∈F

k=1

(w[i] + w0 ) −

zak−1
µk (k − 1)!k!

X

za
µ

)

0

!k 

(w[i] + w )

(

a∈F

k=1

!k 
X

(w[i] + w0 )

,

i∈∂a

(30)
where the last equality is obtained since each sample on
the trace is presented in either F or F0 . The negative loglikelihood function then can be written as
`(w|y, A) = −N kwk1

∞
X
X
zak
−
log 
k
µ (k − 1)!k!
a∈F

k=1

ta
XX
a∈F t=ta

"
At

2




s X
s X
+ I0 (2 za
(w[i] + w0 )) + I2 (2 za
(w[i] + w0 ))
i∈∂a

i∈∂a

−1 #
s X
s X
2 za
(w[i] + w0 )I1 (2 za
(w[i] + w0 ))


i∈∂a

!k 
X

(w[i] + w0 )

Having the gradient of the log-likelihood function the algorithm is similar to the algorithm of section III with the gradient
of the log-likelihood calculated using Eq. (34). However, we
need some additional notations to represent the generalized
algorithm. Let deg(a)
Pn be the number of neighbors of event a,
i.e., deg(a) = i=1 Ata ,i 1 . For j ∈ [deg(a)] let ija be the
j
index of the j th non-zero element of Ata . Similarly, let ia be
j
the index of the j th non-zero element of Ata . Then, [ija , ia ] is
th
the true position of event a on the spectrum if its j neighbor
corresponds to the true scan that caused event a.
The peak picking algorithms are usually sensitive to the
shape of the pulses. Furthermore, the time of arrival of the
ions are noisier than the observation error of the instrument.
Observing many scans enables the instrument to measure the
MCR of the ions with precision significantly better than the
arrival noise level. To overcome the issue of a possible bias in
the estimated MCR in our model, we employ one last trick.
b by
The algorithm constructs an estimate of the spectrum x
assigning each observed event to its most likely neighbor (c.f.
Fig. 3). In other words, let



j

i∈∂a

j∗ = arg max

(31)

j∈[deg(a)]

After some algebra, and keeping only the terms that depend
on w,
`(w|y, A) = − N kwk1
!
X 1
X (w[i] + w0 )
−
log
2
µ
a∈F
i∈∂a
!
s
X (w[i] + w0 ) 
+ log I1 2 za
,
µ
i∈∂a

1
0
i∈∂a (w[i] + w )

P

(34)

i∈∂a

)
X za
= exp −N kwk1 −
µ
a∈F

∞
Y X
zak−1

k
µ (k − 1)!k!

e
∇`(w|y,
A) = −

i∈∂a

a∈F0 i∈∂a

(

e =
Using the same change of variable as before, namely w
1
e w|y,
e
w,
λ
=
N
µ
and
`(
A)
=
`(w|y,
A)
−
N
kwk
we
have
1
µ
e w|y,
e A) + λkwk
e 1 where
`(w|y, A) = `(
!
X 1
X
0
e w|y,
e A) = −
e +w )
`(
log
(w[i]
2
a∈F
i∈∂a



s X
e + w0 ) , (33)
+ log I1 2 za
(w[i]

ia
X

w[i].

(35)

i=ija

Then, given w we reconstruct the estimated spectrum as
∗

∗

b[ija + ∆] = x
b[ija + ∆] + y[t + ∆].
x

(36)

Using this notation the algorithm as as follows.
(32)

1 This definition is slightly inaccurate since an event can potentially fall on
P
Pn
the boundary of an scan resulting in n
i=1 Ata ,i 6=
i=1 Ata ,i but this is
rare and the discrepancy is negligible.

13

Algorithm
Input: trace y, firing times τ , and constants (θ0 , θ1 , λ, γ)
b
Output: estimated spectrum x
1: Calculate the adjacency matrix A as in Eq. (6).
2: Set w(0) = 0, θ = θ0 + θ1 .
2: Repeat until stopping criterion is met:
θ ← θ0 + k12 θ
1

e (k) ; y, A) .
w(k+1) ← ηθ w(k) − γ∇`(w
3:
4:

b=0
Set x
For a ∈ F :

j ∗ = arg max
j∈[deg(a)]

5:

Pija

i=ija

w[i]

For ∆ = 0, . . . , (ta − ta ):
∗
∗
b[ija + ∆] = x
b[ija + ∆] + y[t + ∆]
x
b.
Return x
R EFERENCES

[1] Franz Hillenkamp, Michael Karas, Ronald C. Beavis, and Brian T.
Chait, “Matrix-assisted laser desorption/ionization mass spectrometry of
biopolymers,” Analytical Chemistry, vol. 63, no. 24, pp. 1193A–1203A,
1991, PMID: 1789447.
[2] J.B. Fenn, M. Mann, C.K. Meng, S.F. Wong, C.M. Whitehouse, et al.,
“Electrospray ionization for mass spectrometry of large biomolecules.,”
Science (New York, NY), vol. 246, no. 4926, pp. 64, 1989.
[3] Franz Hillenkamp, Michael Karas, Ronald C. Beavis, and Brian T.
Chait, “Matrix-assisted laser desorption/ionization mass spectrometry of
biopolymers,” Analytical Chemistry, vol. 63, no. 24, pp. 1193A–1203A,
1991, PMID: 1789447.
[4] J. Ragoussis, G.P. Elvidge, K. Kaur, and S. Colella, “Matrix-assisted
laser desorption/ionisation, time-of-flight mass spectrometry in genomics
research,” PLoS genetics, vol. 2, no. 7, pp. e100, 2006.
[5] W.G. Cross, “Two-directional focusing of charged particles with a
sector-shaped, uniform magnetic field,” Review of Scientific Instruments,
vol. 22, no. 10, pp. 717–722, 1951.
[6] J.C. Schwartz, M.W. Senko, and J.E.P. Syka, “A two-dimensional
quadrupole ion trap mass spectrometer,” Journal of the American Society
for Mass Spectrometry, vol. 13, no. 6, pp. 659–669, 2002.
[7] W.E. Stephens, “A pulsed mass spectrometer with time dispersion,”
Phys. Rev, vol. 69, no. 691, pp. 46, 1946.
[8] M.T. Roberts, J.P. Dufour, and A.C. Lewis, “Application of comprehensive multidimensional gas chromatography combined with time-of-flight
mass spectrometry (GC× GC-TOFMS) for high resolution analysis of
hop essential oil,” Journal of separation science, vol. 27, no. 5-6, pp.
473–478, 2004.
[9] J.B. Fenn, “Electrospray wings for molecular elephants (nobel lecture),”
Angewandte Chemie International Edition, vol. 42, no. 33, pp. 3871–
3894, 2003.
[10] A. Shevchenko, I. Chernushevich, W. Ens, K.G. Standing, B. Thomson,
M. Wilm, and M. Mann, “Rapid’de novo’peptide sequencing by a combination of nanoelectrospray, isotopic labeling and a quadrupole/time-offlight mass spectrometer,” Rapid Communications in Mass Spectrometry,
vol. 11, no. 9, pp. 1015–1024, 1997.
[11] John S Janiszewski, Theodore E Liston, and Mark J Cole, “Perspectives
on bioanalytical mass spectrometry and automation in drug discovery,”
Current drug metabolism, vol. 9, no. 9, pp. 986–994, 2008.
[12] O. Trapp, J.R. Kimmel, O.K. Yoon, I.A. Zuleta, F.M. Fernandez, and
R.N. Zare, “Continuous two-channel time-of-flight mass spectrometric
detection of electrosprayed ions,” Angewandte Chemie International
Edition, vol. 43, no. 47, pp. 6541–6544, 2004.
[13] F.J. Knorr, M. Ajami, and D.A. Chatfield, “Fourier transform timeof-flight mass spectrometry,” Analytical Chemistry, vol. 58, no. 4, pp.
690–694, 1986.
[14] A. Brock, N. Rodriguez, and R.N. Zare, “Hadamard transform timeof-flight mass spectrometry,” Analytical Chemistry, vol. 70, no. 18, pp.
3735–3741, 1998.
[15] G.S. Moore, M. Manlove, and A. Hidalgo, “Statistical Inference Algorithm for DitheredMulti-Pulsing Time-of-Flight Mass Spectrometry,”
Submitted, 2012.

[16] Robert Tibshirani, “Regression shrinkage and selection via the lasso,”
Journal of the Royal Statistical Society. Series B (Methodological), pp.
267–288, 1996.
[17] Z. Harmany, R. Marcia, and R. Willett, “This is SPIRAL-TAP: Sparse
Poisson Intensity Reconstruction ALgorithms– Theory and Practice,”
Image Processing, IEEE Transactions on, , no. 99, pp. 1–1, 2010.
[18] J.L. Wiza, “Microchannel plate detectors,” Nucl. Instrum. Methods, vol.
162, no. 1-3, pp. 587–601, 1979.
[19] A. Findling, “A family of logarithmically concave functions defined by
an integral over the modified bessel function of order 1,” Journal of
mathematical analysis and applications, vol. 194, no. 2, pp. 368–376,
1995.
[20] Peng Zhao and Bin Yu, “On model selection consistency of lasso,” The
Journal of Machine Learning Research, vol. 7, pp. 2541–2563, 2006.

14

PLACE
PHOTO
HERE

Morteza Ibrahimi Morteza Ibrahimi received his
PhD in Electrical Engineering from Stanford University in 2013, working with professor Andrea
Montanari. His research interests are in high dimensional statistical signal processing, learning graphical models, optimization through message passing
algorithms, and statistical inference with high dimensional data, specially through fast iterative algorithms. He received his BS from Sharif University
of Technology in 2006, and MSc from University of
Toronto in 2007.

Andrea Montanari Andrea Montanari received a
Laurea degree in Physics in 1997, and a Ph. D.
in Theoretical Physics in 2001 (both from Scuola
Normale Superiore in Pisa, Italy). He has been postPLACE
doctoral fellow at Laboratoire de Physique Thorique
PHOTO
de l’Ecole Normale Suprieure (LPTENS), Paris,
HERE
France, and the Mathematical Sciences Research
Institute, Berkeley, USA. Since 2002 he is Charg
de Recherche (with Centre National de la Recherche
Scientifique, CNRS) at LPTENS. In September 2006
he joined Stanford University as a faculty, and since
2010 he is Associate Professor in the Departments of Electrical Engineering
and Statistics. He was co-awarded the ACM SIGMETRICS best paper award
in 2008. He received the CNRS bronze medal for theoretical physics in 2006
and the National Science Foundation CAREER award in 2008.

PLACE
PHOTO
HERE

George Moore Goerge S. Moore is a research
fellow at Agilent technology. He received his PhD in
1980 from Purdue University and his Bsc/Msc from
Mississippi State University in 1974 all in Electrical
Engineering. He has over 40 years of experience in
the industry.

