DESIGN OF EXPERIMENTS MODEL FOR THE OPTIMIZATION OF SILK FIBROIN BASED NANOPARTICLES

Objective: Silk fibroin based nanoparticles have been utilized extensively in biomedical fields. Amongst many preparation methods, desolvation is a favorable one. However, this method yields nanoparticles with unpredictable parameters. Thus, this investigation aimed to systematically study the effects of three independent variables including fibroin concentration (% w/v, X1), volume ratio between fibroin solution and ethanol (X2), formulation time (h, X3) on three main responses, particle size (nm, Y1), polydispersity index (Y2), zeta potential (mV, Y3). Methods: Fibroin was extracted from degummed Bombyx mori silk. The fibroin calibration curve was constructed by UV-spectrophotometer at 276 nm. The nanoparticles were prepared using the desolvation method of aqueous fibroin solution in ethanol. Design Expert® software was used to design the model. The mean particle size, polydispersity index and zeta potential were determined using ZetaPALS®analyzer. Results: By using D-optimal design with the quadratic model, the results showed that all X1, X2, and X3 variables had significant impacts on the fibroin nanoparticles characteristics Y1, Y2, and Y3. The generated model was also validated and demonstrated to be solid and reliable. The obtained optimal nanoparticles possessed Y1 of 238.1 nm, Y2 of 0.12, and Y3 of -21.78 mV, which were in agreement with the predicted values, 224.8 nm, 0.13 and -19.31 mV, respectively. The optimal actual and theoretical particle characteristics were correlated with a desirable value of R2 = 0.8770. Conclusion: The D-optimal design proved its effectiveness in the prediction and optimization of fibroin nanoparticle properties.


INTRODUCTION
Nanoparticles have gained increasing attention in biomedical fields, especially in the area of cosmeceuticals and drug delivery systems [1][2][3]. Numerous types of nanoparticles have been proposed. They could be classified into three main categories; lipid nanoparticles such as liposomes and solid lipid nanoparticles; polymeric nanoparticles such as natural nanoparticles (i.e., chitosan, alginate, albumin) and synthetic nanoparticles (i.e., polyacrylamide, poly (lactic acid)); and inorganic nanoparticles such as carbon nanotubes and magnetic nanoparticles [2]. Amongst them, polymeric nanoparticles are favorable due to ease of preparation, cheap ingredients and the ability to modify the surface to maximize the desired outcomes [4]. However, the synthetic polymers and their metabolites might cause toxicity to both the human and the environment. Thus, to avoid toxicity issue, natural polymers have been extensively explored.
Fibroin is a natural protein extracted from silk commonly originated from Bombyx mori silkworm. Recently, fibroin has been extensively investigated due to its outstanding properties such as high tensile strength, biocompatible, biodegradable and non-toxic [5]. Fibroin nanoparticles (FNP) could be formulated easily by many methods, such as ball milling, emulsification, spray drying, sonication, desolvation, polyelectrolyte complexation, as well as a supercritical fluid technique [5,6]. The most popular method for preparing FNP in the drug delivery area is desolvation due to its low cost, mild condition, small particle size and ease of operation [7]. Nevertheless, the obtained particles showed significant variation in surface zeta potential, size and homogeneity (polydispersity index, PI) [7]. These variations were mainly dependent upon the type of solvents, formulating conditions, fibroin concentrations, and the volume ratio between fibroin and an organic solvent. Thus, to be able to manifest and predict the FNP properties, it is essential to have a clear understanding of how the preparation conditions affect the FNP characteristics, especially on the size, zeta potential, and PI.
In general, evaluating the effect of a large number of preparation variables demands an appropriate high amount of experiments, which are costly and time-consuming. To minimize the unnecessary waste of resources, the design of experiments (DOE) is used. Mainly, DOE could be utilized for three primary objectives including screening, optimization, and robustness testing. Amongst various DOE types such as full factorial, fractional factorial, D-optimal, Box-Behnken, and composite, the D-optimal design is usually used for the optimization process with mixture factors of qualitative and quantitative (continuous) ones [8]. A D-optimal design is a computer-generated design which composes of the best subset of experiments selected from the variables set [8]. When compared to the traditional full factorial design, D-optimal could significantly reduce the number of required experiments, especially with the continuous factors such as concentrations. Thus, a more suitable model is likely to be established due to the decrease of bias, which possibly comes from the limited chosen levels of the factorial design.
Therefore, in this research, the D-optimal design was used as a tool to systematically study the effects of the preparation conditions on the size, PI, and zeta potential of FNP. The quadratic model was used for the best-fit relationship between variables and responses. Furthermore, based on the optimized equation, the optimized FNP was prepared, and the desirability of the model was determined.

Fibroin extraction
The processes of silk degumming, fibroin extraction and purification were modified from Min et al., method [9]. In brief, the cocoons were

Fibroin calibration curve
The freshly prepared lyophilized fibroin was dissolved in DI water and its absorbance was determined using UV-spectrophotometer (Genesis 10S, Thermo Fisher, USA) at 276 nm. Then, the calibration curve was constructed with concentration ranging from 0.1 to 1.0 mg/ml.

Nanoparticles formulation
FNP was prepared using desolvation method. In brief, lyophilized fibroin was dissolved in DI water to get fibroin solution at various concentrations. Then, it was slowly injected into absolute ethanol at various volume ratio between fibroin solution and ethanol. The spontaneously formed translucent mixture was left to stabilize at 4 °C for two investigate time points of 4 h and 24 h. The particles were then centrifuged at 31514×g (Mikro 220R, Hettich, Germany) for 30 min, re-dispersed in DI water by sonication (40% amplitude, 30 s), washed thrice with DI water and lyophilized at -55 °C for 72 h. The products were stored at 2-8 °C for further characterization.

Experimental design
Using the response surface methodology, the D-optimal design was conducted with two continuous factors and one qualitative factor. The quadratic model was used for the best-fit relationship between variables and responses, and the optimization of FNP. The Design Expert ® software (Version 11.0, Stat-Ease Inc., USA) was utilized for the designation, wherein the lack of fit test, the predicted and adjusted R 2 values were generated. Three independent variables were chosen, based on the formulation process, including fibroin concentration (% w/v, X1), the volume ratio between fibroin solution and ethanol (X2), and formulation time (h, X3). The two variables X1 and X2 were continuously ranged from 0.25% to 2.00% and from 0.33 to 2.00, respectively. Whereas the formulation time X3 was set at 4 h and 24 h. Particles size (nm, Y1), polydispersity index (PI, Y2), and zeta potential (mV, Y3) were selected as the response variables. Statistical polynomial equations were utilized to evaluate the formulation responses, as follows: Where Yi is the response variable, b0 is the intercept, b1 to b8 are the coefficient for the factor X1, X2, X3 and their combinations. 3D surface plots and cube fig. were used to demonstrate these equations.
To analyze and validate the response surface curves and the polynomial equations, the same software was used. The acceptance criteria were set at p-values of less than 0.05, the difference between predicted and adjusted R 2 values of less than 0.2. For the optimal product, the constraints were set at "minimum" for all three responses. Finally, based on the optimized equations, we formulated the formulations using theoretical parameters, and the desirability of the model was determined.

Particles size and charge
The mean particle size and size distribution (PI) were determined by the dynamic light scattering (DLS) method using ZetaPALS ® analyzer (Brookhaven Instrument Corporation, USA). The instrument was equipped with a 35 mW helium-neon laser diode operating at 632.8 nm and a BI-200SM Goniometer connected to a BI-9010AT digital correlator. Samples were diluted with DI water in a disposable polystyrene cuvette with a 10-mm path length. The measurement was performed three times at 25 °C by auto measuring mode at a fixed angle of 90 °.
The zeta potential was determined by phase analysis light scattering method, using the same instrument. The measurement was performed at 14.8 ° to the incident light. All samples were diluted with DI water in a volume ratio of 1:3 in 10-mm cuvette. Data were collected for 10 cycles and each sample was determined in triplicate. The zeta potential was calculated from the electrophoresis mobility based on the Smoluchowski equation included in the system software.

Statistical analysis
All experiments were performed in triplicate to confirm the results.
For quantitative experiments, the mean±SD (standard deviation) was reported. Analysis of variance (ANOVA) was used for statistical purposes, with p<0.05 for significant comparisons.

RESULTS AND DISCUSSION
Fibroin nanoparticles have been extensively utilized in biomedical fields due to their outstanding safety, biocompatibility, and biodegradability [7,[10][11][12][13]. Nevertheless, the systematic study on the designation and optimization of these particles is absent. Thus, in this report, we aimed to use DOE, focused on D-optimal design, to critically investigate the relationships between formulating parameters and the particle properties such as size, homogeneity, and zeta potential. D-optimal design was used due to its benefits in optimizing processing parameters composed of both qualitative and quantitative factors [8], which is superior to traditional factorial designs. Moreover, the model was validated, and the optimal formula was proposed.

Fibroin calibration curve
Previous studies have demonstrated that desolvation method was effective in fabricating FNP by conformational changes from watersoluble amorphous silk I to water-insoluble crystalline silk II [10][11][12][13][14]. However, the most popular available silk form is silk II (i.e., silk fiber) and silk I could only be found naturally in the silkworm gland. Thus, to utilize this preparation method, silk fibers need to undergo a dissolving process in concentrated salt solution, as described in the fibroin extraction method. The resulting silk solution is not stable as they would become a gel, in another word, change back to rigid silk II conformation. Therefore, lyophilization is demanded. We noticed that even keeping lyophilized fibroin at -20 °C; its water solubility was eventually altered leading to some fibroin precipitation.
Thus, to confirm the content of fibroin, the calibration curve was constructed utilizing freshly prepared lyophilized fibroin that could completely be dissolved in water. Its absorbance was determined using UV-Vis spectrophotometer at 276 nm. Then, the calibration curve was constructed from fibroin concentration ranging from 0.1 to 1.0 mg/ml, fig. 1.

Nanoparticles formulation
One of the most popular bottom-up methods for preparing the FNP formulation is desolvation. Upon slowly injected silk I aqueous solution into the water-miscible organic solvent, the FNP was formed due to dehydration leading to the formation of insoluble silk II. However, using this method, the particle properties are heavily depended on the preparation conditions. For example, the particle sizes from different studies were significantly diverse, ranging from 35 nm to 1500 nm [7]. Additionally, in a study combining the desolvation and freezing methods, the fibroin particle size and size distribution could be altered by the amount of ethanol added and the freezing temperature [7]. Thus, the important formulation factors include the fibroin concentration (X1), the volume ratio between the fibroin solution and ethanol (X2), and the formulation time (X3). It is worth to note that the type of polar organic solvents such as methanol, ethanol, propanol, acetone, and isopropanol has no significant effect on the particle size and morphology [13]. Amongst them, methanol and ethanol were utilized earliest and most frequently. To reduce the potential toxicity as well as the environmental harm, ethanol was chosen as a desolvation solvent. Particle size, zeta potential, and PI were chosen as the responses because of their importance in affecting FNP other properties such as dissolution profile and stability. For instance, an increase in PI (>0.5) could enhance the agglomeration tendency in a nanosuspension due to a broad particle size distribution. Table 1 demonstrates 19 formulations, designed by the software, and their respective variables.

Particle size (Y1)
The hydrodynamic diameter of the FNP has been determined using DLS method. All formulations showed the mean particle size in a range of 204.1±7.9 to 595.6±5.1 nm (column Y1, table 1). The Doptimal design was conducted to study the effect of two continuous factors, X1 and X2, and one qualitative factor, X3 on the particle size response (Y1 It was important to note that each coefficient number (i.e., bi) has its own p-value. Thus, only the significant numbers with p<0.05 were selected and put into the equation. The full list of these values and their statistical significance were shown in table 2. The positive and negative values of the coefficients indicate the proportional and inversely proportional relationship, respectively, between the responses and different level combination of independent variables.
In addition, these relations were reliable due to the reasonable concurrence of less than 0.2 between the predicted R 2 (0.8679) and the adjusted R 2 (0.9325). According to the software, a difference of more than 0.2 between these two values indicates a large block effect, which possibly a problem with the model. The lack of fit test revealed a p-value of 0.1375 indicating the model was fit.
From the equation, it was clear that the particle size was significantly affected by the fibroin concentration (p<0.0001), the fibroin solution and ethanol ratio (p<0.0001), and the formulation time (p = 0.0003). For instance, an increase in fibroin concentration from 0.25 to 2.00% resulted in an increase in particle size, fig. 2. This observation was in agreement with the theory, the higher content of fibroin could increase the chance of particle collision due to Brownian motion and eventually leading to particle aggregation. Thus, bigger particles are formed. The same theory could be used to explain the effect of variable X2 on the particle size. When the ratio increased, more fibroin molecules were also subjected to ethanol, leading to bigger particle size.

Homogeneity of the system (Y2)
In the case of particle homogeneity, the PI value was employed. It is a dimensionless number obtained from the fraction of the square of standard deviation and the square of average particle size (SD 2 /Size 2 ). Based on the equation, a high PI indicates a large deviation of particle size in the system.
Thus, the PI value of <0.1 indicates a monodisperse system, while that of<0.5 demonstrates a narrow size distribution. A system with a PI value of >0.5 indicates a broad particle size distribution and is generally unstable due to the Ostwald ripening effect in which the small particles tend to precipitate on the bigger ones. In our data, all designated formulas demonstrated the PI of 0.09±0.01 to 0.36±0.05 (column Y2, table 1), indicating the narrow size distribution leading to system homogeneity. The polynomial equation derived by best fit mathematical quadratic model, with a p-value of <0.0001 was: The coefficient b5 and b7 were cut out from the equation due to their insignificant p-values of 0.2920 and 0.5347, respectively (table 2). The lack of fit test showed a p-value of 0.1630, which was insignificant; therefore, the model was fit. Similar to that of the particles size, the difference between the predicted R 2 (0.7773) and adjusted R 2 (0.9053) was less than 0.2, indicating a reliable relationship between the variables and responses. Visual 3D response surface plots and the representative cube, fig. 3, illustrated the clear effects of all 3 independent factors on the polydispersity index which show the same trend as particle size. The broad particle size distribution was observed with increasing the fibroin content. The reason was based on the mechanism of particle formation in the desolvation process, in which nanoparticles are formed based on the precipitation and aggregation of small fibroin molecules. Thus, when increasing the fibroin content, the system could be composed of a broad range of particle size, consequently, lead to an increase in polydispersity index.

Surface charge (Y3)
The FNP zeta potential was determined utilizing the PALS method. The surface charge of a colloid is measured indirectly based on its electricity on the slipping plane or diffuse layer. Using the DLVO theory, one can predict the aggregation probability of a nanosystem depending on its zeta potential. Generally, there are two forces that contribute to the system stability, the van der Waals attraction and the electrostatic repulsion force. The electrostatic force is related proportionally to the particle zeta potential. Thus, a high zeta potential absolute value indicates a high repulsion force between particles, consequently leads to a more stable system. From the results, the zeta potential of all formulas was ranged from-10.15±0.28 to -21.17±1.74 mV (column Y3, table 1), which was considered appropriate. The effects of variables on the zeta potential response (Y3) were clarified. The polynomial equation derived by the best fit mathematical quadratic model was: With a reasonable p-value of 0.0084. The coefficient b4, b6, and b8 were absent because of their high p-values of 0.3928, 0.6251, and 0.7455, respectively (table 2). The difference between the predicted R 2 (0.6030) and adjusted R 2 (0.7795) was less than 0.2, demonstrating an appropriate prediction power of the model. Additionally, the lack of fit test revealed the insignificant p-value of 0.2830, meaning that the model was fit. The FNP zeta potential was affected by all 3 independent variables as shown by the visual 3D response surface plots and representative cube, fig. 4. From the plots, the formulation time showed remarkable effects on the particle charge. The zeta potential was significantly affected by the variable X1 and X2 at 24 h formulation time, but not at the 4 h. This might be attributed to the duration of the stabilization process at 4 °C. The results suggested that 24 h was an optimal time for completing the particle formation. Four hours were inadequate for this process, thus, the particles appeared to have a similar charge, regardless of other factors such as X1 and X2. This phenomenon also explained the reduced adjusted R 2 in this equation as compared to previous equations of Y1 and Y2.

Formulation optimization
In this study, we aimed to use DOE to demonstrate the relationships between three independent variables; the fibroin concentration, the volume ratio of fibroin solution to ethanol, and the formulation time, and the responses; particles size, homogeneity (PI), and zeta potential. The optimized formulation, in our case, was chosen to be the smallest particle size, the lowest PI, and the highest value of zeta potential. It is worth to note that the higher the zeta potential value, the less chance of particle aggregation due to the repulsion force. As the particles had negatively charge, the highest absolute value means the smallest real number. Thus, amongst the selection criteria of the desired goals (constraints) such as "none", "minimum", "maximum", "target", and "within range", all responses were set at "minimum". However, it is impossible to achieve all desirable goals, as three responses might not go in the same trend. Hence, the desirability number (i.e., the relationship between independent variables and desired response values) was not set at a maximum limit of 1, but at 0.9. The results demonstrated that the optimized formulation with the size of 224.8 nm, PI of 0.13, and zeta potential of -19.31 mV could be theoretically obtained with the fibroin concentration of 0.50%, the volume ratio of fibroin solution and ethanol of 0.33, and the formulation time of 24 h, table 3. The desirability for this formula was 0.8770, which was considered appropriate. 3D cube plot illustrated the relationship between the desirability and the chosen variables, fig. 5.

Model validation
Since the main aim of the DOE and the optimization process is to correspond the effects of independent variables to the predictable values of the response variables. Thus, besides the considerate selection of those variables, method validation is also an important task. To perform this, the optimized formulation was prepared as mentioned in the previous section, following the software parameters.  The results are demonstrated in table 3. Clearly, all actual formulations (n = 3) possessed average response parameters of size (238.1 nm), PI (0.12), and zeta potential (-21.78 mV) within the predicted range, with a confidence level of 95%. Thus, the method was reliable.
The DOE research has been proposed previously in various fields of applications, ranging from engineering to biomedical, as well as in microparticles and nanoparticles development [15][16][17][18][19][20]. Although the D-optimal design is considered to be one of the best design for prediction and optimization of the product properties from a mixture of qualitative and quantitative variables, it still possesses some limitations [21]. The first one is that it requires a model from the software. As "optimal" designs are computer algorithm, and thus, model-dependent, which are different than the "classical" designs (factorial designs). A change in the model such as quadratic or linear could significantly affect the final engineering conclusions [21]. Fortunately, most software like Design Expert ® or Matlab ® could predict and fit each individual model to the design. Thus, users can choose the best fit one. Secondly, although better than the factorial design in term of working with continuous variables, "optimal" designs do not generate the best design points from some continuous region [20]. Therefore, validation is a must to confirm the suitability of the design. Finally, like most mathematical methods, there is no 100% guarantee that the result from the optimal design is the true optimal option. Nevertheless, the results are acceptable and are far superior to any ad hoc designs [21].

CONCLUSION
In this study, we were successful in the development, optimization, and validation of a novel D-optimal design derived model for the formulation of fibroin nanoparticles. Our model suggested that the physicochemical properties of these particles could be controlled and predicted by variation of process parameters. The results of our model also significantly confirmed that the fibroin concentration, the volume ratio between fibroin solution and an organic solvent, and the formulation time were crucial for the fibroin nanoparticles characteristics of size, homogeneity, and zeta potential. The model was also validated and demonstrated to be solid and reliable with well-correlated results between the predicted and actual optimal particles parameters. Consequently, utilizing the model to predict the optimal properties of fibroin-based nanoparticle could reduce cost and time-consuming experiments. Nevertheless, the limitation of our model is that it could only be applied to the fibroin nanoparticles formulation. In conclusion, the mathematical design of experiments, specifically the D-optimal design, can be usefully applied in the prediction and optimization of main nanoparticle characteristics.

ACKNOWLEDGMENT
This work was supported by the Naresuan University Grant number R2561B009. Duy Toan Pham would like to sincerely thank Naresuan University ASEAN Scholarship for financial support during his doctoral study.