Definition of Design Spaces Using Mechanistic Models and Geometric

Article pubs.acs.org/OPRD

Definition of Design Spaces Using Mechanistic Models and Geometric Projections of Probability Maps Salvador García-Muñoz,* Carla V. Luciani, Shankar Vaidyaraman, and Kevin D. Seibert Small Molecule Design and Development, Lilly Research Laboratories, Eli Lilly & Company, 1400 West Raymond Street, Indianapolis, Indiana 46221, United States ABSTRACT: Following a model-centric strategy in the development of a manufacturing process for a new medicine empowers the simultaneous study of a large number of process parameters, which is large enough to exceed the capability of a graphic representation of the interactions across them. This work presents a discussion regarding the identification, description, and communication of multidimensional design spaces of high order. It introduces the reader to mathematical tools developed by the process systems engineering community that become relevant in the challenge to replace graphics as a means to describe and communicate a design space. Concepts like process f lexibility are discussed and illustrated. The paper also introduces geometric projection as a way to capture and describe the shape of the design space in an easier form (than that of the complete mechanistic model) that can be communicated to the regulator. An assessment is presented regarding the key elements communicated by a graphical representation of a design space, and alternate ways of conveying the same information using mathematics are suggested. These ideas are illustrated by applying them to the identification and definition of a design space for a chemical reaction step and the digital risk assessment for a packed bed adsorption step.

■

INTRODUCTION Regulatory guidance for the pharmaceutical industry, contained in documents such as the ICH Q8 (R2) and Q11,1,2 presents the practitioner with the option to adopt an in-depth holistic pharmaceutical process and product development approach, referred to as “Quality by Design” (QbD). The practitioner is encouraged to invest efforts in a more systematic understanding of the complex relationships between materials, processes, and products involved in the manufacture of new medicines. Such understanding can potentially benefit the manufacturer with the regulatory approval to manufacture at any condition in a wide operational envelope; given enough scientific evidence that the process will result in product of acceptable quality. This envelope of operation is referred to as “the design space”. The immediate response of industry at the time is summarized in work by Huang et al.:3 (i) identify the greater “knowledge space”; (ii) carry out a well-designed set of experiments to determine the relationship between process parameters and quality attributes within the “knowledge space” (typically represented by an empirical regression of some form); (iii) use response surface methods to determine the borders of the design space; and (iv) run confirmatory experiments to test and verify the design space. Examples of this empirical approach to design space identification are widely available in literature for drug substance,4 intermediates,5 and drug product.6 The number of publications is of course much larger than this; these selected references are intended to serve as examples of the approach, not as a comprehensive literature review. As this practice matured, this early empirical approach toward the identification of design spaces was replaced with the use of mechanistic models. The equations in these models are derived from conservation laws of chemistry and physics © XXXX American Chemical Society

(not from data) and a set of assumptions about the physical system (e.g., perfectly well-mixed tanks). The validity of these models as they apply to new scenarios relies on how well the assumptions hold. Generally speaking, the concept of extrapolation does not really apply to these models in terms of the numerical inputs but only in terms of the assumptions. Examples of mechanistic models applied to the determination of design spaces are fewer in number but are also found for drug substance7 and drug product.8 An essential component in the successful application of QbD is the need to define the boundaries of the design space and communicate this to the regulatory agencies. Furthermore, the chosen form of communication for the definition and boundaries of the design space has to be such that it can also serve as a tool to audit the process and verify that the practitioner is adhering to its regulatory commitments. With advancements in the practice of QbD came the evolution of the forms of communication for the design space. The use of ranges (upper and lower limits for process parameters) was short-lived, and it did not take long before the use of graphical visualization techniques became perhaps the most popular tool in communicating, defining, and auditing design spaces [see Appendix 2, section C of Q8 (R2)1]. The use of graphical representations of the design space is particularly useful to communicate any necessary conditionality among process parameters. That is, when the allowable value for one parameter depends on the value of the other(s). This is well articulated by Lepore and Spavins9 in the description they provided for the design space: “Use of the design space approach allows more effective dialogue between industry and regulator··· It allows the regulator to more readily see the Received: May 20, 2015

A

DOI: 10.1021/acs.oprd.5b00158 Org. Process Res. Dev. XXXX, XXX, XXX−XXX

Organic Process Research & Development

Article

A Brief Historical Perspective. The use of process simulation in the chemical industry can be tracked back to the 1950s, just after the creation of the FORTRAN language by IBM. Process simulators flourished, and by 1975 Monstanto’s FLOWTRAN simulator had become the standard process simulation tool for industry and academia. Commercial packages to simulate chemical processes became available in the early to mid-1980s with the availability of packages like Aspen Plus (by Aspen Technology), ProII (by Simulation Sciences Inc.), HYSIM (by Hyprotech Ltd.), and gPROMS (by Process Systems Enterprise Ltd.).15 By the early 1990s, the use of computer simulations in lieu of pure statistically driven experimentation had become a pivotal tool in the design of chemical plants.16 The publications and work done by the process systems engineering community is far too extensive to compile, and a review by Stephanopoulos and Reklaitis13 provides a good overview. Other contributions, like the one from Gernaey et al.,17 present a high level perspective of the potential role for a process systems approach to pharmaceutical development, whereas other papers18,19 present proposed applications of systematic product and process development for the development of a manufacturing process for the active pharmaceutical ingredient (API). Relevant Mathematics to QbD and Pharmaceutical Development with Mechanistic Models. Three mathematical exercises are discussed below, as they represent important tools in a model-based approach to process development and QbD. Model Uncertainty and Uncertainty Propagation. Once the model for a system is established, along with the necessary assumptions, the next step is to estimate the model parameters (e.g., kinetic parameters, heat transfer coefficients, etc.). These parameters are needed for the model to be relevant to the specific application and are estimated with experiments designed for this very purpose.20 The parameter estimation exercise for mechanistic models differs from that used for empirical models in the sense that the equations themselves are not in question, just the values of the parameters. Techniques for parameter estimation vary from the simple minimization of the prediction errors (least-squares) to more advanced techniques that account for prior knowledge about uncertainties in the data21 (e.g., Bayesian estimation, errors in variables, and maximum likelihood). Two outcomes result from this exercise: (a) an optimal value for the parameters and (b) their variance-covariance. The variancecovariance of the parameters becomes important as a necessary piece of information to produce the corresponding uncertainty for a prediction by the model. The model parameter uncertainty can then be propagated to the prediction. The resulting prediction uncertainty can be generated either in the form of a confidence interval or an explicit distribution of possible values of the prediction. The latter can be produced with a Monte Carlo simulation,22 accounting for expected variability in model parameters as well as common cause variability. Constraints, Probability of Success, and Design Spaces. An essential component in QbD is a well-defined set of constraints that needs to be fulfilled. These constraints can be product related or process related and are either one or two sided inequality constraints that are easily represented in mathematical terms. A pivotal moment in the advancement of QbD was reached when practitioners recognized that it was possible to mathematically defined a design space as a subspace of the knowledge space, where the estimated probability for a

connection between the design space, criticality and control strategy”. The versatility of plots and graphics resides entirely in the capacity of the human brain to interpret spatial patterns. And while the details of the brain mechanisms that transform vision into knowledge remain very much a mystery,10 it is undeniable that the complexity of the spatial relationships that the human brain can process is vastly superior to that of a computer.11 The human brain is indeed a biological computing machine of great power but with one fundamental limitation: its visual engine only works up to the third dimension. And here resides the critical limitation of graphical representations of the design space: when the practitioner has the need to communicate conditionality among a large number of process parameters to define a design space. In spite of the advances made in the development of visualization techniques,12 the graphical representation of higher dimensional spaces remains one of the greater challenges in computer sciences. This becomes a limitation when mechanistic models are used to understand pharmaceutical processes. Mechanistic modeling technology provides the power to explore the effect of a large number of parameter interactions on product quality. Furthermore, the use of a flow-sheet simulation allows for the simultaneous exploration of the interactions of variables across a manufacturing train. With this approach, the practitioner is empowered to explore the multidimensional interactions across dozens, if not hundreds, of process parameters, which clearly exceed what can be communicated in a graphical format. The purpose of this work is to propose alternatives to define, delineate, and communicate a high-dimensional design space. To do this, we present a discussion of the elements that are communicated with a graphic and suggest potential ways to replace this information with the use of mathematical tools and techniques developed by the process systems engineering community. The overall process is depicted in a flow diagram in Appendix D; the work presented in this paper specifically addresses the last box labeled “Design Space Development”. As such, verification of the model structure and the model parameters is not the scope of this work and is not addressed in full detail. This work is organized as follows: After the introduction, the Discussion introduces the reader to the process systems engineering community and discusses the historical context of the tools that are proposed later to replace graphics. An analysis of the elements we believe are communicated with a graphic, and alternative ways to convey this information using modeling tools are shown in the Replacing Plots with Mathematics section. The Case Studies section shows various examples of these concepts with two different processes in a drug substance manufacturing train, and our final remarks and conclusions are presented in the Conclusions section.

■

DISCUSSION The Process Systems Engineering Community and Relevant Tools to Pharma. Chemical engineersthe process systems engineering community specificallyhave taken a central role in QbD discussions, given that many of the challenges of modern pharmaceutical development and manufacturing have been addressed by this community and implemented for years for other industries.13 The new trends and concepts in pharmaceutical development (e.g., Process Analytical Technology, Quality by Design, Real Time Release) are only new names for old concepts that historically precede the introduction of Q8 and Q11.14 B



Article

the maintenance of satisfactory performance despite adverse conditions, while flexibility is the ability [of the process] to handle alternate (desirable) operating conditions”.27 A review of the current state of research in process flexibility and resiliency is given by Grossmann.32 Replacing Plots with Mathematics. To develop a strategy to replace a graphical representation of a design space (i.e., plots), one first needs to understand the information pieces communicated by such plots. Let us define the following elements as the key pieces of information gathered from a graphical representation of a design space: (i) The boundaries of the design space within the knowledge region: delineation (ii) The location of the operating condition that is furthest away from all constraints (the centroid of the design space): center. (iii) The location of operating conditions of interest that will increase the risk of violating constraints: key boundaries of the design space. (iv) Upper and lower bounds of operation for individual variables: the box. It is assumed that the design space is being defined on a probability map, where the boundaries of the design space already consider the effect of the propagation of uncertainty and expected common cause variability. One way to obtain the information about the first element (design space delineation) is to define the shape of the design space with mathematics; this is purely a geometric exercise. In the most general case, the task is to find a reasonable set of process conditions (set of points) at which the probability function is at its constraint value and then fit an appropriate and conservative geometric approximation to the shape defined by the obtained set of points. This approximation is referred to as the geometric projection of the design space. These geometric models are all that is needed to determine whether a combination of process parameters is part of a design space or not; thus, their equations (specifically the inequalities that arise from them) could be used to represent the process commitments, while the mechanistic model from which the probability map was derived, which led to the geometric projection model, can be discussed in supporting sections of the regulatory filing. A similar concept was presented by Stockdale and Cheng,33 where the design space is subdivided in regular rectangular regions to accommodate the geometric irregularity of the design space. The second element (centroid) can be calculated from the geometric model and represents the safest target operating condition because it is the furthest away from any constraint. While the inclusion of this information as part of the disclosures to the regulator is optional, the information itself is much needed to establish target manufacturing conditions and write master batch records. The third element (boundaries of the design space) is the complement to the first element. For this case, the inverse of the model can be solved to find scenarios that lie in extremes of the design space. These conditions can guide targeted experimentation to confirm or challenge the model framework and provide further confidence that the design space boundaries are well-defined. Careful planning and judicious analysis of the experimental results needs to be exercised due to the probabilistic nature of the boundaries (e.g., if an experiment is carried out under conditions that predict a 90% probability of

system to fulfill all given constraints is greater than a minimum acceptable risk.23−25 This can be written as follows: DS = {x ∈ Kn|h(x) ≥ π )} subject to h(x , σ , Θ , ΣΘ , f , g ) y = f (x , Θ ) g (x , y)

(1)

where DS, design space (a subset of x); f, the model of the system; g, the set of constraints; h, the probability function for each element of x to fulfill the constraints (g); x, all possible combinations of process parameters in the n dimensional knowledge space (Kn); y, model predictions (presumably quality attributes of interest); π, minimum acceptable risk; σ, expected common cause variability; Θ, set of model parameters; and ΣΘ, variance covariance of model parameters Eq 1 is an example of the mathematical formalism that would define the design space irrespective of the number of process and model parameters and constraints. In principle, an expression similar to eq 1, with the proper expressions for f, g, and h, could be an option to define the design space in the commitments section in a regulatory filing. Implementations of this probability calculation have already been reported23,24,26 using a combination of empirical and deterministic models. Additionally, most reported cases involve systems that are small enough that the estimated probability for the entire knowledge space can be represented in a plot. We will refer to these representations as the probability map. A design space established in the probability map has the conservative advantage to be inclusive of all other sources of uncertainty and not just inclusive of the knowledge captured by the model. Flexibility and Resiliency. The concepts of process flexibility and process resiliency were first proposed by Grossman and Morari in 1984.27 Flexibility is presented as “...the flexibility of a design [referring to a chemical process] is determined by its capability to meet constraints and specifications for a range of conditions”.28 The work by Grossman and his collaborators explores this concept and proposes different optimization methods29,30 to calculate the flexibility region for a given process given a model and a set of constraints. In all of these cases, the flexibility region is restricted to be a multidimensional hypercube. The vertices of the hypercube correspond to the maximum and minimum possible values for the process parameters, such that all possible combinations of process parameters within these ranges will meet all the constraints. Very likely one or more of these vertices will hold against one (or more) of the constraints. The problem of calculating the flexibility region can be translated in layman terms as the calculation of the largest multidimensional box in the operational space where all constraints are met. Resiliency31 can be described as “transient flexibility”: flexibility calculates the operating region in which the process is capable of fulfilling all of the constraints in the absence of undesired disturbances. Resiliency, on the other hand, considers the dynamic aspect of the process driven mainly by the design of the automatic control system acting on it. A resilient process is said to be the one that can recover from disturbances in a safe, consistent, and stable manner. The main difference between flexibility and resiliency is that “...resiliency refers to C



Article

Figure 1. Suzuki coupling, main catalytic cycle for the formation of the desired product (P1) and side reaction for the formation of dimeric impurity (Imp1). Reactions numbers and species nomenclature are used within the main text.

success, it is only natural to expect that 1 out of 10 times the quality and/or operational constraints will be violated). The fourth element (the operational box) is easily answered by finding the flexibility region for the process. The mathematical formalization to find the flexibility region for a multidimensional problem (i.e., a feasible set of multivariate process conditions where all constraints are met) is directly applicable to this aspect of QbD. It is especially useful if the design space is being established across multiple unit operations, implying the involvement of a large number of process parameters. The upper and lower ranges for each of the process parameters that are obtained with the flexibility analysis are valuable information in the writeup of manufacturing instructions and master batch records. The following section presents examples of these concepts. Case Studies. In what follows, two case studies are described in which mechanistic models were used to represent the feasible operating region and aid technical decision-making in the development of a pharmaceutical intermediate. To this purpose, the following sequence of unit operations were investigated: (i) Suzuki coupling reaction between a boronic ester (SM1) and an organohalide (SM2) to produce the desired coupling product (pharmaceutical intermediate P1) and an SM1-related dimeric impurity (Imp1) (Figure 1). The investigated Suzuki coupling reaction is biphasic and performed in a fully batch mode. Typically, an aqueous phosphate solution (source of OH− groups, as shown in reactions 10 and 11 in Figure 1) is transferred to the main reaction vessel containing a solution of SM1 and SM2 in THF. After the mixture is intensely degassed, a Pd-based catalyst is transferred to the solution, and the reacting mixture is heated to the desired reaction temperature. During the Suzuki coupling reaction, Pd(II) catalyst is in situ reduced to Pd(0) species, which is the active species in the

main catalytic cycle (reaction 1 in Figure 1). Both the presence of dissolved oxygen and an aqueous phase containing dissolved salts are known to cause Pd speciation in Suzuki coupling reactions.34,35 Pd(II) complexes can promote the (undesired) boronate ester/boronic acid homocoupling pathway (reactions 8 and 9 in Figure 1).36 To prevent excessive Pd speciation, we implemented a thorough oxygen removal protocol as part of the overall control strategy. (ii) Continuous palladium removal in packed bed columns. When the Suzuki coupling reaction is complete (determined by an in-process control), the aqueous phase is separated from the organic phase, and the organic phase is washed with a KCl/ water solution. After this liquid−liquid extraction, the resulting stream is continuously passed through a series of packed bed columns to ensure that the in situ level of residual Pd meets a stringent specification. The column packing material, a silicabased thiol-functionalized scavenger, is used to break down Pd complexes in solution and promote binding to the metal. Recent publications37,38 reported that Pd speciation affects the performance of metal scavenging processes. The detrimental impact of dissolved oxygen to scavenging performance has been confirmed experimentally in our laboratories. In what follows, the concepts outlined in the previous sections are used to aid process development in these two case studies. Suzuki Reaction: Calculation of the Flexibility Region. Although a mathematical representation of the design space comprises all the information required to define the largest feasible region for operation, the ability to define a subset of that region through a “multidimensional box” is appealing. Process parameter ranges provide a useful insight into process parameter criticality and overall process robustness, simplify the communication between R&D and manufacturing organizations, can be easily implemented in control loops, and make the D



Article

determined by running numerous simulations for different combinations of process parameters and sampling model parameters from variance-covariance matrix (i.e., model parameter uncertainty) and accounting for output uncertainty from its common cause variability. For illustration purposes, this particular example was done with 200 realizations of the Monte Carlo simulation using an 85% probability as the acceptable minimum. The outer light gray region in Figure 4 shows the combination of process parameters in which both quality attributes are met (without accounting for uncertainty), whereas the inner darker region is where the probability of meeting both quality attributes was deemed acceptable once uncertainty is accounted for. The lower boundary of the feasible operating region is imposed by the reaction completion constraint, whereas the upper boundary is imposed by the level of Imp1 constraint. As explained before, the flexibility region is the largest multidimensional box that can be fit within the feasible operating region. An explicit enumeration/heuristic approach was used to determine the flexibility region as outlined: (1) Pick a candidate set of center points for the multidimensional box and set an aspect ratio for the box. (2) For each candidate center point, find the largest box size that fits within the feasible operating region. Because the aspect ratio of the box is fixed, this is a single variable optimization and is implemented using a bisection search. (3) Pick the box (from boxes with different center points) that has the largest size. (4) Manually check if certain parameters can be wider for the largest box determined in step 3. This step is needed to address the following limitations of the algorithm used to find the biggest box: (a) While setting the aspect ratio simplifies the search, it also constrains the search space. (b) The input to the algorithm is the set of feasible points that are generated on a discretized space of process parameters. It is necessary to expand or shrink certain parameters of the box to the grid of points used for estimating feasibility. Further details on implementation are described in Appendix C. The flexibility region for the Suzuki reaction is shown in Table 2. Figure 4 also shows a graphical view of the flexible region at an SM1/SM2 molar ratio of 1.03 and initial Pd speciation of 0.1 (see dots). It is proper to disclose that the method followed to estimate the flexibility region was driven by software and the fact that (at the point the work was executed) the model was not implemented in a platform that allowed the use of a formal optimization method. Even though the multidimensional design space was very useful for providing a complete picture of the impact of process parameters to relevant quality attributes, a subsequent risk assessment exercise showed that (a) a validated degassing protocol as well as limiting oxygen in the reactor head space mitigated the risk of both a high level of dissolved oxygen and extensive Pd speciation, (b) an end-of-reaction in-process control ensures reaction completion, mitigating any residual risk imparted by short reaction time, and (c) across the controllability region, the reaction temperature has no impact on quality attributes. As a result of the risk assessment exercise, two process parameters were identified as having some residual risk in achieving acceptable quality attributes: SM1/SM2 molar ratio and catalyst loading. The probability map for the 2-dimensional

detection of potential events or deviations from normal operation straightforward. The mathematical model of the Suzuki coupling reaction is described in Appendix A. The process parameters (model inputs) considered by the model include molar ratio of starting materials, reaction volume, solvent composition, reaction temperature, catalyst loading, initial Pd speciation, O2 level in reactor head space, reaction time, age time, and potassium phosphate charged amount. Table 1 lists Table 1. Model Parameters parameter

value Suzuki Couplinga

ks2,ref (L mol−1 min−1) 1.74 × 102 ks3,ref (L mol−1 min−1) 7.14 × 103 ks4,ref (L mol−1 min−1) 4.37 × 103 −1 ks5,ref (min ) 2.49 × 103 −1 ks6,ref (min ) 2.94 Keq,6 = ks6,ref/ks−6,ref () 2.08 ks7,ref (L0.5 mol−0.5 min−1) 1.04 × 101 ks8,ref (L mol−1 min−1) 7.54 ks9,ref (L mol−1 min−1) 3.84 × 102 ks10,ref (min−1) 5.50 Keq,10 = ks10,ref/ks−10,ref (L−1 mol) 3.96 × 10−2 −1 ks11,ref (min ) 2.34 Keq,11 = ks11,ref/ks−11,ref (L−1 mol) 5.00 × 10−4 ks12,ref 2.74 × 10−3 Packed Bed Columnsb Deff,Pd(0) (m2/min)c 1.12 × 10−13 c 2 Deff,Pd(II) (m /min) 5.30 × 10−16 KL,Pd(0) (g of solution/mg of Pd(0)) exp(−10045/T + 35.743) KL,Pd(II) (g of solution/mg of Pd(II))c 9.1 qm,Pd(0) (mg of Pd(0)/g scavenger) 35.8 qm,Pd(II) (mg of Pd(II)/g scavenger) qm,Pd(0) a ki(T) = ksi,ref exp(−Eai/R(1/T − 1/Tref)) with Tref = 303.15 K; activation energies are as follows (kJ/mol): Ea2 = 27.4, Ea3 = 0.0, Ea4 = 15.3, Ea5 = 22.4, Ea6 = 30.0, Ea,eq6 = −65, Ea,7 = 0.0, Ea,8 = 20.1, Ea,9 = 0.0, Ea,10 = 30, Ea,eq10 = 50, Ea11 = 30.0, Ea,eq11 = 139.0, and Ea,12 = 15.0. b Measured values: Dpart = 66 μm and ε = 0.62; simulations were insensitive to Dax,c, kf,Pd(0), and kf,Pd(II). cNo significant temperature dependence was found within the tested range (48−72 °C).

the model parameters used for simulations, and a few examples are compared to experimental data in Figure 2. Overall model adequacy for this case was assessed partly by using statistical testing on the residuals for the predictions of interest. Figure 3 illustrates the distribution of the model prediction residuals using 54 experimental points for SM2 and 71 points for Imp1. The spread of the residuals and their mean value are then judged as adequate (or not) according to the need. The extent of the specific diligence and efforts toward the verification of the model adequacy (i.e., model validation) depends on the end use of the model per ICH guidance.39 Given the spread and mean of the residuals shown in Figure 3, and the specific needs of our application, the model was deemed to be adequate for the intended use. Two quality attributes limit the feasible operating region or design space: (i) reaction completion, which establishes a maximum amount of unreacted SM2 for the reaction considered to be complete; and (ii) a maximum allowed level of Imp1 below which the stream is known to acceptably process forward. A simplified graphical representation of the multidimensional design space is shown in Figure 4, where four process parameters are represented. The multidimensional design space was E



Article

Figure 2. Comparison of experimental data (symbols) from the extreme limits of oxygen levels and model simulations (traces) for the Suzuki coupling reaction. Experiments were performed at target operating conditions for most process parameters, except for when (a and b) the oxygen level in the reactor head space was varied and when (c and d) the reaction temperature was varied. Error bars represent one standard deviation.

Figure 3. Histogram of the distribution of model residuals for (a) residual amount of SM1 and (b) in situ generated Imp1.

more robust in terms of handling variability in both process parameters and feed stream, which makes it also a more robust option for pharmaceutical processes. The larger the number of columns in series, the greater the ability of the system to handle larger campaign sizes and larger variability in process parameters without observing column breakthrough. On the contrary, operating and capital costs increase significantly as the number of columns in the series increases. An interesting feature of pharmaceutical processes is that campaign size may change significantly over the drug substance life cycle, and the process needs to adequately handle those types of changes also. In this last case study, the mathematical representation of the design space was used to determine the number of columns in series that ensures a robust operation at commercial scale (i.e., the process is capable of delivering high quality material under process parameter-expected variability for large enough campaign sizes).

feasible region is shown in Figure 5. Also shown is the geometric (more conservative) projection of the probability map, which can be simply represented by the following set of equations. 1.001 ≤ [SM1]0 /[SM2]0 < 1.03 ⇒ 0.001 ≤ [Pd]0 /[SM2]0 ≤ 0.003 1.03 ≤ [SM1]0 /[SM2]0 ≤ 1.09 ⇒ 0.001 ≤ [Pd]0 /[SM2]0 ≤ 0.0026

Continuous Pd Removal in Packed Bed Columns: Boundaries for the Design Space to Determine the Number of Columns in the Series. Packed bed adsorption columns inherently operate under nonsteady state conditions. For purity reasons, the scavenger is not regenerated, and fresh resin is repacked after a predefined time, which will determine the maximum campaign size allowed for the selected configuration. Whereas parallel configuration is more efficient in terms of resin usage and pumping requirements, series configuration is F



Article

Figure 4. Probability map projections for Suzuki coupling reaction at an SM1/SM2 molar ratio of 1.03 and initial Pd speciation of 0.1 (molar fraction) within the feasible operating region. Inner axes: Reaction time and equiv of catalyst. Outer axes: O2 level in reactor head space and reaction temperature. The darker region is where the probability of meeting both quality attributes is acceptable. The light gray region is where both quality attributes are met (without accounting for uncertainty). Dots indicate the limits of the flexibility box (inner axes) as they collide with constraints.

Table 2. Flexibility Region for Suzuki Coupling Reaction process parameter

flexibility region

SM1/SM2 (molar ratio) catalyst/SM2 (molar ratio) reaction time (min) reaction temperature (°C) O2 in head space (ppm) initial Pd speciation (molar fraction)

1.01−1.20 0.002−0.0026 110−230 42−64 10−139 0.05−0.14

The model equations and assumptions are presented in Appendix B. Compared to experimental data, the model provided acceptable predictions for the intended use. For instance, simulation results (full manufacturing scale experiments) are compared to experimental data in Figure 6, and a conservative prediction of the Pd levels is provided. Table 1 lists the model parameters used for simulations. Packing material particle size and void fraction are measurable variables rather than model parameters (see Table 1 note). Because the simulation tool is computationally expensive, no probabilistic calculations were performed for this example. Instead, a conservative specification was set as a quality constraint. Model simulations were used to understand the impact of process parameter variability and campaign size for different numbers of columns in series. For this purpose, many combinations of process parameters were investigated, including

Figure 5. Probability map for the process parameters that exhibit some potential risk to quality attributed in Suzuki coupling reactions. Light squares show acceptable probability of meeting specifications. Dashed black line shows the adopted geometric projection of the probability map.

scavenger efficiency (assuming that functional group loading changes the resin capacity from 80 to 100% of the typical observed value); temperature (48 to 72 °C); superficial fluid velocity (0.5 to 0.6 cm/min); total Pd level in the feed stream (from 40 to 55 ppm, accounting for changes in catalyst charge in a Suzuki coupling reaction); and Pd(II) in the feed G



■

Article

CONCLUSIONS

Industry can, and should, leverage state of the art technology available in the modeling and simulation field to aid the development of a manufacturing process for new medicines. The use of mechanistic models of the process train enables analysis of the interactions of a large number of parameters. The advantage of these models is that their applicability resides in the extent to which the assumptions hold and numerical extrapolations are acceptable (as far as the assumptions can be sustained). The method proposed herein assumes that the mathematical abstraction of the system is adequate for the given application. The specific procedure followed for model verification largely depends on the end use of the model. In this respect, the ICH guidance provides a model classification system by their impact (low, medium, and high) and provides guidance as to how the practitioner can implement business processes for model verification and lifecycle.39 Once a model is available, the predictions for the entire knowledge space can be translated into a probabilistic map of meeting all constraints, accounting for the model uncertainty and common cause variability; this can be done either mathematically or explicitly with Monte Carlo simulations. We propose that a second model to be built by finding all points in the probability map (that could be multidimensional) in which the probability of passing all constraints is at a minimum acceptable level and then fitting a geometric model to this shape. This geometric projection can easily be translated into inequality equations that uniquely define the design space. Mathematical tools like flexibility analysis can also be used to determine the largest hypercube (multidimensional box) in the design space to provide manufacturing guidance, write master

Figure 6. Comparison of model result (traces) and experimental data (symbols) for a nominal campaign size run at full manufacturing scale.

stream (from 4.0 to 10.5 ppm, accounting for changes in Suzuki coupling reaction). Desired column geometry (i.e., diameter and length) was determined by space allocation, pressure drop consideration, and so forth. Therefore, column geometry was not changed as part of this analysis. The relative frequency histogram of the total Pd level in the stream leaving the last column is presented in Figure 7 for different campaign sizes and the number of columns in series. On the basis of the simulation results, three columns in series are needed to meet the specification for total Pd for the nominal campaign size (Figure 7a and b). This configuration has been implemented at full manufacturing scale.

Figure 7. Relative frequency histograms of Pd level in the stream, leaving the last column of the series at the end of (a) 0.5 × nominal campaign size; (b) 1 × nominal campaign size; (c) 1.7 × nominal campaign size; and (d) 2.5 × nominal campaign size. H



Article

batch records, and guide operational procedures. Furthermore, as the industry moves to continuous manufacturing, concepts like resiliency analysis represent a well-established methodology to ensure that the necessary automation system will deliver the necessary robustness to keep the process operating as expected in spite of uncertainty. The publication industry has realized that it is imperative to use new forms of communication to not become a bottleneck in the development of science. Publishers have made available a wide range of tools to better describe scientific concepts and disseminate knowledge (e.g., nowadays, one can rotate a 3D plot in online articles and explore live maps). The age of “two dimensional print” is coming to an end in the scientific publication field. We argue that other areas where scientific dialogue takes place (like the pharmaceutical regulatory world) should follow through and embrace new media and new languages (like mathematics) that enable better communication of the complex concepts found in the complex new reality in which we live today.

■

d{[SM1]V } = {− ks6[SM1] + ks − 6[A4]}V dt

(12)

d{[Imp1]V } = {− ks9[A4][A5]}V dt

(13)

(15)

dt

= {ks 4[A4][Z2] − ks11[HBO3−] + ks − 11[K 2B4 O7 ][OH−]}V

(16)

d{[OH−]V } = {ks10[K3PO4 ] − ks − 10[K 2HPO4 ][OH−] dt + ks11[HBO3−] − ks − 11[K 2B4 O7 ][OH−]}V (17)

d{[K 2B4 O7 ]V } = {ks11[HBO3−] − ks − 11[K 2B4O7 ][OH−]}V dt (18)

d{[P1]V } = {ks5[Z3]}V dt

(1)

d{[Pd(II)]V } = {ks7[Pd(0)][O2 ]L1/2 − ks8[Pd(II)][A4]}V dt

(19)

Henry constants as a function of temperature were obtained from Battino et al. (1983).41 The initial conditions to solve the ordinary differential eqs 1−19 are as follows: mPd(0),0 [Pd(0)]0 = McatalystV (20)

(2) (3)

⎛ ⎛ HO ,THF ⎞ ⎟⎟ HO2 = exp⎜⎜ ∑ xj ln HO2, j − x THF ln⎜⎜ 2 ⎝ HO2,water ⎠ ⎝ j = water,THF

d{[Z1]V } = {ks2[Pd(0)][SM2]L − ks3[Z1][OH−]}V dt

(11)

d{[HBO3−]V }

d{[Pd(0)]V } = {− ks2[Pd(0)][SM2] − ks7[Pd(0)][O2 ]L1/2 dt

d{[SM2]V } = {− ks2[Pd(0)][SM2]}V dt

d{[A5]V } = {ks8[Pd(II)][A4] − ks9[A4][A5]}V dt

d{[K 2HPO4 ]V } = {ks10[K3PO4 ] − ks − 10[K 2HPO4 ][OH−]}V dt

On the basis of the assumptions listed above, the following mass balances were used to describe the Suzuki coupling reaction.

d{[O2 ]L V } = {ks12([O2 ]L,sat − [O2 ]L )}V dt

(10)

(14)

Model Equations

⎛ HO ,water ⎞⎞ ⎟⎟ − ln⎜⎜x water + x THF 2 HO2,THF ⎟⎠⎟⎠ ⎝

(9)

d{[K3PO4 ]V } = {− ks10[K3PO4 ] + ks − 10[K 2HPO4 ][OH−]}V dt

(1) Mass transfer between phases is not rate controlling (i.e., equilibrium is adopted to estimate the amount of dissolved oxygen, and pseudo-homogenous assumption was adopted to describe the biphasic system); (2) Pd(II) initial catalyst is mostly instantaneously transformed into active species (see reaction 1 in Figure 1); (3) biphasic system is well-mixed; (4) model inputs are known and certain; and (5) the Wilson model was used to estimate the Henry constant of oxygen in the solvent mixture.40

2

d{[Z3]V } = {ks 4[A4][Z2] − ks5[Z3]}V dt

− ks8[Pd(II)][A4] − ks9[A5][A4]}V

Model Assumptions

yO P = [O2 ]L,sat HO2

(8)

d{[A4]V } = {ks6[SM1] − ks − 6[A4] − ks 4[A4][Z2] dt

APPENDIX A: SUZUKI COUPLING MODEL

+ ks9[A5][A4] + ks5[Z3]}V

d{[Z2]V } = {ks3[Z1][OH−] − ks 4[A4][Z2]}V dt

[Pd(II)]0 =

[SM2]0 = (4)

[SM1]0 = (5)

V= (7)

McatalystV

MSM2V

(22)

mSM1,0 MSM1V

(23)

mK3PO4 ,0 M K3PO4V

∑ j = P1,SM2,water,THF

I

(21)

mSM2,0

[K3PO4 ]0 = (6)

mPd(II),0

(24)

mj ,0 ρj

(25) DOI: 10.1021/acs.oprd.5b00158 Org. Process Res. Dev. XXXX, XXX, XXX−XXX


■

APPENDIX C: FLEXIBILITY REGION DETERMINATION The input to determination of the flexibility region is a set of simulations over a discretized space of process parameters where the probability of passing quality attributes is estimated at each point. The discretizations used are listed in Table 3.

[O2 ]L,0 = [Z1]0 = [Z2]0 = [Z3]0 = [A4]0 = [A5]0 = [Imp1]0 = [P1] = 0

(26)

[K 2HPO4 ]0 = [HBO3−]0 = [OH−]0 = [K 2B4O7 ]0 = 0

■

Article

(27)

APPENDIX B: PALLADIUM ADSORPTION IN A PACKED BED COLUMNS MODEL

Table 3. Discretization of Process Parameter Space where Simulations were Run

Model Assumptions

number process parameter of levels

Adsorption of Pd species in packed bed columns can be described via macroscopic and microscopic mass balances (i.e., at the column and scavenger particle levels, respectively).42 The model assumptions are that (a) the species concentration profile is parabolic within the scavenger particles (microscopic level),38 (b) a linear driving force represents mass transfer of species between the bulk solution and the scavenger particles,38 (c) the equilibrium for multi-component adsorption can be represented by a Langmuir-type adsorption isotherm for competitive adsorption,43 (d) for simplicity, no radial or axial temperature profiles are assumed, (e) there is single liquid phase stream, and (f) the packing material is monodisperse or it exhibits a narrow particle size distribution.


discrete points

9

[0.97,0.99,1.01,1.03,1.05,1.07,1.09,1.12,1.20]

6

[0.001,0.0015,0.00235,0.0025,0.00265,0.003]

5

[60,120,180,240,300]

5

[22,32,40,52,64]

7

[10,20,50,100,200,220,250]

4

[0.05,0.1,0.15,0.2]

Model Equations

Eq 28 describes the mass balance in the packed bed column with boundary conditions given by eqs 29 and 30, and an initial condition given by eq 31. Eqs 32 and 33 represent the mass balance at the scavenger particle level, and eq 34 represents the multi-component adsorption isotherm. ∂C Pd(i)(z , t ) ∂t

=

Dax,C ∂ 2C Pd(i)(z , t ) 2

ε

−

u ∂C Pd(i)(z , t ) ε ∂z

∂z ̅ i)(z , t ) 1 − ε ∂qPd( − ∀ i = 0, II ε ∂t

C Pd(i)(0, t ) = C Pd(i),feed ∀ i = 0, II

(29)

= 0 ∀ i = 0, II

∂z

=

6k f,Pd(i) Dpart

Dpart k f,Pd(i) 10 Deff,Pd(i)

(32)

finer discretization 1

finer discretization 2

1.02−1.12

1.04−1.1

1.01−1.12

0.002−0.0027

0.0021−0.0026

0.002−0.0026

134−226

152−208

110−230

42−62

46−58

42−64

36−164

61−139

10−139

0.06−0.14

0.075−0.125

0.05−0.14

flexibility region

[C Pd(i)(z , t )

* i)(z , t )] ∀ i = 0, II − C Pd(

* (z , t ) = qPd( i)

(31)

[C Pd(i)(z , t ) − C*Pd(i)(z , t )]

∀ i = 0, II * − qPd( ̅ i)(z , t ) = qPd( i)


(30)

C Pd(i)(z , 0) = C Pd(i),initial ∀ i = 0, II ∧ 0 < z ≤ LC

∂t

Table 4. Details on the Determination of the Flexibility Region for Suzuki Coupling Reaction

(28)

process parameter

∂C Pd(i)(LC, t )

∂qPd( ̅ i)(z , t )

A further two-dimensional finer grid was generated for the Pd catalyst/SM2 ratio and reaction time, and the probability of passing quality attributes was interpolated to the finer grid using the grid data function in MATLAB. This was done at every combination of other process parameters. Applying the procedure outlined in the section on calculation of the flexibility region resulted in an estimated box shown in column “finer discretization 1” of Table 4. Because the upper

* i)(z , t ) qm,Pd(i)KL,Pd(i)C Pd( * j)(z , t ) 1 + ∑j KL,Pd(j)C Pd(

bound for O2 in head space was estimated to be 164 ppm and fell between the discretization points 100 and 200 ppm, the same procedure was repeated, but instead of reaction time and catalyst ratio, the finer two-dimensional grid was based on temperature and O2 in head space. The estimated box from this is shown in the column “finer discretization 2” of Table 4. The upper bound of 139 ppm was a better estimate for O2 head space than 164 ppm. These two box ranges were used to manually construct the box that represents the flexibility region. The manual changes consisted of widening the range of each of the process parameters one at a time to verfiy if all points within the box still pass.

(33)

∀ i = 0, II (34)

Except for the void fraction, external mass transfer coefficients, and axial dispersion, all other model parameters depend on intrinsic characteristics of the adsorption process and are inherently scale independent. Moreover, they can be easily determined via well-designed lab-scale batch adsorption experiments. J



■

APPENDIX D

■

AUTHOR INFORMATION

Article

Dpart: HO2:

scavenger average particle size, m Henry constant of oxygen in a mixture of water/ THF, Pa L−1 mol−1 HO2,j: Henry constant of oxygen in pure solvent j, Pa L−1 mol−1 kf,Pd(i): Pd species i external mass transfer coefficient, m/s KL,Pd(i): affinity of Pd species i, g of solution/g of Pd L C: packed bed column length, m mi : mass of species i, g Mi: molecular weight of species i, g/mol P: pressure, Pa qm,Pd(i): maximum adsorbent capacity for Pd species i, g of Pd/g of scavenger qP̅ d(i) (z,t): concentration of the Pd species i adsorbed in the resin, g of Pd/g of scavenger q*Pd(i) (z,t):e Equilibrium concentration of the Pd species i adsorbed in the, g of Pd/g of scavenger t: elapsed time, sec T: temperature, K u: superficial fluid velocity (u = Q/A, where Q is the volumetric flow rate and A is the column crosssectional area)

Corresponding Author

*Phone:+1-317-651-5233. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors would like to acknowledge the help of several colleagues: David Varie and Eric Crick for helpful discussions, and Thomas Wilson, Michael Kobierski, Allison Fields, and Jennifer Groh for their invaluable help gathering experimental data. Also, the authors would like to acknowledge the development of an adsorption column/Suzuki reaction simulation tool by the RES Group Inc. (USA).

■

ABBREVIATIONS CPd(i)(z,t): concentration of Pd species i in solution, g of Pd/g of solution C*Pd(i) (z,t): equilibrium concentration of Pd species i in solution, g of Pd/g of solution Dax,C: axial dispersion in the packed bed column, m2/s Deff,Pd(i): intraparticle diffusivity Pd species i, m2/s K


Organic Process Research & Development V: xi: yO2:

Article

(29) Grossmann, I.; Floudas, C. Comput. Chem. Eng. 1987, 11, 675. (30) Grossmann, I.; Halemane, K.; Swaney, R. Comput. Chem. Eng. 1983, 7, 439. (31) Morari, M. Comput. Chem. Eng. 1983, 7, 423. (32) Grossmann, I.; Calfa, B.; Garcia-Herreros, P. Comput. Chem. Eng. 2014, 70, 22. (33) Stockdale, G.; Cheng, A. Quality Technology & Quantitative Management 2009, 6, 391. (34) Adrio, L.; Nguyen, B.; Guilera, G.; Livingston, A.; Hii, K. Catal. Sci. Technol. 2012, 2, 316. (35) Brazier, V.; Nguyen, B.; Adrio, L.; Barreiro, E.; Leong, W.; Newton, M.; Figueroa, S.; Hellgardt, K.; Hii, K. Catal. Today 2014, 229, 95. (36) Adamo, C.; Amatore, C.; Ciofini, I.; Jutand, A.; Lakmini, H. J. Am. Chem. Soc. 2006, 128, 6829. (37) Mondal, B.; Wilkes, R.; Percy, J.; Tuttle, T.; Black, R.; North, C. Dalton Trans. 2014, 43, 469. (38) Phillips, S.; Kauppinen, P. Platinum Met. Rev. 2010, 54 (1), 69. (39) ICH Q8, Q9, and Q10 Questions and Answers, ICH HARMONISED TRIPARTITE GUIDELINE, 2012. (40) Puri, P.; Ruether, J. Can. J. Chem. Eng. 1974, 52, 636. (41) Battino, R.; Rettich, T.; Tominaga, T. J. Phys. Chem. Ref. Data 1983, 12, 163. (42) Ö zdural, A.; Alkan, A.; Kerkhof, P. J. Chromatogr. A 2004, 1041, 77. (43) Butler, J.; Ockrent, C. J. Phys. Chem. 1929, 34, 2841.

reaction volume, L molar fraction of species i, dimensionless molar fraction of oxygen in head space, dimensionless axial position, m

z:

Greek Symbols

ε: packed bed column void fraction, dimensionless ρi: density of species i, g/L Subscripts

0: K3PO4: Pd(0): Pd(II): SM1: SM2:

■

initial condition (t = 0) potassium phosphate Pd(0) species Pd(II) species boronate ester organohalide

REFERENCES

(1) ICH. Guidance for Industry Q8(R2) Pharmaceutical Development. ICH HARMONISED TRIPARTITE GUIDELINE, 2009. (2) ICH. Development and manufacture of drug substances (chemical entities and biotechnological/biological entities) Q11. ICH HARMONISED TRIPARTITE GUIDELINE, 2012. (3) Huang, J.; Kaul, G.; Cai, C.; Chatlapalli, R.; Hernandez-Abad, P.; Ghosh, K.; Nagi, A. Int. J. Pharm. 2009, 382, 23. (4) Thirunahari, S.; Chow, P.; Tan, R. Cryst. Growth Des. 2011, 11, 3027. (5) Kumar, S.; Gokhale, R.; Burgess, D. J. Int. J. Pharm. 2014, 464, 234. (6) Charoo, N. A.; Shamsher, A. A. A.; Zidan, A. S.; Rahman, Z. Int. J. Pharm. 2012, 423, 167. (7) Brueggemeier, S.; Reiff, E.; Lyngberg, O.; Hobson, L.; Tabora, J. Org. Process Res. Dev. 2012, 16, 567. (8) Prpich, A.; am Ende, M.; Katzschner, T.; Lubczyk, V.; Weyhers, H.; Bernhard, G. Comput. Chem. Eng. 2010, 34, 1092. (9) Lepore, J.; Spavins, J. J. Pharm. Innov. 2008, 3, 79. (10) Freeman, J.; Simoncelli, E. Nat. Neurosci. 2011, 14, 1195. (11) Cox, D.; Dean, T. Curr. Biol. 2014, 24, R921. (12) Liu, S.; Cui, W.; Wu, Y.; Liu, M. Vis. Comput. 2014, 30, 1373. (13) Stephanopoulos, G.; Reklaitis, G. Chem. Eng. Sci. 2011, 66, 4272. (14) Oksanen, C.; Garcia-Munoz, S. Comput. Chem. Eng. 2010, 34, 1007. (15) Spooner, J. A review of computer process simulation in industrial pollution prevention. EPA 600/R-94/128. Aug 1, 1994. (16) Various Editors. CAST Newsletter 1976 to 2000. CAST Division of AICHE. (17) Gernaey, K.; Cervera-Padrell, A.; Woodley, J. M. Comput. Chem. Eng. 2012, 42, 15. (18) Gernaey, K.; Gani, R. Chem. Eng. Sci. 2010, 65, 5757. (19) Cervera-Padrell, A.; Skovby, T.; Kiil, S.; Gani, R.; Gernaey, K. Eur. J. Pharm. Biopharm. 2012, 82, 437. (20) Buzzi-Ferraris, G.; Manenti, F. Chem. Eng. Sci. 2009, 64, 1061. (21) Bard, Y. Nonlinear Parameter Estimation; 1st ed.; Academic Press: New York, 1974. (22) Eckhardt, R. Los Alamos Sci. 1987, 14, 96. (23) Lebrun, P.; Krier, F.; Mantanus, J.; Grohganz, H.; Yang, M.; Rozet, E.; Boulanger, B.; Evrard, B.; Rantanen, J.; Hubert, P. Eur. J. Pharm. Biopharm. 2012, 80, 226. (24) Chatzizacharia, K.; Hatziavramidis, D. Ind. Eng. Chem. Res. 2014, 53, 12003. (25) Castagnoli, C.; Yahyah, M.; Cimarosti, Z.; Peterson, J. J. Org. Process Res. Dev. 2010, 14 (6), 1407. (26) Figueroa, I.; Vaidyaraman, S.; Viswanath, S. Org. Process Res. Dev. 2013, 17, 1300. (27) Grossmann, I.; Morari, M. Proc. 2nd Int. Conf., FOCAPD, CACHE, 1984, p 937. (28) Dimitriadis, V.; Pistikopoulos, E. Ind. Eng. Chem. Res. 1995, 34, 4451. L


Definition of Design Spaces Using Mechanistic Models and Geometric

Recommend Documents