The skeleton of many natural and artificial soft materials can be abstracted as networks of fibers/polymers interacting in a nonlinear fashion. Here, we present a numerical model for networks of nonlinear, elastic polymer chains with rate-dependent crosslinkers similar to what is found in gels. The model combines the worm-like chain (WLC) at the polymer level with the transition state theory for crosslinker bond dynamics. We study the damage evolution and the force—displacement response of these networks under uniaxial stretching for different loading rates, network topology, and crosslinking density. Our results suggest a complex nonmonotonic response as the loading rate or the crosslinking density increases. We discuss this in terms of the microscopic deformation mechanisms and suggest a novel framework for increasing toughness and ductility of polymer networks using a bio-inspired sacrificial bonds and hidden length (SBHL) mechanism. This work highlights the role of local network characteristics on macroscopic mechanical observables and opens new pathways for designing tough polymer networks.

## Introduction

Polymer networks are the basic building blocks of many natural and man-made materials. From cellular cytoskeletons to whole organ tissues, and from rubbers to gels, the load-bearing structures in these materials may be represented as a complex network of polymer chains that are interconnected by cross-linker molecules. Recent studies of soft materials have shown that their mechanical response is rate dependent [1] and the damage may be recoverable [2]. For instance, the crosslinking via ionic bonding leads to formation of sacrificial sites which, when broken, may release additional length of the polymer chains and result in improved toughness and ductility. In the absence of external force or at very low loading rates, new ionic bonds may form leading to time-dependent strength recovery [3]. The sacrificial bonds are generally weaker than the base network's crosslinker bonds. This sacrificial-bonds-and-hidden-length (SBHL) mechanism has been identified in a variety of biological materials including noncollagenous interfaces in bones [4–6], mussel threads [4], and mother of pearls [5,6]. Similar mechanisms were recently introduced for designing tough hydrogels by incorporating recoverable sacrificial bonds into the network [7,8]. Bond damage and formation in these systems depends on the deformation level and loading rates [9]. Furthermore, the mechanical response and fracture properties of the soft polymeric materials highly depend on their topological structures, including local connectivity, bond types, and crosslinking density. However, because the current experimental techniques do not allow for direct visualization of the topology and deformations of the single C–C chain in the polymer networks, the structure-property relation of these materials remains elusive, and the design of the polymeric materials for specific properties is mostly based on trial and error.

Historically, modeling soft materials, like rubber, has been approached using continuum theories. Although these models have been useful in addressing several mechanistic aspects of soft material elasticity in the presence and absence of infused fluids, most classical continuum models from linear elasticity, hyperelasticity [10–15], viscoelasticity [16] to poroelasticity [17–20] lack information about network topology and thus fall short of the capability to guide the network design and may lead to false prediction for fracture susceptibility [1,21]. Understanding the structure-function relations in these networks and in particular the limits on their ductility and rate-dependent responses continues to be one of the long-standing challenges in the field.

To address this challenge of connecting structure to function in polymer networks, several research groups have investigated the mechanics of discrete polymer networks. Some of this prior work has focused on explicit representation of polymer chains using all-atoms molecular dynamics approaches [22,23]. However, these simulations are computationally intensive and are limited to only very short time and length scales. To overcome this limitation, many coarsening approaches were taken. The finite extensible nonlinear elastic (FENE) model [24–26] was proposed as a coarse-graining approach by modeling the polymer chain as a collection of spherical beads (representing monomers) interacting through elastic springs governed by Lennard Jones or other appropriate nonlinear potentials [27]. The FENE is simple and computationally more efficient than the all-atoms representation. However, bending stiffness of polymer chains is not explicitly included and the effect of studying network-function relations remains computationally infeasible. In the current paper, we take one further step in the coarse-graining efforts and use an effective force-extension law on the scale of the individual chain without explicitly modeling the monomer dynamics as in the FENE approach. Our choice of the worm-like chain model also enables us to account for the chain flexibility and bending stiffness.

The eight-chain model of Arruda and Boyce [10] was a significant step forward in encapsulating small-scale network information at a continuum nonlinear elasticity model for polymeric materials. However, the model excludes rate dependence and requires further modification to enable studying damage evolution as it does not contain information about how the network changes structure with loss and formation of polymer chains. One classical model that attempts to incorporate network information into the constitutive equations is due to Lodge [28]. In this work, a loss probability is incorporated to study stress-dependent damage to the network. The work also studies the relation of concentration densities and co-ordination numbers to the constitutive equations. However, that study was limited to rate-independent mechanics of networks composed of units with similar topology. Moreover, affine deformation of junction points is assumed which may not be the case as pointed out later.

Extensive studies in the physics community have been conducted on two-dimensional (2D) and three-dimensional (3D) polymer networks with a primary focus on understanding the role of nonaffine deformation [29] as well as limits on rigidity percolation [30] due to reduction of polymer concentration. These models have been mostly limited to networks with fixed topology and have not considered the consequences of topological evolution due to damage accumulation as a function of stretch or the influence of loading rate [31]. Other network models from the chemical–physics community have investigated effective mechanical properties of randomly crosslinked polymer networks [32–35]. An important result is that the mean field theories fail to explain the variance in the force displacement characteristics of filamentous polymer networks, especially for networks with low concentration [36,37], suggesting a need for explicit incorporation of small-scale networks structure for consistent upscaling. Furthermore, in the works that have considered damage accumulation, the models have largely neglected rate effects and, in most cases, have modeled damage based on a maximum force (stress) criterion for the bonds in the network [38].

Here, we introduce a nonlinear network model with a physics-based rate-dependent damage criterion. The main contribution of this paper is twofold. First, it introduces a physics-based model for the rate-dependent behaviors of polymeric materials in the presence and absence of sacrificial bonds using the transition state theory coupled with bulk viscous dissipation. This framework allows us to investigate how different features of the network topology influence the rate-dependent force-displacement response. Second, it presents a theoretical framework for the design of a SBHL system in the polymeric network that enables significant improvement in toughness and ductility. This framework provides bounds on the tunability of the different network parameters in order to maximize the intended property improvements from the SBHL system.

The paper is organized as follows: In Sec. 2, we review the worm-like chain (WLC) model and introduce the rate-dependent damage model for bond breakage based on the transition state theory. In Sec. 3, we describe the model setup and the governing equations for the quasi-dynamic simulations. In Sec. 4, we summarize our results for the effect of loading rate, chain length disorder, crosslinker density, and internal network architecture (e.g., sacrificial bonds and hidden length) on the system stiffness, peak force, and ductility. In Sec. 5, we discuss the implications of our results for design of tougher networks and for understanding the mechanics of crosslinked polymer chains with and without sacrificial components.

## Theory

where *k _{B}* is the Boltzmann constant,

*T*is the temperature,

*b*is the persistence length,

*x*is the distance between the two ends of the polymer chain, and

*L*is the available contour length. In the limit where the chain is very long compared to this persistence length, it recovers the freely jointed chain behavior, and in the limit where the chain is short compared to its persistence length, it describes rod-like chain behavior. The main distinguishing feature of the WLC model is that the motion of successive segments of the polymer chain is correlated resulting in a finite bending stiffness. This makes the WLC model suitable for modeling a variety of biological materials including DNA, RNA, and several proteins. Given the similarities between hydrogels and many biological tissues, we adopt the WLC in this study. However, we emphasize that the framework presented here is flexible enough to implement other chain-level force-displacement relations including those which may be directly computed based on molecular dynamics simulations.

To model the network of interconnected polymers, we picture a graph $G$ that occupies a bounded domain $\Omega \u2208E2$ (2D Euclidean space). The graph is uniquely determined by a set of nodes, *V* (vertices that represent crosslinking connections) and a set of edges, *E* (the polymer chains). Each of these edges is an unordered two-tuple (*i*, *j*) that uniquely represents a polymer connected across two crosslinker nodes, *i* and *j* (*i*, *j* ∈ *V*). We define the coordination number, *Z* of a node to be the number of connections of that node to other nodes in the graph, $G$, i.e., $Z=\u2211j\u2208V\u2216i1((i,j)\u2208E)$ (where $1$ is the indicator function).

### Damage.

In this model, we consider damage events as the breakage of the bond between the polymer and the crosslinking molecule, which are rate-dependent processes that are influenced by the force at this junction due to the stretch of the polymer chain.

*k*(bond formation) and backward reaction,

_{f}*k*(bond breakage) that depend on the force applied by the polymer chain on the crosslinker—polymer junction as shown in Fig. 1. Locally, the reaction kinetics are assumed to follow Arrhenius rate equations

_{b}*E*is the activation energy for the backward or forward reactions. Following Bell's theory [45], the activation energy consists of two parts: a reference value that represents the intrinsic energy barrier to thermally activated processes in the absence of external forces and an additional force-dependent term that lowers the energy barrier even at constant temperature, i.e.,

_{a}*E*=

_{a}*E*

_{0}–

*F*Δ

*x*. Here, Δ

*x*represents the distance to a transition state (see Fig. 2). The characteristics of bond dynamics, such as the activation energy or the transition distances, may be computed using quantum mechanics or estimated experimentally. They are mainly functions of the bond composition. The

*E*

_{0}term, being a constant, may be incorporated directly in the proportionality constant and finally, we have

*x*and Δ

_{b}*x*are the distances to the transition states and

_{f}*α*and

_{b}*α*are, respectively, the backward and forward rates when no force is applied. For a detailed explanation of these transition distances, cf. [44]. Both of these parameters are assumed to be constant in this study. We assign a binary variable,

_{f}*N*∈ {0, 1} for each edge present in the graph, $G$.

*N*= 0 would represent the crosslinker connection is broken and

*N*= 1 indicates that this bond is intact. To facilitate the computation, we introduce a continuous variable representation of

*N*,

*n*evolving as in Ref. [43]

Thus, by defining $\delta =\u222b0\Delta t(dn/dt)dt$, a bond is considered broken if *δ* ≤ –1 and formed if *δ* ≥ 1.

### Hidden Lengths.

A toughening phenomenon observed in many biological materials [4,46,47] and synthetic hydrogels [48] is the effect of hidden length in the polymer chain. These are essentially self-linked loops within the polymer chain, which reduce the available length to be stretched and stiffen the chain (see Eq. (1), where the chain stiffness increases at a given end-end distance *x*, if the chain length *L* is smaller). However, if the self-links are weaker than the end crosslinker bonds, they would break before the end crosslinker bonds. The breaking of these sacrificial bonds opens additional lengths of the polymer chain, and thus enhances the toughness and/or stretchability of the material [43,49]. We will discuss more on this in Sec. 4.

*m*to be the number of sacrificial bonds in each chain and formulate the rate-dependent damage model as

*m** is the continuous representation of the number of sacrificial bonds and

*k*is the rate of bond breaking. Although, in this work, we do not include healing (or self-recovery) of the sacrificial bonds, it can be easily incorporated using an additional forward term in the rate equation as discussed earlier.

_{sb}*k*is defined as the rate of sacrificial bond breakage per unit time, which is related to the force on the polymer chain as

_{sb}## Setup

All of the simulations reported in this work are in 2D Euclidean space with a rectangular sample chosen in two major topologies. We label the topologies using the coordination number of individual crosslinker molecules. Specifically, we test a network that has an average coordination number of *Z* ≈ 4 (i.e., predominantly quadrilateral mesh) (Fig. 3(b)), and *Z* ≈ 7 (i.e., predominantly triangular mesh) (Fig. 3(a)). We introduce disorder in terms of a distribution of the available chain lengths *L* and the placement of the crosslinker nodes in the network. Other sources of disorder in the network are explained as they appear. We chose our bond parameters as Δ*x _{f}* = 0.25 nm, Δ

*x*= 0.1 nm,

_{b}*α*= 0.1 s

_{b}^{−1}, and

*α*= 0.3 s

_{f}^{−1}as per [43]. We simulate a persistence length of 0.3 nm at a temperature

*T*= 300 K. The initial sample size is 500 × 100 nm.

In this work, we assume a 2D network that is fixed on the bottom and attached to a moving plate on the top. Periodic boundary conditions are applied on the lateral sides. The top plate moves at a constant speed in the *y*-direction as shown in Fig. 4.

*i*, we write

where *C* is the damping coefficient, $Fi=\u2211\u2202iFij$, with *∂i* being the neighbors of node *i* in $G$ and **X*** _{i}* denoting the position vector of the node

*i*. At every time-step, we use predictor-corrector schemes until we reach $L2$-convergence in the positions $Xi\u2009\u2200\u2009i\u2208V$.

We record the force on the top plate as the network is stretched. Since all the chains of the network follow the WLC force-displacement relation, the most stable equilibrium point without any constraints would be if all nodes lump into a point. Since we start with a finite sample size and do not incorporate the excluded volume effects, we have a certain initial force that is required to maintain the initial size. This initial force in the gel system would be the osmotic pressure supported by the solvent. It is also the reason that we do not impose volumetric constraint in the model—as the network of the gel is stretched, the solvent will flow in to fill the extra space. When the top plate is moved, a minimum energy optimization is run over all nodes. Due to the finite bond strength of crosslinker molecules and Dirichlet boundary conditions, we restrict the equilibrium space, and therefore, perform a constrained optimization. Upon stretching the network, appropriate damage characteristics are updated (i.e., Eq. (4) for crosslinker bond damage and Eq. (5) for sacrificial bond damage). Whenever a crosslinker bond is broken, that particular polymer chain is removed from the network, i.e., if (*i*, *j*) represents a polymer chain, a *breaking* event removes the edge (*i*, *j*) from the graph. If all the edges at a node *k* are removed, node *k* is also removed from the optimization routine. The simulation continues until there is no continuous path from the top plate to the bottom plate that bears the load of the moving plate. The sacrificial bonds are also updated in a similar fashion. Furthermore, if the damage integral for the sacrificial bonds crosses −1, a sacrificial bond is broken, a hidden length is added to the parent chain, and the damage integral is reset. This procedure continues until no sacrificial bond remains or the end bond breaks, whichever is first.

With this model, we discuss the effect of the following parameters in the simulation:

Pulling velocity: The top plate is moved at different velocities to investigate the rate effects on the network peak force and ductility.

Chain lengths: All chain lengths

*L*are drawn from a distribution (uniform/Gaussian) with a prespecified mean and a measure of width. This width of the distribution is a measure of the disorder in the system.Sacrificial bonds: In order to observe the effects of hidden lengths in the network, we simulate networks with and without sacrificial bonds. We choose to introduce a uniform random number of hidden lengths,

*m*, drawn from a uniform discrete distribution. The hidden length in each of*m*segments is drawn from a uniform continuous random distribution.Network topology: We consider networks with both quadrilateral and triangular meshes.

## Results and Discussion

In all the results in this section, the force is plotted against the stretch *λ* which is defined as $\lambda =(h\u2212hi/hi)$, with *h _{i}* being the initial height and

*h*the current height (in the

*y*-direction) of the sample. All force-stretch plots have inherent stochasticity in the results due to variability in the chain length distribution in each realization (see Fig. 5). One should note that such uncertainty is more visible around the peak and softening part of the curves (as many chains break in this portion). Moreover, this effect is more pronounced for smaller networks (networks with a few polymer chains, Fig. 5) than for larger networks (>5000 chains or graph edges). For this reason, large networks are simulated in this work.

### Effect of Pulling Velocities.

where *F* is the force magnitude for given chain length *L*, *s* is the distance between the two crosslinkers (i.e., ends of a polymer chain), and *v* is the pulling velocity. The above equation has been written for breaking the crosslinker bond. Therefore, for the same amount of displacement, a higher velocity implies a higher threshold for breaking the bond, i.e., there is a higher stretch allowed. This, consequently, leads to a higher peak force in the worm-like chain paradigm.

Figure 6 shows that the peak force and the maximum stretch ratio vary almost linearly with respect to the log of the pulling velocity. To understand these effects further, we do an order of magnitude analysis as follows. Without inertia, a first-order linear one-dimensional system can be described by $Cx\u02d9+kx=0$. The relaxation time for such a system is given by Θ(*C*/*k*), where *k* is the stiffness of the system and *C* is the damping coefficient as described in Eq. (7). In our analysis, we have a nonlinear system, and hence, this result is not directly applicable, but one may write a linearized version of the system around a given stretch value, *λ*_{0}. This linearization will hold around a small neighborhood of *λ*_{0}, and in this neighborhood, the stiffness may be assumed to be constant. For the simulation results shown in Fig. 6, *C* = 10^{−3} Nsm^{−1}, while the stiffness during loading goes up to 10^{2} N⋅m^{−1} which gives a relaxation time scale of Θ(10^{−5}) s. At the peak load and during the softening phase, the stiffness decreases drastically and often is very close to Θ(10 deg) N⋅m^{−1} (especially during peak). This leads to a relaxation time scale of Θ(*C*) s. Therefore, it is safe to assume that the order of relaxation time scale of the system Θ(*τ*) ≥ Θ(*C*) = 10^{−3} s. However, all the time scales ($1/\lambda \u02d9$) associated with the loading rates considered here are much longer than these relaxation timescales and hence, the bulk viscosity, modeled through *C*, contributes negligibly to the observed rate effects in this case. This suggests that the rate dependence of the peak force and ductility shown in Fig. 6 originates from the bond dynamics that are governed by Arrhenius-like processes. This explains the logarithmic dependence.

In Fig. 7, we probe the time-scales that are comparable to the inherent time scale of the system, i.e., where *C*/*k* is comparable (in order of magnitude sense) to $1/\lambda \u02d9$. Hence, we plot all quantities of interest with respect to the ratio of timescales leaving out the stiffness term. We note a logarithmic dependence in the peak force on loading rate initially which changes to an almost linear relationship and then falls back. During the low $C\lambda \u02d9$ regime, we have the chain force dominating the equilibrium and the response is sensitive to bond dynamics as previously explained. In the high $C\lambda \u02d9$, the viscous force dominates, and hence, the relationship becomes approximately linear. Within this regime, the system is increasingly sensitive to loading rate and initial conditions. In particular, the damage and deformation becomes highly localized in the top portions of the network and there is not enough time for the force to redistribute and percolate downward. This explains why at very high $C\lambda \u02d9$, a very small part of the network takes the entire load and breaks prematurely at a lower load (Fig. 8).

To examine the robustness of these results, we present an estimate for the expected systematic randomness if we conduct these runs for a large number of network realizations. We have done this at lower loading rates and the error bars are shown in Fig. 7, suggesting that the standard deviation is small in this limit due to the engagement of the full network in load resistance. In order to correctly place the error bars at higher loading rates, we account for the reduced engagement and deformation localization by increasing the standard deviation of the results using central limit theorem arguments. Physically, this means the lower the number of chains/edges of graph engaged (stretched out) in the pulling experiment, the larger the variation in different realizations of the same network. As an approximate estimate, if 10% of the edges of the network are engaged, the standard deviation in the simulation is expected to be approximately $1/0.1=3.16$ times the standard deviation in results when the entire graph is engaged (assuming that the variance is normalized by the mean of the distribution). With this estimate, we may conclude that the drop in the peak force and dissipated energy at the highest loading rate is a robust observation and it happens when the loading rate becomes comparable to the intrinsic viscous time scale in the network due to a transition in the damage mode from distributed to localized.

### Effect of Crosslinker Concentration.

We simulate the effect of increasing crosslinker concentration while keeping the polymer concentration constant. This is done by increasing the number of nodes and edges in the graph while keeping the total weight of the network (sum of lengths of all chains) constant. As seen in Fig. 9, the peak force increases up to a maximum value with increasing crosslinker concentration. Further increase in the crosslinker concentration leads to a decrease in the peak force. This may be explained as follows. As the crosslinker concentration increases, the average chain length decreases (to keep the total weight constant). Therefore, the system becomes stiffer and attains higher forces for a given stretch of the network. As the crosslinker concentration further increases beyond a certain threshold, the network becomes critically stressed due to the prevalence of short chains, (low *L*). These short chains are prone to damage accumulation at a faster rate because of higher forces at the scale of individual chains (due to high *x*/*L*). This leads to premature failure at lower global force level. In other words, the forces start off high at the chain level triggering a cascading failure mechanism in the network that lowers the peak force with respect to the initial force.

### Effect of Disorder in Available Chain Lengths.

We simulate disorder in the network by generating available chain lengths of the polymers from a stochastic distribution. We consider both Gaussian and uniform distributions. A higher disorder refers to a larger width (or standard deviation) of the distribution. As seen in Fig. 10(a), the force displacement response weakly depends on the degree of randomized disorder for the range of parameters considered here. The uniform distribution leads to a slightly more pronounced effect than what is observed for the Gaussian distribution case. There is a slight increase in stiffness and slight decrease in peak force. This is due to the nonlinear power law nature of the force curve of a single worm-like chain (i.e., the force rapidly increases for stretch values close to the contour length). The increase in force due to shorter lengths dominates over the decrease in force in larger length chains. This causes an overall increase in the stiffness of the network. However, this also means that damage accumulates faster, and hence, the peak force attained by the network is reduced. The network with higher disorder level also exhibits a larger softening distance. This suggests that disorder may help the network survive longer under damage. We discuss this further in Sec. 5.

### Effect of Coordination Number.

Here, we test the two different topologies. We keep the total weight of the network constant while also maintaining almost the same number of edges in the network. We draw chain lengths from the same random distribution for both cases. As seen in Fig. 11, a network with higher coordination number has higher toughness, peak force and also survives a larger softening (post-peak) distance. This may be explained as follows. As the network topology changes from quadrilateral to triangular, the average chain length decreases (for the same total weight) leading to increased stiffness and higher peak force. Furthermore, due to the increased number of edges per node in the triangular topology, the force redistribution due to the loss of a polymer chain is more efficient. The force increase on adjacent polymer chains during this redistribution is smaller than in the case of quadrilateral and this leads to delayed damage accumulation. On the other hand, a polymer chain detachment in the quadrilateral network leads to a larger force drop and higher force increase in the surrounding polymers. This triggers faster cascading failure and a more brittle behavior. Similar results were obtained in Ref. [51] for nanocomposite gels that increase the coordination number of the network.

### Effect of Sacrificial Bonds.

We consider two scenarios for introducing the sacrificial bonds and hidden length idea into the design of our network. In the first scenario, the total length of the polymers in the network including the hidden lengths is kept equal to the total lengths of the polymers in the network without sacrificial bonds. Figure 12 shows that sacrificial bonds, in this case, increase the peak force and toughness of the network but have negligible effect on the stretchability. In this particular simulation, sacrificial bond parameters were Δ*x _{sf}* = 0.105 nm, while the end bond Δ

*x*= 0.015 nm. This simulates a sacrificial bond that is much weaker than the end crosslinkers. Also, the hidden lengths were kept very small, i.e., a continuous uniform distribution over 5%–10% of the initial available length of the polymer. Therefore, Fig. 12 suggests that even though the two networks have the same (simulated) weight, small-scale changes in the network design, e.g., through introducing small internal loops in the polymer to act as hidden lengths, may allow for extra stiffness, strength, and toughness in the polymeric systems.

_{b}Another possible route is the introduction of additional hidden lengths [7,48] in the network that grant extra stretchability and toughness to the polymeric system. We simulate this by modifying the previous weight constraint. We keep the initial available length of the two networks the same. Hence, any breaking sacrificial bonds add extra length to the polymer chain as in the experiments of [48]. In order to maximize gains in both toughness and stretchability, we investigate the parametric dependence of these gains, within the scope of our model.

*j*=

*x*/

*L*, where

*x*is the node–node distance and

*L*is the available contour length. Instantaneous damage,

*δ*(

*j*), may be written as (

*f*is the force in the chain)

*j*∈ (0, 1), we can ensure that all sacrificial bonds break before the end bonds break, i.e., we use the entire hidden length before breaking end bonds. More formally, if

*m*

_{init}is the initial number of sacrificial bonds, one may write the following condition:

Arguably, this is a very strong condition and thus may rather be regarded as a sufficiency condition. However, it insures that the last surviving sacrificial bond will experience damage at a higher rate than the end bond.

If we choose the strengths and transition state properties such that the right-hand side of the inequality is negative, we satisfy it for all *j*. Therefore, Eq. (13) ensures that we use all the sacrificial bond hidden length before the ends break.

*dj*(

*t*)/

*dt*, we may say that adding a hidden length which changes

*j*to

*rj*,

*r*< 1 such that the damage accumulation rate given by Eq. (10), say

*g*becomes

*g*/

*M*,

*M*> 1∀

*j*>

*j*

_{0}where

*j*

_{0}is the initial stretch in the network. In our simulations,

*j*

_{0}≈ 0.1

If we relate the ratio, *r* to *M*, we may choose a good set of parameters to get the extra stretchability. Examining Eq. (15) closely reveals the minimum is 1 for the left-hand side when *j* = 0. To understand the structure of this inequality, we plot the surfaces parametrized by the tuple (*r*, *j*) for different Δ*x _{e}* values.

From Fig. 13, it may be inferred that as we decrease *r* (i.e., increase the hidden length), we achieve higher *M*. Also, this value keeps getting better as stretch in the network increases. From Fig. 10, we further conclude that in order to get the best stretchability increase between the sacrificial network and the *bare* network, we require a large transition state distance. This result is unexpected because a larger transition state distance increases the damage accumulation rate in the chain, and hence, end-bonds would break faster. However, to compensate for this, we may choose end-bonds that require a higher activation energy, *E _{a}*, to break. This implies a lower rate of breaking at zero force,

*a*. Note that

_{e}*a*values remain the same in both the networks we are testing here. Therefore, a bare network has equally strong end bonds. Also, note that the increase in stretchability does not depend on sacrificial bond parameters (see Eq. (15)). As long as Eq. (13) is met, we can provide a rough estimate of the increase in stretch by plotting the left-hand side of the inequality 15. Moreover, one can see from Fig. 13 that at higher stretch values,

_{e}*M*increases, and hence, we need not have all hidden lengths to be equal. If we allow for gradation in the sacrificial bond parameters and choose the weakest sacrificial bond to have the highest hidden length (

*r*= 0.7) and the other (stronger) sacrificial bonds to have lower hidden length, we can still achieve a significant increase in properties. This is because the other stronger bonds will be opened at a higher

*j*. We leave this strength-gradation study for future work.

We test the above theoretical analysis by simulating a network with *a _{e}* = 0.001 s

^{−1}(i.e., increase

*E*by about 6.9

_{a}*k*) against earlier simulations. We also change Δ

_{B}T*x*to 0.5 nm so as to test whether we get a significant increase as estimated from the aforementioned discussion. Note that the network without sacrificial bonds in Fig. 14 has the same end-bond parameters. We choose an additional length of 40% of the mean length of the network, i.e.,

_{e}*r*= 0.71. Indeed, as shown in Fig. 14, we obtain an increase of ≈ 300% in toughness compared to the reference network with no hidden length. This matches well with our surface plots in Fig. 13 (top-left). Moreover, we observe from Fig. 14 that the initial stiffness of the network with sacrificial bonds is the same as the bare network but it starts softening once the sacrificial bonds start breaking and the hidden length is released. It is interesting to observe that the subsequent breakage of sacrificial bonds and unfolding of the hidden loops lead to saw-tooth-like features in the force-displacement curve that collectively have a similar signature to a yielding plateau. A similar yielding response has been reported for polymer materials as a result of crystallization. Thus, tailoring network topology, using SBHL system, may have similar effects as phase transition.

## Discussion

Soft materials like hydrogels have received considerable interest in recent years. To overcome the intrinsic brittleness of gels and their limited ductility, several strategies for tailoring the network architecture were recently introduced. These gels provide increased fracture resistance. Here, we quantitatively explore the effect of topology on response. More importantly, an alternative mechanism for enhancing ductility and toughness in polymer networks through introduction of hidden lengths via multiple weak sacrificial bonds in the same single network was studied. These bonds break first and upon breakage release additional chain lengths which were previously shielded from direct loading. This SBHL mechanism is biologically inspired and was previously explored by one of the authors in the context of noncollagenous protein interfaces in bone [43,49].

Interestingly, the SBHL system is not limited to bone. Rather, it is widespread in many tough nanocomposite biomaterials such as abalone shell [46], polyprotein dimers [52], and spider capture silk [53]. Recently, Zhu et al. [48] experimentally explored a similar mechanism, inspired by unfolding of protein molecules in mussel threads, to 3D-printed hydrogels with enhanced ductility and toughness. This paper presents a computational paradigm for predicting the rate-dependent response of such artificial networks with sacrificial bonds and suggests a theoretical framework for guiding their design.

The rate-dependent model proposed in this paper has two primary ingredients. The first ingredient originates from the bond breakage and formation at crosslinking sites. This is described using the transition state theory in which each bond is assumed to possess a double well potential with an energy barrier between intact and broken states that depend on the applied force. The second ingredient is a linear viscous force that is hypothesized to simulate the drag force on a polymer network deforming within a fluid domain without explicitly modeling fluid flow. This model is satisfactory for hydrated soft materials with low polymer concentration at low and intermediate loading rates. The model may be extended in the future to account for additional rate-dependent processes that may emerge as the polymer concentration increases such as entanglement effects as well as nonlinear viscous drag forces at extremely high loading rates. A primary result of this paper is that for the range of loading rates and model parameters considered here, the rate-dependent response of the polymer network is not monotonic. At low loading rates, the peak force increases linearly with the logarithm of loading rate. In this limit, bond dynamics dominate the response. As the loading rate increases, the rate dependence of the peak force becomes almost linear suggesting the dominance of the viscous force. At even higher rates, the peak force drops. This nonmonotonic dependence correlates with a change in the deformation mode from being distributed at low and intermediate rates to localized at the highest rates. This nonmonotonic rate dependence represents a testable prediction for future experiments.

Few models exist in the literature for the rate-dependent response of polymer networks with reversible bonding. Krishnan et al. [21] and Hui and Long [54] proposed a 3D finite deformation model for self-healing gels by keeping track of the time evolution of bond breakage and formation. Zhang et al. [55] established a macroscopic model for gel fracture using a cohesive zone law coupled with bulk viscoelasticity idealized using Mullins effect. Our model is distinct from prior work in several aspects. First, our model explicitly accounts for the network topology, an ingredient that has been missing in these previous studies. Second, we use the transition state theory approach to model the bond dynamics which provides a direct link to bond-specific chemical characteristics. Third, none of these models have considered the hidden length mechanism.

Another primary result of this paper is establishing a systematic framework for designing tough polymer networks with sacrificial bonds and hidden length so as to use the entire hidden length efficiently. There are two requirements for this. First, the sacrificial bonds must be weaker than the end bonds so that all the hidden loops unfold before the end bonds break. Second, the released hidden length must be large enough to cause a sufficient force drop on the polymer level that enables it to continue to be stretched before final detachment. We demonstrated one application example in which the toughness of the polymer network is significantly increased (3 times) through the adjustment of the end bonds parameters and the addition of an extra (≈ 40% of initial length) hidden length per sacrificial bond. The outlined procedure complements existing practice in designing tough hydrogels and polymer networks [56] and provides an alternative pathway for enhancing ductility and energy dissipation in polymer networks that does not leverage the double network architecture.

In this paper, we have modeled polymer networks with chain length drawn from both uniform and Gaussian distributions. For the range of length variability explored here, we have found that chain length disorder has a weak effect on the force displacement response. Increased disorder may lead to a slight decrease in the peak force and a slight increase in network ductility as well as a more gradual softening phase. One possible explanation for the lack of strong effect of disorder on network dynamics is the random nature of disorder. It is possible that the effect of disorder to be amplified through introducing an engineered spatial variation in chain lengths rather than through random disorder. For example, one may envision the use of a higher disorder network in crack-susceptible regions and a lower disorder network for other regions. Other possibilities may include introducing patterned blocks of chains with different lengths to localize damage, controlled gradation in chain length or bond strength, as well as strategically placing longer length chains into shorter length networks. These are topics of current investigation in our group.

The results presented in this paper suggest that for networks with simple topology (e.g., clustered networks considered here), the elastic properties are less sensitive to structural randomness in the network and thus may be homogenized using standard procedures [10]. However, the peak force as well as the softening phase exhibits larger variance. We are currently exploring various uncertainty quantification models that are capable of encapsulating this randomness at higher scales.

A future extension of this study will involve investigating the rate-dependent response of 3D networks with several complexities such as multinetwork gels nanoparticle reinforcement, entanglement, and network incompressibility when applicable. These models will ultimately be incorporated in multiscale simulations to explore fracture susceptibility and fatigue resistance due to bond breakage and formation.

## Acknowledgment

We thank the editor and the anonymous reviewers for their comments which improved the manuscript.

## Funding Data

University of Illinois Urbana Champaign, the Campus Research Board of University of Illinois (Grant No. RB17043).

Division of Civil, Mechanical and Manufacturing Innovation (CMMI-Award No. 1435920).

NSF CAREER Award (Award No. 1554326).

Air Force Office of Scientific Research YIP award (Award No. FA9550-17-1-0295).