The hybrid hydrogels formulated in this study were prepared using alginate (alginic acid sodium salt derived from brown algae) and medium viscosity carboxymethyl cellulose (CMC) obtained from Sigma-Aldrich, St. Louis, MO, United States. The G-block of alginate supports gel formation, while the GM and M blocks contribute to increased flexibility. CMC is another water-soluble anionic bio-copolymer composed of β-D-glucose and β-D-glucopyranose-2-O-(carboxymethyl)-mono-sodium salt connected by β-1,4-glucosidic bonds (Han and Wang, 2017). CMC is commonly used as a thickener (Tongdeesoontorn et al., 2011) and is known for its non-toxic and non-allergenic nature. The presence of carboxyl groups in this cellulose derivative enhances its solubility, thickening properties, and stability (Han and Wang, 2017).

Figure 1A provides an overview of the hybrid hydrogel preparation process. Various weight percentages such as 0%–4% of alginate and CMC were mixed into water to prepare a set of Alginate-CMC (AC) compositions. The mixture was stirred overnight to ensure a homogeneous blend of alginate and CMC. Subsequently, the mixture was autoclaved at a temperature of 121°C for 20 min that confirms the sterilization by removing micro-organisms from the bio-inks to reduce the overall probability of contamination. The reduction of contamination of scaffold fabricated with living cells can improve the possibility of higher cell viability and proliferation. In Figure 1A, alginate and CMC are denoted by “A” and “C”, respectively, with the subscripts indicating the corresponding weight percentages.

A rotational rheometer (ARES-LS2, TA instruments, New Castle, DE, United States) equipped with parallel plate geometry (20 mm flat plate) was used to conduct rheological measurements. The measurements were taken at a gap width of 1.0 mm and a temperature of 25°C. Flow curves were generated by assessing the viscosities and shear stresses of various compositions across a range of shear strain rates from 1.0 to 100 s⁻¹. To identify the shear-thinning behavior of the materials considered and determine the coefficients, the Power-Law Equation (Eq. 1) was fitted to the linear region of the shear rate vs. viscosity plot (Therriault et al., 2007). η = K ′ γ ˙ n ′ − 1 (1) Where η is the viscosity, γ ˙ is the shear rate, K ′ and n ′ are shear-thinning coefficients. While the material is extruded through the nozzle, the material experiences shear stress across the cross-sectional area of the nozzle, which is larger near the nozzle wall.

An extrusion-based 3D bio-printer [BioX (CELLINK, Boston, MA)] has been used to fabricate the filaments and scaffolds. The generation of vectorized tool paths for the scaffold was achieved using a Visual Basic-based scripting language, which was then converted into a machine-readable format. The hybrid hydrogels were stored in a disposable barrel reservoir (EFD, Nordson, Westlake, OH, United States) and dispensed by applying air pressure through a dosing nozzle (EFD, Nordson) onto a stationary print bed. The scaffold was fabricated in a grid pattern, with layers of bio-ink extruded in a 0°/90° orientation. Filament width and pore size were measured using ImageJ software, with three random measurements (n = 3) taken for each data point, and the results were presented as the mean ± standard deviation. Table 1 represents various process parameters used in bioprinting experiments for this study.

All process parameters were selected based on the authors’ expertise and widely used in other literature. Process parameters were changed in a systematic way to capture the effects of specific process parameters. To enhance measurement precision, a pixel-based image analysis method was employed to determine filament width as mentioned in our earlier published work (Quigley et al., 2023). In short, microscopic images of deposited filaments were captured within a minute of printing. The scale of the microscope image was determined to calculate the width in pixels. Subsequently, MATLAB’s image batch processing toolbox was utilized for the analysis. The images were loaded into the program, processed, and converted to binary format. The filament width was then determined by identifying the coordinates of the initial and final pixels. These measurements were conducted thrice to obtain an average diameter for the printed filament. The discrepancy between the filament width and the diameter of the extrusion nozzle was defined as the diffusion rate [(actual filament width - nozzle diameter)/nozzle diameter × 100] in this study.

In the first stage, a simple descriptive analysis of data can help to understand the data and its related process behavior. The normality assumption is considered in regression-based modeling, however, for the larger dataset, the regression-based model is still very robust even if data are not perfectly normally distributed.

The Pearson correlation coefficient would be used initially to understand the linear relationship between process variables including the process output. Considering X ₁ and X ₂ represent two variables with having n amount of dataset, the Pearson correlation can be found using the following equation (Weiss and Weiss, 2017): C o r   x 1 , x 2 = ∑ i = 1 n x 1 i − x ¯ 1 x 2 i − x ¯ 2 ∑ i = 1 n x 1 i − x ¯ 1 2   ∑ i = 1 n x 2 i − x ¯ 2 2 (2) Here X ¯ represent the sample mean of each dataset. The Pearson correlation ranged between −1 and 1. The sign represents the relationship between variables whereas the quantitative value represents the strength of the relationship.

As it is intended to apply the machine learning approach to develop our predictive model, it is necessary to split our data into two parts: training and test datasets. The training data would be used to establish different regression-based machine learning models whereas test data would be used to validate and compare the model performances. The data would be randomly split between 75% and 25% of the original dataset approximately for training and test purposes, respectively.

The 3D bio-printed structure’s output and process parameters are a quantitative type of variable that generates level-type data. Therefore, the linear regression-based method will be implemented considering the supervised machine learning modeling. The general least-square approach will be used for regression coefficient estimates. As there exist several process parameters in this study, the multiple-regression modeling will be utilized. The process parameters in engineering applications have higher tendencies to interact and therefore, it is considering all possible interaction terms in the model selection approach. However, higher-order polynomials will not be considered in this work as they have little to no meaning in explaining the 3D bio-printing physical process. Let’s consider a random variable Y represents the filament width of the printed object and a random variable X represents process parameters that influence the filament width. Then the regression model would be, y i = β 0 + ∑ j = 1 q β j x j + ∑ j , k = 1 ,   j ≠ k q β j k x j x k + … … … + ∑ j , k , l … q = 1 ,   i ≠ j ≠ … ≠ q q β j k l … q   x j x k … x q + ϵ (3)

Here q represents the number of process parameters for the 3D bioprinting process, β′s represents the coefficient of the regression equation, and ε represents the error term. The residual sum of squares (RSS) then can be estimated by Eq. 4, R S S = ∑ i = 1 n y i − y ^ i 2 (4)

The total sum of squares (TSS) can be obtained by Eq. 5, T S S = ∑ i = 1 n y i − y ¯ 2 (5)

Both Eqs 4, 5 would be useful for the selection of different statistical model criteria in the later section. The β coefficients of a regression model are estimated considering minimizing the Eq. 4. For a decently larger number of process parameters, the number of total variables will be exponentially increased especially as it is considered all levels of interaction effects between process parameters. To handle this large amount of regression variable it is necessary to have a systematic process to select the best variables which influence the filament width with statistical significance. Three classic stepwise methods are used to select the best variables for the multiple regression model: forward, backward, and mixed stepwise selection methods (Harrell, 2001; Zhang, 2016). The forward selection method started with only intercept and with no variable. Step by step it added one after one variable and stopped the process once the selection criteria were met. The backward selection method starts with all possible variables and step by step it removes each variable following the selection criteria until it reaches an optimal number of variables. The backward selection method requires a larger number of datasets compared to the number of variables and the forward selection method might include some variables very early which become redundant later. The mixed model is a combination of forward and backward section processes where it adds variables to the model but may be removed variables later steps if the variable does not improve the model fit (James et al., 2013). The mixed model method has been adopted in this work for the best set of variable selections for the regression-based machine learning model.

It has been considered three prominent statistical criteria for the variable selection method: adjusted R-square (Adj. R²) Mallow’s Cp, and Bayesian information criterion (BIC) (Zhang, 2016). The traditional R-square value is influenced by the number of variables and hence it is adjusted. Using Eqs 4, 5, the adjusted R-square can be estimated by Eq. 6, A d j   R 2 = 1 − R S S n − 1   T S S n − d − 1 (6)

Here, n is the total number of data points and d represents the total number of all possible variables including interaction terms. As these model selections are implemented using the training data therefore total number of data would be the total number of the training dataset instead of the full dataset. The adjusted R-square ranged between 0 and 1 and it interpreted how much the regression model can explain the variability of the data. The larger adjusted R-square value indicates a stronger model fit.

The second model criterion we adopted in this work is Mallow’s Cp measurement. The Cp can be estimated from the Eq. 7, C p = 1 n R S S + 2 d S 2 (7)

Here S 2 represents the unbiased estimator of residual variance and is calculated as S 2 = R S S / n − d . The Cp is a well-known method for addressing model overfitting. Unlike the adjusted R-square, lower Cp represents a better model fit.

The third model criterion is the Bayesian information criterion (BIC) which considers the Bayesian perspective in model selection. The BIC can be computed using Eq. 8, B I C = 1 n R S S + ⁡ log n d S 2 (8)

Similar to Cp criteria, the model with lower BIC represents the stronger model fit. Due to the logarithmic penalty in BIC estimation, typically it heavily penalizes a larger number of variable inclusions in the model, thus BIC criteria usually have a smaller number of variables in the model (James et al., 2013).

In stepwise method mentioned in the previous section used training data to select the best set of variables for a model fit. As those machine learning models are obtained using training data there are high chances of model overfitting meaning the models perform very well with the training dataset however with the new dataset it does not perform satisfactorily. To address this overfitting issue we will utilize the cross-validation method (James et al., 2013). The test Mean Square Error (MSE) measurement will be used to perform the cross-validation of previously selected models. The test MSE can be obtained by following the equation, M S E = 1 n ∑ i = 1 n y i − y ^ i 2   (9)

In variable selection, we have used Adj. R², Cp, and BIC criteria to choose the best-fit model. All these models will be used to estimate test MSE using the test dataset for selecting the best-performed model. The lower test MSE indicates that the model performs better with unknown or new datasets. The training MSE can be also estimated by Eq. 9 using the training dataset instead of the test dataset.

The cross-validation method would provide the best regression-based machine learning model considering the current dataset. The model coefficient would vary based on training and test data. To address this uncertainty issue and increase the mode coefficient’s estimation accuracy, the bootstrap method is adopted. The bootstrap is a type of nonparametric Monte Carlo method that estimates coefficients using the resampling of data (Efron, 1992). The bootstrap process has been depicted in Figure 2 method framework in detail. In the first step of the bootstrap method, the sample size needs to be determined. In this work, it is considered that the sample size would be the same as the number of original datasets which is n. In the next step, generate a random sample also known as a resample from the original dataset using the random replacement of data points. The usual resample replicates (B) are 50 to 100 in number, however, for higher accuracy of estimation a larger number of replicates is considered such as 10,000 replicates (Limon et al., 2016). Once the resampling replicates are generated then for each replicate the model coefficients (βi) are estimated. Assume the bootstrap replicates for b = 1,2,3, ….B, the bootstrap estimate of the coefficient can be calculated using Eq. 10, β ^ i * ¯ = 1 B ∑ b = 1 B β ^ i b (10)

As a non-parametric type estimation, bootstrap estimates are free from any assumption of distribution fit.

The best-fit model using the test data validation with the bootstrap coefficient estimate would be our final regression-based machine learning model. Several new experiments could be conducted to compare and further validate the experimental scaffold data with machine learning model data. When there are new sets of more experimental data, the total framework can be re-run to update the machine learning model to search for a better model fit with higher accuracy of the coefficient estimates.

Out of 157 data, we have selected 8 of them that were obtained for various process parameters such as 110 kPa and 170 kPa EP, 5 mms⁻¹ and 15 mms⁻¹ PS, 0.21 mm and 0.41 mm ND to discuss the usual effects of them on filament width. All possible combinations are (EP, PS, ND): (110, 5, 0.21), (110, 15, 0.21), (110, 5, 0.41), (110, 15, 0.41), (170, 5, 0.21), (170, 15, 0.21), (170, 5, 0.41), and (170, 15, 0.41). Filament widths for all those combinations are shown in Figure 3. A set of phenomena was observed such as (i) with increasing EP, FW increased for constant ND and PS, (ii) bigger ND extrude filament with larger FW for constant EP and PS, (iii) with increasing PS, FW increased for constant ND and EP. For instance, when the EP was raised from 110 kPa to 170 kPa, a 115% increase in FW was measured during the extrusion process through a 0.21 ND at a PS of 5 mms⁻¹. On the other hand, a similar EP change and ND size, but with a higher PS of 15 mms⁻¹, resulted in a 110% increment in FW. ZH and WF were considered 0.3 mm and 8%, respectively for all combinations. Figure 3C shows filament widths with respect to PS and EP for 51 filaments extruded through 0.41 mm ND maintaining 8% solid content and 522 Pa s viscosity. Figure 3D represents FWs with respect to PS and EP for 38 filaments extruded through 0.21 mm ND maintaining 8% solid content and 522 Pa s viscosity. Both figures even show some trends such as (1) FW increases with increasing EP, (2) EP reduces with increasing PS, and (3) FW increases with increasing ND, the change of FW may not be straightforward to estimate and can be material properties and process parameters dependent (Chung et al., 2013; Diamantides et al., 2017; Habib et al., 2018). Therefore, for a large set of data, ML can be an efficient tool to identify a predictive model to estimate FW.

As mentioned in Section 2, our case study dataset contains five process parameters, and the filament width of the printed object would be the dependent output variable. In total, there are 157 experimental data for this demonstration case study purpose. All analysis of this work has been performed by R Statistical Software (v4.3.1, R Core Team 2023). Figure 4 shows the scatterplot of the relative relation between process parameters including the output variable. The different data density of the scatterplot indicates that all combination of process parameter is not varied in the same proportion during the experimental test.

The scatterplot in Figure 4 provides a visual depiction of the parameters’ relation, however, it does not provide the strength of the relation. The correlation matrix in Figure 5 provides the strength and direction of the variable’s relation using Eq. 2. It is very clear from the correlation matrix that several variables are significantly correlated such as extrusion pressure (EP) with print speed (PS) or material’s weight fraction (WF) with Z-height (ZH) between print nozzle and print bed. This correlation also justifies our initial assumption of possible interaction between process parameters and considers them during the best variable selection method.

As we mentioned in Section 3.2, after the approximately 75/25 random split of the original dataset of 157 data, there are 122 data points for training the machine learning model and 35 data points for testing the model obtained using the training dataset. Considering all possible interactions for five process parameters there are a total of 31 variables. However, while we ran the mixed-type stepwise regression model selection algorithm, 15 of the interaction variables were omitted due to linear dependencies. A higher number of training datasets would help to reduce this issue. Therefore, the mixed stepwise selection algorithm would search for the best model fit among 16 variables using the training dataset.

Table 2 and Figure 6 show the stepwise iteration results considering the three different model criteria. The bold number in Table 2 represents the optimal value of corresponding model criteria whereas the red marks in Figure 4 indicate the optimal number of variables considering each model criteria. Though it shows only the number of variables, the specific variables are not the same in each iteration when it considers different model criteria. It is found that Adj R² criteria reached optimal variable number selection by 9 with the highest value of 0.6491 or 64.91%, the Cp criteria reached optimal variable number by 7 with the lowest value of −12.2301, and the BIC criteria reached optimal variable number by 4 with the lowest value of −97.4532. As mentioned in Section 3, the Cp and BIC methods put extra penalties for additional variables inclusion in the model. Thus, the number of variables is higher in the Adj R² method than in the Cp method and lowest in the BIC method. For our case study, the Adj R² is not in the upper value. However, this model could be useful for predicting the porosity of the printed scaffolds. For the demonstrated case study, we have considered the most critical process parameters in our work for efficient and cost-effective experimentation. Considering a few more environmental variables such as curing temperature and humidity through the Design of the Experiment method can increase the accuracy of the model. In case of a higher number of outliers present, data transformation could further improve the prediction model as well. In our future work, we plan to integrate DoE [showed better R² in our earlier work (Quigley et al., 2023)] and ML, increasing variables, and appropriate data transfer to potentially improve the resulting Adj-R².

The combination of selected variables under each model criterion is listed in Table 3. From all three model criteria, the best-fit model has a heavy presence of interaction terms. This is intuitive because process parameters interact with each other in a 3D printing process. Additionally, this confirms correlation matrix results that indicate a higher correlation between variables.

All model best-fit criteria have been derived using training data. However, the performance of the best model will be cross-validated with the test data. The test MSE has been considered to cross-validate the training best-fit model listed in Table 3. Table 4 tabulated both the training and test MSE of each model obtained using the training dataset. It is observed that for every case, the training MSE is lower than the test MSE. This is very common because the training MSE uses the same data it used for model selection. The test MSE of Adj. R², Cp, and BIC models are 0.0816, 0.0841, and 0.0877, respectively. Thus, the Adj. R² model with 9 variables becomes the best performance machine learning model for our current dataset. The percentage difference of test MSE compared with the best model is also listed in Table 4. Based on the difference it can be concluded that the other two models are also very competitive compared to the best performance model.

The residual analysis of the best performance model is conducted which is the Adj. R² criteria-based model. The residual analysis is conducted using the test dataset. Figure 7 depicts the residual plot and QQ plot based on the test dataset for the best model. Both the figures confirm that data are fairly randomized and follow normality with few outliers both in the upper and lower ends of data.

The presence of an outlier can heavily impact the model coefficient estimate. In machine learning models, data are typically in large numbers with higher chances of being skewed in nature. The Bootstrap estimation method mentioned in Section 3 can significantly help to reduce the uncertainty, outlier issues, and asymmetric nature of data. Table 5 lists the coefficient estimate for the best performance model obtained after the cross-validation method. It is observed that the variation of coefficient estimates among the training and test datasets. In the current case scenario, one of the variables does not have an estimate for test datasets due to the test datasets are typically lower in number. The Bootstrap resampling method generates new sample data using randomized with replacement and estimates the coefficient for larger datasets. Therefore, the Bootstrap estimate is a very robust estimate of the coefficient and very useful especially when datasets are further departed from normality. The last column of Table 5 shows each variable estimate using 10,000 resamples of the Bootstrap method.

Once the coefficient value of the best performance model is estimated, then the final predictive model based on Eq. 3 is established. To further understand and validate our machine learning model, there are conducted 4 more new 3D bio print experiments with different sets of process parameters. The experimental FW (FW-Exp.) would be measured in a similar fashion mentioned in Section 2. The process parameter value of these new experiment samples is used in the best-fit machine learning model and estimates the FW (FW-Mod). The different sets of process parameters, experimental and model-based FW are listed in Table 6. The comparison between experimental and model FW is depicted in Figure 8. The percentage error indicates the model fairly closely predicted the FW within an acceptable range.