A multiscale polymerization framework towards network structure and fracture of double-network hydrogels

Double-network (DN) hydrogels, consisting of two contrasting and interpenetrating polymer networks, are considered as perhaps the toughest soft-wet materials. Current knowledge of DN gels from synthesis methods to toughening mechanisms almost exclusively comes from chemically-linked DN hydrogels by experiments. Molecular modeling and simulations of inhomogeneous DN structure in hydrogels have proved to be extremely challenging. Herein, we developed a new multiscale simulation platform to computationally investigate the early fracture of physically-chemically linked agar/polyacrylamide (agar/PAM) DN hydrogels at a long timescale. A “random walk reactive polymerization” (RWRP) was developed to mimic a radical polymerization process, which enables to construct a physically-chemically linked agar/PAM DN hydrogel from monomers, while conventional and steered MD simulations were conducted to examine the structural-dependent energy dissipation and fracture behaviors at the relax and deformation states. Collective simulation results revealed that energy dissipation of agar/PAM hydrogels was attributed to a combination of the pulling out of agar chains from the DNs, the disruption of massive hydrogen bonds between and within DN structures, and the strong association of water molecules with both networks, thus explaining a different mechanical enhancement of agar/PAM hydrogels. This computational work provided atomic details of network structure, dynamics, solvation, and interactions of a hybrid DN hydrogel, and a different structural-dependent energy dissipation mode and fracture behavior of a hybrid DN hydrogel, which help to design tough hydrogels with new network structures and efficient energy dissipation modes. Additionally, the RWRP algorithm can be generally applied to construct the radical polymerization-produced hydrogels, elastomers, and polymers.


INTRODUCTION
Since the first PAMPS/PAM [poly(2-acrylamido-2-methylpropanesulfonic acid) (PAMPS) as the first network and polyacrylamide (PAM) as the second network] double-network (DN) hydrogel was developed by Gong and co-workers in 2003 1 , DN hydrogels have demonstrated themselves as perhaps the toughest soft-wet materials. In contrast to conventional hydrogels with a singlenetwork (SN) structure that are mechanically weak and brittle (low fracture energies of~10 0-1 J/m 2 and poor elastic moduli of 10 kPa) [2][3][4] , DN hydrogels can readily achieve extremely high mechanical strength (failure tensile stress 1-10 MPa, strain 1000-4000%; failure compressive stress 20-60 MPa, strain 90-95%) and toughness (tearing fracture energy of 1000-4000 J⋅m −2 ) 1,5 , comparable to or even exceed the strength and toughness of some soft-load-bearing tissues and rubbers 6 . Such high mechanical properties of DN hydrogels are mainly stemmed from their unique contrasting but independent network structures and efficient energy dissipation mode. The "sacrificial bond" concept has often been used to interpret the toughening mechanism of DN hydrogels, particularly for chemically-linked DN hydrogels, i.e., upon mechanical deformation, the first, rigid, brittle, and tightly linked network is first fractured into small clusters, which serves as "sacrificial bond" to dissipate energy and protect the second, soft, elastic, loosely linked network from crack propagation. This mechanism well explains toughness enhancement, large hysteresis, and necking phenomena of chemicallylinked DN hydrogels. Based on the "sacrificial bond" design principles, different chemically-linked DN hydrogels were developed, including microgel-reinforced DN gel 7 , void-DN gel 8 , inverse DN gel 9 , jellyfish DN gel 10 , liquid crystalline DN gel 11 , and lamellar bilayers DN gel 12 . The past decade has witnessed the growing interest in these DN hydrogels as driven by their promising applications in diverse fields, including environmental engineering 13,14 , agriculture and food chemistry [15][16][17] , photochemistry 18,19 , wound dressing 20,21 , 3D printing 22,23 , electronic devices 24 , and tissue engineering [25][26][27][28][29][30] .
From a structural viewpoint, the two interpenetrating networks DN hydrogels can be fully chemically-linked, fully physicallylinked, or hybrid chemically-physically-linked. Current knowledge of DN hydrogels, including network structures and toughening mechanisms, mainly comes from fully chemically-linked DN hydrogels. Fundamentally, the fracture of chemically-linked covalent bonds is different from that of physically-linked noncovalent bonds 1,[7][8][9][10] . The chemically-linked DN hydrogels lack selfrecovery and self-healing ability to repair damages and fatigues due to the irreversible break of covalent bonds 31 . The noncovalent bonds (e.g., hydrogen bonds, hydrophobic interactions, π-π interactions, ionic interactions) in physically-linked networks possess intrinsically reversible property to different extents even without any external stimuli. In the past years, a number of physically-linked DN hydrogels have been developed by our and other labs, including agar/PAM 32 , agar/HPAAm 33 , Alginate/ PAAm 34 , PBDT/PAAm 35 , PAM/XG 36 , and CS/PAA 37 . While these DN hydrogels contain either single or double physical (noncovalent) networks, they all exhibit high toughness, strength, extensibility, and self-recoverability to some extents with and without certain external stimuli. These studies demonstrate that the introduction of reversible non-covalent bonds into the DN structure is a general strategy to impart both mechanically strong and self-recovery/self-healing properties in a single material.
While great progress has been made in synthesis methods, design principles, and toughening mechanisms for DN hydrogels 4,5,33 , the research on DN hydrogels is still at a very early stage 4 . Molecular details at different stages of the fracture process are still poorly understood. It still remains challenging to characterize when, where, and how the crack formation and crack propagation occur within and between two networks, how energy dissipate during this crack process, and what similarities and differences are for the cracking details between physically-linked and chemically-linked DN gels and between SN and DN gels. To remedy this issue, computational tools are suitable for providing the atomic-level information about the structure, dynamics, and interaction of polymer networks in DN hydrogels. Notably, only two molecular dynamics (MD) studies were reported on PEO/PAA DN hydrogels by Goddard III and co-workers in 2007 38 and PAMPS/PAAm DN hydrogels by Kubo and co-workers in 2018 28 , which is likely due to the complexity of inhomogeneous DN structures and the diversity of reaction pathways. Using all-atom MD simulations, the stress-strain profiles showed that the tensile stress of PEO/PAA DN hydrogels exhibited~10 times stronger than that of PEO or PAA SN hydrogels alone. However, the short simulation time of~1.2 ns and the pre-constructed DN structure only offer the limited structure-property relationship of this hydrogel. Recently, coarse-grained models were utilized to construct DN hydrogel from monomers and steered MD simulations were further conducted to study the fracture process of PAMPS/PAAm DN hydrogels against stretching. The lack of explicit water molecules in simulations largely neglects the important role of water molecules in determining water-network interactions (particularly hydrogen bonds) on network fracture, and mechanical properties.
Both MD simulations targeted chemically-linked DN gels 5 . From a mechanistic viewpoint, physically-or physically/chemicallylinked DN gels (e.g. agar/PAM) exhibit different mechanical behaviors from chemically-linked DN gels [39][40][41] , including much lower-yielding stress/strain, no flat necking platform, and simultaneous necking 32 , suggesting new energy dissipation modes and toughening mechanisms. To bridge this gap, we developed a multiscale simulation platform to computationally investigate the early fracture of physically-chemically linked agar/PAM DN hydrogels at a long timescale. Specifically, a "random walk reactive polymerization" (RWRP) was developed to mimic a radical polymerization process, which enables to construct a physicallychemically linked agar/PAM DN hydrogel from monomers. The RWRP-constructed agar/PAM DN hydrogels resembled the experimental properties of network structures and dynamics, mechanical stress-strain behavior, and water dynamics to the larger extent. Then, non-equilibrium steered MD simulations were performed to study the structural-dependent mechanical and fractural behaviors of agar/PAM DN hydrogels at the relax and deformation states. Collective simulation results revealed that energy dissipation of agar/PAM hydrogels was attributed to the pulling out of agar chains from the DNs and the disruption of massive hydrogen bonds between and within the DN structures, thus explaining a different mechanical enhancement of agar/PAM hydrogels from chemically-linked DN hydrogels. This computational work provided the atomic details of network structure, dynamics, solvation, and interactions of a hybrid DN hydrogel, as well as a different structural-dependent energy dissipation mode and fracture behavior of a hybrid DN hydrogel, both of which will help to design new tough hydrogels with new network structure and efficient energy dissipation mode. This multiscale polymerization framework can also be generally applied to the radical polymerization-produced hydrogels, elastomers, and polymers with the complicated network structures with and without crosslinkers.

RESULTS AND DISCUSSION
Construction of hybrid agar/PAM DN hydrogel using RWRP Agar/PAM DN hydrogel consists of a hybrid physically-chemicallylinked DN structure, where an agar network as the first physical network linked by hydrogen bonds, while a PAM network as the second chemical network linked by N,N'-methylenebisacrylamide (MBA). To computationally build an agar/PAM DN hydrogel, Fig. 1a shows a step-by-step computational modeling protocol to construct agar and PAM networks sequentially, which mimics a two-step synthesis process of agar/PAM hydrogel in our previous work 32 . Briefly, the first agar network was constructed by selfassembling the agar helixes into helical bundles (Supplementary Movie 1). The resulting agar helical bundle with the typical network morphology was structurally rigid and stabilized by massive hydrogen bonds. Upon the formation of the first agar network, a total of 4389 AM monomers were distributed around the agar network, while 1 mol% of AM monomers were activated and served as initiators. To mimic the free radical polymerization of AM monomers into long PAM chains, we developed a general "random walk reactive polymerization" (RWRP) protocol to construct the second PAM network in the presence of the first agar network, based on a radical reaction mechanism that each PAM chain has only one reactive site (radical) at its end. At the beginning of the RWRP process, 43 reactive sites were randomly created from AM monomers to initiate the polymerization, in line with the initiator concentration (1 mol%) in the experiments. As the RWRP proceeds, one reactive site and one of the surrounding AM monomers are randomly selected to determine whether they could form a bond within a dynamic reaction radius. If successful, a new bond is formed and the corresponding topological connectivity (i.e., bond, angle, and torsion) is created and input into force field parameters. If not, another random selection of reactive site and AM monomer is conducted until all AM monomers are consumed and polymerized into polymer chains. Whenever a new bond is formed, a combination energy minimization and a short MD simulation are performed to further optimize network structures for the next-step polymerization and allow the molecule diffusions in response to the reaction-induced concentration gradients (Supplementary Movie 2). In the way, a random walking reaction between all reactive sites and AM monomers will successively achieve chain-growth polymerization. Figure 1b shows a typical atomic snapshot of agar/PAM DN hydrogel. As shown in Fig. 1c, the RWRP process underwent a twostage chain-growth polymerization, as evidenced by the decrease number of AM monomers. At the early stage of polymerization (1-5 cycles), free MA monomers have the higher accessibility to interact with reactive sites to form low molecular weight (MW) polymers. As the longer polymer chains are formed at the expense of AM monomers, steric effect becomes dominant due to decrease of AM monomers and the large separation distance between free monomers and reactive sites, making radical polymerization less efficient at the later state of polymerization (5-8 cycles). To address this issue, we introduced an acceleration strategy to promote the chain-growth polymerization by adaptively increasing the reaction radius (8-19 cycles) (Fig. 1d). So, the increased reaction radius allowed the activated radicals to be accessible by the remaining AM monomers, all of which were progressively consumed and polymerized into 43 PAM chains, each containing 100 AMs in average. Overall polymerization degree in Fig. 1e was determined to be around~100, consistent with the experimental value. Finally, five MBAs were added and chemically linked to the second PAM network at random crosslinking positions. The resultant agar/PAM DN hydrogel was then optimized, solvated by explicit water molecules, minimized in energy, and equilibrated by 100 ns MD simulations at room temperature of 298 K.
Different from MD simulation that does not involve any bond breakage and formation event, the RWRP allows to form chemical bonds to build long chain polymers from monomers in a progressive and efficient way. This algorithm can not only be generally applicable to the construction of any linear-like polymers from monomers to mimic radical-like chain-growth polymerization, but also bridge a timescale and length-scale gap between simulations and the experiments.
Network structure in agar/PAM DN hydrogel Agar, as a family member of polysaccharide, comprises primarily of alternating D-galactose and 3,6-anhydro-L-galactopyranose linked by α-(1-3) and β-(1-4) glycosidic bonds. Since the gelation of agar network requires agar molecules to adopt α-helical structure and to be associated together to form helical bundles, we generated the triangular agar double helixes as a basic building unit and assembled them into the agar network (Supplementary Movie 1). Figure 2a shows the atomic structure of the first agar and the second PAM network in agar/PAM DN hydrogel, where the first agar network mainly consisted of two structural components: the helical bundles and the soft agar chain connectors ( Supplementary Fig. 1). The agar helical bundles were formed by packing the agar helixes in a triangle way, within which high hydrogen bonding (HB) density offered the strong associate force to assemble three agar helixes together, making it rigid in DN hydrogel. The chain-based root-mean-square-fluctuation (RMSF) values indicate that the first agar network is more stable than the second PAM network (Fig. 2b). The agar chain connectors were relatively soft, facilitating the mesh-like morphology ( Fig. 2a and Supplementary Fig. 1). Differently, PAM chains were more structurally flexible with diverse conformations, thus forming the mesh-like morphology around agar bundles. Hydrophilic crosslinker MBA consists of four hydrogen bonding groups and two double bonds, in which double bonds allow to polymerize with the second network of PAM, while hydrogen bonding groups enable to form physical bonds with both networks. As a result, MBA forms both chemical and physical crosslinks with polymers to create the more branched networks, which help to enhance network strength and dissipate energy efficiently upon deformation ( Supplementary Fig. 2).
We calculated the HB bonding density and the radial paircorrelation functions (g(r)) of HB sites to quantify the HB matrix in the first agar network. The results showed that the HB acceptors are more dominant than the HB donors at the surface of the first agar network. The agar chains contain the massive HB bonding sites, which occupied~50% total agar surface in the DN hydrogel. In comparison, the HB acceptors (~30%) occupied more surface area than the HB donor (~20%) (Fig. 2c). Three g(r) curves of donor-donor, donor-acceptor, and acceptor-acceptor exhibited distinct HB distribution in the agar bundles. Three peaks at~2.5 Å, 4.5,~7.7 Å for HB donor-donor, four peaks at~2.5 Å,~3.7 Å,~4.2, and~5.1 Å for HB donor-acceptor, and five peaks at~2.2 Å,~2.8 Å, 3.6 Å,~4.2 Å, and~4.8 Å for HB acceptor-acceptor, were observed respectively, indicating a sophisticated HB matrix in the first agar network (Supplementary Fig. 3). The combination of these HB peaks offers diverse HB sites not only to associate agar molecules alone in the first network but also to interact with NH 2 and C=O groups in the second PAM network, both of which contribute to the mechanical reinforcement of the agar/PAM DN hydrogel. Figure 2a shows that the computationally polymerized PAM chains in the DN hydrogel displayed more flexible and dynamic structures with different shapes, sizes, and curvatures. The end-to- end distance and the radius of gyration (Rg) of PAM chains in agar/PAM DN hydrogel showed a wide but peak distribution located at~40 and~20 Å, respectively (Supplementary Fig. 4a, b). Considering that AM monomer contains one hydrophilic amide (NH 2 ) and one keto (C=O) groups, individually containing 2 HB donors and 1 HB acceptor, three g(r) curves almost constant HB densities for both donors and acceptors across all pairwise distances ( Supplementary Fig. 3), indicating that PAM chains are evenly and widely distributed in the DN hydrogel. HB donors and acceptors occupied the comparable surface of the second PAM network, indicating that all donors and acceptors contribute equally to form massive and interpenetrating hydrogen bonds within and between the networks (Fig. 2d). In addition to HBs, PAM chain entanglement also contributed to intra-and internetwork interactions for mechanical reinforcement, leading to the inter-PAM chain interaction energy of~−700 kcal/mol (Supplementary Fig. 4c). We also simulated a PAM SN hydrogel in the absence of the agar network. The calculated end-to-end distance, Rg, and inter-chain interaction energy of PAM chains in the PAM SN hydrogel were comparable to those in the hybrid agar/PAM DN hydrogel, suggesting that PAM chains preserve their structural properties in an independent way. However, both SN hydrogels displayed the higher RMSF values than agar/PAM DN hydrogels ( Supplementary Fig. 5), suggesting that the interpenetrating network structures and interactions enhance the overall stability of DN hydrogels.

Network interactions in agar/PAM DN hydrogel
The superior mechanical properties of DN hydrogels mainly come from their unique contrasting, interpenetrating, and independent network and efficient energy dissipation mode, both contributing to network interactions. Visual inspection of MD trajectories showed that the flexible PAM chains penetrated the relatively rigid mesh-like spaces between agar bundles and formed the strong physical interactions with the agar bundles via hydrogen bonds (Fig. 3a). Further structural analysis shows the positional distributions of AM units relative to the agar bundle surfaces (Fig.  3b). While a broad positional distribution of AM around the agar helical bundle was observed across a wide range of separation distances, the highest peak appeared at 3.0-4.0 Å, indicating that PAM chains have close and dominant contacts with agar network. As a result, the strong and favorable interaction between the agar and PAM networks was found to be~5869.9 kcal/mol, a sum of 2636.0 kcal/mol of VDW and~3233.8 kcal/mol of electrostatic interactions (Fig. 3c). At the agar/PAM interface, hydrophilic NH 2 and C=O groups in PAM had more contacts with agar than hydrophobic (CH 2 -CH-) n backbones, confirming the favorable formation of HBs over hydrophobic contacts between two networks. A total of~456 HBs were formed between the two networks (Fig. 3c), in which NH 2 groups of PAMs (~257 HBs) formed almost~1.3 times higher HBs with the agar network than C=O groups (~199 HBs). This was further confirmed by the observation that the dominant HB acceptors at the agar surface tended to attract more HB donors from H atoms in NH 2 groups to form massive HB interactions, thus highlighting the important role of NH 2 groups in promoting network interactions (Fig. 3a).
Our MD simulations have shown that the first agar network is more rigid with the ordered hydrogen bonding matrix, while the second PAM network is more structural flexible with widely distributed hydrogen bonds, which render the two networks to be more accommodate with each other. It appears that HB acceptors in the first agar network and HB donors in the second PAM network contribute more network interactions to stabilize agar/ PAM DN hydrogel. To prove this observation, we eliminated these HB interactions by replacing hydrogen (H) atoms in NH 2 groups in PAM with hydrophobic methyl groups (CH 3 ), and then computationally polymerized the agar/PAM CH3 DN hydrogel using RWRP protocol (Fig. 3b, Supplementary Fig. 6a). Upon modifying the pendant groups of PAM chains, agar/PAM CH3 DN hydrogel displayed large differences in network structure and interactions from agar/PAM DN hydrogel, where hydrophobic (N(CH3) 2 ) groups were more distanced from the agar network (Fig. 3b). Consequently, the nonbonded interaction energy between the agar network and the PAM CH3 network was decreased by~40% (from~−5869.9 kcal/mol to~−3655.5 kcal/mol), and the hydrogen bonds between networks were reduced by~75% (from~456 to~114) (Fig. 3c). Clearly, the decrease of network interactions mainly stems from the loss of hydrogel bonds. From a structural viewpoint, PAM CH3 polymer chains were more extend than PAM chains, as evidenced by the lager end-end distances and Rg values ( Supplementary Fig. 6b, c). The chain-chain entanglement is also substantially weaker, with lower interaction energies of~500 kcal/ mol ( Supplementary Fig. 6d). In parallel, we further experimentally synthesized the agar/PAM and agar/PAM CH3 DN hydrogels and compared their mechanical properties using ensile test. The agar/ PAM DN hydrogel achieved 1.1 MPa of tensile stress at a breaking strain of 20 32 , which was significantly higher than 0.1 MPa of tensile stress at a tensile strain of 13 for agar/PAM CH3 DN hydrogels (Fig. 3d). The computational and experimental data collectively confirm the key components for modulating the network interactions and the important roles of hydrogen bonds in the mechanical reinforcement of agar/PAM DN hydrogel.

Confined water behaviors in agar/PAM DN hydrogel
Water molecules are essential and dominant component (50-90 wt%) for retaining the structural, mechanical, and functional properties of hydrogels. The presence of inhomogeneous network structures and chemistries in hydrogels strongly affect the structure and dynamics of water molecules inside hydrogels, which would, in turn, influence network-network and networkwater interactions. Here, we used multiple parameters of selfdiffusion coefficients (SDC), radial distribution function (RDF), coordination number, and residence time of water molecules to characterize the structure, dynamics, and interactions of water molecules within hydrogels. In agar/PAM DN hydrogel, water molecules formed hydrogen bond chains with both networks for stabilizing network structure and enhancing network interactions.
We first calculated the SDC values of the water molecules in SN and DN hydrogels to characterize their dynamics in confine environment. The results showed that the SDC values of the water molecules in SN and DN hydrogels are much lower than bulk water (~2.3 × 10 −5 cm 2 /s), supporting the role of agar and PAM network in restraining the water molecules inside hydrogels (Fig.  4a). The water molecules are relatively more dynamic in PAM SN hydrogel (~0.29 × 10 −5 cm 2 /s) than agar SN hydrogel (~0.23 × 10 −5 cm 2 /s), implying that agar network has the strong interactions with water molecules than PAM network (Fig. 4a). Notably, the SDC value of agar/PAM DN hydrogel (~0.07 × 10 −5 cm 2 /s) was significantly lower than that of both SN hydrogels, suggesting that the DN structure with densely packed and interpenetrating chains enables to bind and trap water molecules more strongly and efficiently.
To better quantify water distribution and density around the two networks, the radial distribution functions (RDFs) of water molecules around hydrogen bonding sites in both first agar and the second PAM networks were calculated. As shown in Fig. 4b, two networks presented two major, but distinct distribution peaks from water RDFs at 1.5-2.5 Å and 2.5-4.5 Å, respectively. Water molecules with a separation distance of <4.5 Å from both networks formed a strong hydration shell due to strong network-water interactions and network confinement effect. We analyzed the mean residence time (τs) to further evaluate water dynamics in the strong hydration shell (<4.5 Å) (Fig. 4c). For comparison, τs values in agar and PAM SN hydrogels were~450 ps and~508 ps, respectively, which were comparable to each other. But, both values were significantly smaller than τs in agar/PAM hydrogel (~754 ps), again suggesting that DN structure imposes strong structural confinement and interactions on water dynamics in hydrogels. Agar and PAM networks are hydrophilic polymers containing hydrogen-bonding-rich groups, which offer sufficient binding sites to water molecules. As shown in Fig. 4d, in the first agar network, the coordination number of water around hydroxyl (OH−) groups (~5.0) were higher than that of water around carbon oxygen (O−) group (~3.2). Visual inspection of MD trajectories revealed that the O− groups in the agar network were more buried while the hydroxyl (OH−) groups were more exposed at the agar surface, so that OH− groups attracted the more surrounding water molecules via the massive HB interactions. In the second PAM network, the NH2− and C=O group showed the comparable water coordination numbers, i.e., NH2 (~4.8) vs. C=O (~4.2), indicating that both groups contribute equally to water distribution.

Fracture process of agar/PAM DN hydrogel
It is generally accepted that the "sacrificial bonds" mechanism can well explain mechanical behaviors (e.g., toughness enhancement, large hysteresis, and necking phenomena) of chemically-crosslinked DN hydrogels on the basis of irreversible break of covalent bonds. However, this mechanism may or may not work for physical DN gels involving non-covalent bond-induced mechanical enchantment and recovery. To examine the functional role of the first and second networks in the mechanical properties of hybrid agar/PAM DN hydrogel at the deformation state, we performed the steered MD (SMD) simulations that mimic tensile test to study the molecular origin of the fracture and toughening behaviors of agar/PAM DN hydrogel (Fig. 5a and Supplementary  Movie 3). Figure 5a shows the representative snapshots at each stage during the stretching process. By pulling both networks from state ① to ② at a stretching rate of~0.5 Å/ps, i.e., before achieving the fracture point, the overall network organization in DN hydrogel and the morphology of the second PAM network remained almost intact, but the first agar chains were slid down along its helical axis and their helical conformations were mildly elongated, and the corresponding pulling force achieved local maximum at the fracture point ②. The force-distance curve from our uniaxial SMD simulations in Fig. 5b clearly showed a yielding point, as shown by a peak at the low distance. Further increase of the stretch distance led the gel to reach maximal fracture force of~70 nN at ②. The typical elastic zone with a yield point and a plastic zone with a fracture point were in good qualitative agreement with the experimental stress-strain data.
Once achieving the fracture point, the initial fracture occurred in the first agar network through the chain pulling out mechanism and the corresponding energy dissipation is likely contributed by the disruption of strong and rigid hydrogen bonding interactions in the agar bundle, responsible for the yielding behavior of DN hydrogel. The hydrogen bond profiles show that the first agar network breaks first at the stretching distance of~25 Å, followed by the second PAM network at the~35 Å (Fig. 5c). After the fracture point, due to the loss of massive hydrogen bonds in the agar network, continuous pulling from ② → ③ caused a whole agar bundle to be further separated into the small ones, where some agar helixes were almost completely pulled out from the others by unzipping interfacial interactions between two agar helixes, with partial disruption in their helical conformations. Meanwhile, following the disruption of first agar network, PAM chains in the second network experienced the stepwise conformational changes from the entangled coils to the stretched chains, dissipating the energies. Even after the first agar network was completely disrupted, the chain entanglement in the second PAM network remained, leading to the continue consume of the strain force. So, during this deformation process from ② → ③, both networks were largely deformed to dissipate energy in a progressive way, and consequently, the force increased as stretching distance until reaching to maximal force at fracture point. After that, a steep decline at large distance was observed, which was due to the complete dissociation of both networks (see ③ → ⑤). The mechanical reinforcement of agar/PAM DN hydrogel strongly depend on its unique interpenetrating but independent networks with contrasting network structures and different crosslinkers. The agar network is rigid and physically linked, while the PAM network is more flexible and chemically linked, where both networks interact with each other via hydrogen bonds and chain entanglement. Soft and flexible PAM chains tended to bind with rigid agar molecules by acting as a "belt" to wrap around the agar bundle or a "glue" to fill the gap between two agar helixes (Fig. 5a). This indicates that PAM chains are more proactively and readily adapted their conformations to accommodate agar structures and to maximize their intermolecular interactions. Moreover, considering that the agar network is physically linked by hydrogen bonds, at the beginning of the stretching process, the rigid agar network is easier to be fractured than the flexible PAM network by pulling out agar chains from the first network. At this stage, the agar network bears most of stress with a quick decrease of hydron bonds before the fracture, and energy dissipation (Fig. 5b). As applying force increases in the hardening stage, the stress will be transferred from the agar network to the PAM network, due to high HB density and strong chain entanglement between the two networks. Of note, before reaching the fracture point, water molecules in the strong hydration shell (<4.5 Å) diffused to the stretching-induced network voids, resulting in the loss of water molecules in the strong hydration shell (Supplementary Fig. 7). This process helps the stretching-induced network voids to be filled by neighboring water molecules for maintaining the overall integrity of the agar/ PAM DN hydrogel. But, when the network voids become too large to be filled by water molecules, the crack starts to occur and propagate. While the network voids became more exposed to the air, the water molecules showed an opposite behavior. Water molecules that are far away from the network structures diffused back to the network areas. The number of water molecules in the hydration shell increased (Supplementary Fig. 7). These water molecules tend to bind to the stretching-induced exposed HB sites in the agar and PAM networks for enhancing network structures and contributing to the energy dissipation during the fracture process. In addition, PAM chains have more flexible dynamics at the interfaces than agar chains. It would take longer time to stretch PAM chains out from the DN structure and thus give rise to extra stress in the later stage of the tensile elongation process.
In summary, fundamental understanding of the structure-property relationship in DN hydrogels is of importance in rational design and fabrication of high-performance polymer hydrogels with desired properties. While significant efforts and progress have been made experimentally for the design, synthesis, and demonstration of DN hydrogels with superior mechanical properties, molecular-level understanding of network structure-interaction-mechanics relationship remains largely unknown. Herein, we developed a multiscale computational platform, integrating RWRP, MD, and SMD methods, to computationally polymerize the agar/PAM DN hydrogels, and then to examine the structure-induced mechanical enhancement of the agar/PAM DN hydrogels in comparison with those of individual SN hydrogels. The RWRP-constructed agar/PAM DN hydrogels resembled the experimental properties of network structures and dynamics, mechanical stress-strain behavior, and water dynamics to the larger extent. Side-by-side comparison of network structure/ dynamics/interactions between agar SN, PAM SN, and agar/PAM DN hydrogels revealed a distinct fracture behavior of agar/PAM DN hydrogel, during which rigid agar molecules were first pulled out from helical bundles by breaking hydrogen bonds, chain entanglements, and water association, all of which help to dissipate energy efficiently throughout the network and thus to bear large mechanical deformation. This computational work constructs a DN hydrogel from monomers and studies the early fracture behavior at atomic details. The RWRP protocol allows to computational design and predict of new tough materials on the basis of their structureinteraction-mechanics relationship.

Computational construction and modeling of hybrid DN hydrogel
We previously reported a two-step method to synthesize agar/PAM hybrid DN hydrogel by forming the first physically-linked agar network of hydrogen bond associated helical bundles, followed by photopolymerizing and crosslinking the second chemically-linked PAM network by N,N'-methylene-bis-acrylamide (MBA) within the first agar network 32 .
Here, we developed a computational strategy and algorithm to mimic this two-step radical polymerization process and to construct agar/PAM DN hydrogel computationally (Fig. 1a).
Modeling of the first agar network. Agar is a linear polysaccharide chain with alternating D-galactose and 3,6-anhydro-L-galactopyranose. Due to its sol-to-gel property, agar molecules undergo a coil-to-helix structural transition to form a network via helical bundles association by hydrogen bonds. To construct the first agar network in the gel state, we first used three agar double helixes (PDB ID: 1AGA) to construct triangular 3-fold agar bundles as a building unit, and then we assembled ten triangular agar bundles together and minimized them in energy to form the first agar network. The triangular agar bundle contains a total of~60 D-galactoses and~60 3,6-anhydro-L-galactopyranoses ( Supplementary Fig. 1). Two independent trajectories were provided to visualize the formation of the first agar network (Supplementary Movie 1).
Modeling of the second PAM network. Upon the formation of the first agar network, we developed and used a "random walk reactive polymerization" (RWRP) protocol construct the second PAM network crosslinked by MBA around the agar network, which mimics the free radical polymerization. Specifically, a total of 4389 AM (N AM ) monomers, equivalent to AM concentration of~3.4 M in experiments 41 , were randomly added around the first agar network with random orientations and positions. To initiate the polymerization, 1 mol% AM monomers were activated and acted as initiators to react with other AM monomers to form polymer chains. Three index lists of AM monomers (M list ), initiators (R list ), and PAM chains (P list ) were generated at the first polymerization step and then updated for every following polymerization step. For each polymerization attempt, an initiator (r i ) was randomly selected from R list , followed by a distance search of nearby AM monomers within a reaction distance of d i (reaction radius, starting with 10 Å), which can be adaptively increased by a factor of~10% if no AM monomers are available for polymerization. If found, one AM monomer was randomly selected to react with the initiator to form a new bond, to reconstruct a connectivity table (i.e., all bonds, angles, and torsions), and to update M list , R list , and P list lists. Then, the energy minimization was performed to remove any bad contact and optimize the system, followed by a short MD simulation to allow the diffusion and dynamics of AM monomers and PAM chains in the relaxed DN hydrogel. This process was iteratively preceded and then terminated until all monomers in M list were consumed. The formed PAM chains were then chemically linked by MBA. MBA molecules were added into the system by randomly selecting the paired PAM chains in the system using grand canonical Monte Carlo (GCMC) simulations, as determined by an acceptation probability of an addition of P accept where ΔE is the change in configuration energy due to the addition of MBA. We converted chemical potential into the concentration (number) of MBA for controlling insertion trials. Finally, the resultant agar/PAM DN networks were solvated by the TIP3P water molecules. The water molecules within 2.4 Å were removed to eliminate the atom collapse. In the modeled agar/PAM DN hydrogel, the water molecules occupied~68% of the total system, consistent with the 50-90% water content in the experiments. The final system was minimized in energy and equilibrated by MD simulations at the room temperature of 298 K using NPT ensemble and 3D periodical condition. Meanwhile, the similar modeling procedure was applied to agar single-network (SN) hydrogel, PAM SN hydrogel, and the agar/PAM CH3 DN hydrogel for comparison.
Mechanical test of DN hydrogel by steered MD simulations. Steered MD (SMD) simulation, as a non-equilibrium MD method, was used to mimics tensile tests by applying external force on the system of interest to obtain the work on the unfolding pathway of the system. Here, in our SMD simulations, one end of agar/PAM DN hydrogel was fixed by positions, while another was pulled by an external force of~100 kcal/mol/Å 2 at a constant pulling rate of~0.5 Å/ps. The force-related properties as a function of time (i.e., puling distance) were obtained during the pulling process. SMD simulations were performed by the NAMD software with the CHARMM force field for water molecules and the CAHRMM-compatible CGenFF force field for AM, AM CH3, and agar molecules 42,43 . The Langevin thermostat method with a damping coefficient of 1 ps −1 was used to control the room temperature at 298 K. Short-range van der Waals (VdW) interactions were calculated by the switch function, while long-rang electrostatic interactions were calculated using the force-shift methods. The velocity Verlet integration was used to solve the Newtonian motion with a timestep of 0.5-fs with covalent bonds involving hydrogen being constrained by the RATTLE method. We should note that while SMD simulation was performed at high pulling speed (due to the limitation of computing power) as compared to experimental ones, it is generally accepted that the pattern of force-induced changes in physicochemical properties and potential of mean force will remain consistent along the lowest energy trajectories. VMD software was used for simulation and modeling 44 .