Improving the univariate calibration approach with Bayesian modeling for IR reaction monitoring

Abstract

Infrared (IR) spectroscopy is a powerful tool for real-time reaction monitoring in chemical synthesis. For applications that require tracking concentration profiles of reactive species, univariate linear calibration models are commonly used to relate IR signals to analyte concentrations. Despite their simplicity, the accuracy of these models can be limited by spectral overlap and other effects that distort the linear relationship between concentration and signal. To address this limitation, chemometric models are often employed, typically without further examination of opportunities to improve univariate calibration performance itself. Here, we present a novel workflow based on Bayesian statistics to enhance univariate calibration for IR reaction monitoring. The central feature of this workflow is the use of three diagnostic Bayesian probabilistic models, combined with data-preprocessing selection, to screen for IR signals that can potentially improve univariate calibration performance when non-linear effects are present. We applied the workflow to a test reaction system and identified an IR signal in the fingerprint region, along with an uncommon preprocessing strategy, that reduced prediction error by more than 50% compared with the univariate model using the original preprocessing steps. Overall, our workflow aims to improve the usability of univariate calibration approaches and expand the toolbox available to chemists for IR monitoring of complex chemical processes.

---

1 Introduction

Smart sensors integrated with a variety of geometric probes are a suite of tools that enable rapid, automated, and robust measurements of chemical systems in real time.^1,2 These technologies provide continuous and real-time measurements of key species throughout the course of a chemical process. For kinetic studies, the ability to track concentration changes in situ allows researchers to capture transient phenomena, directly observe rate-limiting steps, and construct mechanistic models with higher fidelity.^3,4 In process development and manufacturing, these technologies enable both manual and automated control for product quality and process performance.⁵ This is especially valuable for maintaining product consistency, improving safety, and increasing overall process efficiency.

Infrared (IR) spectroscopy has become a mainstay among analytical techniques for real-time reaction monitoring. The in situ nature of IR spectroscopy allows direct analysis of the reaction medium, avoiding delays and artifacts associated with offline sample handling or transfer. In its simplest use, unique IR signals are first identified, and key process information, such as conversion is obtained by tracking the evolution of these signals throughout the reaction.⁶ For more quantitative analysis, offline samples are collected simultaneously via complementary techniques to obtain analyte concentrations. These concentration values are then linked back to the corresponding IR signals to construct a univariate linear regression model based on Beer–Lambert's law.^7,8

Despite the straightforward nature of the univariate linear regression approach, recent studies^8,9 cautioned against its use, even when a distinct spectroscopic feature can be assigned to the analyte of interest. Chemical reactions often exhibit higher-order effects. The changes in the observed spectroscopic response are influenced not only by analyte concentration but also by additional factors such as pH, temperature, or the presence of other components in the reaction medium. When such effects arise, chemometric models¹⁰ are frequently employed to address the complex dynamics in the data. Very few works, however, have been done to improve the univariate model in the case of higher-order effects. For example, spectral preprocessing optimization is often used to improve chemometric modeling accuracy,^11,12 but such optimization is rarely applied in the development of univariate calibration models.

Bayesian models have gained increasing attention in recent years for diverse applications in chemistry, including reaction and process optimization,^13,14 kinetic parameter estimation, model selection,^15,16 and development of spectroscopic calibration models.^17,18 Their ability to incorporate prior knowledge (e.g., assumptions about measurement noise) and to treat model parameters and predictions probabilistically makes them especially valuable for data analysis in low-data regimes.¹⁶

Here, we present a workflow based on Bayesian modeling to more deeply interrogate the use of the univariate linear regression approach for IR reaction monitoring. The central idea is to perform Bayesian model criticism. Rather than immediately discarding the univariate model following an initial lack of fit, we employ a Bayesian hierarchical linear regression model to assess the potential for performance improvement through optimization of preprocessing steps. When such opportunities are identified, a grid-search optimization campaign is used to locate the optimal preprocessing pipeline. Finally, the improvements in model performance are evaluated using Bayesian posterior predictive checks. The aim of this work is to provide a practical roadmap for improving univariate calibration approaches for IR reaction monitoring, thereby expanding the set of tools available to chemists when working with challenging and dynamic reaction mixtures.

2 Methods

2.1 Materials and reaction details

(S)-(−)-1-Phenylehtanol (1, 97%), isobutyric anhydride (2, 97%), triethylamine (Et₃N, 99.5%) and dichloromethane (DCM) were purchased from Sigma Aldrich and used as received. (−)-Benzotetramisole (btm, 98%) was purchased from TCL America and used as received. The reaction scheme (Fig. 1a) was adapted from the referenced work.¹⁹ The selected concentration ranges for the reactants reflect those typically used in kinetic experiments aimed at elucidating reaction mechanisms.²⁰

Fig. 1 Top: Reaction scheme for the BTM-catalyzed acylation. Bottom: Reaction conditions used in the training and validation experiments.

2.2 In situ IR and online HPLC

IR spectroscopy were performed on a ReactIR 702L Infrared spectrometer with a 6.3 mm AgX DiComp diamond probe (Mettler Toledo). Spectra acquisition was achieved by averaging 64 scans between 4000–650 cm⁻¹ with a resolution of 8 wavenumber per data point (the interval between each spectrum is 1 minute). Automated on-line HPLC analyses were conducted on a 1290 infinity HPLC system (Agilent), a EasySampler 210 probe (Mettler Toledo), and a Directinject-LC system (Telescope Innovations). Details regarding the calibration and operation of this automated online HPLC reaction-monitoring platform are provided in this work²¹ and in the SI (Fig. S3).

2.3 Orthogonal reaction monitoring with IR and online HPLC

All reactions monitored with the in situ Infrared spectrometer and on-line HPLC reaction monitoring platform were performed in a 25 ml 3-neck round-bottom flask (RBF) on a magnetic stir plate. The IR probe was initially inserted into the empty RBF, and an air background spectrum was taken. The flask was then filled with 10 ml of DCM. A solvent background spectrum was recorded, followed by the sequential addition of alcohol 1, triethylamine, and s-benzotetramisole. The IR spectrum acquired at this stage was designated as the T₀ data point for the product concentration in the IR data stream. Simultaneously, an HPLC sample was collected through the on-line sampling setup to set the T₀ data point for the product concentration in the LC data stream. Subsequently, the desired amount of anhydride 2 was added to the flask to initiate the reaction, after which the HPLC and IR sampling sequence were initiated simultaneously to generate the orthogonal reaction monitoring dataset.

A separate nuclear magnetic resonance (NMR) reaction monitoring experiment (Fig. S4) was conducted to confirm the mass balance between starting material 1 and the product (i.e., 1 was converted exclusively to products in this reaction). With this information, the method introduced in this work²² was used to convert the HPLC peak area for 1 and products into concentrations without requiring an external calibration curve (Fig. S4).

2.4 Bayesian inference

All Bayesian analysis were performed in Python using the probabilistic programming package PyMC.²³

3 Results and discussions

3.1 Likelihood function selection for Bayesian inference

The key step in Bayesian modeling is to calculate posterior distributions p(θ|D) of model parameters θ based on observed data D. According to Bayes' theorem:the numerator on the right-hand side comprises the likelihood, p(D|θ), and the prior distribution, p(θ). The likelihood represents the probability of observing data under the assumption of the current model, and the prior encodes existing knowledge or assumptions about the parameters before the data are considered.

The Gaussian likelihood is widely used in chemistry applications,^17,18 based on the assumption that measurement noise is normally distributed and independent of the magnitude of the measurement. In this work, concentration measurements were obtained from an online HPLC reaction monitoring platform. Because the fluidics transfer and dilution steps may introduce noise that scales with signal magnitude, we also evaluated an alternative likelihood function: the lognormal distribution. To assess which likelihood better reflects the measurement process, a double exponential decay function was fitted to the observed concentration profiles of 1 to approximate the mean signal behavior. These mean estimates were then used to infer the noise parameter for each candidate likelihood. The resulting log-likelihood values were similar for both models (Fig. S6), indicating that the current dataset does not provide sufficient evidence to favor one likelihood function over the other. Consequently, we selected the Gaussian likelihood for the remainder of the analysis due to its computational simplicity.

3.2 Diagnostic model and model comparison metric descriptions

Three models were considered in the subsequent Bayesian model comparison. All three models include a noise parameter inferred from the entire training dataset. The first model (M1) is a univariate linear regression model. It contains one additional parameter: the response factor. The second model (M2) is a hierarchical linear regression model. In this model, each experiment in the dataset has a unique response factor. These response factors are related through a common parent distribution, reflecting the assumption that all experiments are influenced by a shared higher-order effect. The third model (M3) is a multivariate linear regression model. This model includes multiple weight parameters corresponding to the number of features in the input IR spectrum. An additional regularization parameter is incorporated to prevent overfitting.

Model performance was assessed using Bayesian leave-one-out (BayesLOO) scores. BayesLOO is calculated by training the model on all but one data point, followed by validation on the left-out point. This procedure is repeated until each data point in the training set has been left out once. The BayesLOO score of M3 was used as a reference (best-case scenario) for subsequent comparisons, as the interpretability of BayesLOO scores is more meaningful when used in relative comparisons.

3.3 Workflow overview

The workflow (Fig. 2) begins at the validation reaction data collection stage. Ideally, more than two validation experiments should be conducted, covering variations in reaction conditions. An initial IR signal is selected according to chemical intuition (e.g., signal around 1600–1900 cm⁻¹ should be selected for analytes with carbonyl function groups) for building M1 and M2. If no such signal can be identified, then the IR signal in the fingerprint region (400–1500 cm⁻¹) that shows a positive correlation with the reactant of interest should be considered next. The full IR spectra is used to build M3. Posterior distributions and BayesLOO scores are then computed for all three models.

Fig. 2 Flowchart illustrating the decision steps in the proposed workflow. The key step is the comparison of BayesLOO scores between the hierarchical linear model (M2) and the multivariate linear model (M3). The outcome of this comparison indicates whether preprocessing optimization is likely to have a meaningful impact on improving the performance of the univariate linear regression model.

In the first scenario, both M1 and M2 exhibit much lower BayesLOO scores than M3. This indicates that higher-order effects are present at the individual level (within each experiment) for the selected IR wavenumber region. In this case, chemometric models should be employed to address these complex effects, or an alternative IR signal can be selected to re-enter the model comparison.

In the second scenario, M1 shows a low BayesLOO score, but M2 performs reasonably well relative to M3. This suggests that the linear relationship between concentration and IR response is maintained within each individual experiment. But higher-order effects, which likely originate from the variation in different experimental conditions, still influence the analyte IR responses. In this scenario, applying data preprocessing optimization will have a considerably higher chance of mitigating the higher-order effects compared to the first scenario. After such a scenario is identified, principal component analysis (PCA) is applied to understand the nature of the variations. This step will inform us which preprocessing categories we should consider for the subsequent optimization. Then, a grid search is performed on all possible combinations of preprocessing steps to identify an optimal pipeline that will enhance the M1 performance. Finally, Bayesian posterior predictive checks are performed to assess the quality of the improvements.

In the third scenario, M1 performs comparably to M3. This suggests that higher-order effects are negligible in the current application, and the univariate model can be reliably used for reaction monitoring.

3.4 Workflow validation on the BTM-catalyzed acylation reaction

For the following analysis, three validation experiments (Fig. 1 experiment 1, 4, and 7) were selected. The reactive species of interest was the acylation product, which contains a carbonyl stretch signal at 1730 cm⁻¹ (Fig. 3a). Using this IR signal, BayesLOO scores were computed for the three diagnostic models.

Fig. 3 (a) Overlay of reaction IR spectra, with greener colors indicating later reaction times. “Valleys” denote the a different set of endpoints used in this work for IR peak area integration. (b) Comparison of Bayesian leave-one-out (BayesLOO) scores for different models and preprocessing strategies. Higher scores indicate better model performance. The BayesLOO score of M3 is used as a best-case scenario to assess the performance of the remaining models.

The results showed that M1 and M2 both yielded substantially lower BayesLOO scores than M3 (Fig. 3b), indicating that this IR region is affected by higher-order effects at the individual experiment level. As a result, the univariate linear regression model is not suitable for calibration in this spectral region. This behavior may arise from the overlap between the product carbonyl stretch and the carbonyl signal of the byproduct isobutyric acid.

Next, two IR signals (1156 cm⁻¹ and 1178 cm⁻¹) that showed strong correlation with the product concentration were identified in the IR fingerprint region. These features likely originate from carbon–oxygen single-bond stretching or bending of the product. Recalculation of the BayesLOO scores using these newly selected IR signals showed substantial improvement in the performance of M2 (Fig. 3b), whereas M1 continued to perform poorly relative to M3 (Table 1, entry 3). These results suggest that the selected fingerprint-region signals, particularly the one at 1156 cm⁻¹, exhibit a strong linear relationship with product concentration. But as the experimental condition changed, a common higher-order emerged and altered the nature of this linear relationship.

Table 1 Overview of the eight preprocessing pipelines evaluated in the optimization campaign. The optimal preprocessing pipeline (bold entry) demonstrated a significant improvement in M1 performance over the original preprocessing pipelines, as indicated by the BayesLOO scores

Index	IR signals	Baseline correction	IR peak range	elpd_loo
1	Peak height	First derivative	Zero	250.7 ± 5.9
2	Peak height	First derivative	Valley	259.7 ± 3.4
3	Peak height	Second derivative	Zero	218.0 ± 4.2
4	Peak height	Second derivative	Valley	223.4 ± 4.1
5	Peak area	First derivative	Zero	84.1 ± 1.5
6	Peak area	First derivative	Valley	292.8 ± 3.6
7	Peak area	Second derivative	Zero	148.5 ± 3.8
8	Peak area	Second derivative	Valley	240.5 ± 3.0

---

3.5 Investigation of higher-order effects and preprocessing step optimization

Principal component analysis (PCA) was applied to the wavenumber region around the peak at 1156 cm⁻¹ to investigate the origin of the higher-order effects. The first principal component accounted for the majority of variance (∼91%). The PC-1 transformed IR data (Fig. 4c, black solid circle), also referred to as the PC-1 scores, closely resembled the product concentration profiles, indicating the largest variation in this IR region is dominated by the product formation.

Fig. 4 Loading plots for (a) the first and (b) the second principal components of the IR spectra. The magnitude of the y-axis denotes the contribution of each wavenumber to the overall variance. (c) PC-1 (dark solid circles) and PC-2 (gray hollow circles) scores plotted as a function of reaction time. Comparing PC scores with measured process variable trends enables identification of the sources of variation captured by each principal component.

PC-2 captured a smaller (∼9%) but still meaningful portion of variance. The PC-2 scores (Fig. 4c, gray hollow cycle) for exp 1 exhibited an upward trend, whereas those for experiment 4 and 7 exhibited downward trends. We therefore reasoned that the variation in the second PC was not directly associated with the change of the major reactive species shown in Fig. 1. This is because such changes would be expected to produce consistent trends across different experiments. We also compared PC-2 scores with the temperature profiles (measured by the in situ IR probe), but did not observe any meaningful correlation. We speculate that the varying amounts of triethylamine used in these experiments may influence the protonation states and hydrogen bonding interactions among reaction species, thereby giving rise to the variations observed in PC-2. In future studies, a pH probe could be incorporated to track pH changes across different experiments. Comparing the pH profiles with the PC scores will provide additional evidence to support or refute the current hypothesis.

PC-2 loading plot showed that the wavenumbers at 1140 and 1164 cm⁻¹ were the two strongest contributors to PC-2 (Fig. 4b). Conventionally, after baseline correction, the maximum peak height of the identified IR peak is used as the IR signal. Based on the observation in PC-2 analysis, we reasoned that incorporating additional information from the surrounding spectral region into the IR signal calculation might help mitigate the higher-order effects. To test this idea, peak area integration was introduced as an alternative way for IR signal calculation, and the integration bounds were extended from a zero-point baseline to the neighboring valley points (Fig. 3a).

To this end, eight different preprocessing pipelines (Table 1) were constructed by permuting three design choices: first vs. second derivative (for baseline correction), peak area vs. maximum peak height (for IR signal calculation), and zero vs. valley (for IR peak bounds definition). Scattering corrections were not included in this optimization. This is because the experimental setup minimized scattering effects through careful control of the IR probe position and the homogeneous nature of the reaction mixture. Normalization and scaling, commonly used to equalize the variation of important and less important variables for multivariate model development,¹¹ were also not considered. BayesLOO scores were calculated for M1 with all eight preprocessing pipelines. The combination of “first derivative + peak area + valley” (Table 1 entry 6) showed substantial improvement relative to the original pipeline (entry 3). These results support our hypothesis that incorporating additional spectral information around the peak maximum in the IR signal calculation helps mitigate higher-order effects, thereby improving the performance of the univariate model (M1).

3.6 Model performance assessment with Bayesian posterior predictive checks

Despite the significant improvements achieved after optimizing the preprocessing steps, a substantial BayesLOO gap remained between M1 and M3 (compare 1 entry 6 and 3 dashed line). We therefore sought to determine whether this gap is detrimental to the practical performance of the univariate model, or whether the current version of M1 is already “good enough”. Bayesian analysis provides an additional tool for assessing model quality: posterior predictive checks (PPCs). One application of PPCs is to visually compare the observed data distribution with the distributions of simulated data generated from the model. This comparison reveals the extent to which model predictions deviate from the observed data, and offers a more nuanced assessment of model adequacy than just a single numerical metric. For better visualization, kernel density estimation was used to approximate probability density functions for all distributions.

The results showed that, as expected, the simulated data distributions from M3 (Fig. 5a, green shaded area) closely resembled the observed data distribution (Fig. 5a, black dashed line). Due to pronounced higher-order effects, the simulated data distributions from M1 using the original preprocessing pipeline (Fig. 5a, blue shaded area) failed to reproduce the observed distribution. However, the simulated data from M1 using the optimized preprocessing pipeline (Fig. 5a, orange shaded area), though displaying slightly greater uncertainty, successfully captured all salient features of the observed distribution.

Fig. 5 (a) Visual comparison of posterior predictive checks among three different models. The model-simulated data distribution from M1 with the optimized preprocessing strategy (orange shaded region) accurately captures the main features of the observed data distribution, indicating good predictive performance. (b) Validation of the univariate model with the optimized preprocessing strategy on four unseen experiments. Accurate predictions were obtained for all experiments except experiment 3.

Based on these observations, we concluded that the univariate linear regression model with the optimized preprocessing pipeline could be used for reaction monitoring, and we applied it to four additional test experiments. The results showed that, except for experiment 3, the new model accurately predicted product concentrations, as indicated by the close agreement between observed and predicted values in Fig. 5b. It also produced an average 70% decrease in the root mean square error of prediction (RMSEP) relative to the original model, and yielded an RMSEP comparable to that of M3 (Fig. S7). We attribute the discrepancy observed in experiment 3 to the inaccurate calculation of the IR signal using the optimized preprocessing pipeline. This is caused by the uncertainty in identifying the correct valley positions in its IR spectra. This limitation could be mitigated in future studies by increasing the spectral resolution of the IR measurements (e.g., from eight to four wavenumbers per data point), which would facilitate more reliable identification of valley points.

4 Conclusions

In this work, we presented a novel workflow based on Bayesian modeling to improve the univariate linear regression approach for IR reaction monitoring. A key feature of the workflow is the use of a Bayesian hierarchical linear model during the model comparison stage. This allowed us to distinguish between individual-level and experiment-level higher-order effects. Identifying the latter is particularly valuable, as these effects are more amenable to correction through preprocessing strategies. Applying the workflow to a BTM-catalyzed acylation reaction system, we identified an IR signal in the fingerprint region that exhibited minimal individual-level but substantial experiment-level higher-order effects. A subsequent grid-search optimization uncovered a combination of less common preprocessing steps that significantly reduced the prediction error by more than 50% for the univariate model. An important consideration during calibration model development is that not all applications demand the same level of accuracy. Recognizing when a model is “good enough” can prevent unnecessary optimization efforts. Bayesian posterior predictive checks informed this decision by enabling direct visual comparison between the observed data distribution and model simulated distributions. These checks indicated that the optimized preprocessing pipeline with the univariate model was sufficiently accurate for reaction monitoring. Together, we demonstrated how the proposed workflow provides a systematic way to improve univariate calibration approaches for reaction monitoring.

The use of Bayesian models is central to the proposed workflow, as they allow measurement noise to be incorporated directly into the models and provide uncertainty estimates for all inferred parameters. These features are particularly valuable when analyzing reaction time-course data, where datasets are typically limited in size and affected by measurement noise. With advances in computing power, programming languages, open-source software, and high-quality online documentation, Bayesian approaches have become increasingly accessible to practitioners without formal training in statistics or computer science. As a result, we anticipate that the proposed workflow can be readily adopted and adapted to other applications with little additional effort.

One limitation of the current study is that only a single reaction system was investigated. For new reaction systems, the current preprocessing optimization strategy may not yield meaningful improvements, even when scenario two is identified. In future work, we aim to apply this workflow to a broader range of reaction systems and to incorporate multivariate curve resolution (MCR) analysis into the workflow. MCR can provide more chemically meaningful interpretations of IR spectra than PCA.¹⁰ This added insight will enable the design of more tailored and effective preprocessing strategies, and ultimately enhance the generalizability of the workflow for improving univariate calibration approaches across a wide range of chemical applications.

Conflicts of interest

There are no conflicts to declare.

Data availability

The Python code implementation of our custom-made IR data processing pipeline and Bayesian analysis are available at: https://github.com/jiayu423/Data-driven-IR-calibration-model-investigations. The entire calibration dataset (7 in total, excel format) that contains the IR raw spectra and corresponding HPLC concentration measurements can also be found in the same GitHub repository.

Acknowledgements

Financial support for this work was provided The University of British Columbia, the Canada Foundation for Innovation (CFI-35883, CFI-44843), the Natural Sciences and Engineering Research Council of Canada (NSERC; RGPIN-2021-03168, Discovery Accelerator Supplement), and the Canada First Research Excellence Fund (CFREF2022-00042).

Notes and references

1. J. Liu, Y. Sato, F. Yang, A. J. Kukor and J. E. Hein, Chem.: Methods, 2022, 2, e202200009.
2. A. M. Kearney, S. G. Collins and A. R. Maguire, React. Chem. Eng., 2024, 9, 990–1013.
3. M. C. Deem, J. S. Derasp, T. C. Malig, K. Legard, C. P. Berlinguette and J. E. Hein, Nat. Commun., 2022, 13, 1–11.
4. Y. Sato, J. Liu, A. J. Kukor, J. C. Culhane, J. L. Tucker, D. J. Kucera, B. M. Cochran and J. E. Hein, J. Org. Chem., 2021, 86, 14069–14078.
5. P. Hamilton, M. J. Sanganee, J. P. Graham, T. Hartwig, A. Ironmonger, C. Priestley, L. A. Senior, D. R. Thompson and M. R. Webb, Org. Process Res. Dev., 2014, 19, 236–243.
6. F. St-Jean, K. A. Piechowicz, L. E. Sirois, R. Angelaud and F. Gosselin, Organometallics, 2018, 38, 119–128.
7. B. C. Haas, N.-K. Lim, J. Jermaks, E. Gaster, M. C. Guo, T. C. Malig, J. Werth, H. Zhang, F. D. Toste, F. Gosselin, S. J. Miller and M. S. Sigman, J. Am. Chem. Soc., 2024, 146, 8536–8546.
8. C. Yang, H. Feng and K. Stone, Org. Process Res. Dev., 2021, 25, 507–515.
9. G. Hutchinson, C. D. M. Welsh and J. Bures, J. Org. Chem., 2021, 86, 2012–2016.
10. S. J. Mazivila and J. L. Santos, TrAC, Trends Anal. Chem., 2022, 157, 116698.
11. J. Engel, J. Gerretzen, E. Szymańska, J. J. Jansen, G. Downey, L. Blanchet and L. M. Buydens, TrAC, Trends Anal. Chem., 2013, 50, 96–106.
12. J. Gerretzen, E. Szymańska, J. J. Jansen, J. Bart, H.-J. van Manen, E. R. van den Heuvel and L. M. C. Buydens, Anal. Chem., 2015, 87, 12096–12103.
13. B. J. Shields, J. Stevens, J. Li, M. Parasram, F. Damani, J. I. M. Alvarado, J. M. Janey, R. P. Adams and A. G. Doyle, Nature, 2021, 590, 89–96.
14. F. Hase, L. M. Roch, C. Kreisbeck and A. Aspuru-Guzik, ACS Cent. Sci., 2018, 4, 1134–1145.
15. M. G. Baltussen, J. van de Wiel, C. L. Fernández Regueiro, M. Jakštaitė and W. T. S. Huck, Anal. Chem., 2022, 94, 7311–7318.
16. X. Li, R. Amirmoshiri, C. R. Davis, I. Muthancheri, A. d. Gombert, S. Moayedpour, S. Jager, A. R. Rotheli and Y. Jangjou, Ind. Eng. Chem. Res., 2025, 6825–6837.
17. M. N. Nounou, B. R. Bakshi, P. K. Goel and X. Shen, AIChE J., 2002, 48, 1775–1793.
18. T. Chen and E. Martin, Anal. Chim. Acta, 2009, 631, 13–21.
19. V. B. Birman and X. Li, Org. Lett., 2006, 8, 1351–1354.
20. D. G. Blackmond, Angew. Chem., Int. Ed., 2005, 44, 4302–4320.
21. T. C. Malig, J. D. B. Koenig, H. Situ, N. K. Chehal, P. G. Hultin and J. E. Hein, React. Chem. Eng., 2017, 2, 309–314.
22. M. C. Deem and J. E. Hein, J. Org. Chem., 2023, 88, 1292–1297.
23. https://www.pymc.io/welcome.html.

---