Melodie
Christensen
*ab,
Yuting
Xu
b,
Eugene E.
Kwan
b,
Michael J.
Di Maso
b,
Yining
Ji
b,
Mikhail
Reibarkh
b,
Alexandra C.
Sun
b,
Andy
Liaw
b,
Patrick S.
Fier
b,
Shane
Grosser
b and
Jason E.
Hein
*acd
aDepartment of Chemistry, University of British Columbia, Vancouver, British Columbia V6T 1Z1, Canada. E-mail: jhein@chem.ubc.ca
bDepartment of Process Research and Development, Merck & Co., Inc, Rahway, NJ 07065, USA. E-mail: melodie.christensen@merck.com
cAcceleration Consortium, University of Toronto, Toronto, ON, Canada
dDepartment of Chemistry, University of Bergen, Bergen, Norway
First published on 11th April 2024
Autonomous process optimization (APO) is a technology that has recently found utility in a multitude of process optimization challenges. In contrast to most APO examples in microflow reactor systems, we recently presented a system capable of optimization in high-throughput batch reactor systems. The drawback of APO in a high-throughput batch reactor system is the reliance on reaction sampling at a predetermined static timepoint rather than a dynamic endpoint. Static timepoint sampling can lead to the inconsistent capture of the process performance under each process parameter permutation. This is important because critical process behaviors such as rate acceleration accompanied by decomposition could be missed entirely. To address this drawback, we implemented a dynamic reaction endpoint determination strategy to capture the product purity once the process stream stabilized. We accomplished this through the incorporation of a real-time plateau detection algorithm into the APO workflow to measure and report the product purity at the dynamically determined reaction endpoint. We then applied this strategy to the autonomous optimization of a photobromination reaction towards the synthesis of a pharmaceutically relevant intermediate. In doing so, we not only uncovered process conditions to access the desired monohalogenation product in 85 UPLC area % purity with minimal decomposition risk, but also measured the effect of each parameter on the process performance. Our results highlight the advantage of incorporating dynamic sampling in APO workflows to drive optimization toward a stable and high-performing process.
Optimization algorithms determine the input parameters that result in the minimum or maximum output parameters for an objective function, f(x). An objective function can be as simple as the quadratic function f(x) = x2, but in reality, objective functions in the chemical sciences are seldom this simple. In fact, mapping out the objective function of a chemical process would require the experimental sampling of each possible input parameter permutation. Due to the impracticality of the large number of experiments required by such an approach, the construction of a surrogate model from a subset of experiments is often preferred. As additional samples are evaluated in each optimization round, the model is refined. This approach to optimization has been coined sequential model-based optimization (SMBO, Fig. 1).8 Common examples in the chemical sciences include Stable Noisy Optimization by Branch and Fit (SNOBFIT),2,9 which relies on linear regression, and Bayesian Optimization (BO),10–15 which commonly relies on Gaussian Process (GP) regression.
We recently applied an autonomous Bayesian Optimization strategy to the optimization of a Suzuki–Miyaura cross-coupling process in high-throughput batch reactors.16 High-throughput batch reactors offer an advantage over the microflow reactors, which have historically been utilized in the APO field, in that they provide broader versatility with respect to the physical processes they tolerate (for example, solid–liquid heterogeneous processes). Batch reactors also offer higher throughput via parallelization in order to generate data for surrogate models more efficiently, however, experiments are typically sampled at predetermined static timepoints rather than at dynamic endpoints.17–20 We devised a novel dynamic sampling strategy in which each high-throughput batch experiment could be sampled over time and terminated upon reaction endpoint detection via a real-time plateau detection algorithm.
The importance of dynamic sampling cannot be overstated in the context of APO, particularly in the optimization of complex chemical reactions through techniques such as SMBO. This crucial technique not only allows for the automatic adjustment of the number of samples taken during run-time, saving precious experimental resources, but also provides an accurate means of capturing the process output upon reaction completion. This is particularly useful in reactions with a propensity for decomposition, where dynamic sampling allows each reaction to run its course.
Dynamic sampling proved particularly useful in the optimization of a pharmaceutically relevant radical photobromination process with the propensity for decomposition through the generation of side product 3 (Fig. 2).21,22 We sought to understand the impacts of process conditions on the process performance and to identify conditions that would avoid decomposition. In addition, we sought to understand the process performance impacts of parameter permutations close to the optimum in the least number of experiments. Without dynamic sampling through a plateau detection algorithm, sampling could potentially terminate prior to the reaction stream reaching stability to the point of reaction completion or decomposition, resulting in the inaccurate capture of the process performance. Herein we describe the implementation of a dynamic sampling-driven SMBO strategy to optimize this photobromination process.
Fig. 2 Photobromination process optimized through dynamic sampling-driven autonomous process optimization. |
Fig. 3 High-throughput reaction profiling evaluating NBS in the presence of 22 acid additives, along with replicate additive-free conditions. |
UPLC area % of starting material 1 (blue), product 2 (green) and side product 3 (red) measured over 9 hours in 1.5 hour time intervals. Conditions: 63 μmol 1, 69.3 μmol NBS and 6.3 μmol acid additive in ACN (0.25 M) irradiated with 405 nm LEDs at 60 mW (level 1) intensity for 9 h at 10–30 °C.
Although we had uncovered a suitable starting point for optimization, we were concerned that several experiments reached a conversion plateau prior to full starting material consumption. We suspected that decomposition of the bromination reagent was responsible for conversion plateaus, but did not have a full grasp of the decomposition mechanism. Lack of understanding around a decomposition mechanism posed a significant risk for the optimization. We therefore embarked upon LED-illuminated NMR spectroscopy studies for monitoring photochemical reactions that had been developed in 2019 and implemented in the monitoring of multiple photochemical reactions,23–25 including a recently reported Wohl–Ziegler bromination.21
UPLC analysis proved to be a suitable method for monitoring starting material 1, product 2 and side product 3, but failed to effectively monitor the bromation reagent (NBS) and succinimide levels due to their lack of strong chromophores. NMR proved to be a complementary method for gaining an understanding of the fate of NBS. Two LED NMR experiments were carried out at 15 °C that monitored starting material 1, product 2, dibrominated side product 3, NBS and succinimide concentrations over time. The first experiment was a light–dark study, where the reaction was irradiated for 10 minutes, aged in the dark for 10 minutes, and then irradiation resumed for 100 minutes (Fig. 4a). No changes in reaction species concentrations were observed over the 10 min dark period, indicating the absence of a dark reaction or decomposition pathway. This observation allowed for us to turn off the lights during autonomous optimization sample analyses. Overall mass balance with respect to 1, 2, and 3 remained steady over time and the mass balance with respect to NBS and succinimide also remained steady over time. The second experiment was a reaction monitoring study under constant irradiation. Here, we noted that 20 mol% of the NBS had converted to succinimide prior to initiation and that the reaction plateaued due to full consumption of the NBS prior to starting material consumption (Fig. 4b and c).
In order to further understand the early formation of succinimide, we prepared a representative mixture of starting material 1 and NBS in anhydrous ACN and observed that 15.3 mol% of succinimide had already formed prior to our first NMR measurement (the level remained constant over 20 hours in the dark). In order to rule out the possibility that the source of succinimide was the reagent bottle, we prepared a 70 mM solution of NBS from the same reagent bottle and observed only 0.6 mol% succinimide. Similarly, the addition of H3PO4 to NBS resulted in a measurement of 0.5 mol% succinimide in solution. In contrast, the addition of starting material 1 to NBS resulted in the immediate formation of 18.6 mol% succinimide and, unsurprisingly, the addition of starting material 1 and H3PO4 to NBS resulted in the formation of 18.4 mol% succinimide.
We speculated that pyridazinone 1 might have formed an N-brominated complex in equilibrium, and were unsure whether this species could serve as an effective bromination reagent. In order to provide evidence for the formation of an N-bromopyridazinone species, we carried out NMR internal standard assay measurements of the concentrations of pyridazinone 1 and succinimide in two solutions, the first containing pyridazinone 1 alone and the second containing pyridazinone 1 along with one molar equivalent of NBS. The first assay measured the single component mixture to be 76.0 mM in pyridazinone 1 while the second assay measured the mixture to be 58.7 mM in pyridazinone 1 and 23 mM in succinimide based on their corresponding N–H signals. Reaction plateau variability could be based on changes in the equilibrium distribution between NBS and the proposed N-bromopyridazinone species under varied process conditions. Nonetheless, the LED NMR studies allowed us to determine that the preparation of stock mixtures containing pyridazinone 1 and NBS could potentially lead to plateau variability issues, and thus, all future stock solutions for autonomous optimization were prepared individually for each reaction component.
Fig. 5 Components of a Chemspeed SWING XL system for high-throughput autonomous photochemistry optimization with dynamic sampling. |
A key feature of this system was a novel photoreactor that could maintain a low temperature while accessing all five light intensities (60–385 mW). For this, we worked with an external vendor to develop a novel temperature-controlled reactor (TCR) with cooling channels that decreased the number of wells from 96 to 48. This reactor achieved excellent temperature control at all five light intensity levels (Fig. 5).
Input parameters | Range | Unit |
---|---|---|
Reagent | NBS, DBDMH | Reagent |
Additive | 8 acids | Acid |
Solvent | ACN, DMC | Solvent |
Reagent equivalents | 1.0–1.5 | Molar equivalents |
Additive loading | 1–25 | mol% |
Rxn temperature | 5–35 | °C |
Light intensity stage | 1, 2, 3, 4, 5 | Stage |
Output parameters | Objective | Unit |
---|---|---|
LC area % product | Maximize | LCAP |
Acid | pKa1 (H2O) |
---|---|
Hydrochloric acid (HCl) | −8.0 |
Sulfuric acid (H2SO4) | −3.0 |
2-Picolinic acid | 1.0 |
Phenylphosphonic acid (PPA) | 1.9 |
Phosphoric acid (H3PO4) | 2.1 |
DL-lactic acid | 3.9 |
Acetic acid (HOAc) | 4.8 |
Water (H2O) | 15.7 |
Finally, the search space was focused on a single optimization output, the UPLC area % of monobrominated product 2, measured at 210 nm wavelength. Our goal was to maximize the amount of desired product 2 while minimizing the amounts of starting material 1 and side product 3. Prior experiments revealed that the UPLC area % of product 2 correlated very well with the solution yield of product 2; thus, the analysis was limited to this output parameter measurement for simplicity (see ESI† for details).
The first optimization campaign explored the implementation of a linear regression (LM) model based sequential optimization strategy with a predicted mean acquisition function and was executed for 46 iterations (Fig. 6). The second optimization campaign explored the implementation of a Gaussian process (GP) model based Bayesian optimization (BO) strategy with an expected improvement (EI) acquisition function and was executed for 34 iterations (Fig. 7). The third optimization campaign explored the implementation of a Gaussian process (GP) model based Bayesian optimization strategy (BO) with three alternating acquisition functions, including expected improvement (EI), probability of improvement (PI) and upper confidence bound (UCB). The final campaign was executed for 48 iterations (Fig. 8).
The time course data reveals three general categories of kinetic profiles. In the first kinetic profile, excellent rate acceleration is observed, but this acceleration is accompanied by rapid decomposition after reaching a maximum product level (for example, plot 5 in Fig. 6). In the second kinetic profile, moderate and controlled rate acceleration is observed, resulting in moderate to high product formation (for example, plot 28 in Fig. 7). In the third kinetic profile, the reaction is slow and the plateau is reached at very low product levels (for example, plot 26 in Fig. 8). The ideal process would display the second kinetic profile, reaching a high level of product in a controlled fashion, with a low risk of decomposition. We observed that both BO campaigns converged to optima exhibiting the second kinetic profile category (Fig. 7 and 8), while the LM campaign did not appear to converge at all, but this can be mostly attributed to the optimizer code focusing on random sampling past iteration 21 (Fig. 6). These results, along with the virtual benchmarking studies provided in the ESI section,† indicate the implementation of a Gaussian process (GP) model-based Bayesian optimization approach (BO) as a superior optimization strategy.
Moreover, monitoring the reaction plateau and reporting the reaction outcome at the plateau point resulted in a much more meaningful comparison among the varied process conditions. Take, for example, the cases of iteration 19 and 39 in the first optimization campaign (plots 19 and 39 in Fig. 6, with their 4 minute and final samples colored red). If both iterations were sampled at 4 minutes, it would have appeared that 19 outperformed 39. Instead, sampling both iterations at their plateau points (10 minutes for 19 and 12 minutes for 39) revealed that, in fact, 39 significantly outperformed 19 upon reaction completion, because decomposition was observed in 19 after 4 minutes of reaction.
Although the time course data presented in Fig. 6 through Fig. 8 demonstrated the various kinetic profiles that could be observed in the photobromination reaction under study, additional visualizations that would provide deeper insights around local and global optima were still needed. Given that three categorical parameters were under evaluation, including reagent, additive and solvent, we hypothesized that multiple local optima were likely to exist. This certainly made the optimization more challenging. For multivariate data, we found the combined bar chart, line and scatter plot format to be the most informative data visualization technique (Fig. 9–11).
Visualization of the multivariate data from first LM based optimization campaign demonstrates the importance of reagent, solvent and additive selection (Fig. 9). NBS outperformed DBDMH under a majority of conditions and DMC appeared to yield optimal results in combination with a broader selection of acid additives. The optimizer appeared to focus on HCl and H2SO4, which promoted significant rate acceleration, however the kinetic profiles with these additives aligned with the first category (decomposition after a maximum product level) at higher additive loadings. Visualization of the multivariate data from the second and third BO campaigns answered some of the questions that arose from the first campaign (Fig. 10 and 11). The second campaign also focused on HCl and H2O4 with NBS in DMC, with similar observations around decomposition under high loadings of these two additives. Here, acetic acid with NBS in DMC was sampled in higher detail, and this additive appeared to promote moderate rate acceleration to generate high product levels. The decomposition observed with the lower pKa acids (HCl and H2SO4) at higher loadings did not appear to be an issue with acetic acid. The third campaign shifted the focus from the lower to higher pKa acids, such as phosphoric, lactic, and acetic acid. This is likely because of the implementation of the UCB acquisition function, which was designed to explore unsampled regions of the parameter space, thus, it is not surprising that additional local optima were revealed with the implementation of this acquisition function. Although lactic acid promoted moderate rate acceleration, the reactions were still quite fast, reaching their plateaus within 8 and 20 minutes. What is preferred about lactic acid, from a process chemistry perspective, is the robustness in product purity levels observed across the entire sampled range, with minimal decomposition risk.
The autonomous process optimization experiments ultimately identified two optimal conditions for this process: (1) 1.5 equivalents NBS with 1.0 mol% H2SO4 in DMC, under a light intensity stage of 2 at 5 °C for 8 minutes (1) 1.5 equivalents NBS with 8.5 mol% lactic acid in DMC, under a light intensity stage of 1 at 15 °C for 20 minutes and, both resulting in the generation of 85 UPLC area % of product 2. The algorithmic optimization of a multivariate parameter space in tandem allowed for a broad variety of parameter combinations to be explored, unveiling two local optima, the latter displaying behavior more amenable to large scale processing.
When applied to the combined APO data, RF modeling was used to rank the influence of each parameter on the process outcome, as well as model the partial dependence of continuous parameters such as additive pKa and additive loading on predicted product 2 UPLC area % (Fig. 12). The ranking of parameter impacts aligned with our qualitative observations, where reagent, solvent, additive pKa and additive loading were determined to be most critical, while light intensity stage and temperature were determined to be less critical. The minimal impact of light intensity on reaction performance is not surprising because the Wohl–Ziegler bromination is proposed to proceed through a radical chain mechanism. In radical reactions with low quantum yields, light intensity can be a critical optimization parameter and should be investigated.28 The partial dependence plots of additive pKa and loading on predicted product 2 UPLC area % also aligned very well with our qualitative observations, where optimal performance with lower pKa additives was observed at lower loadings.
As multivariate optimization data sets become more complex, parameter importance modeling will become critical for the interpretation of APO data. The random forest feature importance scores are computed by permuting each feature and calculating the percent increase in mean squared error on the out-of-bag data. Although it can provide valuable insights on which features are more informative for the model prediction, we need to be cautious about the limitations that the feature importance ranking may be biased by many factors, such as overfitting, correlated features, imbalanced data, and categorical variables with more levels. The partial dependence plot provides a straightforward summary of the marginal effect of a feature on the outcome, but this average effect may also be biased in the presence of correlated features. Therefore, it is important to consider multiple model interpretation techniques in conjunction with domain knowledge to gain a comprehensive understanding of the input feature impacts. If resources allow, it is more rigorous to evaluate model generalizability or validate any scientific observation/hypothesis on prospective out-of-sample experimental data. In addition, we could have opted to rank parameter importance through the underlying linear and Gaussian process models, but decided a uniform modeling technique across all optimization campaigns, which leverage different modeling methods, would be more consistent.
Our dynamic sampling strategy was key to the identification of three kinetic profiles associated with the process, which were largely influenced by reagent, solvent, additive, additive pKa, and additive loading, leading to the identification of decomposition-free conditions under high weak acid loadings or low strong acid loadings. Finally, the plateau detection capability allowed for reporting the process outcomes once process stream stability was reached, capturing accurate purity readings. As a bonus, the implementation of random forest (RF) modeling unveiled valuable process insights.
Future optimizations around processes with the propensity for decomposition will be more successful the development of algorithms capable of modeling time course data. The conclusions reached through visual inspection of the reaction profiles proved highly valuable and the incorporation of this valuable data in the automated decision-making process would significantly enhance future algorithmic process optimization efforts.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3sc06884f |
This journal is © The Royal Society of Chemistry 2024 |