Design Space for Polymorphic Co-crystallization: Incorporating

Jun 17, 2014 - ... and Research), 1 Pesek Road, Jurong Island, Singapore 627833 ... National University of Singapore, 10 Kent Ridge Crescent, Singapor...
1 downloads 0 Views 2MB Size
Article pubs.acs.org/crystal

Design Space for Polymorphic Co-crystallization: Incorporating Process Model Uncertainty and Operational Variability Zai Qun Yu,*,† Pui Shan Chow,† and Reginald B. H. Tan*,†,‡ †

Institute of Chemical & Engineering Sciences, A*STAR (Agency for Science, Technology and Research), 1 Pesek Road, Jurong Island, Singapore 627833 ‡ Department of Chemical & Biomolecular Engineering, National University of Singapore, 10 Kent Ridge Crescent, Singapore 119260 ABSTRACT: A methodology has been developed to determine the design space of caffeine−glutaric acid co-crystallization for process development and scale-up. Polymorphic purity is chosen as a critical quality attribute, and its functional relationship with five process parameters is described by a first-principles process model and an empirical supersaturation threshold. After model parameters are estimated and the supersaturation threshold is obtained from a few experiments in a 1 L crystallizer, the risks associated with model parameter uncertainty and operational variability are quantified using Monte Carlo modeling. These risks, while assessed to be significant, could be effectively mitigated by sufficient aging after seeding. Operating ranges of individual process parameters that are robust against operational variability are then modeled sequentially and verified experimentally by scale-up to a 10 L crystallizer.

1. INTRODUCTION First-principles process models describing conservation laws and rate phenomena have been increasingly utilized in design space determination for pharmaceutical processes, e.g., drying,1 coating,2 chemical synthesis,3−5 reactive distillation,6 etc. They capture the functional relationships between quality attributes of intermediates or products, process parameters, and input variables in the form of ordinary differential equations. The prior mechanistic knowledge embedded in them allows efficient designs of experiments and prospectively setting of process parameters for optimal outcome. Consequently, much time and resources can be saved by integrating these models into process development. As an integral unit of pharmaceutical manufacture, crystallization processes define the purity, polymorphism, morphology, particle size distribution, and other solid-state properties of intermediates and products. Crystallization development usually relies on factorial design and empirical models to correlate quality attributes with process parameters and input variables.7 At the same time, advances in process modeling and optimization,8−10 model regression,10,11 and process analytical technology (PAT)12−14 have laid the foundation for a transition to model-assisted development methodologies. In a previous study, Yu et al. demonstrated an efficient procedure to capitalize on first-principles process modeling and PAT in defining the design space of seeded cooling crystallization on a bench scale.15 The design space was totally model-based, and a procedure to map it into operating ranges of individual process parameters for scale-up is yet to be © 2014 American Chemical Society

established. For the time being, there lacks clear guidance from regulatory authorities concerning model verification. Model maintenance and execution may not be easy to implement in routine operation. Therefore, it is of practical interest to map model-based design space into fixed operating ranges of individual process parameters that are easy to implement.16,17 Two factors must be accounted for in design space mapping. One is uncertainty in model parameters and the other operational variability of process parameters upon scale-up. The former can be represented in the covariance matrix resulting from nonlinear model regression using standard algorithms. The latter is represented by a standard deviation from its set point when systematic variation is absent. It can be estimated from historical data of related equipment and facilities.18 In this study, the same procedure as that described in Yu et al.15 will be extended to co-crystallization that makes use of first-principles process modeling. Caffeine−glutaric acid− acetonitrile is chosen as the model system. Caffeine and glutaric acid can form a 1:1 co-crystal in two polymorphs when their molar ratio is within a specific range.19 Form I is acicular in shape and metastable, whereas Form II is prismatic and stable. Form II is the desired product. For this study, polymorphic purity is chosen as the critical quality attribute. Five input variables and process parameters have impacts on supersaturation dynamics and, in turn, on polymorphism of coReceived: April 21, 2014 Revised: May 27, 2014 Published: June 17, 2014 3949

dx.doi.org/10.1021/cg500547m | Cryst. Growth Des. 2014, 14, 3949−3957

Crystal Growth & Design

Article

moment transformation.8 μj denotes the jth moment of PSD. mT is related to μ3 by mT = ρS φv μ3 (4)

crystals, including cooling profile, seed loading, seeding temperature, seed particle size distribution (PSD), and starting concentration. After the design space is determined on a bench scale, the impacts of uncertainty in model parameters and operational variability in some process parameters will be evaluated using Monte Carlo simulations that were adopted by Figueroa et al. in design space delineation for a reactive distillation.20 Model parameter and process parameters are assumed to obey normal distributions. Then, fixed operating ranges of individual process parameters will be modeled sequentially via Monte Carlo simulations as well. Following that, a scale-up experiment will be conducted for one combination of process parameters in the resulting operating ranges for verification purpose. It is worthwhile to mention that there exist different methods other than Monte Carlo simulation to address uncertainty in model parameter and operational variability. For example, Hermanto et al.21 and Nagy et al.22 addressed the impacts of model parameter uncertainty in robust control of crystallization based on worst-case performance. Simone et al.23 and Yu et al.24 dealt with operational variability by implementing feedback control of solute concentration. As there are difficulties to overcome before applying feedback control of solute concentration to commercial production, this study provides a universal and straightforward alternative method to handle uncertainty and variability.

where ρS denotes crystal density, equal to 1486 kg/m3. φv denotes volume shape factor, with a measured value of 0.4003. GL stands for linear crystal growth rate, and it is expressed empirically as G L = k GΔc g

where Δc is supersaturation, defined as the difference between actual concentration, c, and solubility, c*. Equation 3 defines a temperature profile for cooling crystallization. T is slurry temperature with p(t) being a predetermined temperature profile. Temperature dependence of solubility, c*, is given as c* = 0.0068 exp(0.0777T )

dμj dt

= jμj − 1G L , j = 1, 2, 3

dT = p(t ) dt

(6)

Nucleation is negligible when relevant process factors are coordinated in such a way that supersaturation is maintained at relatively low levels. Other assumptions for this process model include ideal mixing, size-independent crystal growth rate, no particle agglomeration and breakage, and uniform volume shape factor. As a result, there are two model parameters, kG and g, to be estimated. Online measurements of solute concentration afford adequate information for regression of these two parameters. However, solid-side information is necessary for kinetics study when particle creating events are occurring, such as nucleation and agglomeration. Curve fittings of online solute concentration measurements were performed with a Matlab routine, lsqnonlin, which adopts standard algorithms for nonlinear regression. The standard methods make a few assumptions about errors, such as being uncorrelated, following a normal distribution, constant variance, etc. Parameter estimates are also assumed to obey normal distributions with their uncertainty being contained in a covariance matrix. In comparison, alternative regression methods such as those based on Bayesian statistics can be used to quantify parameter uncertainty with fewer assumptions.26

2. A MATHEMATICAL MODEL OF CAFFEINE−GLUTARIC ACID CO-CRYSTALLIZATION It has been revealed that polymorphic purity is determined by supersaturation level, and there exists a supersaturation threshold below which the stable form can be obtained consistently in seeded co-crystallization.24,25 Therefore, control over polymorphic purity can be translated into control over supersaturation. The five process factors, including cooling profile, seed loading, seeding temperature, seed particle size distribution, and starting concentration, can be accounted for in the following ordinary differential equation set that describes supersaturation dynamics when nucleation is negligible: dmT d(Mc) + =0 dt dt

(5)

3. EXPERIMENTAL SECTION

(1)

3.1. Experimental Setup. Bench-scale experiments were conducted in a 1 L jacketed crystallizer with a flat bottom, as shown in Figure 1. The crystallizer is 110 mm in inner diameter, and there are four built-in baffles. It is equipped with an overhead stirrer and a marine-type impeller with a diameter of 42 mm. The stirrer rotated at 500 rpm for all bench-scale experiments, and it was found that solids could be well suspended and there was no significant breakage of cocrystals. Crystallization temperature was controlled by a circulator (Julabo FP50). For the scale-up experiment, a 10 L crystallizer was used with a round bottom and an inner diameter of 240 mm. A marine-type impeller as in the bench-scale experiment was in place for agitation (96 mm in diameter). The power number of the marine-type impeller was 0.34.27 The stirrer speed was set at 340 rpm based on the criterion of equivalent power input per unit volume. The temperature was controlled by another circulator (Julabo LH46). Attenuated total reflection Fourier transform infrared spectroscopy (ATR-FTIR) was used to measure solute concentration during crystallization (Nicolet 4700 with an ATR probe made by Axiom Analytical INC). The metastable form of the caffeine−glutaric acid cocrystal has a distinctive needle-like habit, and hence, its appearance can be detected by particle vision measurement (PVM, V700S-5-C22-K, Mettler-Toledo).

(2) (3)

Equation 1 represents the mass balance between liquid and solid phases. t is time and mT, solid mass in slurry. c is dimensionless solute concentration defined as the mass of solute per unit mass of solvent. M is the mass of solvent. A definition of caffeine−glutaric co-crystal concentration and solubility was given in a previous study.24 In co-crystallization, the concentrations of caffeine and glutaric acid are changing in accordance with stoichiometry. Therefore, it can be imagined that molecular pairs of caffeine and glutaric acid exist in liquid phase, and it is these pairs that crystallize out as a co-crystal. With this definition, a co-crystal can be treated as a single compound as in crystallization of single components. Equation 2 describes the dynamics of PSD moments, and it is derived from a corresponding population balance equation by 3950

dx.doi.org/10.1021/cg500547m | Cryst. Growth Des. 2014, 14, 3949−3957

Crystal Growth & Design

Article

Starting solutions were saturated at 35 °C. After appropriate amounts of caffeine, glutaric acid, and acetonitrile were added into the crystallizer, the temperature was raised to 38 °C to dissolve all solids under agitation. After that, the temperature was cooled in a preset ramp to 10 °C. Seeds were applied at a chosen temperature. Cocrystals were harvested at 10 °C and dried in air at ambient temperature.

4. RESULTS AND DISCUSSION 4.1. Experiment Summary. Six bench-scale experiments and one scale-up experiment, as shown in Table 1, were conducted. The purpose of bench-scale experiments was to probe the supersaturation threshold for the emergence of the metastable form, and to estimate model parameters. The scaleup experiment was for design space verification. The solvent mass varied between 383.8 and 413.4 g while the initial concentration of co-crystal was kept constant in benchscale experiments. The seed loading and cooling profile was adjusted from experiment to experiment based on the results of previous experiments and the target of the current experiment. Seeds were applied between 34 and 35 °C to avoid dissolution. There was no aging after seeding in these six experiments. After the fourth experiment, the supersaturation threshold was determined and model parameters were estimated by curve fittings. The settings for the fifth and sixth experiments were then guided by a series of process modeling with the constraint of peak supersaturation being lower than the supersaturation threshold. Supersaturation trajectories in the six experiments are shown in Figure 3. Supersaturation peaked in the first half of each

Figure 1. Experimental setup. 3.2. Chemicals and Procedures. Anhydrous caffeine of 99% purity was supplied by Fluka, glutaric acid of 99% purity was supplied by Alfa Aesar, and HPLC grade acetonitrile was acquired from Fisher Chemical. Two batches of stable form seeds (Form II) had been prepared from unseeded crystallization (labeled as Seed 1 and Seed 2). Metastable co-crystals (Form I) nucleated first and then converted to stable form in unseeded crystallization. The undersize cumulative PSDs of seeds are presented in Figure 2

Figure 2. Undersize cumulative mass percentage of seeds. Equation 4 is used to obtain various moments of PSD from sieve analysis. Seed 1: μ0 = 3.33 × 106/g; μ1 = 2.13 × 102 m/g; μ2 = 1.72 × 10−2 m2/g; μ3 = 1.68 × 10−6 m3/g. Seed 2: μ0 = 1.30 × 106/g; μ1 = 1.25 × 102 m/g; μ2 = 1.39 × 10−2 m2/g; μ3 = 1.68 × 10−6 m3/g.

Figure 3. Supersaturation trajectories in experiments 1−6.

Table 1. Summary of Bench-Scale Experimentsa

a b

no.

M, g

C0, g/g

Ms, wt %

1 2 3 4 5 6 7

383.8 410.3 387.4 413.4 399.8 411.1 5904

0.1075 0.1070 0.1054 0.1049 0.1053 0.1057 0.1077

5.6 5.3 2.8 1.4 4.1 4.0 5.6

seed PSD Seed Seed Seed Seed Seed Seed Seed

1 1 1 1 1 2 2

Ts, °C

linear cooling time, min

S, g/g

34.0 34.0 34.5 34.8 34.5 34.8 34.6

120 80 82 124 82 94b 123

0.012 0.015 0.019 0.022 0.016 0.014 0.011

polymorph Form Form Form Form Form Form Form

II II I I II II II

M stands for solvent mass, C0 for initial concentration of solute, Ms for seed loading, Ts for seeding temperature, and S for peak supersaturation. Solution was cooled to 28 °C at 0.2 °C/min and then cooled to 10 °C at 0.3 °C/min. 3951

dx.doi.org/10.1021/cg500547m | Cryst. Growth Des. 2014, 14, 3949−3957

Crystal Growth & Design

Article

the fifth experiment. The model parameters were updated by including online data in experiment 5. Now, the supersaturation trajectories can be simulated for any process parameter combinations with the mathematical model (eqs 1−6). All of those combinations resulting in a peak supersaturation lower than the threshold (0.015 g/g) constitute the bench-scale design space of the co-crystallization, if the model is assumed to be deterministic. One of those combinations was verified in experiment 6 wherein the temperature was cooled from 35 to 28 °C at a cooling rate of 0.2 °C/min. Then, another cooling rate of 0.3 °C/min was used until the final temperature of 10 °C (total cooling time was 94 min). Simulated peak supersaturation was 0.013 g/g, quite close to the actual value of 0.014 g/g, as shown in Table 1. Co-crystals of the metastable form were not found in PVM images. 4.3. Model Parameter Estimation and Residuals Analysis. Curve fittings for experiments 1, 2, 5, and 6 are shown in Figure 5. In general, the process model agrees very well with experimental data. Estimates of model parameters are updated as kG × 104 = 1.09 m/min, g = 0.700. The variances of kG × 104 and g are 0.00524 and 0.000202, respectively. Their covariance is 0.00102. In model regression, kG was multiplied by a scaling factor of 104 to avoid ill-conditioned matrices. kG × 104 is of the same order of magnitude with g. The residual scatter plot of the mathematical model and their histogram are shown in Figure 6a,b, respectively. Figure 6a indicates that the residuals are correlated to some degree. The magnitude of residuals is larger in the first half of experiments 5 and 6 than in the second half. Figure 6b discloses that the residuals are not of a normal distribution around a zero mean. A mode appears at 0.004 g/g, and the distribution shape is slightly skewed. These observations suggest that the assumptions in standard methods for nonlinear regression are not violated seriously. The estimates of model parameters and their covariance matrix will be used in process modeling later on. The contour plot of probability density function of kG and g is illustrated in Figure 7. Contours on the kG−g plane are narrow ellipses, suggesting that there is strong correlation between these two parameters. The impacts of model parameter uncertainty are evaluated using Monte Carlo simulation; 103 combinations of kG and g randomly sampled are plugged into the process model, resulting in 10 3 supersaturation trajectories (green curves in Figure 8). Peak values of some green curves exceed the supersaturation threshold (horizontal bold red line). The blue curve is from a simulation without model parameter uncertainty. It can be seen that uncertainty in model parameters leads to diffusion of supersaturation trajectories and poses a risk to crystallizations whose peak supersaturation is already close to the threshold. 4.4. Operational Variability in Process Parameters and Its Impacts. As mentioned earlier, five process parameters are considered in caffeine−glutaric acid co-crystallization, including starting concentration, cooling profile, seed loading, seeding temperature, and seed PSD. Their operational variability differs. Variability in actual starting concentration is described by a normal distribution around a mean value of 0.105 g/g with a standard deviation of 0.0025 g/g (±2.5%). The combined impacts of variability in starting concentration and model parameter uncertainty are assessed by Monte Carlo simulation, as shown in Figure 9. The blue curve is from a process simulation with neither model parameter uncertainty nor

batch. By the end of experiments 1, 2, 5, and 6, supersaturation dropped to zero and it indicated that solute concentration was approaching equilibrium. However, supersaturation measured by ATR-FTIR did not drop to zero by the end of experiments 3 and 4 due to encrustation of the ATR probe at relatively high supersaturation. 4.2. Supersaturation Threshold for Metastable Form. In experiment 1, a seed loading of 5.6 wt % and a linear cooling time of 120 min (corresponding to a cooling rate of 0.2 °C/ min) were chosen based on experience. Peak supersaturation was 0.012 g/g. No co-crystals of the metastable form were found in PVM images in experiment 1. In order to push up peak supersaturation and to probe the supersaturation threshold for the metastable form, the linear cooling time was shortened to 80 min (corresponding to a cooling rate of 0.3 °C/min) in experiment 2 while seed loading was 5.3 wt %. Peak supersaturation rose from 0.012 to 0.015 g/g. Co-crystals of the metastable form did not appear either in experiment 2. To further stress crystallization condition, seed loading was cut down to 2.8 wt % and a cooling rate of 0.3 °C/min was employed in experiment 3. The metastable form appeared, and the peak supersaturation reached 0.019 g/g. Stressed conditions were also used in experiment 4 by cutting seed loading to 1.4 wt % and the linear cooling time was 124 min. The metastable form appeared again, and peak supersaturation was 0.022 g/g. The correlation between the appearance of metastable form and the supersaturation level was confirmed. On the basis of the results of the first four experiments, it was decided that the supersaturation threshold can be set at 0.015 g/g (the peak supersaturation in experiment 2). Model parameters were estimated from online measurements of solute concentration in experiments 1 and 2. Nucleation was negligible in these two experiments as undersize cumulative mass percentage of products in Figure 4 depicts. The mass

Figure 4. Undersize cumulative mass percentage of co-crystals from experiments 1 and 2.

fraction of co-crystals smaller than 250 μm was less than 0.5 wt % in both experiments. Data from experiments 3 and 4 were not included for parameter estimation because nucleation of the metastable form occurred therein. With estimates of model parameters, experiment 5 was conducted to test the threshold. Process simulation showed that the settings in experiment 5 would lead to a supersaturation of 0.015 g/g, while actual peak supersaturation was 0.016 g/g, slightly higher than simulated results. Co-crystals of the metastable form did not emerge in 3952

dx.doi.org/10.1021/cg500547m | Cryst. Growth Des. 2014, 14, 3949−3957

Crystal Growth & Design

Article

Figure 5. Curve fittings for model regression using measured solute concentration in four experiments: (a) experiment 1, (b) experiment 2, (c) experiment 5, and (d) experiment 6.

with a standard deviation of 1 °C. Likely, the standard deviation of the actual seed loading from its set point is assumed to be ±2.5%. The impacts of operational variability in these two process parameters are not depicted here individually for brevity. Their variability will be taken into account in the same way as that of starting concentration. In summary, operational variability in starting concentration, seeding temperature, and seed loading will be considered in this study, while variability in seed PSD and time (cooling profile) is ignored. 4.5. Operating Ranges of Individual Process Parameters. The choice of starting concentration is dictated by technical and economic considerations. It usually has a single set point. In this study, it is set at 0.105 g/g. As mentioned in the previous section, seed PSD can be specified by D-values. The D-values of Seed 2 will be used as the upper limit for seeds PSD at larger scale. Simulations have revealed that an aging time not less than 30 min after seeding can effectively consume the initial supersaturation and thus lower the risk of exceeding the supersaturation threshold. For example, with an aging time of 30 min and the settings of other process parameters being kept the same as those in Figure 9, 97.3% of simulated supersaturation trajectories in Figure 10 have a peak value below the threshold when model parameter uncertainty and variability in starting concentration are taken into account. In contrast, the probability is 63.4% when there is no aging period, as shown

variability in starting concentration (the settings of other process parameters are the same as those in Figure 8). Green curves are from 105 simulations wherein model parameters and starting concentration being randomly sampled in their respective distributions. Compared with Figure 8, it can be seen that the risk of exceeding the supersaturation threshold becomes much higher by the addition of variability in starting concentration. Only 63.4% of the green curves have a peak value below the supersaturation threshold. A cooling profile for crystallization often contains an aging period after seeding, one or more cooling ramps and hold periods between them if necessary, and another aging period at the final temperature. For illustrative purposes, three-stage profiles are adopted herein: an aging period after seeding, a linear cooling period, and another aging period at the final temperature. To simplify the situation, operational variability in the time length of each period is ignored. Seed PSD is usually specified by D-values. For example, D50 = 110 μm means that 50% of a sample’s mass is smaller than 110 μm. In-process control during seed preparation can ensure the actual D-values of seeds be equal to or lower than a certain upper limit. It is well-known that smaller seeds lead to lower supersaturation peak when seed loading is kept constant. Therefore, the risk associated with operational variability in seed PSD is trivial. Operational variability in actual seeding temperature is also assumed to follow a normal distribution around its set point 3953

dx.doi.org/10.1021/cg500547m | Cryst. Growth Des. 2014, 14, 3949−3957

Crystal Growth & Design

Article

Figure 8. Probable supersaturation trajectories in co-crystallization. The blue curve is simulated supersaturation trajectory without considering uncertainty in model parameters. Green curves are obtained from 1000 process simulations wherein the values of kG and g are randomly sampled in the distribution depicted in Figure 6. The bold red horizontal line stands for the supersaturation threshold. Settings of process parameters are initial concentration: 0.105 g/g; seeding temperature: 34 °C; seed loading: 5 wt %; no aging time after seeding; linear cooling time: 120 min; and aging at the final temperature: 10 min. Seed 2 was used.

Figure 6. (a) Residuals scatter of the mathematic model and (b) histogram of residuals.

Figure 9. Simulated supersaturation trajectories without (blue curve) and with (green curves) model parameter uncertainty or operational variability in starting concentration. 63.4% of the green curves have a peak value below the threshold (horizontal bold red line). The settings of process parameters are the same as those in Figure 7.

have to be decoupled from each other in order to derive a fixed operating range for each of them. Two risks are associated with seeding temperature. First, the initial supersaturation upon seeding may exceed the threshold or be negative (seed dissolution) if seeding temperature is too low or too high. Second, the supersaturation peak appearing after seeding may exceed the threshold if seed mass is too low or linear cooling time is too short. The operating ranges of seeding temperature can be determined through Monte Carlo simulation when the second risk is ignored. The initial supersaturation upon seeding depends on starting concentration and seeding temperature. At a certain set point of seeding temperature, 105 combinations of possible actual starting concentration and possible actual seeding temperature were obtained by random sampling from their respective distributions, which resulted in 105 possible values of initial supersaturation. Then, the probability that initial supersaturation falls between zero and the threshold can be calculated for this set point of seeding temperature. Figure 11 shows the

Figure 7. Contour plot of probability density function of kG and g.

in Figure 9. As a result, the lower limit for aging time after seeding is set at 30 min. In addition, the lower limit for aging time at the final temperature is set at 10 min, which has been found to be sufficient to deplete the remaining supersaturation for the model system. Now, we have three process parameters whose operating ranges are to be sought out, including seeding temperature, seed loading, and linear cooling time. There are strong interactions between them. When one of them is changed, the settings of the other two must be adjusted accordingly to have a peak supersaturation lower than the threshold. They 3954

dx.doi.org/10.1021/cg500547m | Cryst. Growth Des. 2014, 14, 3949−3957

Crystal Growth & Design

Article

Figure 10. Simulated supersaturation trajectories with an age time of 30 min after seeding. The blue curve is obtained from a simulation without either model parameter uncertainty or operational variability in starting concentration. 97.3% of green curves have a peak value below the supersaturation threshold (horizontal bold red line). The settings of other process parameters are the same as those in Figures 7 and 8.

Figure 12. Probability plot of peak supersaturation not exceeding the threshold versus linear cooling time when seed loading is set at 5 wt %. Uncertainty in model parameters and operational variability in starting concentration, seeding temperature, and seed loading are taken into account.

each point in Figure 12 takes 105 process simulations that account for uncertainty in model parameters, variability in starting concentration, seeding temperature, and seed loading. Likewise, the lower limits of linear cooling time for three more seed loadings are calculated as 180, 120, and 100 min for a seed loading of 3, 4, and 6 wt %, respectively. Any one of the pairs can be used as operating ranges on scale. It seems that the lower limit of linear cooling time does not change much when seed loading is higher than 5 wt %. Therefore, the lower limit for seed loading is chosen as 5 wt %. The operating ranges of all process parameters for co-crystallization of caffeine−glutaric acid are summarized in Table 2. It should be mentioned that these ranges are predicated on the assumed operational variability indicated in Table 2 and the model parameter uncertainty depicted in Figure 7.

Figure 11. Probability plot of initial supersaturation falling between 0 and 0.015 g/g (threshold) at different set points of seeding temperature. When seeding temperature is set between 33.9 and 34.6 °C, the probability will be equal to or higher than 70%.

Table 2. Operating Ranges and Assumed Variability of Process Parameters process parameters starting concentration, g/g D-values of seeds D10, μm D50, μm D90, μm seed loading, wt % seeding temperature, °C temperature profile age time after seeding, min linear cooling time, min age time at final temperature, min

probability of initial supersaturation falling between zero and the threshold at different set points of seeding temperatures. The probability curve has a convex shape. When seeding temperature is set between 33.9 and 34.6 °C, the probability will be equal to or above 70%. A probability target of 70% was chosen for illustrative purposes in this study. The operating ranges of seed loading and linear cooling time are then calculated at the boundaries of the seeding temperature range, i.e., 33.9 and 34.6 °C, respectively. At a certain different seed loading, the probability of peak supersaturation being lower than the threshold changes with linear cooling time and seeding temperature. For example, the probability is plotted against linear cooling time in Figure 12 when seed loading is set at 5 wt %. With a seeding temperature of 33.9 °C (the curve marked by squares), the probability is 49.7% when the linear cooling time is 60 min. In comparison, the probability is 27.2% when seeds are applied at 34.6 °C with a linear cooling time of 60 min (the curve marked by diamonds). Both curves cross the horizontal line of 70% when the linear cooling time is lengthened to 100 min. As a result, the lower limit of linear cooling time is chosen as 100 min for seeding between 33.9 and 34.6 °C. Probability calculation at

operating range

standard deviation

0.105

0.00250

≤67.0 ≤113 ≤151 ≥5.00 33.9−34.6

0 0 0 0.125 1.00

≥30.0 ≥100 ≥10.0

0 0 0

4.6. Verification of the Design Space on Scale-Up. The design space represented in Table 2 was verified with a scale-up experiment in a 10 L crystallizer for one combination of process parameters (experiment 7 in Table 1). Simulated and experimental supersaturation trajectories are shown in Figure 13 along with the temperature profile. Seeds were applied at time zero, and there was no aging period after seeding because the variability of process parameters on this scale is still negligible. It can be seen that the model agrees well with scale3955

dx.doi.org/10.1021/cg500547m | Cryst. Growth Des. 2014, 14, 3949−3957

Crystal Growth & Design



up results. Metastable co-crystals did not appear throughout the scale-up batch.

5. CONCLUSION First-principles process modeling was employed to explore the design space for co-crystallization of caffeine−glutaric acid on a bench scale and then to map the model-based design space into fixed operating ranges of individual process parameters for larger scales. Uncertainty in model parameters and operational variability in process parameters were demonstrated to cause dispersion of supersaturation trajectories and pose risks to polymorphic purity. These risks, meanwhile, could be mitigated effectively by sufficient aging after seeding. The operating ranges of process parameters on larger scales were determined sequentially via Monte Carlo simulation whereby uncertainty in model parameters and variability in starting concentration, seeding temperature, and seed loading were accounted for. The operating ranges were verified with one scale-up experiment. The uncertainty and operational variability are assumed to obey normal distributions in this study for simplicity, which does not exclude the use of real distributions in industrial applications. AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected] (Z.Q.Y.). *E-mail: [email protected] (R.B.H.T.). Notes

The authors declare no competing financial interest.



c g t M T kG mT GL c* Δc μj ρS φv

REFERENCES

(1) Dobry, D. E.; Settell, D. M.; Baumann, J. M.; Ray, R. J.; Graham, L. J.; Beyerinck, R. A. A Model-Based Methodology for Spray-Drying Process Development. J. Pharm. Innovation 2009, 4, 10. (2) Chen, W.; Chang, S.-Y.; Kiang, S.; Marchut, A.; Lyngberg, O.; Wang, J.; Rao, V.; Desai, D.; Stamato, H.; Early, W. Modeling of Pan Coating Processes: Prediction of Tablet Content Uniformity and Determination of Critical Process Parameters. J. Pharm. Sci. 2010, 99 (7), 3213−3225. (3) Hallow, D.; Mudryk, B.; Braem, A.; Tabora, J.; Lyngberg, O.; Bergum, J.; Rossano, L.; Tummala, S. An Example of Utilizing Mechanistic and Empirical Modeling in Quality by Design. J. Pharm. Innovation 2010, 5 (4), 193−203. (4) Brueggemeier, S. B.; Reiff, E. A.; Lyngberg, O. K.; Hobson, L. A.; Tabora, J. E. Modeling-Based Approach towards Quality by Design for the Ibipinabant API Step. Org. Process Res. Dev. 2011, 16 (4), 567− 576. (5) Burt, J. L.; Braem, A. D.; Ramirez, A.; Mudryk, B.; Rossano, L.; Tummala, S. Model-Guided Design Space Development for a Drug Substance Manufacturing Process. J. Pharm. Innovation 2011, 6 (2), 181−192. (6) Figueroa, I.; Vaidyaraman, S.; Viswanath, S. Model-Based ScaleUp and Design Space Determination for a Batch Reactive Distillation with a Dean-Stark Trap. Org. Process Res. Dev. 2013, 17 (10), 1300− 1310. (7) Castagnoli, C.; Yahyah, M.; Cimarosti, Z.; Peterson, J. J. Application of Quality by Design Principles for the Definition of a Robust Crystallization Process for Casopitant Mesylate. Org. Process Res. Dev. 2010, 14 (6), 1407−1419. (8) Randolph, A. D.; Larson, M. A. Theory of Particulate Processes Analysis and Techniques of Continuous Crystallization, 2nd ed.; Academic Press: New York, 1988. (9) Nagy, Z. K.; Fevotte, G.; Kramer, H.; Simon, L. L. Recent Advances in the Monitoring, Modelling and Control of Crystallization Systems. Chem. Eng. Res. Des. 2013, 91 (10), 1903−1922. (10) Worlitschek, J.; Mazzotti, M. Model-Based Optimization of Particle Size Distribution in Batch-Cooling Crystallization of Paracetamol. Cryst. Growth Des. 2004, 4 (5), 891−903. (11) Togkalidou, T.; Tung, H. H.; Sun, Y.; Andrews, A. T.; Braatz, R. D. Parameter Estimation and Optimization of a Loosely Bound Aggregating Pharmaceutical Crystallization Using in Situ Infrared and Laser Backscattering Measurements. Ind. Eng. Chem. Res. 2004, 43 (19), 6168−6181. (12) Barrett, M.; McNamara, M.; Hao, H. X.; Barrett, P.; Glennon, B. Supersaturation Tracking for the Development, Optimization and Control of Crystallization Processes. Chem. Eng. Res. Des. 2008, 88 (8A), 1108−1119. (13) De Calderon Anda, J.; Wang, X. Z.; Lai, X.; Roberts, K. J.; Jennings, K. H.; Wilkinson, M. J.; Watson, D.; Roberts, D. Real-Time Product Morphology Monitoring in Crystallization Using Imaging Technique. AIChE J. 2005, 51 (5), 1406−1414. (14) Chen, J.; Sarma, B.; Evans, J. M. B.; Myerson, A. S. Pharmaceutical Crystallization. Cryst. Growth Des. 2011, 11 (4), 887−895. (15) Yu, Z. Q.; Chow, P. S.; Tan, R. B. H.; Ang, W. H. PAT-Enabled Determination of Design Space for Seeded Cooling Crystallization. Org. Process Res. Dev. 2013, 17 (3), 549−556. (16) Garcia, T.; McCurdy, V.; Watson, T. N. J.; am Ende, M.; Butterell, P.; Vukovinsky, K.; Chueh, A. Verification of Design Space Developed at Subscale. J. Pharm. Innovation 2012, 7 (1), 13−18. (17) Lepore, J.; Spavins, J. PQLI Design Space. J. Pharm. Innovation 2008, 3 (2), 79−87. (18) Mitchell, J. D.; Abhinava, K.; Griffiths, K. L.; McGarvey, B.; Seibert, K. D.; Sethuraman, S. Unit Operations Characterization Using Historical Manufacturing Performance. Ind. Eng. Chem. Res. 2008, 47 (17), 6612−6621. (19) Yu, Z. Q.; Chow, P. S.; Tan, R. B. H. Operating Regions in Cooling Cocrystallization of Caffeine and Glutaric Acid in Acetonitrile. Cryst. Growth Des. 2010, 10 (5), 2382−2387.

Figure 13. Simulated and experimental supersaturation trajectories and temperature profile in scale-up experiment. Seeds were applied at time zero.



Article

NOMENCLATURE solute concentration, g/g crystal growth order appearing in eq 5, dimensionless time, min solvent mass, kg temperature, °C crystal growth rate constant appearing in eq 5, m/min solid mass in slurry, kg linear crystal growth rate, m/min solubility, g/g supersaturation, g/g the jth moment of number-based particle size distribution crystal density, kg/m3 volume shape factor, dimensionless 3956

dx.doi.org/10.1021/cg500547m | Cryst. Growth Des. 2014, 14, 3949−3957

Crystal Growth & Design

Article

(20) Figueroa, I.; Vaidyaraman, S.; Viswanath, S. Model-Based Scaleup and Design Space Determination for a Batch Reactive Distillation with a Dean−Stark Trap. Org. Process Res. Dev. 2013, 17 (10), 1300− 1310. (21) Hermanto, M. W.; Chiu, M.-S.; Woo, X.-Y.; Braatz, R. D. Robust Optimal Control of Polymorphic Transformation in Batch Crystallization. AIChE J. 2007, 53 (10), 2643−2650. (22) Nagy, Z. K. Model Based Robust Control Approach for Batch Crystallization Product Design. Comput. Chem. Eng. 2009, 33 (10), 1685−1691. (23) Simone, E.; Saleemi, A. N.; Tonnon, N.; Nagy, Z. K. Active Polymorphic Feedback Control of Crystallization Processes Using a Combined Raman and ATR-UV/Vis Spectroscopy Approach. Cryst. Growth Des. 2014, 14 (4), 1839−1850. (24) Yu, Z. Q.; Chow, P. S.; Tan, R. B. H.; Ang, W. H. Supersaturation Control in Cooling Polymorphic Co-Crystallization of Caffeine and Glutaric Acid. Cryst. Growth Des. 2011, 11 (10), 4525−4532. (25) Myerson, A. S. Handbook of Industrial Crystallization, 2nd ed.; Butterworth-Heinemann: New York, 2002. (26) Hermanto, M. W.; Kee, N. C.; Tan, R. B. H.; Chiu, M.-S.; Braatz, R. D. Robust Bayesian Estimation of Kinetics for the Polymorphic Transformation of L-Glutamic Acid Crystals. AIChE J. 2008, 54 (12), 3248−3259. (27) Paul, E. L.; Atiemo-Obeng, V.; Kresta, S. M. Handbook of Industrial Mixing: Science and Practice; John Wiley & Sons: New York, 2003.

3957

dx.doi.org/10.1021/cg500547m | Cryst. Growth Des. 2014, 14, 3949−3957