Gelatin methacryloyl (GelMA) hydrogels have attracted considerable attention because of their exceptional biocompatibility, tunable mechanical properties, and capacity to support cell adhesion and proliferation, making them highly suitable for applications in tissue engineering, drug delivery, and regenerative medicine [1–3]. Numerous strategies have been proposed to modulate the mechanical properties of GelMA hydrogels for their wide application. Traditional chemical modifications, such as adjusting the degree of methacrylation or employing dual crosslinking mechanisms, have enabled precise tuning of the stiffness and elasticity [4–6]. Similarly, hybridization with other polymers, such as hyaluronic acid or chitosan, has been investigated to increase elasticity and tensile strength while preserving biological compatibility [7–9]. Other reinforcement methods, including the incorporation of nanomaterials such as carbon nanotubes, clay nanosheets, and metallic nanoparticles, offer alternative means to strengthen GelMA hydrogels without extensive chemical alteration [10,11]. Although various strategies reflect ongoing efforts to optimize the balance between mechanical robustness and functional performance in GelMA hydrogels, they frequently introduce trade-offs that compromise biocompatibility and degradation behavior.

Recent breakthroughs in physical entanglement engineering offer a promising alternative pathway. Chansoria et al. [12] utilized the synergy between thermal crosslinking and photopolymerization to achieve a wide range of mechanical properties within a fixed GelMA composition. This thermal crosslinking mechanism leverages the intrinsic temperature responsiveness of the polypeptide structure of gelatin, where cooling drives helix formation through restored hydrogen bonding [13,14]. Notably, this approach alters the microstructure of the polymer network without affecting the chemical or biological properties of the hydrogel. On the other hand, the transient nature of these physical interactions introduces time-dependent mechanical behavior—a feature that is absent in covalently crosslinked systems [15–17]. However, the lack of quantitative understanding regarding how the cooling process regulates entanglement limits the predictive design of such systems. This knowledge gap becomes particularly critical when engineering hydrogels for dynamic physiological environments where both static stiffness and time-dependent viscoelasticity must be precisely coordinated [18–21].

Our investigation bridges this gap by establishing a direct correlation between thermal history and mechanical response through multiscale analysis. We first investigated the effects of cooling-induced entanglement on the mechanical behavior of GelMA hydrogels through uniaxial tensile and stress relaxation tests. These experiments demonstrate that extended cooling durations are correlated with enhanced mechanical properties. Microstructural characterization further revealed that cooling induces a more compact polymer network structure, with MD simulations confirming the pivotal role of entanglement in these improvements. On the basis of these findings, we developed a theoretical model to characterize the cooling-induced entanglement mechanisms and predict the mechanical responses of GelMA hydrogels under various cooling durations. This work not only advances the fundamental understanding of the relationships among cooling conditions, microstructure, and mechanical performance but also establishes a framework for the rational design of GelMA hydrogels with optimized properties for biomedical applications.

In this study, dumbbell-shaped GelMA hydrogels (5 wt%, 50 % degree of substitution) with a gauge length of 12 mm, width of 2 mm and thickness of 1 mm were fabricated via an LAP photoinitiator (0.5 wt%). The samples in the precooling groups underwent thermal crosslinking at 4 ºC for 5, 10, 15, or 20 min before UV curing (405 nm, 10 mW·cm^–2, 90 s), whereas those in the control groups were immediately photopolymerized at 25 ºC. All the samples were equilibrated for 1 h postfabrication to relieve residual stresses. The mechanical response of the GelMA hydrogels was evaluated through uniaxial tensile and stress relaxation tests via an Instron 5944 universal testing system equipped with a 10 N load cell. For monotonic tensile characterization, samples were subjected to constant deformation rates of 0.014 and 0.14 s^–1. For stress relaxation analysis, samples were first stretched to target strains (50 %, 100 %, 150 %) at 0.14 s^–1 and then held at constant strain for 20 min while the stress decay was recorded.

As demonstrated in Fig. 1, the GelMA hydrogels exhibited significant cooling time-dependent mechanical reinforcement under monotonic tensile loading. Noncooled control samples (0 min, black curves) displayed baseline mechanical behavior with an ultimate tensile stress of 4 kPa at strains below 200 %. Extending the cooling duration to 15 min amplified this stress by 10-fold (Fig. 1(a)), highlighting the critical role of thermal processing in network reorganization. This cooling-driven mechanical enhancement was accompanied by a progressive increase in the Young’s modulus. While the control samples exhibited linear elasticity until rupture, the cooled hydrogels exhibited significant nonlinear hardening. Notably, longer cooling durations resulted in both higher Young’s moduli and more substantial strain hardening. A strain rate sensitivity analysis comparing Fig. 1(a) and (b) revealed consistent viscoelastic responses across all the cooling conditions. For example, at the higher strain rate of 0.14 s^–1 (Fig. 1(b)), the reinforcement effect intensified further, with the ultimate tensile stress increasing >15-fold from the control (4 kPa) to the 15 min cooled samples (> 60 kPa).

To further explore the viscoelastic behavior of the cooled GelMA hydrogels, Fig. 2 systematically characterizes the stress relaxation dynamics of the GelMA hydrogels under various cooling durations and applied strains. The original experimental data are plotted in Fig. 2(a), (c) and (e), which show that the uncooled control hydrogels exhibited negligible stress relaxation. The stress decay plateaus beyond 1200 s, but prolonged testing was intentionally limited to avoid dehydration artifacts that could compromise measurement reliability under ambient conditions. For comparative purposes, the relaxation stress (Fig. 2(b), (d) and (f)) was normalized relative to the initial stress at the onset of relaxation. At 50 % strain (Fig. 2(b)), the hydrogels cooled for 20 min retained 80 % of the initial stress after 1200 s, whereas the initial stress in the 5 min cooled group was 71 %. This trend persisted at higher strains, with 20 min cooled samples maintaining 76 % (100 % strain, Fig. 2(d)) and 68 % (150 % strain, Fig. 2(f)) residual stress, consistently outperforming shorter cooling durations. Critically, while the positive correlation between the cooling time and residual stress held for all strains, the magnitude of this effect diminished with increasing deformation. At 50 % strain, the residual stress difference between the 20 min and 5 min cooled hydrogels reached 10 %, whereas at 150 % strain, this difference narrowed to 5 %. This strain-dependent convergence suggests that large deformations partially override cooling-induced structural advantages, likely because increased chain mobility overwhelms topological constraints at high strains.

We systematically characterized the microarchitecture of GelMA hydrogels by scanning electron microscopy (SEM) to elucidate the structural basis of the cooling-induced mechanical enhancement. Following standardized preparation protocols, including freeze-drying (DYYB-10, Deyangyibang), gold sputter coating (Supro Instrument), and triplicate imaging (Phenom Pro SEM), we observed distinct cooling-dependent architectural transformations (Fig. 3). The microstructural evolution exhibited clear temporal progression: noncooled controls (0 min) displayed a large pore size network that progressively transformed into a densely interconnected architecture with extended cooling. This structural transformation occurs progressively with increasing cooling duration, demonstrating a clear time dependence of network reorganization. These microstructural modifications provide direct physical explanations for the mechanical behavior of GelMA hydrogels: the reduced pore size corresponds to the increased elastic modulus in tensile tests, whereas the improved connectivity accounts for the superior stress retention during relaxation.

To bridge the research gap in understanding the effect of cooling on the physical crosslinking of GelMA hydrogels, MD simulations were employed to characterize the microstructural evolution as well as the mechanical behavior during deformation, overcoming the limitations of experimental methods in revealing the molecular interactions with and without cooling. To capture the qualitative impact of cooling, we did not simulate different cooling durations explicitly. Instead, a sufficiently long equilibration at 4 ºC was applied to ensure that key cooling-induced structural features were well established within the simulation window. The polymer chain with 6 structural units of GelMA hydrogels (shown in Fig. 4(a)) was initially constructed. The 50 GelMA hydrogel chains and 500 water molecules were packed into the simulation box with the size of 7.74 × 7.74 × 7.74 nm³, α = β = γ = 90º and a density of 1 g·cm^–3, containing 38,500 atoms. The crosslinking reactions were performed via a crosslinking algorithm [22] at temperatures of 277 K and 300 K, respectively. The process for the crosslinking reaction is shown in Fig. 4(b), which involves the new formation of C-C bond and therefore chemical linking between two GelMA hydrogel chains. During the crosslinking process, the reactive C atom on one GelMA chain breaks its double bond and forms a single bond with another reactive C atom on the adjacent chain, and this process continues until the desired degree of crosslinking (50 %) is achieved. The constructed model was geometrically optimized, and the energy was minimized via the conjugate gradient approach for a stable structure. Equilibrium was performed under the NPT ensemble at 300 K and a pressure of 1 atm for 1 ns to achieve an equilibrated state. The condensed-phase optimized molecular potentials for atomistic simulation studies (COMPASS) forcefield [23] was used to describe the intra- and intermolecular interactions within the system, and the Ewald method was used to handle long-range electrostatic interactions. Figure 4(c) shows the chemical and physical crosslinking networks in the cooled and uncooled GelMA hydrogels. Since directly identifying topological entanglements from dynamic simulations remains challenging, we employed the fractional free volume (FFV) as a widely accepted proxy for assessing void space and capturing structural changes associated with entanglements, particularly during molecular-level deformation [24]. To characterize the FFV, the Connolly surface was calculated using a probe radius of 1.5 Å. The calculated initial FFV was lower in the cooled GelMA hydrogel than in the uncooled hydrogel, indicating a denser network architecture due to enhanced physical crosslinking induced by the cooling process

To investigate the mechanical behaviors of these two material systems, uniaxial tensile simulations along the x-axis direction were performed by progressively increasing the applied stress in uniform increments, with a fixed duration of 30 ps at each loading step. Before every tensile loading, a 15 ps equilibration was performed to ensure the accuracy of each measured stress-strain value. The microstructural evolution during the tensile deformation process was captured to characterize the state of the cooling-induced entanglements. Figure 5(a) and (b) show snapshots of the structural evolution of GelMA hydrogels with and without cooling during tensile deformation, revealing how their thermal history governs their mechanical behavior. At a strain of 0.5, the cooled GelMA hydrogel shows tensile behavior characterized by the gradual activation of physical crosslinks under loading, as shown in Fig. 5(a). At a strain of 3, the cooled GelMA hydrogel has a more entangled network, with sliding occurring at some localized regions. This seemingly compact network of GelMA hydrogels with local elongated cavities along the tensile direction is the result of cooling-induced entanglement. These induced physical crosslinks act as temporary reinforcement nodes distributed throughout the matrix, contributing to the enhanced elastic modulus measured experimentally. However, the network of noncooled GelMA hydrogels dominated by permanent chemical crosslinks with limited physical interactions shows different tensile behaviors during loading. At a strain of 0.5, the sliding movement of chains with the generation and nucleation of cavities in the matrix was clearly observed in Fig. 5(b). The cavities in the noncooled GelMA hydrogel continued to grow both longitudinally and laterally under tensile loading through chain sliding movement with increasing strain. At a strain of 3, the noncooled GelMA hydrogel undergoes progressive network degradation characterized by uncontrolled cavity coalescence and lateral expansion. Chain sliding becomes predominant when chemical crosslinks act as fixed points, resulting in an elongated network along the loading direction. The structural evolution observed in these two systems underscore the critical role of cooled-induced entanglement in improving the mechanical properties of GelMA hydrogels. Correspondingly, we calculated the fractional free volume (FFV) during the deformation process. As shown in Fig. 5(c), the FFV increases noticeably with strain, which is consistent with the observed network evolution. Notably, the increase in FFV is significantly smaller in the cooled GelMA hydrogel than in its uncooled counterpart, further suggesting that cooling promotes a denser and more stable network structure under mechanical loading. We further extracted the stress–strain curves for the cooled GelMA and uncooled GelMA systems, as shown in Fig. 5(d). The simulation results exhibit a similar trend to the experimental data, capturing the overall increase in network stiffness. Notably, the difference in stress magnitude between the cooled and uncooled systems is less pronounced in the simulations than in the experiments. This discrepancy likely arises from the idealized nature of the simulated networks: the polymer architectures are constructed via predefined crosslinking protocols, which may underestimate the structural differences induced by cooling in real samples.

To quantitatively understand the impact of cooling-induced chain entanglements on the mechanical properties of GelMA, we developed a constitutive model to analyze the contribution of entanglements to the mechanical response of the hydrogel. Our experimental observations led to three key findings: (1) noncooled gels exhibit purely elastic behavior without strain hardening effects; (2) cooling processing induces polymer chain entanglements that significantly enhance mechanical performance, including the emergence of strain hardening; and (3) the entanglement contribution displays viscoelastic characteristics with obvious stress relaxation behavior. On the basis of these observations, we decompose the true stress into two components:where σ_en and σ_cr represent the stresses contributed by entanglements and chemical crosslinks, respectively.

The elastic response of the chemically crosslinked network under uniaxial deformation is described by a Neo-Hookean modelwhere μ_cr is the shear modulus of the chemically crosslinked network, and λ is the stretch ratio in the loading direction.

The reinforcing effects of chain entanglements on polymer stiffness have been well characterized through various theoretical models, including both elastic and hyperelastic formulations [25–27]. Recent studies have further identified their significant contribution to viscoelastic behavior [28,29]. In this work, we adopt a viscohyperelastic formulation based on the eight-chain modelwhere μ_en is the shear modulus arising from the entanglement; n represents the number of chain segments; I₁ = λ²(t) + 2λ^–1(t) denotes the first invariant of the right Cauchy-Green tensors; β signifies the inverse Langevin function written as ; and ф(t) stands for the normalized density of entanglement as a function of time, which is governed bywhere k_{d_0} is the reference value of the disentanglement rate, and α is a scaling exponent reflecting the entanglement decay. This equation combines a baseline disentanglement rate with a power-law decay to simplify the multiscale relaxation dynamics of the microstructure [30]. The exponential term e^σ_en/σ₀ captures the stress-dependent acceleration of disentanglement, with σ₀ representing the applied stress and a characteristic stress scale [15]. The initial condition is ф(t) = 1, corresponding to a fully entangled state at the onset of loading.

In total, the model contains six parameters: the shear modulus of the chemically crosslinked network μ_cr, the shear modulus of the entangled network μ_en, the number of effective chain segments in the entangled network n, the basic disentanglement rate k_{d_0}, the stress scaling factor of disentanglement rate σ₀, and the decay factor α. These parameters were systematically calibrated against experimental data from Figs. 1 and 2 through the following protocol: μ_cr was directly calculated from the stress-strain relation of non-cooled hydrogels; μ_en and n were optimized to capture the hyperelastic behavior of the GelMA hydrogel at a strain rate of 0.14 s^–1, where viscous effects are minimized; k_{d_0}, σ₀, and α were simultaneously fitted to stress relaxation profiles across multiple strain levels.

Table 1 summarizes the parameter values obtained for the GelMA hydrogels processed for different durations of cooling. The entanglement shear modulus μ_en monotonically increases with prolonged cooling, quantitatively reflecting the progressive formation of topological constraints. Moreover, the basic rate k_{d_0} decreases markedly with increasing cooling time, indicating a slowdown in chain disentanglement. In contrast, both the stress sensitivity factor and the scale factor tend to increase. This inverse relationship suggests that mature entanglements formed during extended cooling possess enhanced kinetic stability. The coupled evolution of these parameters implies a self-reinforcing mechanism: extended cooling not only promotes the formation of more entanglements but also improves their mechanical integrity.

The model with incorporated parameters in Table 1 successfully predicts both the rate-dependent tensile response and stress relaxation behavior of GelMA hydrogels. As shown in Figs. 6 and 7, the numerical solutions demonstrate excellent agreement with the experimental measurements across all the cooling conditions. Three key mechanistic insights emerge from this validation:

(1) Rate dependence origin — The strain rate sensitivity predominantly arises from time-dependent disentanglement processes. Lower deformation rates permit more extensive disentanglement during loading, reducing the instantaneous modulus.

(2) Stress relaxation dynamics — accelerated relaxation at higher strains stems from a stress-activated disentanglement mechanism. The exponential dependence of k_d on σ_en creates positive feedback where increased stress promotes faster entanglement loss.

(3) Cooling time effects — Hydrogels with extended cooling durations exhibit reduced rate sensitivity and slower relaxation. This originates from their elevated k_{d_0} values, which reflect more stable entanglements requiring higher activation energies for dissociation. Remarkably, at large strains (150 %), the relaxation trajectories of long-cooled (20 min) and short-cooled (5 min) hydrogels are similar, suggesting that the stress levels eventually overcome the initial entanglement stability differences.

This study elucidates the critical role of cooling-induced topological entanglements in enhancing the mechanical performance of GelMA hydrogels. Through experimental and computational analyses, we demonstrate that extended cooling durations promote the formation of compact polymer networks with reduced mesh sizes, leading to significant enhancements in tensile strength, elevated Young’s modulus, and superior stress retention during relaxation. The proposed viscohyperelastic constitutive model quantitatively correlates the evolution of physical entanglements with mechanical outcomes, providing a mechanistic understanding of how cooling-induced structural reorganization governs hydrogel viscoelasticity. These findings advance the fundamental understanding of noncovalent engineering in soft matter mechanics, providing a predictive framework to optimize GelMA hydrogels for biomedical applications requiring tunable viscoelasticity, such as cell scaffolding materials [31].

While the current study establishes a clear structure-property relationship under controlled conditions, several limitations remain that warrant further investigation:

(1) This study specifically focused on the effect of cooling duration at a fixed temperature (4 ºC) to isolate the temporal aspect of entanglement formation. We acknowledge that broader thermal histories—including variations in the cooling rate, subzero temperatures, and thermal cycling—may further influence entanglement dynamics and network stability. These factors represent promising directions for future optimization and generalization of the proposed framework.

(2) To quantitatively capture these phenomena, we developed a visco-hyperelastic constitutive model that decomposes the stress contributions from chemical crosslinks and physical entanglements. The model effectively links the evolution of microstructural entanglements to macroscopic mechanical outcomes, enabling accurate prediction of rate-dependent tensile and stress relaxation behavior. However, it assumes additive and decoupled contributions from entanglement and chemical crosslinks. While this simplification is consistent with our experimental design and allows for robust parameter calibration, it may overlook synergistic or competitive interactions between these networks under more complex physiological conditions, such as dynamic loading, hydration, and biochemical remodeling [32,33]. Future model extensions incorporating such couplings would be valuable for improving predictive fidelity in biological environments.

(3) While mechanical enhancement by cooling-induced entanglement is clearly established, the study does not include postcooling biocompatibility or degradation assays. Given the importance of cytocompatibility and degradability for biomedical applications, future investigations should systematically evaluate cell adhesion, proliferation, and enzymatic degradation to ensure that cooling treatments preserve GelMA’s biofunctionality.