Federico
Florit
^{ab},
Anirudh M. K.
Nambiar
^{a},
Christopher P.
Breen
^{c},
Timothy F.
Jamison
^{c} and
Klavs F.
Jensen
*^{a}
^{a}Massachusetts Institute of Technology, Department of Chemical Engineering, Cambridge, MA 02139, USA. E-mail: kfjensen@mit.edu
^{b}Politecnico di Milano, Dipartimento di Chimica, Materiali e Ingegneria Chimica “G. Natta”, Piazza Leonardo da Vinci, 32, 20133, Milano, Italy
^{c}Massachusetts Institute of Technology, Department of Chemistry, Cambridge, MA 02139, USA
First published on 21st September 2021
Batch and continuous reactors both enable exploration of a chemical design space. The former rely on transient experiments, thus experiencing a wide variety of operating conditions over time, whereas the latter are usually operated at steady state and are representative of only one set of conditions. Operating a continuous reactor under dynamic conditions allows more efficient exploration of the underlying reaction space for extraction of kinetics and optimization of performance. We present a methodology to efficiently explore a design space using a tubular flow reactor installed on an automatic platform (equipped with FTIR and HPLC analysis) operated in a transient regime using sinusoidal variations of the parameters. This data-dense method proves to be quicker with respect to steady-state operations because of the larger amount of information collected during a single experiment. A computational analysis provides a simple criterion for the design of dynamic experiments in order for them to be representative of steady-state conditions. The methodology is applied experimentally to the synthesis of a pharmaceutical intermediate via an esterification reaction in the presence of base. In the experiments, up to three parameters (reaction time, base equivalents, and temperature) are changed simultaneously. Proper design of the trajectories in the design space allows verification of the consistency of the results by exploiting the self-crossings within each trajectory and crossings between different trajectories. The experiments further validate the developed criterion for dynamic operations.
One way to reduce the number of experiments required to find an optimal condition is the use of design of experiment (DoE) methods,^{1,2} which provide reaction conditions aimed at describing the variations of some objective function with respect to the optimization parameters while using a small number of experiments. These approaches employ a simple model, the response surface, to represent the objective function (e.g., yield) by parameters that link the objective function to the optimization variables (e.g., reaction temperature, contact time, and reagent concentration).^{3,4} Thus, these methods rely on the response surface providing a good approximation the actual system.
Tubular reactors used in flow chemistry typically have very low dispersion, and thus, operate as plug-flow reactors (PFRs).^{5,6} For such systems the residence time distribution is narrow and represented by the average residence time, while fluid properties (e.g., composition and temperature) are uniform across the section of the reactor. The continuous reactor can then be represented as a sequence of batch reactors (BRs) moving in space, the segregated model.^{5} A batch reactor enters the system at some time and leaves the reactor after a time interval equal to the average residence time, thus the residence time corresponds to the reaction time of the batch reactor. The initial conditions for each BRs correspond to the inlet boundary conditions of the continuous reactor. When the continuous reactor is operated dynamically, i.e., residence time and inlet conditions change in time, the batch reactors experience different reaction times and initial conditions. In contrast, a continuous system at steady state represents only one set of initial conditions (one equivalent batch reactor). Different reaction conditions are attained along the reactor, but it is difficult to obtain information from these intermediate points of the steady system due to geometrical, mechanical, and optical difficulties in accessing the flow path.^{7}
Steady-state conditions are generally used in the DoE-based methods, which discard the information gathered during the change of reaction conditions (the transient behavior of the reactor). The dynamic (transient) regime of a reactor is therefore richer in information provided to the user compared to steady-state methods. The larger amount of information provided by a single dynamic experiment has been exploited in several kinetic studies in which the reaction outcomes were monitored for varying residence times.^{7–9} Hone et al.^{10} have shown that the approach can be extended to account for the effect of axial dispersion, and Aroh and Jensen^{11} have explored effects of two simultaneous variations of parameters with time exemplified by linear variations of flowrate and temperature. Most recently, Taylor et al.^{12} have used linear ramps of flowrate at discrete temperatures combined with machine learning to discriminate among reaction mechanisms. All the aforementioned experimental studies required inline/online analysis of the reactor effluent, such as FTIR,^{8,11} Raman,^{7,9} HPLC,^{10,12} and NMR.^{13}
A tubular reactor with low dispersion operated in a dynamic regime allows exploration of the effects of different variables on reactor performance as they change over time, resulting in a powerful tool for chemical design space exploration and eventually (dynamic) optimization. Wyvratt et al.^{14} used dynamic experiments to explore a design space for a chemical reaction (Knoevenagel condensation) by linearly changing (one at a time) the residence time and the catalyst equivalents with time, using an online FTIR. The parameters were also changed together in a linear orthogonal fashion, as suggested by other researchers,^{11} or using sinusoidal variations for one parameter and linear variations for the other. A response surface was fitted to the obtained dynamic data which are in good agreement with the steady-state results. Haas et al.^{15} have adopted a similar concept to study the chemical design space of a photochemical cycloaddition reaction, where the parameters investigated were reaction time and one reactant concentration. In this case the residence time was changed by using multiple step changes of the inlet flowrate. An online HPLC/MS system was used for the analysis. The concept of transient flow experiments is an extension of batch dynamic experiments, for which the aim is to obtain a trajectory of the optimization parameters in order to optimize some performance function.^{16–18} While batch reactors optimized in such a way are then operated with dynamic variations even after obtaining the optimal conditions, continuous reactors optimized using dynamic experiments are run at constant conditions (those obtained dynamically at the optimal conditions).
Herein we develop experimental design rules for dynamic exploration of the chemical design space for a reaction performed in a continuous reactor. These experimental design rules are model-free, thus not require estimating a response surface or any information regarding the kinetics of the system. Using the synthesis of a pharmaceutical intermediate as a case study, the resulting methods demonstrate the feasibility of using dynamic experiments, specifically sinusoidal variations, to explore chemical reaction space on an automated flow chemistry platform.^{19}
(1) |
(2) |
(3) |
For changing variables that are not pointwise (determined solely at the inlet), but involve the entire reactor (such as the reactor temperature in general or the light intensity/wavelength for photochemical reactions), the BRs traveling inside the continuous system experience a range of such variables in time. If these variables do not change widely over the effective residence/reaction time, then they can be represented by averaged values.
By means of mass balance equations, the concentration in the tubular reactor can be modeled using dimensionless quantities (details in the ESI†) describing how the dimensionless concentration of species i, ĉ_{i} (concentration divided by total initial inlet concentration), changes along the axial dimensionless coordinate, (position divided by reactor length), and over the dimensionless time, (time divided by initial residence time). In order for the dynamic experiment to be representative of steady states, ∂ĉ_{i}/∂ ≈ 0, i.e. the system is in a pseudo-steady state (PSS). The dimensionless concentration is determined by the intrinsic properties of the system (reaction kinetics and reactor properties) and the dynamics of the dimensionless parameters, ŷ_{k} (value of the parameter divided by its initial value), to be explored. ŷ_{k} thus represents how the considered parameter varies with respect to its initial value. Consequently the PSS condition is reached when:
(4) |
(5) |
In this work sinusoidal changes are considered for all variables, which implies:
(6) |
(7) |
Reagent solutions were prepared using volumetric glassware and N,N-dimethylacetamide (DMA) as the solvent, except for N,N-diisopropylethylamine which was used neat. Infrared spectroscopy data was collected using a Mettler Toledo FlowIR instrument in 15 second sampling intervals and processed using the iC IR 7.0 software by subtraction of the pure solvent spectrum. Calibration of 3 by was performed by injecting standard solutions into the flow cell such that the observed signal heights correlated linearly with concentration. A single point baseline with 1834 cm^{−1} as reference was applied to the in situ IR data. The height of the newly formed ester bond in 3 was observed in the range of 1783–1746 cm^{−1}. Concentration of 3 was then computed by referencing the calibration curve.
High-performance liquid chromatography data was acquired using an Agilent 1260 Infinity series HPLC instrument equipped with an Agilent Poroshell 120 SB-C18 2.7 μm (3.0 × 50 mm) column (Agilent part # 689975-302T). The analytical method used in this study utilised a water (+0.1 v/v% formic acid) and acetonitrile (+0.1 v/v% formic acid) gradient. From 0–3.3 minutes, acetonitrile was increased from 5 to 100 v/v% followed by a hold at 100 v/v% acetonitrile from 3.3–5 minutes. At 5 minutes, acetonitrile was dropped back down to 5 v/v% and the column was re-equilibrated from 5–8 minutes. Starting material 1 and product 3 eluted at 1.7 and 2 minutes, respectively.
The reaction was run in a 1 mL tubular reactor module of a plug-and-play automated flow chemistry platform previously reported by our lab^{19} (Fig. 1B and C). Inlet streams were pumped into the reactor using positive displacement pumps (VICI Valco Instruments, model M6HP, part # CP2-4111-DHP) and all lines had check valves (IDEX Health & Science, part # CV-3301). The reactor module contains a 0.125 mL pre-mixing chamber made of blow formed perfluoroalkoxy alkane (PFA) polymer film equipped with a magnetic stir bar. The mixed streams flow into a 1m tubular reactor (PFA tubing, 1/32′′ ID, 1/16′′ OD) that is pressed into channels in a heated aluminum housing which is kept at a set temperature with controlled electrical heating. The reactor effluent then passes through two inline analytical modules in series: an FTIR module connected to the inline Mettler Toledo FlowIR and an HPLC module connected to an online sample injection valve (VICI Valco Instruments, part # EUHB-CI6WHC.06) that injects a 60 nL sample of the reaction mixture into the Agilent HPLC every 10 minutes. Finally, the outlet module directs the effluent to a back-pressure regulator (Zaiput Flow Technologies, part # BPR-10) and then to waste. A graphical user interface built using Python's Tkinter package served three purposes: (i) interfacing with the reactor, pumps, and injection valve using serial communication, (ii) allowing the user to enter the parameters in eqn (6) that fully define the dynamic trajectory for each variable, and (iii) automatically generating flow rate and temperature trajectories and dynamically updating setpoints once the experiment is started.
The variables to be explored were the residence time inside the reactor, the mole equivalents of base with respect to 1, and temperature. The design space studied was: instantaneous residence time, τ_{I}, between 0.5 and 10 min, base equivalents (ratio between Hünig's base and 1 inlet concentration) between 1:0 and 1:3, and temperature, T, between 40 and 100 °C. 1 and 2 were always fed in stoichiometric quantities. Base equivalents were lowered by decreasing the base flow rate and increasing the DMA diluent flow rate by the same amount, allowing us to modulate base equivalents without altering the concentrations of 1 and 2.
Initially, two runs with constant temperature were performed by changing residence time and base equivalents. Afterwards, a run changing all variables at the same time was performed. All variations were sinusoidal in time as from eqn (6) and the parameters of average value, variation amplitude, and period are reported in Table 1, which also reports the experiment duration, t_{exp}. In order to ensure that the reaction mixture remained a liquid, the back-pressure regulator was set to 50 psi. Changes in residence time and base equivalents over time were realized by manipulating the flow rate of reactants and solvent with time.
Run | t _{exp} [min] | y _{ k } | y ^{0}_{ k } | δ _{ y k } | T _{ y k } [min] | K _{ y k } |
---|---|---|---|---|---|---|
1 | 180 | τ _{I} | 0.5 min | 19 | 360 | 0.17 |
B.E. | 1.15 | 0.13 | 60 | 0.0068 | ||
T | 333.15 K | 0 | — | 0 | ||
2 | 180 | τ _{I} | 0.5 min | 19 | 360 | 0.17 |
B.E. | 1.15 | 0.13 | 60 | 0.0068 | ||
T | 353.15 K | 0 | — | 0 | ||
3 | 300 | τ _{I} | 5.25 min | −0.9048 | 75 | 0.40 |
B.E. | 1.15 | 0.13 | 60 | 0.072 | ||
T | 313.15 K | 0.1916 | 600 | 0.011 |
The system was primed with solvent, brought to the initial reaction temperature, and a blank IR spectrum was acquired. A steady state condition was obtained with the parameters computed at t = 0, by waiting a time equal to three initial residence times (3τ_{I}(0)), and verifying the constant output through the IR spectra. Afterwards, the dynamic (transient) experiment was started for the duration reported in Table 1. The FTIR and HPLC data were collected starting from t = 0 at constant time intervals. Afterwards a new steady state condition was reached at the final values of the explored parameters. Further details on the experimental procedure are reported in the ESI.†
Sinusoidal variations were chosen in order to scan the parameters space in a smooth way allowing many combination of variables (high, low, intermediate values). Additionally, by adopting proper values of T_{yk}, one can design the variations in such a way that the same conditions (the same parameters) are obtained at different times (creating a crossing in the parameters space). If the same value of the objective function is obtained at the crossing, then the value obtained is path-independent (as it is approached from different directions in the parameters space) and this can only happen at a steady state. Thus the check of the objective function at a parameter crossing is an a posteriori check of the proper design of dynamic experiments.
The CDF becomes skewed towards larger values of K_{T}, while it broadens for smaller values of K_{T}. This means, as expected, that slower variations (small K_{yk}) lead to results of the dynamic experiment that coincide with steady states. The threshold value can be selected by choosing a limiting value of 2σ_{Δ} and a value of the CDF K_{yk}. The limit value of 2σ_{Δ} is selected to be equal to 0.01, while the value of the CDF to 0.90, meaning that guarantees yield deviations (between dynamic and steady simulations) below 1% in at least 90% of the cases possible in practice. As a consequence, the values of to be used in eqn (7) for residence time, concentration, and temperature can be computed for sinusoidal variations:
K^{max}_{τ} = 0.2, K^{max}_{c} = 0.2, K^{max}_{T} = 0.04 | (8) |
Fig. 3 reports the time values of the studied parameters of run 1 from Table 1. The instantaneous residence time was varied with half a sinusoidal oscillation, while base equivalents with three full oscillations (see Fig. 3a), consequently the flowrate of each reactant and solvent changes non-sinusoidally in time (Fig. 3b). The yield measured with time is reported in Fig. 3c, which outlines the effect of both residence time and base equivalents variations.
Fig. 3 Time variation of the input and output parameters for run 1 of Table 1: (a) instantaneous residence time and base equivalents; (b) species flowrate; (c) measured yield of the desired product. |
Fig. 4a reports trajectory in design space. The space is bi-dimensional as only two parameters were changed sinusoidally in time, thus yield of 3 can be easily observed in Fig. 4b as a function of the two input parameters. The use of sinusoidal variations allowed two crossings in the parameter space (black circles in Fig. 4): the yield at these locations is the same on both path portions (as seen from the value of the yields), therefore all the points plotted represent steady-states of the system at the corresponding values of τ and base equivalents.
Fig. 4 Results for run 1 of Table 1: (a) design space; (b) yield as a function of the design space parameters (circles: FTIR; squares: HPLC). The crossings of the trajectory are circled. |
The best reaction conditions obtained, for this particular temperature, are located at the extreme of the design space, namely at large τ and base equivalents. The large number of data collected via IR allowed checking the progress of the dynamic experiment and the HPLC data confirmed the location of the best conditions. It must be noted that the IR yields were systematically higher than the HPLC yields by approximately 5% due to overlap between the peaks of reactant 1 and product 3 in the IR spectrum which increased the product height. Nevertheless, the same yield trends were observed in both the IR and HPLC data. While IR enabled densely mapping out the design space along the trajectory, HPLC was a crucial complementary tool for providing the ground truth yield.
Fig. 5 reports the time behavior of the parameters of run 3 from Table 1. The instantaneous residence time undergoes four full oscillations, while base equivalents five (Fig. 5a), thus again, the flowrates do not follow sinusoidal variations (Fig. 5b). Temperature is varied by half an oscillation (Fig. 5c). The collected yield is depicted in Fig. 5c, where larger variations are observed mainly as a consequence of the residence time variations (lower residence time correspond to lower yields). The phase space is three-dimensional, thus it can be plotted in Fig. 6 as a 3D graph where yield can be reported as color of the various points. Note that the temperature reported is the integral average of the temperatures observed by the fluid during the residence time. Again, the large number of IR data collected allows to find a good indication of the optimal values of the parameters for this synthesis and the location of the best conditions studied is at the extrema of the analyzed space, namely, large residence time, base equivalents and temperature.
Fig. 5 Time variation of the input and output parameters for run 3 of Table 1: (a) instantaneous residence time and base equivalents; (b) species flowrate; (c) measured temperature; (d) measured yield of the desired product. |
Fig. 6 Results for run 3 of Table 1: design space where yield is represented by the color (circles: FTIR; squares: HPLC). This trajectory has no crossings. |
The sinusoidal variations allowed efficient exploration of the parameters space. Variations were chosen in order to collect more data closer to the expected location of the optimal conditions (in fact the region at low temperature and residence time was excluded), as inferred from runs 1 and 2. In the studied case, yield increases monotonically at increasing values of the three parameters in the analyzed space, but base equivalents do not change yield significantly at large values of residence time and temperature.
Results from the three experiments are combined in Fig. 7. The three experiments cross each other at three locations (marked on the figure), where the measured yield differs by approximately 1%, confirming that the dynamic experiment is able to describe steady states (as parameter crossings lead to the same result). Therefore, an a posteriori check of the choice of the dynamic variations can also be performed between intersecting experiments with different input parameters.
Fig. 7 Design space where yield is represented by the color (data from FTIR) where all experiments are reported. The crossings between different experiments are circled. |
The developed methodology reduces dramatically the experimental effort and especially time to explore the design space. In steady-state experimental campaign, the same amount of information (same number of data-points for all three experiments) would have been obtained in approximately 270 h of experiments (considering just one residence time to reach steady-state and no cleaning times) compared to 11 h for the dynamic approach. In practice, at least three residence times are needed to ensure steady state and flushing procedures might also be performed in-between experiments, further increasing the experimental time. Thus, the dynamic experiments with the proposed methodology provide more than 75 times information than steady-state methods in the same time period. Finally, the task of manually preparing the exploration trajectory in the design space takes approximately the same time as manually preparing a DoE, thus becoming a negligible time compared to experimental time (both methods could be automatized). An additional advantage is that the trajectory designed in one design space can be easily transposed to any other design space with the same dimensions, and the trajectory will have roughly the same shape of the original one.
It is worth noting that the parameters chosen for experiments 1 and 2 lead to values of K_{yk} (see Table 1) satisfying the criteria obtained numerically, namely eqn (8). This confirms that slow variations satisfying eqn (7) lead to a dynamic response that is representative of steady states with very small deviations in yield. Run 3 does not satisfy the criterion for τ, nonetheless the experiment led to acceptable results (dynamic experiment coincident with steady states). The criteria developed are quite stringent as they allow a maximum yield deviation of 1%. If instead a 2% yield deviation is deemed acceptable, the values to be used in eqn (7) are:
K^{max}_{τ} = 0.6, K^{max}_{c} = 0.5, K^{max}_{T} = 0.09 | (9) |
The method was validated experimentally by studying the yield of an esterification reaction between a carboxylic acid and an alkyl halide. Small deviations (below 1%) from the steady-state yield of the desired product were obtained with the dynamic experiment if the proposed criteria are followed. Notably, no detailed information on the system studied (such as the kinetic mechanism and the physical–chemical properties) was required to apply the proposed method. The data-rich exploration of the design space further enabled identification of the best reaction conditions with a great reduction of experimental time in comparison with methods using steady states. This method could become a useful tool for automated optimization by allowing for efficient collection of initialization data. The dynamic approach could also be advantageous in discerning between reagent candidates for a reaction, or for extracting chemical kinetics.
Footnote |
† Electronic supplementary information (ESI) available: Model and experimental details. See DOI: 10.1039/d1re00350j |
This journal is © The Royal Society of Chemistry 2021 |