Modeling the Additive Effects of Nanoparticles and Polymers on Hydrogel Mechanical Properties Using Multifactor Analysis

Interpenetrating networks (IPN)s have been conceived as a biomimetic tool to tune hydrogel mechanical properties to the desired target formulations. In this study, the rheological behavior of acrylamide (AAm) [2.5–10%] hydrogels crosslinked with N,N′-methylenebis(acrylamide) (Bis) [0.0625–0.25%] was characterized in terms of the saturation modulus affected by the interaction of silica nanoparticle (SiNP) nanofillers [0–5%] and dextran [0–2%] at a frequency of 1 Hz and strain rate of 1% after a gelation period of 90 min. For single-network hydrogels, a prominent transition was observed at 0.125% Bis for 2.5% AAm and 0.25% Bis for 5% AAm across the SiNP concentrations and was validated by retrospective 3-level factorial design models, as characterized by deviation from linearity in the saturation region (R2 = 0.86). IPN hydrogels resulting from the addition of dextran to the single network in the pre-saturation region, as outlined by the strong goodness of fit (R2= 0.99), exhibited a correlated increase in the elastic (G’) and viscous moduli (G”). While increasing the dextran concentrations [0–2%] and MW [100 kDa and 500 kDa] regulated the increase in G’, saturation in G” or the loss tangent (tan(δ)) was not recorded within the observed operating windows. Results of multifactor analysis conducted on Han plots in terms of the elastic gains indicate that amongst the factors modulating the viscoelasticity of the IPN hydrogels, dextran concentration is the most important (RDex = 35.3 dB), followed by nanoparticle concentration (RSiNP = 7.7 dB) and dextran molecular weight (RMW = 2.9 dB). The results demonstrate how the Han plot may be systematically used to quantify the main effects of intensive thermodynamic properties on rheological phase transition in interpenetrating networks where traditional multifactor analyses cannot resolve statistical significance.


Introduction
Hydrogels are highly porous three-dimensional polymer networks with physical and biochemical properties that make them ideal for many biomedical applications such as tissue scaffolds, drug delivery, cell supports, and corneal implants [1][2][3][4][5]. However, unmodified hydrogels typically have low mechanical strength, which limits their applicability such as use for load bearing tissue scaffolds, for example, bone regeneration and artificial cartilage [1,3,6]. Recent studies have developed a number of ways to improve the mechanical properties of hydrogels including through the incorporation of additives such as polymers and nanoparticles [1,3,[6][7][8][9][10][11][12]. The use of polymer additives, either to create interpenetrating polymer networks (IPN) in which the hydrogel is made of two different crosslinked polymers, or semi-IPNs in which only one of the polymers is crosslinked and the other polymer is trapped within this crosslinked network, is a common way to improve the mechanical properties [1,[13][14][15][16][17][18][19][20]. For example, one study found that chitosan (CS)-poly(acrylic acid) (PAA) IPNs resulted in a hydrogel with higher tensile strength and adsorption capacity than either the CS or PAA hydrogels, resulting in a wider range of applications for the hydrogel [21]. Incorporating nanoparticles such as clay nanoparticles and graphene oxide

Polymerization Reaction
Single network hydrogel samples (pAAm) were prepared as previously reported [24]. Briefly, AAm and Bis stocks were diluted to their desired concentrations in pH 7.2, 250 mM Tris-HCl buffer, followed by the addition of TEMED (0.1% of the final reaction volume) and 10% w/v APS solution (1% of the final reaction volume). For the pAAm-dextran IPN hydrogel samples, 100 kDa and 500 kDa dextran was first dissolved in Tris-HCl buffer by stirring at room temperature, and then diluted to the desired concentrations before adding it to the reaction mixture prior to the addition of APS and TEMED. Similarly, for the nanocomposite hydrogels, various amounts of SiNPs were added to the single network or IPN hydrogel reaction mixtures prior to the addition of APS and TEMED.

Rheological Measurements
Rheological measurements of the single network and IPN hydrogels with or without nanoparticles were carried out using the MCR302 rotational rheometer (Anton Paar, Graz, Austria) using parallel plate geometry, as described previously [27]. Briefly, 500 µL of a well-mixed reaction mixture was pipetted onto the lower plate of the rheometer. The upper plate was lowered until the desired gap distance (1 mm) was achieved. Amplitude sweeps at a constant frequency of 1 Hz were then carried out to ensure measurements were carried out in the linear viscoelastic regime of the hydrogels. Next, dynamic sweep tests over frequencies ranging from 0.1 to 100 Hz were recorded in the linear viscoelastic regimes (strain amplitude of 0.01) to determine the shear elastic modulus. Elastic and viscous moduli were determined by following the gelation for 90 min (to ensure complete gelation; gelation usually occurs within 20 min) at 1 Hz and 1% strain for all samples. Relative elastic (G') and viscous moduli (G") were calculated by normalizing the values for the pAAm-SiNP hydrogels to the corresponding values for the control pAAm gels. The moduli are related by the tangent of the loss factor (delta) given by Equation (1).
All measurements were carried out at room temperature and the final parameters were reported as an average of three independent measurements.

Factorial Design
A retrospective design of experiments was carried out for the single network hydrogels using the elastic shear modulus (G') as a response. The effect of curvature was examined by combinations of 3 level factorial designs generated using JMP-SAS v.16 (JMP Statistical Discovery LLC, Cary, NC, USA) and conducted over the experimental ranges outlined in Tables 1 and 2. Results were analyzed according to the methodology proposed by Box et al. [43] and subsequently validated by software. The levels of coded and uncoded variables as well as respective responses and associated hydrogel network types are presented throughout the study.  Table 2. Interpenetrating network factor levels for the MANOVA analyses for investigation into the G', G", and tan(δ) responses.

Multivariate Analysis of Variance
Multivariate analysis of variance (MANOVA) was conducted for the IPN hydrogels using MATLAB version r2021b with a significance threshold of 0.05 according to the variable and levels outlined in Table 2. The concentration of the acrylamide (AAm) was fixed at 5%. The responses were G', G", and tan(δ).

Discriminant Analysis Using Han Plots
Han plots were generated to discriminate between the elastic and viscous rheological behaviors of the polymer networks. This method is used to investigate the effect of intensive thermodynamic properties on rheological phase transition [44,45]. As polymer, cross-linker, and filler concentrations are also intensive properties governing the viscoelastic state, the use of the method was justified. For each plot, the Euclidean distance between the discrimination axis and datapoints, denoted as elastic gain, was calculated and averaged across levels for the IPN hydrogel design variables outlined in Table 2 and ranked in order of importance using Equations (2) and (3).
Difference, D i , was given by Equation (2), where X is the Euclidean distance at a given level averaged over the number of levels (L). Rank, R, expressing the significance of the variable in terms of the range of difference D i , is given by Equation (3).
Equations (2) and (3) will also be used to evaluate the change in elasticity in terms of G', from the single network initial condition to the final IPN hydrogel configuration expressed by the variable ∆G'. A positive difference indicates an increase in G' due to dextran incorporation.

Saturation in Nanoparticle Mediated Enhancement of Hydrogel Mechanical Properties
First, the impact of the incorporation of nanoparticles as well as different polymer and crosslinker concentrations on the elastic modulus (G') was investigated. Various concentrations of silica nanoparticles, AAm polymers, and Bis crosslinkers were used to create single network hydrogels, and the G' values were investigated for various formulations of single network hydrogels. Upon preliminary investigation, the hydrogels exhibited a saturation at 4% SiNPs, which can be visualized in the plot of G' values versus SiNP percentage (Figure 1a-c). The hydrogels also exhibited a saturation trend with respect to Bis concentration, as shown in Figure 1d-f. The saturation behaviors for the hydrogels with respect to Bis and SiNP concentrations were less pronounced for the highest AAm concentration of 10%.

Factorial Design Models Identify the Significance of Hydrogel Nanocomposite Parameters
To further model the curvature observed at higher nanoparticle concentrations using the design parameters, 3 level factorial designs (N = 27 runs) using either high or low Bis concentration ranges were analyzed using the values presented in Tables S1-S3 at higher SiNP values. The regression equation for the 3 3 design is given by Equation (4). Parameters and their respective significance levels are reported in Table 3. Ŷ = β 0 + β 1 X 1 + β 2 X 2 + β 3 X 3 + β 12 X 1 X 2 + β 23 X 2 X 3 + β 13 X 1 X 3 + β 123 X 1 X 2 X 3 At lower Bis concentrations, the regression model was dominated by the effects of AAm (X1) and SiNP (X3) at a significance threshold of 0.1 (p = 0.1), as deviation from linearity is characterized by an R 2 value of 0.86 (Figure 2a). At higher Bis concentration, the regression model was dominated by the effects of AAm (X1), Bis (X2), and their respective interaction (X1 × X2) (p = 0.1), characterized by an adequate coefficient of determination of

Factorial Design Models Identify the Significance of Hydrogel Nanocomposite Parameters
To further model the curvature observed at higher nanoparticle concentrations using the design parameters, 3 level factorial designs (N = 27 runs) using either high or low Bis concentration ranges were analyzed using the values presented in Tables S1-S3 at higher SiNP values. The regression equation for the 3 3 design is given by Equation (4). Parameters and their respective significance levels are reported in Table 3.
At lower Bis concentrations, the regression model was dominated by the effects of AAm (X1) and SiNP (X3) at a significance threshold of 0.1 (p = 0.1), as deviation from linearity is characterized by an R 2 value of 0.86 (Figure 2a). At higher Bis concentration, the regression model was dominated by the effects of AAm (X1), Bis (X2), and their respective interaction (X1 × X2) (p = 0.1), characterized by an adequate coefficient of determination of 0.97. For the former, the effect of curvature could be observed in the mid-section of the surface plot due to the gap between the experimental data and predicted model ( Figure   : 4% SiNP, and : 5% SiNP and (d) 2.5% AAm, (e) 5% AAm, or (f) 10% AAm. N = 3 for each hydrogel preparation.

Factorial Design Models Identify the Significance of Hydrogel Nanocomposite Parameters
To further model the curvature observed at higher nanoparticle concentrations using the design parameters, 3 level factorial designs (N = 27 runs) using either high or low Bis concentration ranges were analyzed using the values presented in Tables S1-S3 at higher SiNP values. The regression equation for the 3 3 design is given by Equation (4). Parameters and their respective significance levels are reported in Table 3. At lower Bis concentrations, the regression model was dominated by the effects of AAm (X 1 ) and SiNP (X 3 ) at a significance threshold of 0.1 (p = 0.1), as deviation from linearity is characterized by an R 2 value of 0.86 (Figure 2a). At higher Bis concentration, the regression model was dominated by the effects of AAm (X 1 ), Bis (X 2 ), and their respective interaction (X 1 × X 2 ) (p = 0.1), characterized by an adequate coefficient of determination of 0.97. For the former, the effect of curvature could be observed in the mid-section of the surface plot due to the gap between the experimental data and predicted model (Figure 2a). For the latter, the superimposition of the empirical and theoretical was seamless, except for the lower limit region of the AAm (2.5%) (Figure 2b). depicted in Figure 1, with 2a corresponding to low Bis, high silica nanoparticles, 2b corr to high Bis, high silica nanoparticles, and 2c corresponding to low Bis, low silica nanopa els, respectively. In order to operate in a region where the rheological behavior of the system the saturation region, as presented in Figure 2a (G' vs. SiNP), another 3 3 factori outlined in Table S3, was conducted in the lower SiNP range. This model is dom the main effects (X1, X2, X3) and AAm/Bis (X1 × X2) as well AAm/SiNP (X1 × X3) int The fact that the significance threshold has the same magnitude (p < 0.0001) for effects as well as the interactions is indicative of the onset of the competitive me for crosslinking between the Bis and SiNP. The strength of the model is reflected value of 0.9912 and the superimposition of the theoretical and experimental fi the surface response plots except for the lower limit region of the AAm (2.5%) (F Combining the analytical results from the models above, the operating windo IPN hydrogels at the mid AAm level (5%) and below saturation was proposed.  Table 3, where X1 is represents the AAm concentration.  (4) is shown in grey to which the actual experimental surface (Y) is superimposed in black for all three models in Table 3, where X 1 is represents the AAm concentration.
In order to operate in a region where the rheological behavior of the system is below the saturation region, as presented in Figure 2a (G' vs. SiNP), another 3 3 factorial design, outlined in Table S3, was conducted in the lower SiNP range. This model is dominated by the main effects (X 1 , X 2 , X 3 ) and AAm/Bis (X 1 × X 2 ) as well AAm/SiNP (X 1 × X 3 ) interactions. The fact that the significance threshold has the same magnitude (p < 0.0001) for the main effects as well as the interactions is indicative of the onset of the competitive mechanisms for crosslinking between the Bis and SiNP. The strength of the model is reflected by an R 2 value of 0.9912 and the superimposition of the theoretical and experimental findings in the surface response plots except for the lower limit region of the AAm (2.5%) (Figure 2c). Combining the analytical results from the models above, the operating window for the IPN hydrogels at the mid AAm level (5%) and below saturation was proposed.

Design of the Interpenetrating Network Hydrogel Nanocomposite Studies
Dextran was used to make IPN hydrogels to study the combined effects of dextran and SiNPs on the hydrogel elastic moduli, as previous studies have demonstrated additive effects of different fillers on hydrogel properties [1,3,6-12]. As described previously, one of the goals of the work was to investigate the individual and additive impacts of incorporating different fillers. It was therefore important to use SiNP values in the linear range of the single network. In a similar vein, Bis concentrations below saturation were used-in the previous study, the G' values saturated at 0.25% Bis, so Bis concentrations at and below 0.25% were selected [29]. Dextran concentrations and molecular weights were chosen based on the literature [4]. Additionally, a single AAm concentration was selected to study the interpenetrating networks because AAm (X 1 ) was a highly significant factor for the single network data throughout the three models proposed by Equation (4), as the goal was to focus more on the impact of dextran and the interactions between the variables for the IPN hydrogel studies.

Effect of Dextran on the Elastic and Viscous Properties of Hydrogel Nanocomposites
Using the parameter selection outlined in Section 3.3., pAAm-Dex-SiNP IPN hydrogel samples were prepared to investigate the role of dextran on hydrogel nanocomposites. Although there was no evidence of the saturation of G' for dextran with a molecular weight of 100 kDa (Figure 3a), an onset of saturation was observed at 1% dextran for the higher molecular weight (500 kDa) (Figure 3b). This saturation was more prominent for hydrogel nanocomposites prepared with a higher concentration of Bis ( Figure S1 and Table S4). Furthermore, the effect of dextran on G' was more pronounced for pAAm hydrogels with SiNP (2% SiNP samples) compared to those without SiNP (0% SiNP). Collectively, the results indicate an interaction between dextran and SiNP for the IPN hydrogel nanocomposites that was more pronounced at higher Bis concentrations.
Leveraging previous results that demonstrated the independent effect of nanoparticles and polymer fillers on hydrogel viscous modulus, this study also investigated the combined impact of these fillers on the viscous modulus of pAAm hydrogels [28,39,40]. As the dextran concentration and molecular weight increased, there was also an increase in G" (Figure 3c,d). However, in contrast with G', there was no saturation for G", and correspondingly for tan(δ), at higher concentrations of dextran. These results suggest that the maximum potential for fluidity of the pAAm-Dex IPN hydrogel nanocomposites has not been reached (Figure 3e,f)  Leveraging previous results that demonstrated the independent effect of nanoparticles and polymer fillers on hydrogel viscous modulus, this study also investigated the combined impact of these fillers on the viscous modulus of pAAm hydrogels [28,39,40]. As the dextran concentration and molecular weight increased, there was also an increase in G" (Figure 3c,d). However, in contrast with G', there was no saturation for G", and correspondingly for tan(δ), at higher concentrations of dextran. These results suggest that the maximum potential for fluidity of the pAAm-Dex IPN hydrogel nanocomposites has not been reached (Figure 3e,f).

Statistical Analyses of the Interpenetrating Network Hydrogel Nanocomposites
Multifactor statistical methods aimed at characterizing the additive effect of dextran enabling the formation of the IPN hydrogels are summarized below. Results of the MANOVA analysis investigating the main effects in terms of the rheological parameters for the interpenetrating network are presented in Table S4. All factors were significant (p < 0.05) for the related rheological parameters examined, the exception being the MW of the dextran for the viscous modulus and loss factor.
As this preliminary result does not lend itself to a quantitative pareto effect, the analysis was extended to Han plots. Han plots displaying the discrimination between elastic

Statistical Analyses of the Interpenetrating Network Hydrogel Nanocomposites
Multifactor statistical methods aimed at characterizing the additive effect of dextran enabling the formation of the IPN hydrogels are summarized below. Results of the MANOVA analysis investigating the main effects in terms of the rheological parameters for the interpenetrating network are presented in Table S4. All factors were significant (p < 0.05) for the related rheological parameters examined, the exception being the MW of the dextran for the viscous modulus and loss factor.
As this preliminary result does not lend itself to a quantitative pareto effect, the analysis was extended to Han plots. Han plots displaying the discrimination between elastic and viscous behavior as a function of dextran concentration are presented across the SiNP, Bis, and molecular weights of dextran in Figure 4a-d. All datapoints reside over the 45 • discrimination axis, indicating that the system is cross-linked. Furthermore, the data demonstrates the effect of dextran on the viscous modulus of the IPN hydrogel nanocomposites, as exhibited by the higher elastic than viscous behavior at lower dextran concentrations. and viscous behavior as a function of dextran concentration are presented across the SiNP, Bis, and molecular weights of dextran in Figure 4a-d. All datapoints reside over the 45° discrimination axis, indicating that the system is cross-linked. Furthermore, the data demonstrates the effect of dextran on the viscous modulus of the IPN hydrogel nanocomposites, as exhibited by the higher elastic than viscous behavior at lower dextran concentrations.    The increase in the elastic moduli computed from the comparison of single and interpenetrating network formulations due to dextran incorporation, in terms of G' values, are presented in Table 4. In terms of ranking the factors modulating the elastic behavior, dextran concentration was the most important (RDex = 4.46 dB), followed by the effects of the dextran MW (RMW = 1.53 dB) and nanoparticles (RSiNP = 1.21 dB).  The increase in the elastic moduli computed from the comparison of single and interpenetrating network formulations due to dextran incorporation, in terms of ∆G' values, are presented in Table 4. In terms of ranking the factors modulating the elastic behavior, dextran concentration was the most important (R Dex = 4.46 dB), followed by the effects of the dextran MW (R MW = 1.53 dB) and nanoparticles (R SiNP = 1.21 dB).

Discussion and Future Work
The goal of this work was to investigate the impact of incorporating two different additives on the mechanical properties of hydrogels. Other groups have previously explored combining two different additives: (a) nanoparticles to form nanocomposites, and (b) other polymers to form an IPN hydrogel to increase the mechanical strength of hydrogels [29,39]. To combine these approaches, dextran was added to the pAAm-SiNP nanocomposite hydrogels to form an IPN hydrogel nanocomposite. By combining SiNPs and dextran, the mechanical properties of the hydrogels increased beyond the saturation point for SiNPs, consistent with studies reported by other groups [24,46,47]. Results from these experiments also indicated a saturation behavior in the increases in elastic modulus for the IPN hydrogel nanocomposites, which indicates the existence of a 'global' saturation point, beyond the 'local' saturation observed for either Bis-or SiNP-mediated enhancements. Interestingly, the viscous modulus did not saturate, and, in fact, increased with increasing dextran concentrations, as confirmed by the Han plots ( Figure 4). As stated in previous literature, the incorporation of additives that do not lead to chemical crosslinking into hydrogels can lead to increases in G" [29,39]. These results demonstrate that the use of different additives may allow us to independently control the elastic and viscous modulus of hydrogels.
Statistical modeling was used to interpret the data and plan subsequent experiments. For the single network hydrogel, a factorial design was used to model the data, thereby narrowing down the experimental space for the IPN hydrogel tests, a by-product of which may be the reduction of environmental waste. First, the impact of incorporating nanoparticles (SiNPs) into pAAm hydrogels was analyzed. The G' values for pAAm hydrogels saturated at higher concentrations of SiNPs, indicating a maximum impact of SiNPs on the elastic modulus. This saturation behavior was observed using a 3 3 factorial design, which illustrated the curvature of the system in a 3D space, along with 2D plots of the data, which also present the curvature (Figure 2a-c). The experimental space bound by the ranges of the variables in Equation (4) highlight the saturation of G' with increasing concentrations of Bis and Aam, and, as a result, subsequent experiments focused on low concentrations of Bis, and one concentration of AAm, to better see the impact of dextran on the hydrogels. As the use of factorial ANOVA did not discriminate between the main effects of the IPN hydrogel design matrix outlined in Table 2, the Euclidean distance to the 45 • discrimination axis in the Han plot was used as an optimization parameter to rank the factors ( Figure 5).
Future studies will encompass the effect of frequency on complex moduli explored through the Carreau Yasuda model used for the characterization of shear-thinning behavior of the modeled polymer nanocomposite formulations [48,49]. Frequency sweep analysis will in turn shed light on whether the rheological transition/saturation as a function of SiNP rheology modifiers is accompanied by the previously recorded conductive percolation limit [27] to expand the extrusion-based 3D printing processing capabilities of the custom AAm networks as well as expand the applications to SMART conductive gels [50]. The applications of this work fall into two major categories: (1) increasing the mechanical strength and (2) independent control of viscous and elastic modulus. First, the addition of dextran enabled us to increase the mechanical properties of the single network nanocomposites, exceeding the saturation point due to the addition of SiNPs, thus indicating an additive effect between the SiNPs and dextran. This allows for applications that require high mechanical strength such as its use as bone or cartilage scaffolds for tissue engineering [1,3]. Second, adding dextran allows researchers to tailor hydrogels with specific mechanical properties because the addition of dextran regulated increases in G' and G" to different extents. For instance, different concentrations of dextran would allow for preparations of hydrogel nanocomposites with similar values of G', but significantly different values of G". These efforts will entail frequency dependent studies to confirm the correlated spectra of the elastic versus viscous modulus across custom multivariate nanocomposites. Future applications will include the development of hydrogels for drug delivery and 3D bioprinting, where tunable viscoelastic properties may increase the ease of printing and extrusion [51][52][53].

Funding:
The research received no external funding. The work was supported by an internal grant from the School of Engineering at Santa Clara University. The funders had no role in the study design, data collection, or analysis. Data Availability Statement: Data will be made available upon request to the corresponding author.