Design of Experiments for Modeling of Fermentation Process Characterization in Biological Drug Production

: Biological products are increasingly important, and therefore the industry has begun to adopt quality by design, as recommended by the ICH and the U.S. FDA. Smaller companies, however, have faced difﬁculties in employing full-scale experiments or the quality by design strategy. Thus, this study provides an alternative way to build a model from existing data with experimental software that does not require full-scale experiments. This empirical study hopes to provide a practical way to improve the efﬁciency of smaller biopharmaceutical companies and researchers. Moreover, the models provided here can be applied to process characterization in recombinant protein production.


Introduction
Recently, the pharmaceutical industry has been advised by the ICH [1] and U.S. FDA to adopt a new quality control strategy. The U.S. FDA provides new guidelines to assist the industry under the new strategy, known as quality by design (QBD, ICH Q8 R2) [1][2][3][4]. An increasing number of new drug applications are based on QBD [3].
During the past decade, the number of biological products has been increasing. The target product of this study is a recombinant biological drug from Green AbioTechnology Co., Ltd., in Taiwan. The application for this drug will be sent to the U.S. FDA. The traditional manufacturing process of a biological product often utilizes organisms in the first step to overexpress the product. In the follow-up downstream process, the product is purified. There are several developmental stages before the final commercial manufacturing process can be determined [2]. Often, lab-scale, process characterization, qualification, and final validation studies are necessary.
In the early stages before process validation, manufacturers in the biopharmaceutical industry often engaged in troubleshooting with an existing process that was built by a traditional approach, known as one-factor-at-a-time (OFAT). Few researchers are aware that the recently introduced quality by design (QBD) approach can help in this. In the QBD approach, a useful tool, the design of experiments (DOE), can find the direction and optimized parameters. DOE and model building can potentially help [5,6], not only in terms of time-efficiency, but also in terms of the completeness of the science (design space) [7,8].
A central concept of the principle is to systematically correlate the relationships between the process parameters and the quality attributes [8]. This also implies that the empirical relationships and formulas can be built to assist researchers throughout the process development stages [2,8,9]. In the process of development, continuous improvement of the process parameters is crucial [2]. Model building (during the development stages) can help in determining the next direction or experiments in the beginning, and eventually help with optimizing the current parameters [4]. In the late developmental stages, good control of the quality attributes ensures an acceptable final quality [8]. Finally, an overall picture of the whole design will ensure filing and yield a better quantity or quality of the product [2,8]. The U.S. FDA has granted regulatory flexibility for analytical methods applying QBD [4].
In the early developmental stages and the scale-up stages, however, using a full-scale DOE for the process parameters is usually not a suitable plan. The fermentation process may take 1.5 weeks [10], and the follow-up purification process may take another 1.5 weeks (a recombinant protein produced by Escherichia coli). Thus, it is likely infeasible to fully investigate a fermentation process with five critical process parameters, especially if the initial operation space is not well-defined (due to the possibility of out-of-specification (OOS) events, which waste more time). Thus, the minimum number of experiments is assumed to be on the order of 10, and the full set will take about 15 to 25 weeks (3 to 6 months). Often, such an experiment involves about 10 steps, dozens of quality attributes, and some 50 process parameters, or even dozens of critical process parameters across the steps. For biopharmaceutical researchers in smaller companies or those with very limited budgets and time, full-scale DOE can still be time-consuming and resource-intensive to implement. Moreover, in reality, during this transition stage (to QBD), the first version of the process parameters is built with the OFAT approach. A concept of quality by design and the DOE built on the existing data would be useful for researchers. In this study, DOE researchers can only develop as few as three additional datasets to optimize (at least locally) the volumetric yield and total yield. This approach helps the entire pharmaceutical industry to make more drugs, which have been eagerly requested by patients. This paper explores the suitability and provides a case study on a biological drug developed using this new concept.
In DOE analysis, researchers provide input and collect responses. Usually, the volumetric yield and the final yield are prioritized. Thus, they were chosen as the responses here.
In lab-scale studies (presented below), the research and development studies on quality attributes followed ICH Q8 (R2) guidelines, and the analytical methods for determining purity and activity in the final lot-release specification were pre-validated for at least linearity, specificity, and repeatability. All the quality attributes in the quality target product profile were fully qualified within the lot-release specification after the scale-up batch was finally confirmed. The analytical methods for the lot-release specification were also validated.

Chemicals, Bacterial Strains, and Vectors
The chemicals listed here are the ones used in the lab-scale study. The suppliers of the chemicals used on a larger scale (such as in the pilot plant) could be different. Isopropyl β-D-thiogalactoside (IPTG) was purchased from MDBio Inc. (Taipei, Taiwan). The E. coli strain BL21 (DE3) was purchased from Novagen (Madison, WI). Expression vector pET19b carrying cDNA encoding a biological drug, MI001-S, was constructed in the lab.

Protein Expression, Fermentation, and Purification
The fermentation was performed in a Firstek fermenter (FB-10B or an equivalent that works at 5-L, Firstek Scientific, Taiwan). The cells were harvested by centrifugation at 9000× g for 30 min at 4 • C (pre-cooled). If the cells were not lysed within 12 h, the pellets were refrigerated at −20 • C. Cell lysis was performed by a sonicator (for a smaller volume at or below 0.3 L) or a high-pressure homogenizer at around 85 MPa with the temperature controlled at 10 • C. The process was repeated three times. The initial solution temperature before each lysis process was around 10 • C. The lysis buffer was 0.1 M Tris, 0.3 M NaCl, and 8 M urea at pH 8. The lysate was centrifuged at 30,000× g for 30 min before being filtered through a 0.2 um filter. The purification process required three columns. The filtrate was then dialyzed against the same buffer without 8 M urea. Alternatively, for a volume larger than 1 L, the filtrate was directly loaded into the column, and then three column volumes of buffer were used to slowly elute the buffer without the urea for 3 h. Buffers were temperature controlled to less than 10 • C. The first column was the Ni-NTA affinity column with buffer. The equilibrium buffer contained 0.02 M imidazole. The elution buffer contained 0.3 M imidazole. The next column was an ion-exchange column with an initial NaCl salt concentration of 0.2 M and elution buffer with a salt concentration of 1 M. The follow-up purification was performed with an RP-HPLC machine equipped with a column packed with C18 resin. The equilibrium buffer was water mixed with 0.1% trifluoroacetic acid. The elution buffer contained water, trifluoroacetic, and a high concentration of acetonitrile. The slow linear gradient required three column volumes. The purity and protein concentration were calculated from the UV-VIS absorbance ratio of A280 nm [11] The extinction coefficient was predicted by the ExPASy ProtParam program [11,12].

Scaling Up Attempts
The current fermentation parameter settings were transferred and scaled up to a 250 L scale fermenter, and the final products of the batches were successful in terms of quantity and quality. The results indicated that the full set of fermentation process parameters could be scaled up. In addition, the scaled-up purification process was also successful.

Model Building in Design Expert
The predictive model developed during the optimization process was built in the Design Expert software (Minneapolis, MN) version 13.0.8.0. There were two types of models. One was a coded model that normalized the maximum of each variable as +1, the minimum as −1, and the middle point as zero. This model was suitable for visualizing the relative effects across the variables. For the actual model, all the variables were in units. The criteria to accept the most suitable model were the corrected Akaike information criterion (AICc) and the Bayesian information criterion (BIC, priority in this study) [13]. The procedure is shown below, including the steps of model building, selection, and validation:
Use the default regular two-level factorial design to create the table; 3.
Choose the number of variables (insert additional runs if needed); 4.
Model building: generate both the factorial design model and the polynomial design models; 5.
Select candidate models aiming at low AICc and BIC values; 6.
Validation: compare the models by high R 2 , high adjusted R 2 , low predicted residual error sum of squares (PRESS), low p-values in variables, high p-values in lack of fit, low Std. Dev. (also called root mean square error, RMSE) and low (relative) coefficient of variation (C.V.%), monitor summary of fit and ANOVA; 7.
Perform model diagnostic analysis to study the residuals.

Model Building in JMP
The additional predictive model was built in the JMP software (SAS Institute Inc., Cary, NC, USA) trial version 16.1.0. The most suitable model was selected by comparing all the models with the intercept, the main effect, and the primary interactions. The lowest or the next lowest BIC (priority) and AIC models were selected. Furthermore, the root mean square error (RMSE) and the coefficient of determination (R 2 , and adjusted R 2 , typically larger than 0.85) were used to select the best model. We also considered the analyses listed in the previous section, including number 6, validation. After the selection of the best model, the prediction module from the profiler function of the JMP software could be used later, after the model reached an optimized and stable region. For the simulation study, the Monte Carlo simulation approach was conducted in 100,000 simulation runs. The simulation provided a tolerance interval of the process parameters. This function could be useful for researchers.

The Existing Data and the Analyses
This work aimed to bridge the gap between traditional OFAT and the (full-scale) DOE approach. We started with the fermentation record from the first 8 runs (batches), which is presented in Table 1. Table 1. Experimental data based on OFAT approach (runs 1 to 8) and three consecutive confirmation runs (runs 9, 10, and 11). The optical density at 600 nm is defined as OD. All the data were collected using a fed-batch fermenter described in the Materials and Methods section. Thus, the first eight data points were initial scans across three sets of temperature and optical density (OD) values. The latter three were then executed as the final batches. On the basis of the initial search of eight points, we found that the preliminary models suffered from noise and the volumetric yield model was slightly better. The ANOVA analyses indicated that the volumetric yield model had a low p-value at 0.043, and the total yield had a moderate p-value at 0.06. The analyses indicated that the volumetric yield had a significant variable of A = [IPTG] at p-value = 0.003. The RMSE was 0.033, average was 0.226, and C.V.% was 14.68. In contrast, the total yield model returned a RMSE of 106.2, average of 575.8, and C.V.% of 18.45, and the B = temperature variable was moderately significant with p = 0.05. Thus, we looked at the numerical optimization analysis and Figure 1. Both models were in agreement that the optimization direction was (red) color-coded to be at a lower temperature and have a higher inducer concentration. The coded (normalized) equation of the volumetric yield is shown in Equation (1), where A = [IPTG] and B = temperature.

Optimization and Verification
On the basis of the analysis with the first eight data points, we found that there were two optimization attempts. The verifications are described as follows. In the first, condition 1, the temperature was 18 • C, and the inducer concentration was 1 mM. In the second, condition 2, the temperature was 17 • C, and the inducer concentration was 2 mM. Condition 1 was executed as runs 9, 10, and 11. The quality data confirmed that these conditions worked (Figure 2).

Optimization and Verification
On the basis of the analysis with the first eight data points, we found that there were two optimization attempts. The verifications are described as follows. In the first, condition 1, the temperature was 18 °C, and the inducer concentration was 1 mM. In the second, condition 2, the temperature was 17 °C, and the inducer concentration was 2 mM. Condition 1 was executed as runs 9, 10, and 11. The quality data confirmed that these conditions worked ( Figure 2).
On the basis of the preliminary data, we found that condition 2 returned a higher quantity than condition 1 on the lab scale (a smaller fermenter). Then, both were examined in a CMO plant with a larger fermenter. Both conditions yielded better total product quantities, with 25% increases over run 5 (Table 1). However, the total yield of condition 2 was similar to that of condition 1 in the CMO fermenter. In summary, condition 2 returned a higher yield than condition 1 in a smaller fermenter but returned a similar quantity in a large fermenter.
The results may indicate differences between large and smaller fermenters (i.e., the scale-up parameters might not be suitable in this condition). In addition, the temperature of around 17 °C was a few degrees above the limit of the cold-water system, and therefore the process took more time. In the end, condition 2 was not used because of the higher cost of the inducer and longer time before induction (lower temperature and the limit of the system).

Final Set of Data and the Characterization of the Model
On the basis of the full data set of 11 runs, we built both of the models. The ANOVA analyses indicated that the total yield model had a low p-value of 0.0013. The RMSE was 89, average was 695, and C.V.% was 12.8. The R 2 was 0.93, and the adjusted R 2 was 0.88. The analyses indicated that the total yield had a significant variable of B × C = temperature On the basis of the preliminary data, we found that condition 2 returned a higher quantity than condition 1 on the lab scale (a smaller fermenter). Then, both were examined in a CMO plant with a larger fermenter. Both conditions yielded better total product quantities, with 25% increases over run 5 (Table 1). However, the total yield of condition 2 was similar to that of condition 1 in the CMO fermenter. In summary, condition 2 returned a higher yield than condition 1 in a smaller fermenter but returned a similar quantity in a large fermenter.
The results may indicate differences between large and smaller fermenters (i.e., the scale-up parameters might not be suitable in this condition). In addition, the temperature of around 17 • C was a few degrees above the limit of the cold-water system, and therefore the process took more time. In the end, condition 2 was not used because of the higher cost of the inducer and longer time before induction (lower temperature and the limit of the system).

Final Set of Data and the Characterization of the Model
On the basis of the full data set of 11 runs, we built both of the models. The ANOVA analyses indicated that the total yield model had a low p-value of 0.0013. The RMSE was 89, average was 695, and C.V.% was 12.8. The R 2 was 0.93, and the adjusted R 2 was 0.88. The analyses indicated that the total yield had a significant variable of B × C = temperature × OD at p = 0.034 ( Table 2). The coded (normalized) equation of the volumetric yield is shown in Equation (2), where A is [IPTG], B is temperature, and C is OD.
The volumetric yield model returned a low p-value of 0.0002. The RMSE was 0.036, average was 0.282, and C.V.% was 12.6. The R 2 was 0.93, and the adjusted R 2 was 0.90. The analyses indicated that the volumetric yield had significant variables of A = [IPTG] at p = 0.003, B = temperature at p = 0.017, and A × B at p = 0.013 (Table 3). The coded (normalized) equation of the volumetric yield is shown in Equation (3), where A is [IPTG] and B is temperature.
Because the highest total yield and the volumetric yield were achieved in runs 9, 10, and 11, then using the Point Prediction function in the "Post Analysis" directory, we were able to predict the total yield to be on average 1002 ± 89 (885 to 1120, at 95% C.I., Figure 2a). Similarly, the volumetric yield was predicted to be 0.426 ± 0.036 (0.382 to 0.470 at 95% C.I., Figure 2b).

Final Set of Data Executed by JMP
On the basis of the full dataset of 11 runs, we also built both the total yield and volumetric yield models with JMP software. After similar procedures, the analyses indicated that both JMP and Design Expert software returned the same final total yield and volumetric yield models. The ANOVA analyses of the total yield are listed in Table 2, and those of the volumetric yield are listed in Table 3. The coded formula for total yield was the same as Equation (2). The coded formula for volumetric yield was the same as Equation (3). The plot of actual versus predicted values is shown in Figure 3.   The highest total yields and volumetric yields were those from runs 9, 10, and 11. Thus, the optimization focused on condition 1, with an inducer concentration of 1 mM and a temperature of 18 • C. The average of the simulations of the total yield was 1035.3 (927.7 to 1142.9). The average of the simulations of the volumetric yield was 0.425 (0.390 to 0.461). Data and distribution are presented in Figure 4.

Discussion
Software comparison: This study was analyzed in Design Expert and in JMP. Both programs returned the same models. The differences were in the optimization in the final stage. In Design Expert, a single average with a standard deviation was calculated. In JMP, we set the conditions such that the S.D. of inducer concentration, induction temperature, and induction OD were 1.0 ± 0.1 mM, 18 ± 1 °C, and 50 ± 10, and the additional distribution can be seen in Figure 4. We found these two kinds of software to be similar to the point of equivalence on the basis of this study.
After the models were built on the basis of eight data points, the Design Expert simulations indicated that, for an inducer concentration of 1 mM, induction temperature of 18 °C, and induction OD of 50, the predicted total yield would be 998 ± 110 and the predicted volumetric yield would be 0.40 ± 0.03. The later runs validated the predictions.
In this study, quadratic terms were forced in both software programs. None of the attempts returned significant results, and the ANOVA results were no better than those of these simpler models. Thus, the models were not included.
The final conditions for this drug were found to be 1.0 mM, 18 °C, and 50 induction OD. The expected total yield was simulated to be around 1020, and the volumetric yield, around 0.43. According to experience, the operation space of the inducer concentration

Discussion
Software comparison: This study was analyzed in Design Expert and in JMP. Both programs returned the same models. The differences were in the optimization in the final stage. In Design Expert, a single average with a standard deviation was calculated. In JMP, we set the conditions such that the S.D. of inducer concentration, induction temperature, and induction OD were 1.0 ± 0.1 mM, 18 ± 1 • C, and 50 ± 10, and the additional distribution can be seen in Figure 4. We found these two kinds of software to be similar to the point of equivalence on the basis of this study.
After the models were built on the basis of eight data points, the Design Expert simulations indicated that, for an inducer concentration of 1 mM, induction temperature of 18 • C, and induction OD of 50, the predicted total yield would be 998 ± 110 and the predicted volumetric yield would be 0.40 ± 0.03. The later runs validated the predictions.
In this study, quadratic terms were forced in both software programs. None of the attempts returned significant results, and the ANOVA results were no better than those of these simpler models. Thus, the models were not included.
The final conditions for this drug were found to be 1.0 mM, 18 • C, and 50 induction OD. The expected total yield was simulated to be around 1020, and the volumetric yield, around 0.43. According to experience, the operation space of the inducer concentration was 1.0 ± 0.1 mM, the induction temperature was 18 ± 1 • C, and induction OD was 50 ± 10 ( Figure 4). Studies including simulations may serve as supporting data for the determination of the (putative) critical process parameters and development of the fermentation process.
In the downstream process, after three purification steps, the lot-release specification of the final drug substance was around 98 ± 1 (%) in purity. After the first step of purification, the purity was low, with very large uncertainty. On the basis of the preliminary results of the scaling-up batches (data not shown), we found that the purity with a fixed condition after the first column was around 40% to 50%. Thus, the contribution of the fermentation process parameters to the quality of the product seems difficult to detect. Models will suffer from noise. After three steps of purification, the contributions were minimized. Thus, the quality attributes were assumed not to be correlated to the current three process parameters in the fermentation stage.
Provided that more data are analyzed, we can then know whether the simplest models can be applied to other systems. In addition, we hope that the above strategy or models can be applied to other recombinant protein systems.

Conclusions
This case study shows that linear regression models built on OFAT data can help researchers easily achieve optimization. From the addition of a few points, the volumetric yield model gained low p-values of 0.0002, with only main effects and primary interaction terms. The variables all significantly contributed to the model. The analysis can serve as supporting data for the determination of the operating space of the process parameters in future drug applications.