Model Adequacy in Assessing the Predictive Performance of Regression Models in Pharmaceutical Product Optimization: The Bedaquiline Solid Lipid Nanoparticle Example

Abstract

This study aimed to assess the predictive performance of first- and second-order regression models in optimizing bedaquiline (BQ) solid lipid nanoparticle (SLN) formulations. A three-step central composite design and graphical optimization process was employed. A design of experiments method was used to evaluate the impact of BQ, Tween 80 (T80), polyethylene glycol (PEG), and lecithin on the formulations’ response variables, including Z-average (PSD), polydispersibility index (PdI), and Zeta potential (ZP). Secondly, we quantified the relationship between experimental variables using the regression model coefficients. Lastly, we predicted the responses and verified the models’ adequacies to ensure accurate representation and effective optimization. The first-order polynomial showed poor model adequacy and required further refinement due to its lack of explanatory power and significant predictors. Conversely, the second-order models provided superior fitness, sensitivity to variability, complexity, and prediction consistency. The optimized formulation achieved a desirability value of 0.9998, indicating alignment with the desired criteria. Specifically, the levels of BQ (19.4 mg), T80 (25.2 mg), PEG (39.2 mg), and lecithin (200 mg) corresponded to PdI (0.41), PSD (250.99 nm), and ZP (−25.95 mV). Maintaining a BQ concentration between 10 and 20% and T80 between 15 and 18% is vital for maximizing ZP and minimizing PdI and PSD, ensuring stable SLN formulations. This study underscores the significance of precise model selection and statistical analysis in pharmaceutical formulation optimization for enhanced drug delivery systems.

1. Introduction

Developing efficient drug delivery systems in the pharmaceutical industry is crucial for enhancing therapeutic outcomes and patient compliance. Solid lipid nanoparticles (SLNs) have emerged as promising carriers for drug delivery due to their ability to improve drug solubility, stability, and bioavailability. Prior research demonstrated the potential of nanotechnology in the treatment of pulmonary inflammatory diseases, highlighting the benefit of focused therapeutics with lower toxicity and dosage requirements. Some advantages of nano-drug delivery systems include improved drug product stability, solubility, and controlled release. SLNs have been used to deliver poorly water-soluble molecules and have the potential to achieve sustained drug release or targeted delivery to the site of interest [1,2]. One similar drug of interest with poor aqueous solubility is bedaquiline (BQ), which is used in combination with other drugs in treating multi-drug-resistant tuberculosis (MDRTB). In earlier research, our group detailed the formulation of BQ as SLNs for oral delivery, presenting an opportunity to enhance the drug product performance [3]. In the current study, we have progressed to comparing the adequacy of the first- and second-order regression models in accessing predictive performance regarding product optimization for the BQ SLN formulation. The various regimens for treating tuberculosis (TB) require using multiple drug combinations for at least six months [4]; therefore, the prior research aimed to enhance patient adherence by formulating BQ as SLNs over an existing parenteral approach [5]. However, the current study focuses on comparing first-order and second-order regression models for optimizing BQ SLN formulation. The regression equations were methodically employed to investigate how levels of the drug and excipients impact critical quality attributes (CQAs) such as Z-average (PSD), polydispersibility index (PdI), and Zeta potential (ZP). Utilizing a full factorial design, the group systematically investigated the impact of these key formulation variables on the CQAs of the BQ SLN formulation.

After identifying a potential nanocarrier system, a critical phase in product development involves designing and optimizing the formulation [6,7]. Various statistical experimental design methods have been employed in pharmaceutical manufacturing to optimize process parameters. Traditional optimization techniques typically involve altering one factor at a time, necessitating numerous experiments, and failing to capture the combined effects of multiple variables. These methods also demand extensive data and time to identify optimal levels, making them unreliable [8]. In contrast, the design of the experiment (DoE) approach aims to understand interactions among variables, facilitating the optimization process and providing relevant statistical models [9,10]. In line with this, multiple techniques and approaches can be utilized to achieve this optimal nanoparticle formulation; one such method is the central composite design (CCD) method [6,10]. Therefore, this study examined the combined effects of input factors such as BQ, polysorbate or tween 80 (T80), polyethylene glycol (PEG), and lecithin on sonicated and unsonicated SLN formulations. Furthermore, the research optimized the formulation process parameters using CCD and the response surface methodology (RSM). This approach seeks to provide contour and response surface plots to further evaluate the effects of input factors on response variables, allowing for a more comprehensive understanding of the formulation optimization process. Therefore, in product optimization, a typical first step is to use a first-order statistical analysis to screen the variables that will be further investigated under a multifactorial experimental design.

The first-order model offers an initial exploration of the response surface and identifies regions of interest for further investigation. The results obtained can guide the design of experiments for a second-order model or other multifactorial experimental design. For instance, if certain factors are identified as significant in the first-order model, they can be prioritized in the design of the second-order experiments [11]. These steps were used to ensure a more efficient and targeted approach to experimentation, used in the current study.

1.1. First-Order Statistical Analysis

Generating a first-order model using statistical software like JMP involves several critical steps. Initially, data collection is paramount, requiring the identification of input variables (factors) and the response variable, collected systematically for accuracy. Various metrics such as R-squared (R²), adjusted R², p-values, and residual plots are provided to evaluate the model’s fit and explain the variability in the response variable [12].

1.2. Multifactored Experimental Design

The parameters were analyzed using an RSM design called CCD. RSM is a systematic approach that uses statistical and mathematical tools to model and analyze the relationships between several variables and the responses they produce in a process [11]. It is useful for fitting a quadratic surface to the data, which helps to identify the optimal settings for process parameters. By minimizing the number of required experiments, RSM efficiently explores the relationships and interactions between multiple variables, enabling a thorough analysis and optimization of the process [13,14]. It employs quantitative data from the experiments to construct regression models and optimize an outcome (dependent variable) affected by the independent factors (input variables).

Equations (1) and (2) demonstrate linear and quadratic functions, respectively.

Y = β₀ + β₁ x₁ + β₂ x₂ + … + β_q x_q + ϵ (i = 1, 2, …, N)

(1)

$$\mathrm{Y}={\beta }_0+{\sum }_{j=1 }^k{\beta }_j{x}_i+{\sum }_{j=1}^k{\beta }_{jj}{x}_j^2+\dots +\mathsf{\Sigma }{\sum }_{i<j=2}^k{\beta }_{ij}{x}_i{x}_j+ϵ\dots$$

(2)

In these equations, Y represents the response, β₀ is a constant, and β_i, β_ii, and β_ij are coefficients representing linear effects and interaction terms. The terms x_i and x_j are the independent variables, and ϵ denotes the error.

1.3. Central Composite Design (CCD)

CCD is an experimental approach commonly used in RSM [12]. It offers an efficient way to explore the influence of multiple variables (n) on a response. CCD combines several elements:

• Factorial runs (2n): these experiments investigate each variable at two pre-defined levels.
• Axial points (2n): these additional runs explore values beyond the initial two levels for each variable, often denoted as ±α from the center point.
• Center points (nc): replicates are performed at a central point, where all variables are set at their mid-point values.

This design allows researchers to assess both linear and quadratic effects of the variables on the response. With a growing number of variables (n), the total number of runs in a complete CCD increases slower than in a full factorial design, making it more efficient. CCD is useful for modeling quadratic relationships, where the individual second-order effects cannot be isolated using only factorial runs. This makes it an invaluable tool for process optimization, as it can identify the combination of variables that achieve the desired outcome, particularly in exploring the potential of BQ SLNs for oral delivery. The combined effects of these factors are optimized to achieve the desired outcomes, either by minimizing or maximizing specific conditions [10,12,15]. This study aimed to assess the predictive performance of first- and second-order regression models in optimizing the BQ SLNs formulations.

2. Materials and Methods

2.1. Experimental Design

The method for the preparation of the BQ SLNs was comprehensively described in the first paper [3]. Summarily, 10–20 mg of BQ and a target mass of 212 mg lecithin were used to form the organic phase, while the aqueous phase was prepared using a 1:2 ratio of T80:PEG. The optimization process for the BQ SLN formulation used CCD involves three fundamental steps. The first was to design statistically based experiments (DOE) using a selection and control of four independent variables (BQ, T80, PEG, and lecithin) and systematically studying their effects on the CQAs of the BQ SLN formulations. Secondly, the coefficients within a mathematical model were estimated to quantify the relationship between the independent variables and the formulation’s CQAs (PSD, PdI, and ZP). Lastly, we predicted the response and verified the model’s adequacy, to ensure that the model accurately represents the experimental data and can be used to optimize the process effectively. Four independent variables were chosen for the DOE: BQ (X₁), T80 (X₂), PEG (X₃), and lecithin (X₄); see Equation (3). The range and levels of these factors were varied according to the experimental design. These four variables, along with their respective ranges, were identified as critical parameters for improving the drug-release profile of the solid lipid nanoparticle formulation [3].

Y = f (X₁, X₂, X₃ … X_n)

(3)

where Y is the response of the system and X_n is the independent variables (input factors). It is assumed in the design that the independent variables in the experiments are continuous and regulated by experiments with negligible errors. The CQAs were identified as PdI (Y₁), PSD (Y₂), and ZP (Y₃).

2.2. Generating Experimental Data from DOE Trial Runs

The various ranges and levels of the input variables were used to generate twenty-four trial samples of BQ SLNs. These samples were then evaluated to obtain the response values (observed data) needed to generate predicted data to optimize the formulation. A full factorial design was created to consider all factor-level combinations before fitting the model. The collected data were formatted and entered into JMP Pro version 16 statistical software [16] to fit the model. This software program performs statistical analysis, modeling, and data visualization. It is used to support DOE and design generation, where it interactively builds and refines graphs and tables with the graph builder, and tabulates tools, respectively. The software also employs the least squares method for fitting linear regressions such as the first-order model.

2.3. Model Development and Analyses

This section was divided into two parts. The first involved model fitting using a first-order polynomial, and then statistical analysis and optimization was carried out using a second-order model.

2.3.1. Model Fitting Using a First-Order Polynomial

The first-order model was fitted using the JMP software, followed by regression analysis to estimate coefficients that minimize the difference between observed and predicted response values. This first-order regression model was applied to the data using ordinary least squares (OLS) regression methods. The resulting coefficients (β₀, β₁, β₂, β₃, etc.) were analyzed to understudy the linear relationships between the input variables and the responses. The first-order regression model is expressed as

Y = β₀ + β₁X₁ + β₂X₂ + … + β_pX_p + ϵ

(4)

Y_PdI, Y_PSD, and Y_ZP are the dependent variables for unsonicated and sonicated BQ SLN formulation, while X₁, X₂, X₃, and X₄ are the independent variables. β₀, β₁, β₂, … β_p are the coefficients, and ϵ is the error term (residual).

2.3.2. Statistical Analysis and Optimization Using a Second-Order Model

The JMP software also performed parameter combinations for the CCD. The design evaluated the effect of BQ, T80, PEG, and lecithin as independent variables, on PSD, PdI, and ZP as response variables. Polynomial coefficients were computed from experimental data to predict these response variables. The model adequacy was tested using the analysis of variance (ANOVA) test. By applying a multiple regression analysis, the results were adjusted to a second-order polynomial equation. The equation for responses in both the unsonicated BQ formulation and the sonicated versions were obtained from the group’s prior study, shown in

Table 1. The ranges and levels for the independent variables used in the DOE to target achieving CQAs for the BQ SLN formulation.

Input Variables (mg)			Ranges and Levels
			−β	−1	0	+1	+β
X1	Bedaquiline (BQ)		11.1	13.5	15.8	18.2	20.5
X2	Polysorbate 80 (T80)		11.3	14.6	17.9	21.1	24.4
X3	Polyethylene Glycol (PEG)		24.1	30.8	37.5	44.2	50.9
X4	Lecithin (Lec)		202.6	215.5	228.4	241.2	254.1
Response Variables		Goals	Acceptable Ranges
Y1	Polydispersity Index (PdI)	Minimize	≤0.4
Y2	Z Average (PSD)	Minimize	≤500 nm
Y3	Zeta Potential (ZP)	Match target	−30 mV > X > + 30 mV

.

The experimental runs were carried out in a randomized sequence to minimize errors and the impact of uncontrolled factors. The research optimized the response variables (Y) by closely approximating the relationship between the independent variables and the response surfaces. The response variables were utilized to develop an empirical model that correlates experimental variables through a second-degree polynomial equation, as represented in Equation (5) [11,17].

$$\mathrm{Y}={\beta }_0+{\sum }_{j=1 }^k{\beta }_j{x}_i+{\sum }_{j=1}^k{\beta }_{jj}{x}_j^2+\dots +\mathsf{\Sigma }{\sum }_{i<j=2}^k{\beta }_{ij}{x}_i{x}_j+ϵ\dots$$

(5)

In this equation, Y represents the response predicted by the model. β₀ is the constant coefficient, β_j is the linear coefficient, β_jj is the quadratic coefficient, and β_ij is the interaction coefficient. The variable parameters for the BQ SLN formulation process are denoted by x_i and x_j, with n representing the number of factors being studied and optimized. The term ϵ signifies the random error component in the model.

Table 2. The relationship between the coded and actual values of the variables.

Code	Actual Values of Variables
−β	Xmin
−1	[(Xmax + Xmin)/2] − [(Xmax − Xmin)/2α]
0	[(Xmax + Xmin)/2]
+1	[(Xmax + Xmin)/2] + [(Xmax − Xmin)/2α]
+β	Xmax

Note: α is 2^k/4; K is the number of independent variables, which is 4. The contents of the table are adapted from Ref. [18].

was used to calculate the function of the range of interest of each input factor shown in Table 1. This study utilized a face-centered (FC) CCD, where the number of tests (N) includes the standard 2^k factorial points with the origin at the center. It placed axial points at the face centers of the factorial cube to simplify the design and avoid issues with extreme points. In this design, the quadratic terms were generated from the center at a distance of α from 2k axial points with no replicates or center points. Here, the number of independent variables is K = 4; therefore, α = 2.00, obtained from 2^4/4. The axial points (2k) are used for screening analysis and maintaining readability, ensuring the constant variance of model predictions at equidistant points from the design center [18].

The total number of experimental trials (N) used in the DOE is given by Equation (5), where C is the number of center points. The study design has four independent variables (k = 4), and the number of center points is assumed to be 0; therefore, the total number of trial experiments conducted was N = 16 + 8 + 0 = 24.

N = 2^k + 2k + C

(6)

Thus, for these experiments with four variables, the design includes sixteen factorial points, eight axial points, and zero replicated center points, totaling 24 runs. The variables are coded to ±1 for factorial points, 0 for center points, and ±β for axial points, facilitating a straightforward and efficient experimental setup.

2.3.3. Model Fitting Parameters

In this study, various fit statistics were used to assess the adequacy and performance of the regression models. The primary parameters included the coefficient of determination (R²), adjusted R², and predicted R², which indicate the proportion of variance in the response variable explained by the model. A higher R² value signifies better model fit. Adjusted R² accounts for the number of predictors in the model, offering a more accurate measure of fit when multiple variables are involved. Predicted R^2k evaluates the model’s ability to predict new data, with higher values indicating strong predictive performance. Additional metrics such as the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) were used to balance model fit with complexity, where lower values suggest a more parsimonious model. Finally, the Root Mean Square Error (RMSE) and Coefficient of Variation (CV) were calculated to assess the model’s precision and variability, with lower values reflecting higher accuracy and consistency [11].

2.4. Stepwise Evaluation of the Regression Model’s Adequacy

To evaluate the adequacy of regression models applied to the BQ formulations, the model adequacy check compared the performance of first-order regression models between the unsonicated and sonicated BQ formulations. The analysis involved fitting first-order polynomial equations to the data to reveal which formulations demonstrate acceptable model adequacy. If initial analyses suggest that the models need further refinement to enhance their predictive accuracy and overall reliability, second-order models will be utilized.

2.5. Contour Plots

The JMP statistical software was employed to perform second-order polynomial regression analyses on experimental data generated from the BQ SLNs sonicated formulations. The process involved generating critical diagnostic plots, such as the actual vs. predicted plots, contour plots, and response surface plots, which are essential for evaluating model performance and understanding the relationships between variables.

3. Results

3.1. Model Fitting Using a First-Order Polynomial

The following first-order regression models were obtained for the sonicated and unsonicated BQ SLN formulations. The model equations listed below were derived by executing parameter estimations in the JMP software. The linear equations describe the overall relationship between predictors and expression levels (response variables) within the chosen statistical model. The resulting coefficients illustrate the linear relationships between the input variables and the responses.

For unsonicated bedaqualine formulations:

Y_PdI = 0.2854 + 0.0131(X₁) − 0.0004(X₂) + 0.0011(X₃)

(7)

Y_PSD = 339.979 + 8.354(X₁) − 0.4992(X₂) − 1.2104(X₃) − 0.0107(X₄)

(8)

For sonicated bedaqualine formulations:

Y_PdI = 0.0290 + 0.0213(X₁) − 0.0004(X₂) + 0.0004(X₃)

(9)

Y_PSD = 70.8680 + 7.1795(X₁) + 0.7935(X₂) − 0.4162(X₃) − 0.0113(X₄)

(10)

Y_ZP = −42.0184 + 0.1594(X₁) + 0.0907(X₂) − 0.1958(X₃) − 0.0009(X₄)

(11)

3.1.1. ANOVA Test Results for the BQ SLN Formulations

Full details of the ANOVA results for the sonicated and unsonicated BQ SLNs formulations were communicated in an earlier paper [3]. Summarily, the ANOVA results from the DOE data indicated that the model lacked significance for all dependent variables in the unsonicated formulations; the p-values for PdI and PSD were estimated as 0.5877 and 0.0814, respectively. However, in the sonicated formulations, significant effects were more pronounced, particularly with PdI, compared to other response variables such as PSD and ZP. The p-value of PdI was found to be statistically significant (p = 0.0003) for a level of significance of α = 0.05.

3.1.2. Summary of the First-Order Regression Models

The optimization of pharmaceutical formulations, such as BQ, requires precise and reliable models to predict key response variables. The adequacy of first-order models for predicting response variables in the sonicated and unsonicated BQ formulations is summarized in

Table 3. Model summary for unsonicated and sonicated bedaqualine formulations using first-order linear regression models.

Model	R2	Adj R2	Predicted R2	Precision Ratio	Std Dev.	CV (%)	AIC	BIC
Unsonicated Formulation
PdI	0.1319	−0.0584	−0.2568	0.8809	0.16	31.77	−5.97	−3.85
PSD	0.3403	0.2014	−8.03	0.5158	67.68	14.20	276.4	278.5
Sonicated Formulation
PdI	0.6537	0.5808	−12.2718	0.9748	0.01691	32.60	−45.68	−43.55
PSD	0.3403	0.2014	−15.8241	0.00008	11.4569	28.04	267.16	269.29
ZP	0.07806	−0.1160	−11.3790	0.0025	2.0795	28.83	185.25	187.38

.

The study compared the optimization of BQ formulations using first-order models for sonicated and unsonicated variants. The unsonicated formulation’s PdI model showed poor fit and prediction capabilities (R² = 0.1319, adjusted R² = −0.0584, predicted R² = −0.2568), with a Precision Ratio of 0.8809, a standard deviation of 0.0331, and a Coefficient of Variation (CV) of 20.41%. The Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) were −5.97 and −3.85, respectively. The combined metrics of low R², high variability, negative AIC and BIC values, negative adjusted and predicted R², and low precision collectively highlight the limitations in the first-order polynomial to accurately predict the PdI model in the unsonicated formulations. The PSD model for unsonicated formulations had a moderate fit (R² = 0.3403, adjusted R² = 0.2014, predicted R² = −15.8241), but poor precision (Precision Ratio = 0.00008, CV = 28.04%), with AIC and BIC values of 267.16 and 269.29. The ZP for unsonicated formulation was not included because the ZP values in sonicated formulations provide more meaningful insights into the formulation’s stability and performance. Sonication, unlike the unsonicated counterpart, can alter particle size and surface charge, making ZP a critical factor in evaluating and optimizing the formulation [19].

For the sonicated formulation, the PdI model exhibited better fit (R² = 0.6537, adjusted R² = 0.5808), but poor prediction (predicted R² = −12.2718), with a Precision Ratio of 0.9748, a standard deviation of 0.01691, and a CV of 32.60%. It also had lower AIC and BIC values, −45.68 and −43.55. The PSD model for sonicated formulations had similar fit and prediction issues to the unsonicated ones, but the ZP model showed poor fit and prediction (R² = 0.07806, adjusted R² = −0.1160, predicted R² = −11.3790), with AIC and BIC values of 185.25 and 187.38. The later suggests that the first-order polynomial is not appropriate for accurately predicting the ZP model, as evidenced by the combination of a very low R² value, negatively adjusted R², significantly negative forecasted R² and high AIC and BIC values. A lower AIC or BIC value indicates a better fit [20]. Overall, the sonicated formulation’s PdI model demonstrated better fit and precision, making it more adequate despite prediction challenges, highlighting areas for improvement in both formulations. Considering the poor indices obtained using the first-order modeling, there was a need to employ a second-order model approach.

3.2. Model Fitting Using a Second-Order Polynomial (CCD)

The software was also used to fit the second-order model equations that evaluated the effect of BQ, lecithin, T80, and PEG as independent variables, on PSD, PdI, and ZP as response variables. Second-order polynomial coefficients were computed from experimental data to predict these response variables. The model adequacy was tested using the ANOVA test. Multiple regression techniques were used to analyze the data and fit a model. Which, unlike the first order, includes linear terms, quadratic terms, and interaction terms. By this, second-order polynomial equations for the sonicated BQ formulations that captured more complex relationships and curvature in the data were generated as shown below:

Y_PdI = 0.411 + 0.00 + 0.028X₂ + 0.001X₃ + 0.001(X₂)²

(12)

Y_PSD = −34.309 − 44.184X₁ + 38.489X₂ − 5.840X₃ + 2.637X₄ + 0.927(X₁)² + 0.562X₁X₂ − 0.662(X₂)² + 0.226X₁X₃ + 0.234X₂X₃ + 0.063(X₃)² − 0.012X₁X₄ − 0.101X₂X₄ − 0.030X₃X₄

(13)

Y_ZP = 53.184 − 5.253X₁ − 1.943X₂ + 2.593X₃ + 49.309X₄ − 0.051(X₁)² + 0.135X₁X₂ − 0.180(X₂)² + 0.027X₁X₃ + 0.024X₂X₃     0.031(X₃)² + 0.01X₁X₄ − 0.003X₂X₄ − 0.004X₃X₄

(14)

Model-Adequacy Checks for the First- and Second-Order Regression Models Used to Evaluate Sonicated BQ SLN Formulations

The earlier step, using first-order models in optimizing both formulations, shows that in sonicated formulations, statistical significance, higher precision, and improved predictive model fits the polydispersity index (PdI) and a better model fitness and complexity was achieved in the formulation’s particle size distribution (PSD). This implies that sonication drastically changes the physical characteristics of the nanoparticles, such as particle size and surface charge, which are essential for their stability and functionality, hence the need for greater emphasis on sonicated formulations in subsequent investigations. To evaluate the model adequacy of the first-order and second-order models in predicting the behavior of sonicated BQ formulations, the main effects from both models for sonicated BQ SLN formulation are compared in

Table 4. ANOVA table comparing main effects from first- and second-order model sonicated BQ SLN formulation.

	Sonicated Formulation (First-Order Model)			Sonicated Formulation (Second-Order Model)
Variables	PdI	Z-Average (PSD)	ZP	PdI	Z-Average (PSD)	ZP
Main Effects	p-Value	p-Value	p-Value	p-Value	p-Value	p-Value
Model	0.0003 *	0.0814	0.8046	0.0073 *	0.0918	0.6389
BQ	0.0001 *	0.0091 *	0.7264	0.00006 *	0.0064 *	0.5856
T80	0.6259	0.6803	0.7950	0.5000	0.2728	0.4226
PEG	0.7577	0.6689	0.2744	0.9695	0.5925	0.4344
Lecithin	0.6971	0.6434	0.8293	0.5088	0.3806	0.3697

* Significant p-values.

. To further compare fitness, variability in response factors, complexity, and prediction consistency in the two models, other model evaluation metrics were computed in

Table 5. A comparative analysis of sonicated bedaqualine formulation parameters in the linear and quadratic model.

Model	R2	Adj R2	Predicted R2	Precision Ratio	Std Dev.	CV(%)	AIC	BIC
Sonicated BQ Formulation for First-order Model
PdI	0.6537	0.5808	−12.2718	0.9748	0.01691	32.60	−45.68	−43.55
PSD	0.3403	0.2014	−15.8241	0.00008	11.4569	28.04	267.16	269.29
ZP	0.07806	−0.1160	−11.3790	0.0025	2.0795	28.83	185.25	187.38
Sonicated BQ Formulation for Second-order Model
PdI	0.8950	0.7317	0.0111	8.5101	0.0591	16.36	18.44	−40.42
PSD	0.7907	0.4652	−0.5806	3.7789	40.8672	22.94	332.38	273.51
ZP	0.5623	−0.1184	−2.0822	1.2850	9.0742	28.87	260.14	201.28

.

The ANOVA in Table 4 for the sonicated BQ formulations provides crucial insights into the statistical significance of the main effects in both the first-order and second-order models. The primary (input) variables under analysis were BQ concentration, T80, PEG, and lecithin, while the response variables examined included PdI, PSD, and ZP. The p-values for each variable are highlighted, indicating the statistical significance of their effects. In the first-order model, the analysis reveals several key findings. For the PdI, the model is highly significant (p = 0.0003), indicating that the main effects significantly explain the variability in PdI. BQ concentration is highly significant (p = 0.0001), suggesting that BQ is a major factor influencing PdI. In contrast, T80 (p = 0.6259), PEG (p = 0.7577), and lecithin (p = 0.6971) are not significant, indicating that these variables do not significantly affect PdI in the first-order model.

For the PSD, the model is marginally significant (p = 0.0814), implying the main effects explain some of the variability in particle size. BQ concentration is significant (p = 0.0091), indicating it has a notable effect on particle size. However, T80 (p = 0.6803), PEG (p = 0.6689), and lecithin (p = 0.6434) are not significant, indicating that these variables do not significantly impact particle size. For ZP, the model is not significant (p = 0.8046), indicating that the main effects do not significantly explain the variability in the ZP. BQ concentration (p = 0.7264), T80 (p = 0.7950), PEG (p = 0.2744), and lecithin (p = 0.8293) are not significant, suggesting these variables do not significantly affect ZP.

In the second-order model, the analysis reveals more nuanced findings. For the PdI, the model is significant (p = 0.0073), indicating that the second-order model significantly explains the variability in PdI. BQ concentration is highly significant (p = 0.00006), indicating that it is a major factor influencing PdI. T80 (p = 0.5000), PEG (p = 0.9695), and lecithin (p = 0.5088) are not significant, indicating that these variables do not significantly affect PdI. For the PSD, the model is marginally significant (p = 0.0918), implying the second-order model explains some of the variability in PSD. BQ concentration is significant (p = 0.0064), indicating that BQ has a notable effect on particle size. However, T80 (p = 0.2728), PEG (p = 0.5925), and lecithin (p = 0.3806) are not significant, showing that these variables do not significantly impact particle size. The ZP model is not significant (p = 0.6389), indicating that the second-order model does not significantly explain the variability in ZP. BQ concentration (p = 0.5856), T80 (p = 0.4226), PEG (p = 0.4344), and lecithin (p = 0.3697) are not significant, suggesting these variables do not significantly affect ZP.

3.3. Model Summary for First-Order Regression for the Sonicated Samples

The analysis of the sonicated BQ formulation parameters, derived from the linear and quadratic models, reveals significant differences in model performance and prediction accuracy. The linear model for PdI demonstrates moderate explanatory power with an R² of 0.6537 and an adjusted R² of 0.5808. However, the predicted R² is highly negative (−12.2718), indicating poor predictive capability. The Precision Ratio is close to 1, suggesting low precision. The model’s standard deviation is low (0.01691), but the high CV of 32.60% indicates substantial variability relative to the mean. For PSD, the linear model shows low explanatory power with an R² of 0.3403 and an adjusted R² of 0.2014.

3.4. Model Summary for Second-Order Regression for the Sonicated Samples

The quadratic model for PdI exhibits high explanatory power with an R² of 0.8950 and an adjusted R² of 0.7317. The predicted R² is positive (0.0111), indicating reasonable predictive capability, compared to the first-order model approach. Similarly, the Precision Ratio is higher (8.5101), reflecting good precision. Furthermore, the standard deviation (0.0591) and CV (16.36%) are both lower than those in the first-order model, indicating less variability. These results suggest the second-order model gives better predictions for the formulations.

3.5. Optimization for BQ SLNs Formulations

3.5.1. Contour and Response Surface Plots (2D and 3D RSM)

These visual tools are essential for diagnosing model fit, understanding variable interactions, and optimizing conditions to achieve the best possible outcomes. The statistics obtained from the second-order modeling were used to generate 2D and 3D RSM.

3.5.2. Effect of BQ and T80 on PdI

In Figure 1a, the actual versus predicted residual plot shows that the observed relationship appears reasonable, supported by the high significance of the quadratic model (p < 0.0073). The dispersion of data (dots) around the trend line indicates a good model fit, evidenced by an R² value of 0.89. This suggests that changes in the input factors explain 89% of the variability in the PdI. Additionally, only two actual observations fall outside the shaded confidence interval region.

The response surface plot in Figure 1b appears to be elliptical, suggesting a quadratic relationship between the input factors (BQ and T80) and the response variable (PdI). This implies the presence of an optimal region within the design space. As BQ increases, PdI initially increases to a maximum value and then decreases, indicating a quadratic effect. As T80 increases, PdI shows a similar quadratic trend, first increasing to a peak and then decreasing. The gradients show that the system is sensitive to changes in both factors, especially in the central region. This plot is useful for identifying the precise levels of BQ and T80 that maximize PdI, guiding the optimization process in formulation development. The legends in the response surface plots are critical components that provide a color-coded guide to the values of the response variable across the surface. It also visually maps the response values and aids in trend identification, optimal region detection, and data comparison.

The contour plot in Figure 1c shows a range of PdI values, from approximately 0.229 to 0.458. There is a curvature in the contour lines, suggesting a nonlinear relationship between the variables and the response. This indicates that the effect of one variable on PdI depends on the level of the other variable. The central region of the contour plot, where the PdI values are lower (around 0.229 to 0.3206), appears to be an optimal range for both BQ and T80 concentrations. This region represents conditions where the PdI is minimized, indicating more uniform particle sizes. As one moves away from this central region, the PdI values increase, indicating less optimal conditions for the formulation. Increasing the concentration of BQ tends to increase the PdI initially, as indicated by the rising contour lines moving from left to right. Similarly, increasing the concentration of T80 also tends to increase the PdI initially, as indicated by the rising contour lines moving from bottom to top. Furthermore, the contour lines are closely packed, indicating a more sensitive relationship between levels of BQ and T80, and the PdI. This suggests that small changes in the level of these input factors in the sonicated formulation have a substantial impact on PdI.

3.5.3. Effect of BQ and T80 on Z Average (PSD)

The plot in Figure 2a demonstrates that the model provides a moderate fit (R² = 0.79) for predicting the PSD of the sonicated BQ formulation, explaining a substantial portion of the variability. However, the prediction error (RMSE = 40.867) and the p-value (0.0918) suggest that while the model is fairly accurate, there is room for improvement in precision and statistical significance. The 3D and contour plots for the PSD model are shown in Figure 2b,c.

For the response surface plot in Figure 2b, as BQ increases, the PSD decreases sharply initially, indicating a strong influence on particle size reduction at lower concentrations. However, beyond a certain point, the decrease in PSD starts to plateau, showing diminishing returns with higher BQ concentrations. The influence of T80 on the PSD is also significant. At lower T80 concentrations, the PSD decreases, but this trend reverses at higher concentrations, likely due to micelle formation and increased surfactant layer thickness, which could increase particle size. The optimal region for the PSD is at the lower front corner of the plot, where both BQ and T80 concentrations are lower. This area corresponds to the smallest PSD values, indicating the most favorable conditions for achieving minimal particle size.

In the contour plot in Figure 2c, the highest PSD values (greater than 281.089) occur at the top of the plot, particularly where T80 is high and BQ is low. The lowest PSD values (less than 129.7) are at the bottom left, where BQ and T80 are low. The optimal region, where PSD is maximized, is indicated by the outermost contour lines at the top of the plot. This suggests that the highest PSD values occur when BQ is approximately 18–20 and T80 is approximately 20–25.

3.5.4. Effect of BQ and T80 on ZP

The R² value of approximately 0.5765 indicates that the second-order model predicts only about 58% of the variability in the actual ZP of the formulation. This suggests a moderate level of predictive power. The RMSE of approximately 9.1 indicates that on average, the model’s predictions deviate from the actual ZP values by about 9.1 mV. Figure 3 shows the 3D response surface and 2D contour plots that explain the effect of BQ and T80 on the ZP at fixed levels of lecithin and PEG concentration.

The response surface plot in Figure 3b shows a circular pattern around the peak region, indicating a local maximum for ZP. As BQ concentration increases from 10 to around 20 units, ZP values initially increase, reaching a peak. However, beyond this point, further increases in BQ concentration lead to a slight decrease in ZP. Similarly, an increase in T80 concentration from 10 to around 20 units increases ZP; reaching a peak value. The optimal region for ZP is around 20 units for both BQ and T80 concentrations, where ZP reaches its peak value. This region is marked by the highest contour levels (dark brown).

The contour plot in Figure 3c shows that higher values of ZP (less negative) are observed in the regions where T80 is between 20 and 25, and BQ is between 14 and 18. Lower ZP values (more negative) are seen where T80 is lower and BQ is between 10 and 12. The optimal region, where ZP is less negative (around −11.059), is located where T80 is high (above 20) and BQ is around 18. Conversely, the most negative ZP values (around −39.995) are seen where T80 is low (around 10) and BQ is moderate (between 12 and 16).

3.6. Graphical Optimization

Graphical tools such as contour plots and response surface plots are essential in visualizing the relationship between factors like BQ and T80 concentrations and their impact on the desired response. The use of these tools to guide the optimization of BQ (10–20%) and T80 (15–18%) can be explained as follows: first, the contour plots were designed with the use of JMP software to show regions of constant responses based on varying concentrations of BQ (10–20%) and T80 (15–18%). When the concentration of the variables are consistently adjusted, it can identify the optimal regions where the CQAs meet the desirable product target profiles, such as minimal PSD. For example, the range of BQ (12–18%) and T80 (16–17%) concentrations, where a PSD of ≤ 500 nm is achieved, forms the acceptable region considered as the optimal conditions for achieving that CQA. When the concentration of either input variable deviates from the specified limits, the formulation’s PSD increases, moving outside the desired target profile. The interplay in varying the BQ and T80 concentrations to identify the region of optimal conditions for the response variable is used to showcase the systematic optimization of the formulation parameters. On the other hand, the 3D plot provides a three-dimensional view of the response, showing how it varies across the entire range of BQ and T80 concentrations. With this, the optimal region where the desired response is maximized or minimized is identified.

The optimization of the developed sonicated BQ SLN formulation was performed to find the levels of concentration of BQ (X₁), T80 (X₂), PEG (X₃), and lecithin (Lec) (X₄). After plotting the desired values of dependent and independent variables on JMP software, it suggested that for the factors X₁, X₂, X₃, and X₄, the values were 19.43, 25.21, 39.20, and 200.00, and PdI (Y₁), PSD (Y₂), and ZP (Y₃) were 0.41, 250.99 µm, and −25.95 mV, respectively; see Figure 4.

A statistical profiler tool predicts the optimal levels of input factors required to achieve output variables that provide a design space that maintains the CQAs of a target formulation. The JMP software model profiler for the second-order approach generated a high desirability level for the sonicated BQ formulations. The desirability value of 0.9998 indicates that the conditions under which this value is obtained are almost perfect in meeting the desired criteria. This suggests that the BQ SLN formulation or process conditions are nearly optimal, with the levels of the input factors set to produce response variables that almost perfectly align with the desired targets. With such a desirability value, further optimization may be minimally needed, as the current conditions yield optimal performance. The high level of desirability reflects strong confidence in the process or formulation, indicating excellent model fit and predictive capability. It also implies that the trade-offs between different response variables have been effectively balanced, resulting in a robust design space where small variations in input factors will not significantly affect the outcomes.

In designing a preliminary design space for the BQ SLN formulations, the overlay plot shown in Figure 5 depicted a white region, representing the area where optimal conditions for all three response variables overlap. The design space region is bounded by PSD, PdI, and ZP. The respective values are 110.6 µm (indicating smaller and more desirable particle sizes), 0.39 (showing a more uniform particle size distribution), and −19.87 mV (indicating moderate stability), which are generally acceptable for the formulation. Moreover, it also illustrates that the optimal formulation is achieved at a T80 concentration of 15–20% and a BQ concentration of 14–18%. This range exhibits high ZP and optimal PSD, indicating stable and appropriately sized BQ SLN formulations. As T80 increases beyond 18%, ZP decreases, signifying reduced particle stability.

The design space is robust within these concentrations, but exhibits sensitivity to changes, with small deviations potentially leading to significant impacts on PSD and stability.

Therefore, maintaining concentrations within the optimal range defined in the design space is crucial to ensure the desired stability and performance of the BQ SLN formulations. Future steps will validate and standardize the optimized conditions to ensure consistent high-quality outcomes.

4. Discussion

BQ SLNs were generated from sonicated and unsonicated formulations. Based on the results obtained, the final regression model analyses focused on the sonicated formulations. The first-order model (linear regression) captures only linear relationships with fewer parameters, making it easier to interpret. In contrast, the second-order model (quadratic regression) includes linear, quadratic, and interaction terms, offering more flexibility to capture complex relationships and curvature. However, this requires more data and is therefore more complex. This study compared the performance matrices of these two regression equations to evaluate the model performance of the individual response variables (CQAs) for the BQ SLN formulations. This comparison aimed to identify which regression analysis—first or second—yielded models with superior characteristics in terms of fitness, sensitivity to variability, complexity, and prediction consistency. It also helps determine which model better explains the variability in the data and provides a more accurate representation of the underlying relationship between the variables. CCD was used to fit the quadratic models in this study. It is useful in RSM, where the goal is to optimize a response influenced by several variables. Furthermore, CCD helps explore the design space and fit a quadratic model [9,15].

4.1. Model Fitting Using First-Order Polynomials for Sonicated and Unsonicated Samples

Several metrics were used in evaluating the adequacy of the statistical models used in this study; including R², adjusted R², predicted R², Precision Ratios, standard deviation, CV, AIC, and BIC. These metrics assessed the model fit, complexity, and predictive performance. Each metric plays a crucial role in determining regression model adequacy. CV provides a measure of relative variability and helps compare model consistency. AIC and BIC are used for model selection, with AIC balancing fit and complexity and BIC imposing a stronger penalty for additional parameters [11,19]. For the unsonicated formulations, the R² values were relatively low, indicating limited explanatory power. High RMSE and CV values also suggest substantial prediction errors and high variability [19,21].

4.2. Model Summary for First- and Second-Order Regression for the Sonicated Samples

Although the R² value for PdI was lower in the first compared to the second (0.6537 vs. 0.8950), both regression models displayed a fair prediction; but showed varied results for the other two output variables. The highly negative predicted R² (−15.8) from the first-order model for PSD indicates poor predictive performance. The Precision Ratio is extremely low (8 × 10⁻⁵), and the standard deviation is high (11.5), reflecting significant dispersion. The high CV shows considerable variability; similarly, high AIC and BIC values indicate a poor fit [22]. The linear model for ZP performs poorly with a very low R² (0.08) and a negative adjusted R² (−0.12); indicating first-order modeling does not explain the variability well. Likewise, this negative predicted R² also reflects poor predictive ability. The low Precision Ratio (0.003) and high standard deviation (2.1) suggest low precision and high variability. The AIC and BIC values are relatively high, indicating a poor model fit.

In contrast, the second-order models show improved performance. For the PSD response, the quadratic model has R² = 0.7907, and a predicted R² of −0.5806; which, while still negative, represents an improvement over the linear model. The Precision Ratio is higher (3.7789), and the lower CV indicated improved precision. The ZP model shows moderate improvement in R² = 0.5623 over the value from the first order. Although the predicted R² is still negative (−2.0822), the Precision Ratio is higher (1.2850), and the standard deviation is reduced (9.1). Although the quadratic model for PdI demonstrates high explanatory power with R² (0.8950) suggesting a clear and reliable signal [20], the presence of high Precision Ratios indicates that while the model performs well, further validation is required to ensure its predictive reliability [23].

4.3. Effects of Input Variables in First- and Second-Order Model Sonicated BQ SLN Formulations

The analysis of variance for both first-order and second-order models provides essential insights into sonicated BQ formulations. In both models, BQ concentration consistently shows significant effects on the PdI and PSD. This indicates that BQ plays a crucial role in determining the quality and characteristics of the nanoparticles. However, T80, PEG, and lecithin do not show significance, suggesting their roles are less critical in these specific formulations or that their optimal concentrations might have already been achieved. Furthermore, the second-order model provided a better fit for the PdI, highlighting the importance of considering nonlinear effects and interactions in the formulation process. Overall, optimizing BQ concentration will yield significant improvements in PdI and PSD and a higher likelihood of achieving the desired quality attributes for the BQ SLN formulations.

4.4. Significance of the Second-Order Polynomial Model Equations

The regression analysis of sonicated BQ formulations reveals crucial insights into how input variables affect key responses. T80 and BQ concentrations significantly influence the PdI and PSD in the sonicated formulations. Specifically, T80 decreases PdI due to its surfactant properties, and bedaqualine reduces particle size. Significant positive interaction terms between BQ and T80, and between T80 and PEG, indicate synergistic effects in reducing particle size. In the sonicated formulation, increasing concentrations of both BQ and T80 decrease ZP, potentially affecting nanoparticle stability.

4.5. Justifying the Use of the Quadratic Model in RSM Optimization

The study compared linear and quadratic models in analyzing sonicated BQ formulations. The linear model showed weaknesses in predicting the behavior and explaining the variability in key characteristics like PSD. An important aspect of model adequacy is the statistical significance of its predictors. In the first-order models evaluated, only BQ showed significant effects on response variables in either of the formulations, while other predictors did not show statistically significant effects. This is evident from the high p-values associated with other predictors, which are all greater than 0.05. This lack of significance further underscores the inadequacy of the first-order models in explaining the variability of the response factors. The drawback in the linear model’s poor predictive performance was evident by low R² values (0.3404 for PSD) and negative predicted R² values across all response variables. Also, prominent is the model’s inability to adequately account for the variability among these variables (as in adjusted R² values of 0.5808 for PdI), thereby allowing for significant variation, and inaccurate predictions, which raises the possibility of inconsistent product quality. These result in suboptimal formulation, longer and more expensive development times due to repeated cycles of formulation testing and adjustments, a lack of reliable data for regulatory submission, and detrimental effects on customer confidence and business reputation. Overall, the quadratic (second-order) model is superior to the linear (first-order) model for the sonicated BQ formulation. It provides higher explanatory power (adjusted R² values of 0.7317 for PdI), better predictive capability, though low (predicted R² values of 0.0111 for PdI), improved precision (a Precision Ratio of 8.5101 compared to 0.9748 for the linear model), and lower variability across all response variables (PdI, PSD, and ZP). All these results suggest that it is the better model for optimization purposes. The result agrees with findings made in a study by Kura Beyene, where the quadratic model equation exhibited an RMSE of 0.0747 while the linear model equation demonstrated a significantly higher RMSE of 0.1596. The results reported the quadratic model provides superior accuracy and reliability in its predictions compared to the linear model [24]. To overcome these limitations, the study opted for a second-order model, which was believed to ensure better efficacy and stability of the BQ formulations by providing a more accurate and reliable framework for optimization [16]. The model performed significantly better. It captured nonlinear relationships between variables, allowing for a more precise optimization of the formulation. This ultimately leads to better drug delivery through achieving the desired particle size and stability.

4.6. Optimization Using 2D and 3D RSM

The 2D and 3D plots create visual representations to illustrate the relationship between input and response variables. The 2D RSM or contour-plot two-dimensional graphs typically show the effect of two independent variables on a response variable while keeping other variables constant. On the other hand, the 3D or response-surface-plot plots show the interaction between two or more independent variables and their combined effect on the response variable [10].

4.6.1. Effect of BQ and T80 on PdI

The sonication of the BQ SLN formulation process enhances particle size uniformity, thereby improving the overall quality of the formulation. The response surface analysis identifies a central region with lower PdI values, indicating uniform particle sizes crucial for consistent drug delivery. Sonication breaks down aggregates, resulting in a homogeneous mixture. The process shows a nonlinear, quadratic relationship between BQ, T80, and PdI, highlighting the need to optimize these concentrations. The quadratic model, with a high R-squared value of 0.65 and a significant p-value (<0.0073), indicates precise control is essential. The process’s sensitivity to BQ and T80 changes ensures stability and efficacy. The strong model fit and predictive accuracy supports the reproducibility, reliability, and scalability of the sonicated BQ SLN formulations.

4.6.2. Effect of BQ and T80 on Z Average (PSD)

The 2D contour plot and 3D response surface explains the maximum PSD as a function of BQ and T80 concentration at a constant level of lecithin and PEG for BQ SLN formulations, considering that at a lower concentration of T80, it was seen that the PSD initially increases significantly, attains maximum effect at values 16, and then levels off at a lower concentration. As the BQ concentration rises from 10 to 20, the quadratic effect for BQ is positive, indicating a convex curvature, meaning the decrease in PSD with increasing levels of BQ starts to level off at higher concentrations. Furthermore, the increase in PSD with a higher T80 concentration post sonication can be attributed to micelle formation, increased surfactant layer thickness, and potential reduction in sonication efficiency.

4.6.3. Effect of BQ and T80 on ZP

The reduction in ZP at a T80 medial concentration value (e.g., 18) can be attributed to the effective adsorption of surfactant molecules onto the lipid particle surface, enhancing stability through steric and possibly electrostatic mechanisms [25]. However, beyond a certain concentration, the formation of micelles and the saturation of the particle surface with surfactant molecules result in a decrease in ZP due to reduced surface charge, destabilizing interactions, and changes in the electrical double layer around the particles. Higher stability is generally believed to be associated with absolute ZP values >30 mV [26]. The optimal region for ZP is around 20 units for both BQ and T80 concentrations, where ZP reaches its peak value. This region is marked by the highest contour levels (dark brown in plot b of Figure 3), indicating the maximum ZP. This region indicates the levels of the input variables (T80 and BQ) that would yield maximum stability for the BQ SLN formulations [25].

4.7. Graphical Optimization to Define a Design Space for the BQ SLN Formulations

An overlay plot in the design space is a powerful visualization tool in DOE and RSM for optimizing processes involving multiple responses. It helps to identify the feasible region where all the specified criteria are met, enabling more informed and efficient decision making in experimental design and process optimization [27,28,29]. Statistical tool profilers can be used to relate the feasible region in an overlay plot to optimal settings. Such prediction profilers use a systematic approach to identify feasible regions, interactively adjust variables, and fine-tune settings to generate acceptable ranges to create a design space. This integration ensures a comprehensive understanding and optimization of the design space, leading to robust and efficient process optimization [29]. In this study, a preliminary design space was created for the BQ SLN formulations from estimates from the prediction profiler obtained from the JMP software. By linking the feasible region in an overlay plot with the optimal settings identified in the prediction profiler, one can generate acceptable ranges within the design space where the BQ SLN formulations can consistently maintain their desirable CQAs. This integration ensures a comprehensive understanding and optimization of the design space, leading to robust and efficient process optimization.

4.8. Future Studies

It is envisaged that further validation using the Monte Carlo simulation will be performed when larger-scale BQ SLN formulations are made. The Stochastic optimization (Monte Carlo simulation) involves defining the design space, developing a DOE model, and running simulations to assess the variability and reliability of the experimental outcomes. This approach helps identify risk areas and informs decisions for process optimization by comparing simulation results with experimental data to ensure consistency and robustness [30].

4.9. Study Limitations

This study represents our ongoing effort to discover the most effective approach for optimizing the process of formulating the BQ SLNs. We recognize that the initial models applied to optimize these formulation batches may not fully capture the complexities of the data. For example, it was observed that a reduced model was generated for the y variables when the study dataset was simulated using the k-fold cross-validation technique. This infers stronger predictive performance and optimal fit for the various fixed statistics or optimality criteria. In future studies, particularly during the scale-up stage, we plan to implement regularization techniques such as k-fold cross-validation or Lasso regression [31,32]. These methods will help refine our models and achieve more effective outcomes with optimal fit statistics.

5. Conclusions

The study highlights the importance of model adequacy in assessing the predictive performance of regression models. The findings indicate that while first-order models may have limited explanatory power, the quadratic model significantly captures nonlinear relationships between variables, allowing for the more precise optimization of formulations. To overcome these limitations, the study opted for a second-order model, which was believed to ensure better efficacy and stability of the BQ formulations by providing a more accurate and reliable framework for optimization. The study demonstrates that for optimal solid lipid nanoparticle formulation, maintaining BQ concentration between 10% and 20% and T80 concentration between 15% and 18% is crucial. This balance maximizes ZP and minimizes PdI, ensuring particle stability and appropriate size distribution. The high Precision Ratio (8.5101) and significant model fit (p < 0.0073 for PdI) validate the reliability of these findings. The results highlight the importance of the careful adjustment of BQ and T80 concentrations to achieve desired nanoparticle characteristics and optimize the formulation process for effective drug delivery systems. Finally, the broader significance of this study can be applied to the importance of precise model selection in pharmaceutical product optimization. Identifying and using the right regression model during product development is essential in improving formulations in drug delivery systems.