Addressing the stiffness–toughness conflict in hybrid double-network hydrogels through a design of experiments approach

Abstract

An open challenge in soft matter science is the ability to create hydrogels that are soft but also have high fracture energy. A possible solution to this stiffness–toughness conflict has arisen through a class of material known as hybrid double-network hydrogels, which combine two polymeric networks with diametrically opposed chain stiffness and with both covalent and physical crosslinking. The vast parameter space inherent to such a system means it is difficult to identify the precise compositional parameters that lead to both high toughness and low stiffness. In this work, we address this challenge through a design of experiments (DoE) framework used to establish the statistical relationship between factors and mechanical properties of a hybrid double-network hydrogel. The crosslinking density of the networks is noted to play a prominent role in determining the stiffness of the hydrogel, while the network characteristics of the ductile network determine the toughness of the hydrogel. We also report that contrary to observations in current literature, it is possible to toughen the hydrogel without stiffening it. Therefore, the present experimentation and optimization exercise provides a hands-on guide for the use of DoE to determine the conditions for optimised mechanical properties of thin hybrid double-network hydrogels for various applications.

---

1 Introduction

Hydrogels are water-swollen polymer networks with a multitude of practical uses across medical devices and wound dressings and as structural materials in food and personal care products. Hydrogels constructed from a single polymeric network are typically mechanically weak with low fracture energies¹ and low elastic moduli.² Many hydrogels have been shown to exhibit intriguing mechanical behaviours such as stress relaxation that mimics soft tissues, such as the skin, tendons, and ligaments^3–7 making them good candidates for applications in biomedical engineering,^8,9 and tissue engineering.¹⁰ However, prevailing issues with mechanical weakness mean that single-component hydrogels are typically used in environments that involve low mechanical stress. To solve this shortcoming, a class of material known as a double-network hydrogel has been designed with both high stiffness and high toughness.¹¹ Conventionally, double-network hydrogels are constructed from two polymeric networks with diametrically opposed physical properties, namely a rigid, highly crosslinked polyelectrolyte network and a sparsely crosslinked neutral polymer network that is highly extensible.¹² Double-network hydrogels have fascinating mechanical properties in that they fracture at high energies (up to 10³ J m²) and have large tensile stress (up to 10 MPa) and high tensile strain (up to 2000%).¹³ A refinement on this initial design principle is the hybrid double-network, in which the charged network is ionically crosslinked, which provides transient crosslinking, allowing fracture and re-healing of the rigid network.¹⁴ A schematic of this design is shown in Fig. 1(a), in which the rigid network polymer is alginate, crosslinked by divalent calcium ions, and the extensible network is polyacrylamide (PAAm), crosslinked by N,N′-methylene bisacrylamide (MBA). Numerous factors, such as the stretchability of the stiffer network,¹⁵ the relative strength of the two networks,¹⁶ and the sacrificial nature of the stiffer network,¹² are thought to contribute to making a hybrid double-network hydrogel characteristically tough. In general, the mechanical properties of a hybrid double-network hydrogel can be affected by other factors, such as entanglement¹⁷ and heterogeneity.¹⁸ Crosslinking between the two networks in the form of hydrogen bonding¹⁴ is also a significant factor since it allows load sharing, which we have recently reported, through shear rheology experiments that identify the precise role of these crosslinks on the initial stage of yielding.¹⁹ However, despite this seemingly extensive investigation, the independent control of stiffness and toughness hasn't been achieved. As such, the two properties tend to be mutually exclusive.²⁰ Therefore, a major stiffness–toughness conflict²¹ has arisen in double-network design, meaning the resistance of the materials to fracture is typically regarded to represent a combination of high strength and large deformability. As a result, the design of a hydrogel with low stiffness but high toughness remains a challenge.

Fig. 1 (a) Schematic of the Ca²⁺-alginate/PAAm hybrid double-network hydrogel with a 3D rendering of a glass slide and silicon spacer mould designed to synthesise the hydrogel samples; (b) setup for tensile testing of one of the Ca²⁺-alginate/PAAm hybrid double-network hydrogel samples; (c) stress vs. strain curves for two different hybrid double-network hydrogel samples – sample (1) (in red) has composition: C_MBA = 1.75 mM, C_APS = 1.1 mM, C_CaCO3 = 0.045 M, C_AAm = 1.8 M, and C_Alg = 1.6% w/v and sample (2) (in blue) has composition: C_MBA = 0.25 mM, C_APS = 0.14 mM, C_CaCO3 = 0.028 M, C_AAm = 1.8 M, and C_Alg = 1.6% w/v; (d) representation of input parameters and output parameters for the DoE in this study.

How then, can we address this challenge in the context of hybrid double-network hydrogels? The main challenge is posed by an evident lack of a direct, clear correlation between various mechanical properties of this hydrogel and its combined constituent networks parameters, such as crosslinking density and monomeric concentration.

There is often a complex correlation between properties and networks parameters, making it difficult to grasp the relation holistically. The mechanical properties of the hydrogel can be influenced by more than one network parameter. To illustrate this point, one may consider a relatively well-studied Ca²⁺-alginate/polyacrylamide (PAAm) hybrid double-network hydrogel.¹⁴ Here, the PAAm network is polymerized through free radical polymerization such that the concentration of the crosslinker affects the effective length of the polymer chain between the crosslinks. By contrast, the concentration of initiator affects the size of the resultant polymer chain itself. Since the fracture energy of the hydrogel depends on the polymer chain length, the relation of the FE of the hybrid double-network hydrogel can be surmised to depend both on the initiator concentration and the crosslink concentration. Another example is illustrated in recent work from our group in which we show that changes in the elastic modulus depend more on the alginate crosslink concentration than on the PAAm crosslink concentration.¹⁹ Furthermore, it has been observed that a compositional parameter of a hybrid double-network hydrogel can affect different mechanical properties in different ways. In Ca²⁺-alginate/polyacrylamide (PAAm) hybrid double-network hydrogel, increasing the crosslinking density of either of the networks has been reported to increase the overall elastic stiffness of the hydrogel; however, the toughness characterized by fracture energy of the hydrogel increases at first but decreases when the crosslinking density is high.¹⁴

An additional challenge is presented by the mode of investigation, which typically entails changing one factor at a time. When the response is complex, with few factors interacting, this traditional method falls short as the parameter space is not explored enough in depth. However, extensively exploring parameter space could involve a larger number of experiments and thus infeasible experimental timescales, particularly in industrial or biomedical applications. At the same time, if the parameter space is not explored systematically, the analysis of data result in misleading findings. This seemingly impassable problem is often sidestepped in industrial research, through a statistical method called design of experiments (DoE). Using DoE, it is possible to systematically explore a parameter space more extensively with fewer experiments using careful statistical analysis that investigates the relationship between control factors and the responses mathematically. The independence of parameters in a Ca²⁺-alginate/PAAm double-network hydrogel also makes it a good candidate for DoE because the hydrogel is a combination of two independent polymer networks with distinct physical properties (stiff alginate and flexible PAAm) and distinct crosslinking chemistries (covalently crosslinked PAAm and ionically crosslinked alginate). It is therefore entirely expected that compositional factors relating to crosslinking and composition would lead to changes in the mechanical properties. However, it is unclear how interactions between these compositional factors would influence these mechanical properties. The conventionally used one-factor-at-a-time approach cannot test these interactions, while the DoE approach presents itself as an ideal alternative candidate as it analyses effects that not only arise from independent single parameters but also from interactions between the parameters. Recently, DoE has been demonstrated to be capable of correlating the parameters related to the synthesis process of Ca²⁺-alginate/PAAm hybrid double-network hydrogel and the mechanical properties.²² In this work, we now aim to gather an indicative understanding how this approach can be used to address the stiffness–toughness conflict. We, therefore, focus on two parameters: Young's modulus (E) and fracture energy (FE). The use of DoE enables a significant reduction in the number of experiments required while retaining maximum certainty in the effects of the experimental parameters. Additionally, using the statistical models, DoE finds the best set of parameters to achieve optimum properties within the range of experiments.

2 Experimental section

2.1 Materials

Acrylamide (AAm), alginate, N,N′-methylene bisacrylamide (MBA), ammonium persulfate (APS), N,N,N′,N′-tetramethylethylenediamine (TEMED), calcium carbonate (CaCO₃) and glucono-δ-lactone (GdL) were purchased from Sigma-Aldrich. All chemicals were used as received.

2.2 Synthesis of the hybrid double-network hydrogel

The hydrogel synthesis was motivated by the protocol reported previously.¹⁹ AAm and alginate stock solutions were prepared with concentrations of 5 M and 5.4% w/v, respectively. Solutions of CaCO₃ with concentrations specified in Table 1 were prepared with MilliQ water. The volume of water was varied to maintain a constant total sample volume of 5 ml. The solution was sonicated at 25 °C for 20 minutes to reduce the particle size of the CaCO₃ so that it would react more readily. 0.04 M APS stock (initiator) solution was prepared by mixing 45.6 mg of APS with 5 ml of deionised water and was vortexed at 2000 rpm in the Fisherbrand ZX4 Vortex Mixer. Once sonication was complete, 2.7 mM of TEMED (cross-linking accelerator) and AAm, alginate, APS, and MBA (cross-linker), at molar concentrations also specified in Table 1, was added to the CaCO₃ solution and vortexed at 2000 rpm. Nitrogen gas was then purged through the solution for 10 seconds to prevent exposure to oxygen. Following this, GdL was added from a freshly prepared 1 M stock solution. The hydrolysis of GdL is known to proceed gradually through a series of intermediate steps, which liberates calcium ions from the insoluble CaCO₃ salt over time. Finally, the reaction mixture was vortexed at 2000 rpm for 5 seconds. The molar concentration of GdL was kept at a 2 : 1 constant ratio with the molar concentration of CaCO₃ to ensure full dissociation of CaCO₃. Each sample was formulated to produce 5 ml of Ca²⁺-alginate/PAAm hybrid double-network hydrogel, which translated to two samples.

Table 1 13 unique formulations of hybrid double-network hydrogel produced by the Minitab software for initial screening and corresponding mechanical properties. Additional replicates were synthesized and tested for each of these formulations to include in the DoE's statistical analysis

Run unit	MBA, mM	AAm, M	Alginate, % w/v	APS, mM	CaCO3, M	YM, kPa	FE, J m^−2
A	0.25	0.89	0.99	0.14	0.045	1.5 ± 0.1	27 ± 2
B	3.0	1.8	0.39	0.14	0.045	39.7 ± 3.0	377 ± 28
C	0.25	1.8	0.39	0.62	0.011	6.4 ± 0.5	1750 ± 130
D	3.0	0.89	1.6	0.62	0.045	55.3 ± 4.1	1020 ± 77
E	1.6	0.89	0.39	0.14	0.011	7.0 ± 0.5	54 ± 4
F	1.6	1.3	0.99	0.62	0.028	45.8 ± 3.4	852 ± 64
G	3.0	1.8	0.99	1.1	0.011	34.8 ± 2.6	512 ± 38
H	0.25	1.3	0.39	1.1	0.045	53.3 ± 4.0	876 ± 66
I	3.0	1.3	1.6	0.14	0.011	6.3 ± 0.5	456 ± 34
J	3.0	0.89	0.39	1.1	0.028	18.5 ± 1.4	388 ± 29
K	1.6	1.8	1.6	1.1	0.045	44.9 ± 3.4	1130 ± 85
L	0.25	1.8	1.6	0.14	0.028	14.2 ± 1.1	8040 ± 600
M	0.25	0.89	1.6	1.1	0.011	6.3 ± 0.5	456 ± 34

---

2.3 Synthesising thin films of the hybrid double-network hydrogel

To produce the thin film samples for tensile testing, a mould was designed (Fig. 1(a)). The mould contained two glass slides with a 700-micron thick silicon spacer in between them. The silicon spacer was precision cut using the OMTECH 80 W CO₂ laser engraving machine and cutter. To prepare the mould, the silicon spacer was placed on one of the glass slides, and 1.6 ml of the final reaction mixture was pipetted into the mould created by the singular glass slide and silicon. This mould was then placed in the vacuum chamber for one minute to remove air bubbles. The second glass slide was then placed on top and held in place by bulldog clips. The specimen was placed upright in the binder oven at 50 °C for 3 hours to allow for complete polymerization of the AAm chains.

2.4 Mechanical testing of the thin films

Mechanical specimens were tested using the LRK-500 Plus universal tensile testing machine (Fig. 1(b)). A 50 N load cell was used at a rate of 10 mm min⁻¹ per the ASTM F2150 guidelines.²³ Grade 60 sandpaper was used as end tabs to prevent slippage during the experiment and provide a barrier between the hand-tightened grips and the sample, minimizing damage. Every specimen was formed to 70 × 15 × 0.7 mm dimensions and measured with a micrometer to ensure consistency between experiments. Each sample had a tolerance of ±10% of the desired dimensions. The thickness was measured at three distinct points, with the median value being taken in accordance with ASTM D412.²⁴ The LRK-500 Plus recorded data on the resistive force exerted by the specimen and the extension. From this, the nominal tensile stress (Pa), σ_t (= force/cross-sectional area), was calculated. The tensile strain, ε_t, was calculated as h/h₀, where h denotes the extension and h₀ denotes the original gauge length. h₀ was calculated by subtracting the length of the end tabs from the total length of the gel and was kept constant at 0.04 m. The specimens were tested until fracture. To characterise the stiffness, the E of the specimens was calculated in the initial linear elastic region of each curve, at 20% strain per ASTM D638²⁵ as E = σ_t/ε_t while FE was calculated by integrating the area under each stress–strain curve using the trapezium method, where n denotes the number of subintervals, set to be equal to the number of data points, specifically 1000 (Fig. 1(c)). It produced the FE in J m⁻³, which was converted into standard units of J m⁻² by multiplying by h₀.

2.5 Design of experiments (DoE)

The DoE methodology builds a predictive and phenomenological model to relate input variables (factors) to output variables (response variables) through a process of sequential optimization.²⁶

For this work, E and FE were identified as the response variables for the hybrid double-network hydrogel as they characterize the stiffness and toughness of the material, respectively. Compositional variables were selected as the factors, specifically the concentrations of the alginate and AAm and their crosslinkers (MBA and CaCO₃) and the initiator for AAm polymerization (APS) (Fig. 1(d)).

The Minitab 21 software guided the DoE process. The DoE framework required limits for each factor to be established. The molar concentration of each factor was varied individually to extremes while adjusting the water content to maintain a constant sample volume of 5 ml. The initial experimental limits were defined as broadly as feasibly possible to maximise the parameter space explored and ensure that the resulting model would apply to future work. As such, limits were defined solely by practical handling constraints. Specifically, the lower limit was defined by samples that either did not appear to form gels or formed gels that were too fragile to insert into the load cell. The upper limit was determined as being reached when gelation happened too rapidly, preventing the sample from being pipetted into the mould.

Experiments were devised using fractional factorial design principles, in which a carefully chosen subset of the possible combination of factors are selected. This approach is designed to maximise the potential to reveal the most important features of the experimental problem while minimising the number of experiments.²⁷ We adopted a specific approach known as a definitive screening design (DSD) that has been shown to minimise the possibility of confounding factors.²⁸ An additional advantage of the DSD is that it provides accurate insight into two-way interactions in which two factors simultaneously and cooperatively influence a response variable. The determined limits (Table 2) were input, and this produced a fractional factorial design, resulting in 13 hydrogel formulations, as seen in Table 1. Each formulation had two replicates, each with two repeats. Here, a replicate refers to conducting an entire experiment multiple times. A repeat refers to performing the same experimental procedure multiple times on the same sample. Each formulation was synthesized and tensile tested. From the tensile data, values for E and FE for each sample were derived and input into the Minitab software.

Table 2 Upper and lower limits of concentrations for factors used in DoE

Limit unit	C MBA, mM	C AAm, M	C Alg, % w/v	C APS, mM	C CaCO3, mM
High	3.0	1.8	1.6	1.1	0.045
Low	0.25	0.89	0.39	0.14	0.011

---

Once acceptable screening results had been established, custom three-level factorial designs were conducted on both E and FE, using only experimental data of significant factors identified by the screening as inputs. A 3-level factorial design investigates each factor at three levels: a minimum, a maximum (equivalent to the determined limits), and a midpoint between the two levels. This design was chosen as it efficiently assesses the main interaction effects in DoE. A statistical model was produced for each of the responses. Instead of the coefficient of determination, R², we used outputs of R²_adj, and R²_pred as follows:

1. R²_adj adjusts R² for the number of predictors penalizing the addition of unnecessary variables and is defined as 1 − ((1 − R²)(n − 1)/(n − p − 1)) where n = number of datapoints and p = number of predictors or independent variables.²⁹ Using R²_adj thus avoids the effects of overfitting, providing a more accurate measure of the models goodness of fit.

2. R²_pred is calculated with a formula that is equivalent to systematically leaving each observation from the data set and determining how well the regression model predicts the removed data point. This strategy provides a more reliable measure of model accuracy when the sample size is small by eliminating the effects of overfitting.³⁰

3 Results and discussion

3.1 Selection of factors

The main aim of this work is to understand the conflict between two response variables: toughness and stiffness. We therefore begin our study by selecting the factors that are known to influence toughness and stiffness individually. Fundamentally, hydrogel mechanical properties are influenced by the polymer concentration and the crosslink concentration. For example, in single-network PAAm hydrogels, increasing the cross-linking density by adjusting the cross-linker concentration leads to an increase in the elastic modulus³¹ and a reduction in the loss tangent, reflecting a stiffer, more elastic hydrogel.³²

For the Ca²⁺-alginate/PAAm hybrid double-network hydrogel that forms the basis of this work, the highest FE is reported to be achieved at the intermediate concentrations of crosslinkers.¹⁴ The reason for this can be surmised from the classic Lake–Thomas theory of fracture mechanics.³³ According to the theory, the ideal fracture energy (G_c) of a crosslinked polymer network depends on the area density of load-bearing polymer chains on fracture surfaces (ρ), the average number of monomer units between the crosslinks, and the bond dissociation energy as G_c = ρN_eU. The highest FE is achieved at an intermediate crosslinker concentration due to its two contrasting effects on G_c. An increase in the crosslinker concentration decreases N_e. However, if the crosslinker concentration is very low, the number of load-bearing chains and, in turn, ρ also remain low, and hydrogels can lack structural integrity, and failure might occur due to chain pullout rather than due to bond scission.³⁴ As a result, FE starts low at a low crosslinker concentration as low ρ dominates over the high N_e. FE increases till ρ increases, dominated by the effect of crosslinker concentration until an optimal value of FE is reached. After this, the effect of the reduction in N_e plays a dominant role in decreasing FE with increasing crosslinker concentration. Independently, the stiffness of such networks is also known to increase with an increase in crosslinker concentration in either of the two networks.^35,36 The crosslinking density of the alginate network is determined solely by the Ca²⁺ concentration, however, the crosslinking density of the PAAm polymer network is also dependent on the concentration of the initiator that drives the free radical polymerization of PAAm monomers. For example, a low initiator concentration enables extended chain propagation, as the initiator radicalizes only a few monomeric chains.³⁷ This leads to the formation of long-chain polymers, which tend to be more flexible, resulting in a lower overall elastic modulus. Based on this, we limit our factor space to five compositional parameters. For the PAAm network, we include monomer concentration of acrylamide (C_AAm), initial crosslinker concentration (C_MBA), and initiator concentration (C_APS) as factors for investigation in this study, as they are likely to affect the PAAm network's mechanical properties. Similarly, for the alginate network, we include the polymer concentration of alginate (C_Alg) and crosslinking density (C_CaCO3).

3.2 Screening for significant factors

Having defined our factors, we now seek to identify those factors that significantly affect the response variables of the films, by using the DSD framework. This involves analyzing data from 13 sets of formulations with 2 replicates, as shown in Table 1. This screening framework provides information on both the linear interactions that relate one factor and one response variable and two-way interactions in which the response variable depends on simultaneous changes in two factors. Two-way interactions are believed to be of central importance in predicting the mechanical response of hybrid double network hydrogels due to the known co-dependence of the network stiffness on the concentrations of the different crosslinks^35,36 and because the FE and the strain at fracture have been indicated to have a dependence on two-way interactions.²²

For all 26 independently synthesized samples, the value of E ranges from 1.5 to 55.3 kPa, while the FE ranges from 27.2 to 8040 J m⁻². The values of the FE are comparable to previously reported values for Ca²⁺-alginate/PAAm hybrid double-network hydrogels.^14,22 We note that the values of E are lower than those reported in previous studies,¹⁴ which can be attributed to the lower thickness of samples used in our experiments. To reveal which of the experimental factors are significant in determining our response variables, we use a statistical approach known as analysis of variance (ANOVA). ANOVA is a statistical test that examines the correlations between factors and response variables. The ANOVA analysis provides two outputs of significance:

1. The p-value, which is the probability of measuring the given response variable if there is no correlation between the response variable and the factor. A smaller p-value, therefore, indicates that the response variable would be unlikely to be measured unless a correlation were present.

2. The standardized effect, a T-value, measures the influence of a factor in relation to variation in the sample data. A larger standardised effect indicates a greater probability that a factor affects the response variable that exceeds the sample variance. A standardised effect of 0 indicates no statistical significance. We define significance between factors as occurring when a 95% confidence interval is exceeded.

The results of the ANOVA approach are shown in bar chart form in Fig. 2, with factors exceeding the 95% confidence interval shown in bold. The individual term of a factor A is represented in the form of A to represent C_A, and the interaction terms between factors are denoted in the form A*B to represent C_A × C_B, where A and B denote the factors. To summarise, we find that C_MBA, C_APS, and C_CaCO3 are significant factors in determining the response variable E. C_CaCO3 and C_MBA control the crosslinking densities in alginate and PAAm networks respectively while C_APS controls the degree of polymerization of AAm monomers. Meanwhile, the FE is significantly influenced by C_APS, C_MBA, C_AAm, and C_Alg. This indicates that the hydrogel failure is driven predominantly by the covalent PAAm network as characterized by its network parameters at synthesis: concentration of AAm monomers (C_AAm) and crosslinking (C_MBA) and the degree of polymerization (C_APS) of the network, in addition to C_Alg. This agrees well with the prevailing hypothesis on the fracture mechanism of such gels, which posits that the ionically crosslinked alginate provides a mechanism of energy dissipation¹⁹ at low and intermediate strains but that the eventual failure of the hydrogel is determined by the fracture of the covalently crosslinked acrylamide network.¹⁴

Fig. 2 Analysis of variance (ANOVA) for (a) Young's modulus (E) and (b) fracture energy (FE) presented through p-values and standardized effect for different factors and their interactions estimated from DSD.

Remarkably, within the concentration limits used, no two-way interactions have any significant effect on E. By direct contrast, for FE, the interactions of MBA with both AAm and alginate are significant. In other words, the influence of MBA-mediated covalent crosslinking on the FE is enhanced by an increase in the concentration of either polymer network. A similar interaction between MBA and AAm has also been reported previously for Ca²⁺-alginate/PAAm hybrid double-network hydrogels.¹⁴

3.3 Secondary screening DoE – prediction of Young's modulus (E)

From the initial screening, a number of factors were eliminated from the analysis since they were shown to have a negligible influence on the response variables. We now utilise the remaining significant factors and use multivariate regression to construct a preliminary model to predict the first of our response variables, E. This preliminary model was characterised through goodness of fit, yielding an adjusted goodness of fit R²_adj of 65.9% and a predicted goodness of fit R²_pred of 42.6%. Standard industrial protocols define an acceptable upper limit of acceptability of 65.0%, therefore, it is necessary to refine this model.

To improve the model, we designed a secondary screening stage through a series of further experimental tests that systematically vary the remaining significant factors, C_CaCO3, C_MBA, and C_APS while keeping the insignificant factors, C_AAm and C_Alg, constant at 1.8 M and 1.6% w/v respectively. The compositions used in these further tests are shown in Table 3. Upon preparing these compositions, the lower limit of C_MBA was increased to 0.5 mM from 0.25 mM as it was observed that this higher concentration was required to ensure sample integrity.

Table 3 Details of the 13 unique formulations of hybrid double-network hydrogel produced by for secondary DoE and corresponding mechanical properties. Additional replicates were synthesized and tested for each of these formulations to include in the DoE's statistical analysis

Run unit	MBA, mM	AAm, M	Alginate, % w/v	APS, mM	CaCO3, M	YM, kPa
A′	3.0	1.8	1.6	0.14	0.045	57.1 ± 4.3
B′	0.5	1.8	1.6	1.1	0.011	14.4 ± 1.1
C′	3.0	1.8	1.6	0.62	0.045	76.5 ± 5.7
D′	3.0	1.8	1.6	1.1	0.028	49.3 ± 5.7
E′	3.0	1.8	1.6	1.1	0.011	30.5 ± 2.3
F′	1.8	1.8	1.6	1.1	0.045	68.7 ± 5.2
G′	0.5	1.8	1.6	0.62	0.011	11.0 ± 0.8
H′	0.5	1.8	1.6	0.14	0.045	33.0 ± 2.5
I′	1.8	1.8	1.6	0.62	0.028	50.1 ± 3.8
J′	0.5	1.8	1.6	0.14	0.028	32.7 ± 2.5
K′	0.5	1.8	1.6	1.1	0.045	34.7 ± 2.6
L′	1.8	1.8	1.6	0.14	0.011	14.1 ± 1.1
M′	3.0	1.8	1.6	0.14	0.011	19.1 ± 1.4

---

Fig. 3(a) shows the experimental design of the secondary stage, and Table 3 shows the E and FE values. A significant improvement in R²_adj (= 88.6%) and R²_pred (= 80.4%) is observed. The significance of the regression coefficients for each of the values is shown in Table 4. Among the three factors that were retained from the initial screen, the value of C_APS was found to be non-significant and was excluded from the model. Therefore, E was found to be dependent on the two remaining factors, C_CaCO3 and C_MBA. Furthermore, and in contrast to the preliminary screening, the two-way interaction between MBA and CaCO₃ was also found to be significant, which is clearly shown by the curved contours in Fig. 3(b). This is an interesting outcome of the model since it suggests that the crosslinking of the two networks does not have a simple additive effect on E. In other words, an increase in crosslinking of either of the two constituent networks increases the contribution of the crosslinking of the other network to the network stiffness. One possible explanation for this is that the two polymers are directly coupled through hydrogen bonding, leading to load-sharing between the two networks.¹⁴ Recent rheological experiments from our group have confirmed that the absence of hydrogen bonding leads to quantifiably different viscoelastic responses.¹⁹

Fig. 3 (a) A schematic representation of the fractional factorial design approach used for secondary screening; (b) contour plot showing the effect of MBA and CaCO₃ on E; (c) contour plots of the predicted values of FE, derived from secondary screening for different concentrations of MBA and AAm.

Table 4 Significance of the regression coefficients fitted for E utilizing secondary screening. The concentrations in the equation have the units corresponding to the units described in Table 3. The p-values marked ‘n.s.’ represent the terms that were found not significant

Unit	a 0	a MBA	a APS	a CaCO3	a MBA,APS	a MBA,CaCO3	a APS,CaCO3
p-Value	<10^−4	7 × 10^−4	n.s.	<10^−4	n.s.	0.018	n.s.

---

Combining all regression coefficients allows us to devise the following equation for E:

E = 6314 − (507 × C_MBA) + (532 345 × C_CaCO3) + (391 638 × C_MBA × C_CaCO3) (1)

Using the model, we re-examine and assure that increasing C_MBA limit from 0.25 mM to 0.5 mM does not inadvertently affect other factors by predicting the values of E for 2 sets of concentrations (L in Table 1 and J′ in Table 3) where only C_MBA is increased (see S1 for further details, ESI†).

3.4 Secondary screening DoE – prediction of the fracture energy (FE)

To build the model for describing FE, we perform fractional factorial analysis on initial formulations (Table 1) using only the significant factors as noted in Fig. 2(b). Although the developed model exhibited high R²_adj (= 87.2%) and R²_pred (= 80.40%), determining 7 parameters from 13 unique data points likely meant there was a high scope for overfitting. As a result, we drop the interaction term between MBA and alginate, i.e., C_MBA × C_Alg, and the term for the degree of polymerization of AAm chains, i.e., C_APS, from the model as they had the highest p-values amongst the significant terms. The 5-parameter model is observed to have increased resistance to overfitting (Fig. S2, ESI†). The model also exhibits improved R²_adj = 88.9% and R²_pred = 83.2% and is represented by the following equation:

FE = −89 105 + (31 310 × C_MBA) + (75 771 × C_AAm) + (20 831 × C_Alg) − (26 459 × C_MBA × C_AAm) (2)

We observe that the constituents related to the second network, i.e., AAm and MBA, dominate the contribution to the FE. Compared to E, the role of MBA in determining the FE is more convoluted. In general, although MBA is negatively correlated with the FE, i.e., FE increases as C_MBA decreases, the regression coefficient is positive (eqn (2)). This shows that MBA interacts strongly with AAm, as shown by the curved contours in Fig. 3(c). These interactions have inverse effects on the FE such that FE increases with increasing C_AAm and decreasing C_MBA. This can be seen as we track FE's evolution along the diagonal starting from the highest C_MBA and the lowest monomeric concentration to the lowest C_MBA and the highest monomeric concentration in Fig. 3(c).

The inverse effect of MBA and AAm effect was also observed between certain weight percentage ratios of MBA/AAm previously.³⁸ This can be explained by considering the movement of polymer chains. A denser network, resulting from higher cross-linking density, leads to shorter average molecular lengths between cross-link points. Consequently, this restricts the movement of polymer chains. When the hybrid double-network hydrogel is stretched, the network with shorter chains fractures because it cannot effectively dissipate the energy stored between crack tips. The network with the longer chains then provides a platform of decreased resistance for crack propagation. So, as C_MBA continues to increase, the influence of shorter chain length outweighs the increased density of cross-linking, and so FE decreases.³⁹ Applying this same logic; the increase in chain length caused by increased C_AAm has the inverse effect, resulting in increased network elasticity and thus a larger FE (Fig. S3, ESI†). Although the MBA and alginate interaction term is eliminated from the model here, owing to the concern of overfitting, similar inverse effects of MBA and alginate have been observed for this class of hydrogel.²²

3.5 Predictive power of the models

To test the predictive power of our models independently, we construct hydrogels from 9 entirely new formulations (Table S5, ESI†) and probe whether the resulting response variables are well-approximated by the model summarized in eqn (1) and (2). The predicted values for E agree with the measured values within a tolerance of 30% barring 2 outliers (Fig. 4(a)) while for the predicted values of FE, they agree within the tolerance 40% with 2 outliers (Fig. 4(b)). We surmise that tolerance values highlight potential errors caused by differences in the kinetics of hydrogel formation for different formulations. These differences can induce heterogenities within hydrogel samples that are not accounted for by the statistical models. Additionally, an inadequate representation and an uneven spread of E and FE values measured in the data used for building the model themselves, as presented in Fig. 4(c), could be contributing adversely to the tolerance values of the test formulations. One example of such nature can be inferred from Table 1 where the maximum E value measured is 55.3 ± 4.1 kPa (Run D) which is much less than the highest E measured among the test formulations in Fig. 4(a) which is nearly 100 kPa. Therefore, we emphasize that these models (eqn (1) and (2)) represent the capability of a highly efficient design of experiments that is necessarily restricted to a small number of experiments. Thus, the predictive accuracy and tolerance values could likely be improved by switching to a more detailed analysis through the full factorial design, where all possible combinations of levels are tested, which could lead to more evenly spread responses across the range.

Fig. 4 Final output of the model of FE (a), and E (b) in relation to measured values for new independent formulations; (c) presents the distribution of E and FE values in the initial DoE data presented in Table 1 and the predictions for 3⁵ = 243 set of concentrations corresponding to all possible combinations for 3 levels of the 5 factors.

With the models, we also optimized and found the optimum set of concentrations within the range of the factors to achieve minimum and maximum E and maximum FE (Fig. S4, ESI†). As C_MBA and C_CaCO3 positively correlate with E individually, the minimum and maximum E values coincide with the lowest and highest concentrations of both constituents, respectively. E_max (∼81 kPa) is about 5.5 times higher than E_min. FE decreases with increasing C_MBA, FE_max (∼1856 J m⁻²) corresponds to the lowest C_MBA. Furthermore, since FE does not depend on C_CaCO3, we can formulate an optimum set of concentrations (C_MBA = 0.5 mM, C_AAm = 1.78 M, C_Alg = 1.6% w/v, C_APS = 0.1375 mM, C_CaCO3 = 0.0112 M) such that E_min and FE_max could be achieved, demonstrating that stiffness and toughness can be tuned independently and simultaneously.

4 Conclusions

This study helps elucidate the significant factors that define network properties and various interactions among them, which together play a major role in determining the E and FE. Importantly, we also establishes that E and FE can be tuned independently within the limits of the concentrations used. Specifically, C_CaCO3 and C_MBA, which define the crosslinking densities of their respective polymer networks, are shown to influence the E, in agreement with the literature (Fig. 2(a) and Table 4). The FE, representing the strength of the hybrid double-network hydrogel, is found to be heavily dependent on the PAAm network (Fig. 2(b) and Table 5). The interaction of C_MBA with C_AAm in determining the FE, and the cross-network interaction of C_MBA and C_CaCO3 in determining the E are observed to be significant. Since the DoE only reveals the correlations between the parameters and the properties and not the causations, the fundamental origins of interaction terms, particularly those between constituents of different networks, require further investigation. The cross-network interaction in determining the E could be hypothesized to arise from diffusion-controlled processes during hydrogel formation. Additionally, although the study has succeeded in elucidating the constituent-property relationships, we believe that it will benefit from more data-intensive DoE approaches such as full factorial design. It also mitigates the potential risk of overfitting and assists with a more robust training dataset of the E and FE values, overcoming the issue of uneven spread, as presented in Fig. 4(c).

Table 5 Significance of the regression coefficients fitted for FE utilizing fractional factorial design. The concentrations in the equation have the units corresponding to the units described in Table 1. The p-values marked ‘r.t.c.o.’ represent the terms that were removed to counter overfitting while those marked ‘n.s.’ represent the terms that were found not significant

Unit	b 0	b MBA	b AAm	b Alg	b APS	b MBA,AAm	b MBA,Alg	b MBA,APS	b Alg,APS
p-Value	<10^−4	0.0716	<10^−4	<10^−4	r.t.c.o.	<10^−4	r.t.c.o.	n.s.	n.s.

---

The study also provides a general approach to solve complex formulation problems. To our knowledge, efforts to increase the toughness of hybrid double-network hydrogels have always resulted in stiffer hydrogels. The study challenges this misconception as it enables the exploration of a broad parameter space, using a minimal set of experiments and achieving optimal mechanical properties through DoE. To illustrate this, Fig. 4(c) plots the values of FE vs. E for the original set of experiments shown Table 1 in relation to the output data of the model for all five factors at three-levels each which amounts to 3⁵ = 243 distinct formulations, far more than would be feasible using standard one-factor-at-a-time approaches.

Additionally, the study highlights that the mechanical properties of bulk hydrogels are retained in sub-millimeter thin films with minimal loss, which could open avenues for the use of hybrid double-network hydrogels as films and substrates in practical applications. Although the DoE approach does not reveal causation, we believe that the approach could be fundamentally useful to gain firsthand knowledge of any new multicomponent material and reveal any prevalent complex relationships between parameters and properties of interest.

Author contributions

VK, NS, and AA-R conceptualised the study. VK, LC, and EL carried out the experimental work and analysed the data. SK contributed to the DoE protocol design. VK and LC prepared the original draft. NS, AO, and AA-R supervised the work. NS and AA-R reviewed and edited the article.

Data availability

Raw tensile testing data and minitab output data for initial screening and mechanical properties data are available at the University of Liverpool Repository at https://doi.org/10.17638/datacat.liverpool.ac.uk/2972.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

VK is supported by an Engineering and Physical Sciences Research Council (EPSRC) PhD studentship through the Soft Matter for Formulation and Industrial Innovation Centre for Doctoral Training (SOFI2 CDT) (EP/L015536/1), co-funded by Reckitt, UK. VK thanks the SOFI2 CDT for training on statistics and DoE. The authors thank Akram Karimian and Craig Hansford for their guidance with the equipment and machinery. Finally, AA-R is grateful to the late Professor Stefan Egelhaaf, to whom this special issue is dedicated, for his support and encouragement in the early stages of his academic career.

Notes and references

1. G. J. Lake and A. G. Thomas, Proc. R. Soc. London, Ser. A, 1967, 300, 108–119.
2. S. Naficy, H. R. Brown, J. M. Razal, G. M. Spinks and P. G. Whitten, Aust. J. Chem., 2011, 64, 1007–1025.
3. J. M. Korde and B. Kandasubramanian, Chem. Eng. J., 2020, 379, 122430.
4. M. R. Islam and M. L. Oyen, Exp. Mech., 2021, 61, 939–949.
5. Y. Jian, S. Handschuh-Wang, J. Zhang, W. Lu, X. Zhou and T. Chen, Mater. Horiz., 2021, 8, 351–369.
6. S. Mondal, S. Das and A. K. Nandi, Soft Matter, 2020, 16, 1404–1454.
7. H. L. Lim, Y. Hwang, M. Kar and S. Varghese, Biomater. Sci., 2014, 2, 603–618.
8. F. Ullah, M. B. H. Othman, F. Javed, Z. Ahmad and H. M. Akil, Mater. Sci. Eng., C, 2015, 57, 414–433.
9. P. Lu, D. Ruan, M. Huang, M. Tian, K. Zhu, Z. Gan and Z. Xiao, Signal Transduction Targeted Ther., 2024, 9, 166.
10. M. Suzuki, Images of the Twenty-First Century, Proceedings of the Annual International Engineering in Medicine and Biology Society 1989, pp. 916–917.
11. J. P. Gong, Y. Katsuyama, T. Kurokawa and Y. Osada, Adv. Mater., 2003, 15, 1155–1158.
12. J. P. Gong, Soft Matter, 2010, 6, 2583–2590.
13. Q. Chen, H. Chen, L. Zhu and J. Zheng, J. Mater. Chem. B, 2015, 3, 3654–3676.
14. J.-Y. Sun, X. Zhao, W. R. K. Illeperuma, O. Chaudhuri, K. H. Oh, D. J. Mooney, J. J. Vlassak and Z. Suo, Nature, 2012, 489, 133–136.
15. T. Matsuda, T. Nakajima, Y. Fukuda, W. Hong, T. Sakai, T. Kurokawa, U.-I. Chung and J. P. Gong, Macromolecules, 2016, 49, 1865–1872.
16. K. Fukao, T. Nakajima, T. Nonoyama, T. Kurokawa, T. Kawai and J. P. Gong, Macromolecules, 2020, 53, 1154–1163.
17. M. Huang, H. Furukawa, Y. Tanaka, T. Nakajima, Y. Osada and J. P. Gong, Macromolecules, 2007, 40, 6658–6664.
18. Y. Kawauchi, Y. Tanaka, H. Furukawa, T. Kurokawa, T. Nakajima, Y. Osada and J. P. Gong, J. Phys.: Conf. Ser., 2009, 184, 012016.
19. V. Kopnar, A. O’Connell, N. Shirshova and A. Aufderhorst-Roberts, Mechanistic Origins of Yielding in Hybrid Double Network Hydrogels, 2024, https://arxiv.org/abs/2409.05765.
20. R. O. Ritchie, Nat. Mater., 2011, 10, 817–822.
21. W. Fan, T. Du, A. Droce, L. R. Jensen, R. E. Youngman, X. Ren, L. Gurevich, M. Bauchy, P. Kristensen and B. Xing, et al. , ACS Nano, 2022, 16, 9748–9761.
22. O. Gorke, M. Stuhlmüller, G. E. M. Tovar and A. Southan, Mater. Adv., 2024, 5, 2851–2859.
23. ASTM, Characterization of biomaterial scaffolds used in tissue-engineered medical products, Accessed: 01 April, 2024, https://www.instron.com/en/testing-solutions/astm-standards/polymer-hydrogels-astm-f2150.
24. ASTM, The Definitive Guide to ASTM D412 Tensile Testing of Elastomers, Accessed: 01 April, 2024, https://www.instron.com/en/testing-solutions/astm-standards/astm-d412.
25. ASTM, How to perform a tensile strength test on plastics according to ASTM D638, Accessed: 01 April, 2024, https://www.instron.com/en-gb/testing-solutions/astm-standards/astm-d638.
26. N. M. Bozorg, M. Leclercq, T. Lescot, M. Bazin, N. Gaudreault, A. Dikpati, M.-A. Fortin, A. Droit and N. Bertrand, Biomater. Adv., 2023, 153, 213533.
27. G. Box, S. Hunter and W. Hunter, Statistics for experimenters. Design, innovation, and discovery, 2005, vol. 2, pp. 235–274.
28. B. Jones and C. J. Nachtsheim, J. Quality Technol., 2011, 43, 1–15.
29. J. Miles, R Squared, Adjusted R Squared, 2014 DOI:10.1002/9781118445112.stat06627.
30. Q. Chen and J. Qi, J. Appl. Econ., 2023, 26, 2207326.
31. L. B. Hazeltine, C. S. Simmons, M. R. Salick, X. Lian, M. G. Badur, W. Han, S. M. Delgado, T. Wakatsuki, W. C. Crone, B. L. Pruitt and S. P. Palecek, Int. J. Cell Biol., 2012, 2012, 1–13.
32. D. Calvet, J. Y. Wong and S. Giasson, Macromolecules, 2004, 37, 7762–7771.
33. G. J. Lake and A. G. Thomas, Proc. R. Soc. London, Ser. A, 1967, 300, 108–119.
34. H. R. Brown, Macromolecules, 2007, 40, 3815–3818.
35. M. M. Fitzgerald, K. Bootsma, J. A. Berberich and J. L. Sparks, Biomacromolecules, 2015, 16, 1497–1505.
36. S. Naficy, S. Kawakami, S. Sadegholvaad, M. Wakisaka and G. M. Spinks, J. Appl. Polym. Sci., 2013, 130, 2504–2513.
37. W.-F. Su, Radical Chain Polymerization, Springer Berlin Heidelberg, 2013, pp. 137–183.
38. J.-Y. Sun, X. Zhao, W. Illeperuma, O. Chaudhuri, K. Oh, D. Mooney, J. Vlassak and Z. Suo, Nature, 2012, 133–136.
39. Y. Li, C. Wang, W. Zhang, Y. Yin and Q. Rao, J. Appl. Polym. Sci., 2015, 132 DOI:10.1002/app.41342.

---

Footnote

† Electronic supplementary information (ESI) available: Information S1: verification of absence of inadvertent effects on changing C_MBA limits; Fig. S2: 5-parameter FE model performance summary; Fig. S3: stress–strain curves to demonstrate the effect of increasing MBA and AAm concentration on FE; Fig. S4: optimization Plots for E and FE; Table S5: composition of 9 formulations explored for validation of the models. See DOI: https://doi.org/10.1039/d4sm01470g

---