First, we model the mass spectrometry system based on parameters like the mass and charge of each constituent, assuming a certain number of constituents in the sample. Here, a constituent is defined as a substance with a specific monoisotopic mass. We then perform a MAP (Maximum A Posteriori) estimation of these parameters from the observed spectrum. By comparing the maximum posterior probability in models with different numbers of constituents, we determine the model with the most appropriate number of constituents.

However, this model has a large dimensionality of the number of constituents multiplied by 6. Moreover, the posterior probability for one of the parameters, the monoisotopic mass, is flat over a large portion of the search space and has several sharp peaks locally. Hence, gradient-based methods are not suitable for this case due to anticipated gradient vanishing. Therefore, to estimate the parameters, we combine the No-U-Turn Sampler (NUTS (Hoffman and Gelman, 2014), a type of Markov Chain Monte Carlo (MCMC), with Simulated Annealing (Kirkpatrick, Gelatt, and Vecchi, 1983).

The purpose of using Simulated Annealing is to introduce a temperature parameter. By selecting a high-temperature exploration parameter distribution, we can actively explore parameters even in areas where the posterior probability is flat or has sharp peaks. This ensures a broader search across the parameter space, reducing the chance of overlooking the global solution and getting trapped in local minima.

Furthermore, NUTS can explore parameters sparsely in areas with small gradients and can explore parameters in detail in areas with large gradients. Thus, introducing NUTS allows efficient exploration of the vast, high-dimensional parameter space.

On the other hand, while MCMC is good at searching for global solutions, it does not always reach the optimal solution within a certain number of search steps. Therefore, we use the parameters with the highest posterior probabilities obtained from NUTS and Simulated Annealing as initial values and apply stochastic variational inference. By doing this, we search for the optimal parameter where the posterior probability is maximized in the vicinity of that initial value, aiming to improve the accuracy of parameter estimation.

However, simultaneously searching for parameters for all possible numbers of constituents leads to a curse of dimensionality, where the search space explosively expands as the number of constituents increases, potentially reducing search efficiency and accuracy. To avoid this problem, we sequentially increase the number of constituents from k = 1 to the maximum conceivable number k = k max , and for k constituents calculate the optimal parameters and their posterior probabilities, and efficiently focus the parameter search areas for the k+1 constituents by these posterior probabilities.

To balance the complexity of the model (number of constituents) and its fit (loss against the data), in addition to the prior distribution of each parameter, we introduce a prior distribution for the number of constituents. We also incorporate a prior distribution on the differences between the monoisotopic masses of multiple constituents. For analytical purposes, we have defined a single constituent as a substance with a distinct monoisotopic mass, thereby ensuring that their masses do not mutually take the same value. When seeking to separate isomers, it is essential to integrate other techniques such as fragmentation, ion mobility spectrometry, and chromatography, in addition to the proposed method. We first construct a model with k = 1 constituent, obtain the optimal parameters and the maximum posterior probability based on the above prior distributions and observed data.

Next, we construct a model with k = 2 constituents. For one of the two constituents, we use a prior distribution centered on the optimal parameters already estimated for k = 1 , narrowing its range. This suppresses the significant increase in the parameter search space. Based on this new prior distribution, we estimate the optimal parameters and obtain the maximum posterior probability.

Subsequently, we seek the maximum posterior probability for each model with constituent numbers up to the upper limit k max   by efficiently exploring the optimal parameters in the same manner.

Finally, we compare the maximum posterior probabilities corresponding to each model with different numbers of constituents. We select the model with the highest probability and obtain the estimates for the monoisotopic masses and ion counts.

In a time-of-flight mass spectrometer, the relationship between time of filght t t o f _ t h e o r y , mass m , and charge z can ideally be represented as in Eq. 1 (Boesl, 2017). Here, V represents the acceleration voltage to eject ions, and y is the flight distance of ions. t t o f _ t h e o r y = 0.72   y m z ∙ 1 V (1)

The actual signal obtained is a convolution of the delta function δ t t o f with the detector’s response waveform R t . Preliminary experiments have shown that t t o f has a stochastic ToF jitter ε following the logarithmic normal distribution Λ t o f μ t o f , σ t o f , E ε = ⁡ exp μ t o f + σ t o f 2 2 = 0 ( E : e x p e c t e d   v a l u e ). Additionally, the height (intensity) of the response waveform also has a stochastic response factor α following the logarithmic normal distribution Λ r e s μ r e s , σ r e s .

Reflecting these variances in time of flight and detector response intensity, the response S t when a single ion with mass m and charge z enters the detector according to Eq. 1 can be written as: S t , m , z = α ∙ R t * δ t − t t o f _ t h e o r y − ε = α ∙ R t * δ t − 0.72   y m z ∙ 1 V − ε (2)

The mass m and charge z of constituent j follow a distribution as described below. Since the natural isotopic abundances differ by element, and each functional group has different charge rates, the probability mass functions of mass and charge are represented by a multinomial distribution as: p j m = m j ′ + m ! ∏ a ∊ H , C , N , O , S 1 n a ! u a ξ m (3) q j z = l j ! ∏ b = 1 l j v b z (4) j : c o n s t i t u e n t   I D s   j = 1,2 , ⋯ k k : n u m b e r   o f   c o n s t i t u e n t s   i n   t h e   s a m p l e m j ′ : m o n o i s o t o p i c   m a s s   o f   c o n s t i t u t e n t   j n a : n u m b e r   o f   h y d r o g e n , c a r b o n , n i t r o g e n , o x y g e n , p h o s p h o r u s , s u l f u r   a t o m s   o f   c o n s t i t u e n t   j   a = H , C , N , O , S u a : n a t u r a l   i s o t o p i c   a b u n d a n c e   o f   H 2 , C 13 , N 15 , O 18 , S 34   o f   c o n s t i t u t e n t   j ξ a : i n c r e a s e d   n u m b e r   o f   n e u t r o n s   a = H , C , N , O , S l j : n u m b e r   o f   f u n c t i o n a l   g r o u p   o f   c o n s t i t u t e n t   j v b : e a c h   c h a r g e   r a t e   o f   f u n c t i o n a l   g r o u p   o f   c o n s t i t u t e n t   j   b = 1,2 , ⋯ l j

However, this model becomes too complex for Bayesian inference due to its large number of parameters ∑ j = 1 k l j + 10 . Therefore, we introduce an approximation that assumes equal isotopic abundances for all elements, an increase of one neutron in isotopes, and equal charge probabilities for all functional groups.

Consequently, the mass distribution of constituent j can be represented by a binomial distribution based on the monoisotopic mass as a reference, where n j atoms contained in a molecule are replaced by isotopes with an increase of +1Da at a probability of u j . Similarly, for the charge distribution, l j functional groups can be represented by a binomial distribution where they acquire a charge of +1 (in the case the mass spectrometry system is in positive mode) at a rate of v j . This reduces the number of parameters down to six dimensions. p ∼ j m = m j ′ + n j m u j m 1 − u j n j − m (5) q ∼ j m = l j z v j z 1 − v j l j − z (6) n j : r e p r e s e n t a t i v e   n u m b e r   o f   a t o m s   o f   c o n s t i t u t e n t   j u j : r e p r e s e n t a t i v e   n a t u r a l   i s o t o p i c   a b u n d a n c e   o f   c o n s t i t u t e n t   j v j : r e p r e s e n t a t i v e   c h a r g e   r a t e   o f   f u n c t i o n a l   g r o u p   o f   c o n s t i t u t e n t   j

Considering that constituent j contains up to i j ions, its spectrum is represented as: S ¯ j t = ∑ i = 1 i j S t , m i , z i (7) m i ∼ p ∼ j z i ∼ q ∼ j i j : n u m b e r   o f   i o n s   i n   c o n s t i t u e n t   j

If the ion counts are sufficiently large, it can be approximated as: S ¯ j t ∼ i j ∑ m ∈ Q + ∑ z ∈ Z + p ∼ j m ∙ q ∼ j z ∙ S t , m , z (8)

Noise in the time-of-flight mass spectrometer is known to be stationary and follows a normal distribution based on preliminary experiments. It is known that the thermal noise of detection circuits in such as mass spectrometers follows a normal distribution (Johnson, 1928), suggesting that in this case, thermal noise is the dominant factor in overall noise. Therefore, the noise to be added to the entire spectrum is represented as r t ∼ N μ n o i s e , σ n o i s e .

Considering these, the conclusive spectrum combined with the multiple spectra of single constituent is represented as: S ¯ t = ∑ j = 1 k S ¯ j t + r t   (9)

When the observation data from the mass spectrometer S ¯ o b s t is obtained, assuming the number of constituents as k , the posterior probability distribution P θ | S ¯ o b s t for parameters θ : m 1 ′ , i 1 , n 1 , u 1 , l 1 , v 1 , ⋯ , m k ′ , i k , n k , u k , l k , v k is defined as per Bayes’ theorem. Note that P S ¯ o b s t | θ represents the likelihood of parameters θ when S ¯ o b s t is provided, and P θ denotes the prior distribution. P θ | S ¯ o b s t ∝ P S ¯ o b s t | θ P θ (10)

Here, in addition to the prior distribution of each parameter (uniform distribution), we incorporate a regularization term, w b i c to achieve a suitable balance between model complexity (number of constituents) and model fit (loss with respect to data). We also introduce a regularization term, w e x , to prevent multiple constituents within the same model from assuming the same monoisotopic mass. Hence, we introduce the following logarithmic prior distribution: log P θ ∝ − w b i c + w e x (11)

To determine the appropriate number of constituents k , we define the regularization term w b i c representing the complexity of the model with k based on the Bayesian Information Criterion (BIC). The BIC is a statistical measure that balances the fit to the data and model complexity (Schwarz, 1978; Neath and Cavanaugh, 2012). Here, n represents the dimension of the observation data S ¯ o b s t , which in this study is the number of data points in the time direction. w b i c = λ ⋅ k 2 ⋅ ⁡ log ⁡ n λ : 300 h y p e r   p a r a m e t e r

Furthermore, we define a constituent by its unique monoisotopic mass. Therefore, if the estimated values of the monoisotopic mass parameters of multiple constituents are the same in the algorithm, the count of constituents will not be accurate. Here, we define the logarithmic prior distribution (regularization term) w e x as follows, using a penalty that increases exponentially according to the difference in estimated monoisotopic mass values, as shown in Figure 1. The integral of the spectrum ∫ 0 ∞ S ¯ o b s t   d t is also multiplied as a coefficient to ensure that the impact of the penalty does not change depending on the scale of the observed data. w e x = ∑ i = 1 k − 1 ∑ j = i + 1 k a ∫ 0 ∞ S ¯ o b s t   d t 2 b exp − m i − m j b (12) a = 0.001 , b = 0.1   h y p e r   p a r a m e t e r

Here, by substituting the parameter θ generated from MCMC into model (9), we obtain the spectrum as S ^ p r o p t . We assume a normal distribution for the noise. The standard deviation of the noise denoted as σ is set to 2,000. Consequently, the logarithm of the posterior probability distribution is as follows. log P θ | S ¯ o b s t ∝ ⁡ log P S ¯ o b s t | θ P θ = ⁡ log P S ¯ o b s t | θ log P θ ∝ L : = β ⁡ log 1 2 π σ 2 exp − S ^ p r o p t − S ¯ o b s t 2 2 σ 2 − β w b i c + w e x = − β 1 2 σ 2 S ^ p r o p t − S ¯ o b s t 2 − n   log σ − n 2   log 2 π − w b i c − w e x (13)

Here, β represents the inverse temperature. As described in 2.1, we use Simulated Annealing to ensure active parameter exploration in flat areas or sharp peaks of posterior probability. This is achieved by multiplying the inverse temperature β   < 1 to the posterior probability. Initially starting from a low inverse temperature value (i.e., high temperature) and gradually increasing to a higher value (i.e., low temperature). At low inverse temperatures (high temperatures), the system explores a wide parameter space. Conversely, at high inverse temperatures (low temperatures), the system converges to the optimal solution. This time, we set the temperature change in three stages: β = 0.2 5 → 0.2 4 → 0.2 3 .

To obtain the maximum posterior probability and parameters θ that maximize the posterior probability (formula (13)), we conduct sampling from this posterior probability distribution using MCMC.

From the posterior probability distribution P θ | S ¯ o b s t , we sample the parameter θ to select the one that maximizes the posterior probability. We use the No-U-Turn Sampler (NUTS) for sampling, a recent and popular variant of the Markov Chain Monte Carlo (MCMC) method. NUTS is a type of MCMC, especially a derivative of the Hamiltonian Monte Carlo method (HMC) (Neal, 2011; Betancourt, 2017).

After executing MCMC, the parameters of the maximum posterior probability obtained are inherited as initial values, and optimization of the parameters is performed using Stochastic Variational Inference [SVI (Kingma and Max, 2013; Wingate and Weber, 2013; Ranganath et al., 2014)]. For more details, please refer to the Supporting Material.

The overall picture of the workflow to determine the optimal parameters and posterior probability for each assumed number of constituents from the observational data of the mass spectrometer is as shown in Figure 2.

First, as described in 2.3.1, (i) input the observational data of the mass spectrometer with dimensions of flight time and ion counts. Then (ii) assume that the number of constituents, k , contained in the sample is 1. (iii) Set the inverse temperature to 0.2 5 .

Next, as described in 2.3.2, (iv) sample 1,000 × 4 times from the posterior probability distribution, (v) set the MAP solution obtained by MCMC as the initial value for the next MCMC. Then (vi) divide the inverse temperature by 0.2. Repeat steps (iv) to (vi) three times.

As described in 2.3.3, (vii) set the parameter of the maximum posterior probability obtained by MCMC as the initial value for SVI. Then (viii) optimize the parameters with SVI, and (ix) set the parameter of the maximum posterior probability obtained by SVI as the initial value for the next MCMC.

Following 2.3.1, increase the number of constituents, k , by 1. Repeat steps (iii) to (vi). (x) Continue this until the maximum possible number of constituents, k max . Finally, (xi) compare the maximum posterior probabilities of models from constituents k = 1 to k max , and (xii) select the model with the largest posterior probability. Also, obtain the optimal parameters at that time.

The specifications of the PC used for verifying the proposed method, as well as the software versions, are detailed in Table 1. The proposed method handles data with 1 million dimensions along the time axis, requiring a large memory size. Additionally, to rapidly explore a wide 6-dimensional parameter space using MCMC, the high-speed probabilistic programming library, NumPyro, along with its compatible CUDA and GPU, were used.

Based on the nucleic acid drug Fomivirsen (Perry and Balfour, 1999) (ID: A), four impurity constituents with modified base sequences were added, and spectra for a total of five constituents were generated via simulation. Specific values are as per Table 2. This enables the replication of a system where the principal constituent’s isotopic distribution and the impurity spectra are mixed.

Ion counts for each constituent were set at 20,000. To facilitate the interpretation of results and to ensure that the algorithm treats each constituent fairly, we will conduct evaluations using a 1:1 concentration ratio for each component in the proposed method. The atomic counts n a for H, C, N, O, and S were obtained from the molecular formula of each constituent. Natural isotopic abundance ratios u j followed the NIST Atomic Weights and Isotopic Compositions for All Elements. The representative functional group number l j and the representative charge rate v j were set to 224 and 0.035, respectively, to ensure that the generated spectra resembled real data.

The procedure involved sampling from the multinomial distribution represented by Eqs 3 and 4 20,000 times (total incoming ion counts) for each constituent. Subsequently, spectra were formed following the procedures in Eqs 2 and 7.

The mutation from C (Cytosine) to U (Uracil) is called deamination and is generated in the synthesis process due to solvent conditions and thermal stress (Stavnezer, 2011; Gao, Choudhry, and Cao, 2018).

The spectra of the generated single constituents A to E were combined according to Eq. 9 in the 15 combinations listed in Table 3. This allows for a comprehensive combination of 2-3 constituents based on constituent A, as well as an evaluation of each individual constituent. We use these as test data.

The results of estimating the number of constituents in the spectra of the test data (Mixture No.1∼15) using our proposed method are as shown in Table 4. The values within the table represent the negative logarithm of the maximum posterior probability in the model of constituent count k. Therefore, the smallest value should be selected.

By choosing the most suitable number of constituents based on this criterion, the success rate for estimating the true number of constituents was 80% (12/15). Additionally, the presence or absence of impurities (distinguishing between k = 1 and k≧2) could be determined with 100% accuracy. We believe this is sufficient as a standard for recognizing the presence and number of impurities in pharmaceuticals and taking appropriate measures.

The optimal monoisotopic masses and ion counts estimated in the model where the posterior probability is maximum for each test data are shown in Table 5. The monoisotopic mass had an average error of 1.348 Da and a maximum error of 4.931 Da. This is insufficient to determine how many mutations have occurred, making it unsuitable for examining the cause of impurity generation with a difference of 1 Da. Regarding the ion counts, there was an average error of 4% and a maximum error of 82%. For instance, the standards for total desamido impurity and total impurities in injectable glucagon are 14% or less and 31% or less, respectively (Bao et al., 2022). Therefore, the accuracy of the ion count estimation in the proposed method is insufficient to estimate the impact of impurities.

For reference, a comparison between the spectra reconstructed from the estimated parameters and the original signal is shown in Figure 3. The overall view in (a) represents the charge distribution, and the enlarged view in (b) represents the isotopic distribution. From these results, it is clear that the spectrum we generated closely matches the observed data. Despite the spectra matching, errors in parameter estimation occurred because of the high degree of freedom in isotopic parameters that trade-off with monoisotopic mass. Even if the monoisotopic mass was lower than the true value, by increasing the representative atomic number n j or the representative isotopic natural abundance u j , it is possible to make it fit the observed data to some extent.

Also, the estimated ion counts of each constituent showed errors of up to 82% from the true values. This is presumed to be due to the trade-off relationship between the ion counts of each constituent, with a decrease in the ion count of one constituent being compensated by an increase in another. This is further supported by the fact that the average error in ion counts settles at 4%.

Deconvolution of the test data was performed using the existing method, UniDec as well. The results of deconvolution for each observed spectrum by UniDec are shown in Table 6. According to these results, the accuracy for the correct number of constituents was 13% (2/15). This is presumed to be because the UniDec algorithm, which obtains the number of constituents after multiple iterations of deconvolution, does not necessarily guarantee the number of constituents. Please note that this use of UniDec to determine the number of constituents is not its intended application.

For the verification above, we used UniDec (Version 6.0.2). The particularly set parameters during this verification are shown in Table 7. The Mass Range was set to the same range as the proposed method, and Sample Mass Every (Da) was set to 0.1 to sufficiently detect impurities with a difference of 1 Da. For parameters not mentioned, default values were used.