A central question in soft matter is understanding how several individual, weak bonds act together to produce collective interactions. Here, gel-forming telechelic polymers with multiple stickers at each chain end are studied through Brownian dynamics simulations to understand how collective interaction of the bonds affects mechanical response of the gels. These polymers are modeled as finitely extensible dumbbells using an explicit tau-leap algorithm and the binding energy of these associations was kept constant regardless of the number of stickers. The addition of multiple bonds to the associating ends of telechelic polymers increases or decreases the network relaxation time depending on the relative kinetics of association but increases both shear stress and extensional viscosity. The relationship between the rate of association and the Rouse time of dangling chains results in two different regimes for the equilibrium stress relaxation of associating physical networks. In case I, a dissociated dangling chain is able to fully relax before re-associating to the network, resulting in two characteristic relaxation times and a non-monotonic terminal relaxation time with increasing number of bonds per polymer endgroup. In case II, the dissociated dangling chain is only able to relax a fraction of the way before it re-attaches to the network, and increasing the number of bonds per endgroup monotonically increases the terminal relaxation time. In flow, increasing the number of stickers increases the steady-state shear and extensional viscosities even though the overall bond kinetics and equilibrium constant remain unchanged. Increased dissipation in the simulations is primarily due to higher average chain extension with increasing bond number. These results indicate that toughness and dissipation in physically associating networks can both be increased by breaking single, strong bonds into smaller components.

## INTRODUCTION

Polymeric gels are found in both natural and synthetic materials as soft solids that perform important mechanical functions. In natural materials, gels are used as a building block for a number of components found in the human body. The extracellular matrix provides a supportive scaffold for cells while also regulating intercellular communication.^{1,2} Cartilage acts as a protective barrier in joints that experience high levels of impact.^{3–5} Additionally, protein gels assist in the adhesion of mussel byssi to underwater rock formations, the formation of bacterial biofilms on ships,^{6,7} and the protective responses of animals like whales and hagfish.^{8} In synthetic materials, gels are used in many common household items. Gels are used in pressure-sensitive adhesives as a way to form a molecular contact with rough or non-wetting surfaces and dissipate energy during peeling.^{9–12} Soft contact lenses utilize the ability to keep eyes hydrated and minimize irritation while also remaining optically clear and mechanically robust.^{13,14} Gels are also used in drilling fluids,^{15} diapers,^{16} food products,^{17} and paints.^{18} In all of these cases, better understanding of how specific components affect the mechanical properties is important both for the development of bioinspired mimics and better material design. In particular, many applications benefit from increased material toughness.

Increasing toughness or dissipation in polymer networks is frequently achieved by engineering structures with multiple associating functionalities into a single polymer system. It is possible to couple multiple networks via the formation of interpenetrating or double networks, where different types of associating systems are combined in order to capitalize on specific properties like strength or extensibility.^{19–23} In these and other network systems, methods for energy dissipation include the breaking of bonds,^{24} storing energy in network junctions through additional chain length,^{25} or incorporating multiple associating bonds along the polymeric backbone.^{26} Recent advances have shown the importance of sacrificial bonds in increasing material toughness, where pre-stressed chains are incorporated into the network in order to redirect energy away from the primary network.^{27} Theoretical insight into how to design the arrangement of bonds in physical networks could provide valuable guidance for achieving the next generation of these toughened gels.

Transient network theory is the state-of-the-art theoretical approach for characterizing polymer networks under deformation when the networks have reversible bonds. Transient network theory was initially developed to explain the stress relaxation behavior of rubbers by modifying rubber elasticity to incorporate the idea that bonds can transiently associate and dissociate.^{28} Since then, it has been applied to associative polymer systems, using chemical kinetics to model bond association and dissociation under stress.^{29} Further work has shown that topological characteristics can affect the stress response of the system and decrease the elastic effectiveness of multiple network junction points.^{30,31} The presence of multiple bonds along the backbone of the chain results in deviation from traditional linear Maxwell behavior,^{32,33} and activated chain pullout helps theoretical predictions better capture non-monotonicities in stress.^{34,35}

With the exception of theories that incorporate multiple bonds along the chain backbone, transient network theory has focused on understanding telechelic polymers, where there is just one associating group at each end of the polymer chain. These theories predict Maxwellian behavior during small amplitude oscillatory shear for all telechelic polymers. However, many block copolymer gels, such as pressure sensitive adhesives, tend to have larger values of dissipation and viscous modulus.^{9,10} One key shortcoming of the current approaches is that telechelic materials are treated much the same regardless of whether their associating groups comprise coiled-coils^{36–38} or long-chain hydrophobes.^{39} In particular, the existing approaches to transient network theory are incapable of capturing the behavior of systems with complex endgroups such as block copolymers which may disengage from the network in stages.

This work presents a method for modeling telechelic molecules with associative endgroups where each endgroup is composed of multiple small bonds that are capable of linearly associating into an aggregate or dissociating from an aggregate. This model is formulated with telechelic block copolymers in mind, where the associating endblock domain is composed of many monomers that associate with the aggregate, enabling step-wise disengagement of the chain from the network junction as it is pulled out. This model is then used to study chain relaxation in the linear regime, illustrating the dependence of chain relaxation on the relative kinetics of association and dissociation and the number of bonds per endgroup. Application of the model to nonlinear behavior in shear and extension then provides insight into the effect of increasing the number of bonds on network dissipation during flow.

## METHODS

### Kinetics

The telechelic polymers modeled in this study are composed of end-groups of a defined size that are capable of being fully or partially dissociated. A useful molecular example to visualize the concept is a long associating group composed of many monomers which can dissociate in series from an associative cluster (Scheme 1). Both partial association and full association result in a network of associated domains in solution.

Under the assumption that these polymer chains obey unentangled dynamics that can be approximated as Rouse dynamics under equilibrium and undergo homogeneous flow (a strong assumption in certain types of nonlinear flow where shear banding^{40,41} may occur, but reasonable in quiescent conditions or at low flow rates on the order of or below the chain relaxation time), the distribution of chain configurations as a function of time and degree of association can be modeled by solving the Smoluchowski equation for the probability density of a single chain of length *N*. As in previous work, the deformation processes are modeled in the frame of reference of one chain end.^{34} While one chain end is fixed at the origin, the other chain end is free to react with the surrounding network, transitioning between the dangling state and both the partially and fully associated bridged and looped states. To enable modeling of the system with Brownian dynamics, the equations are formulated in the Langevin form of the Smoluchowski equation for convection, reaction, and diffusion such that the position of the chain end for species *i* (*i* = *D*, *Lk*, *Bk* for dangling, looped, and bridged, where *k* denotes the number of attached stickers from *k* = 1 to *nmax*) can be tracked as a function time ($r\u2192i*t*$),

where the spring constant *H* can be rewritten in three dimensions as *H* = 3*k*_{B}*T*/*Nb*^{2}, $dW\u2192*$ corresponds to a three-dimensional Wiener process (each component is a random variable drawn from normal distribution with zero mean and a variance of *dt*^{*}), and $\u2207\u2192*u\u2192*$ is the velocity gradient, which in shear and extension is (respectively)

Equation (1) in its entirety corresponds to the stochastic differential equation that governs the evolution of dangling chains. Only the first term on the right-hand side of Eq. (1) is included for bridged chains, which are assumed to deform affinely with the network. Looped chains are assumed to remain fixed at $r\u2192*=0$.

Equations (1) and (2) can be made dimensionless in time and space using *t*^{*} = (*ζ*/2*H*)*t* = *τ*_{R}*t* and $r*=Nb2/3r$ [in previous work, space was non-dimensionalized using $r*=Nb2r$].^{34} The dimensionless unit for time was chosen to be half the Rouse time in order to match the terminal time in *G*(*t*) for a system of dangling chains. The dimensionless unit for space was chosen to be the mean square length of a Rouse chain in *x*, *y*, and *z*. Both dimensionless units were chosen in this way because a majority of this work will be studying stress relaxation. The non-dimensional equations have the form

where the Weissenberg number (*Wi*) is the dimensionless rate of shear ($\gamma \u0307\tau R$) or extension ($\epsilon \u0307\tau R$) in non-linear start-up flows. The non-dimensional equations for looped, dangling, and bridged species are

The different equations are coupled via first order association and dissociation reactions where $ka*$ is the rate of association from dangling to bridged or looped conformations, and $kd*$ is the rate of dissociation from bridged or looped to dangling conformations. These rates can be non-dimensionalized by the dangling chain relaxation time such that $ka*=ka/\tau R$ and $kd*=kd/\tau R$, and these relative rates can be related by the equilibrium constant *K*_{eq}, which is defined as the dimensionless rate of dissociation divided by the rate of association (=*k*_{d}/*k*_{a}). To model bonds with many states of partial association, it is necessary to determine relationships for the kinetic rates of attachment and detachment for the individual bonds relative to the total rates.

The dimensionless intermediate rates of association and dissociation can be written as $ka\u2032$ and $kd\u2032$. The relation between $ka\u2032$ and *k*_{a} and between $kd\u2032$ and *k*_{d} can be determined using the reaction kinetics. A detailed derivation for this relationship can be found in Sec. 1 of the supplementary material. Assuming pseudo-first order reaction kinetics and a fixed total rate of association and dissociation, the relationship between the total rate of dissociation and the intermediate rate of dissociation becomes

Similarly, for association,

Assuming that the equilibrium between each intermediate step ($Keq\u2032$) is constant, *K*_{eq} and $Keq\u2032$ can simply be related as

where both the equilibrium constant and the equilibrium between each intermediate step are defined as the rate of dissociation over the rate of association,

Equations (6) and (7) can be plugged into Eq. (9), showing that they are all consistent. It is important to note that the relationships in Eqs. (6)–(8) represent the equivalent kinetics for a system with the same total bond free energy, regardless of the value of *n*. Another reasonable assumption would be to linearly increase the total bond free energy with increasing *n* (i.e., *E*_{total} = *nE*_{bond state} instead of *E*_{total} = *E*_{bond state}/*n*). This is discussed briefly in Sec. 10 of the supplementary material.

The detachment of bridged chains to form either dangling chains or bridged chains with one less bond is mediated by Bell’s law for the force-activated detachment of bonds. It is also necessary to adjust this Bell’s law term in order to reflect the number of intermediate bonds. The extension-dependent Bell’s law term for chain-end dissociation can be easily modified as a function of the number of intermediate bonds, the force on the chain end, the length scale for binding, and the thermal energy.^{42} When there are *n* intermediate bonds, the rate constant for dissociation from a bridged chain with *k* bonds to a less associated bridged chain with *k* − 1 bonds (*R*_{Bk→Bk −1}) or the reaction rate constant for dissociation from a bridged chain with *k* = 1 bond to a dangling chain (*R*_{B1→D}) can be expressed as

Here, *β* is the dimensionless scale for binding.

By imparting a specific association size to the chain end, it is necessary to adjust the value of *N* used to determine the finitely extensible nonlinear elastic (FENE) force on a particular conformation. The adjusted chain length is instead determined using *N*_{adj} = *N* − *nx*, where n is the number of stickers and *x* is the number of Kuhn segments per sticker (here *x* = 1, so minimal stored length is released by bond dissociation) such that the resultant value for the chain’s finite extension is

Similarly, the rate of association for dangling chains is multiplied by an additional term to make the associated conformation (looped or bridged) location-dependent. However, this probability of associating from a dangling chain to a looped or bridged chain does not need to be adjusted for the number of bonds. Both rates are governed by a Gaussian probability of forming a looped chain, expressed in dimensionless units as

Thus, the rate constant for forming a looped chain can be expressed as

and the rate constant for forming a bridged chain can be expressed as

For purposes of this work, *N* = 1000, α = 0.25, and β = 0.35. These values were chosen to stay consistent with previous work.^{34} The value for α represents the range of interaction between two chain ends such that the dimensional version of α is $\alpha *\u2009=\u20090.5Nb22\u2009=\u20090.25Nb2$. The value for β corresponds to the dimensionless characteristic bond length corresponding to a fraction of *R*_{g}, approximately corresponding to a larger associating domain, all of which form aggregates with characteristic dimensions ranging from 2 to 40 nm.^{43–46}

### Numerical integration

The specific numerical integration used depended upon the state of the chain end. For bridged chains, integration was performed using a first-order explicit Euler algorithm. This gives evolution equations for *r* in *x*, *y*, and *z* in shear as

and in extension as

For dangling chains, a semi-implicit predictor-corrector algorithm as described by Öttinger^{47} or Herrchen was used.^{48} The algorithm is semi-implicit because only the spring-force law is treated implicitly. For the predictor step, the equation for the dangling chain is as follows:

The corrector step using the same random number (Δ*W*) is

The direction of *r*(*t* + Δ*t*) is given by the direction of the right-hand-side of Eq. (18), and the length is determined by solving the cubic equation that results from substituting $r\xaft+\Delta t$ calculated in Eq. (17) into (18) and taking the absolute value of the equation $r2$. Equation (18) has a unique solution between 0 and the maximum extensibility of the chain (*Nb*).

These stochastic and deterministic differential equations were coupled with a stochastic method to implement the attachment/detachment processes. In order to minimize the time per simulation without losing important physics, two methods were used: the Gillespie algorithm, which is a standard stochastic algorithm for chemical kinetic systems, and the more efficient tau-leap method.^{49} In both the Gillespie and tau-leap algorithms, the total number of reactions between different populations of chemical species is calculated using the propensity function *a*_{j} defined such that *a*_{j}Δ*t* is the probability that reaction *j* will occur during the next time interval Δ*t* for a particular state. All reactions are assumed to be unimolecular. The propensity is calculated by multiplying the number of molecules in a particular state by their conformation dependent rates of dissociation or association. The total propensity is the sum of propensities of all possible reaction channels [i.e., for looped to dangling, the total propensity is *a*_{tot,L →D} = *n*_{L}*k*_{d}, and for dangling to looped, the total propensity would be *a*_{tot,D→L} = $\u2211i\u2208D$*k*_{a}*f*(*r*_{i})]. It is necessary to calculate the propensities individually for each molecule as a function of its particular conformation. Specifics pertaining to the tau-leap algorithm are discussed in greater detail in the supplementary material.

For this work, the number of reactions per time step along each reaction pathway is a random number drawn from a Poisson distribution with a mean equal to the propensity multiplied by the time step. However, care had to be taken so that the total number of calculated reactions does not exceed the number of available molecules. This was done by choosing the time step for each simulation on a case-by-case basis through analysis of the number of reactions occurring during short simulations. The use of a Poisson distribution instead of a Gaussian distribution to generate the random number makes it unnecessary for the statistics to require a minimum number of reactions per pathway.

### Special case: Looped to dangling

When a chain end goes from the looped to the dangling state, it is necessary to not just switch the identity of the chain but also adjust the spatial location of the chain end. Since looped chain ends are assumed to be at *r* = 0, simply switching the identity to dangling will affect the physics of the system and result in too many chains returning to the looped state [*f*(*r* = 0) = 1]. Instead, the chain end is re-distributed about *r* = 0 in order to obey the detailed balance between the looped and dangling conformations. Thus, the molecules that undergo the reaction from a looped to a dangling chain are recreated with a random conformation drawn from the probability distribution of chains that undergo the reaction from dangling to looped chains. This probability distribution can be calculated as the average of *f*(*r*) integrated over the equilibrium distribution of Gaussian dangling chains. This results in a Gaussian with a mean at zero and a variance of 3*α*/(2 + 3*α*).

## RESULTS AND DISCUSSION

### Effect of bond number and reaction kinetics on quiescent network relaxation

The primary determinant of the network relaxation is the rate of dissociation and re-association of network junctions. For a physical network, the average time a chain spends detached is always less than the time a chain spends attached. This can be expressed using the constants described above as

where *τ*_{a}, the time between association events, represents the time a chain spends detached and *τ*_{d}, the time between dissociation events, represents the time a chain spends attached. Both times are a function of the intermediate rates of association and dissociation. Network relaxation is dependent upon the presence of dissociated dangling chains and has been shown experimentally for triblock systems to be dependent upon the time scale for pullout of one chain end from an associative junction.^{50} This suggests that network relaxation is dependent upon *τ*_{d} as well as the fraction of detached chains, where this fraction is dependent upon the value of *K*_{eq}.

When a chain is allowed to relax fully, the time scale for that relaxation is moderated by the terminal relaxation time for dangling Rouse chains (*τ*_{R}) assumed to be fixed at one end. For a coarse-grained dumbbell of length *N*, this time corresponds to

This terminal relaxation time was used to non-dimensionalize the various kinetic constants. As such, the value of *τ*_{a} (and by definition $ka\u2032$) relative to 1 can be used to divide the sample space into two possible cases for the rate of association:

$\tau a\u22651ka\u2032\u22641$—the chain is able to fully relax while it is dissociated.

$\tau a<1ka\u2032>1$—the chain is not able to fully relax while it is dissociated.

These definitions use the intermediate rates of association, which ultimately define the lifetime of a dangling chain. For context, polystyrene systems exhibit Rouse times on the order of 10^{−5} to 10^{−1} s for *N* = 10–1000.^{50–52} In order to get case I behavior where the chains are able to fully relax before re-associating, $ka\u2032$ would need to be less than τ_{R}^{−1}. The rules demarcating the boundary between cases I and II can be re-written in terms of the total rate of association and the total equilibrium constant in order to better compare systems with multiple bonds such that the boundary between case I and case II is

where *n* represents the number of bonds per association, and *k*_{a} and *K*_{eq} reflect the total rate of association and the equilibrium constant for an equivalent network with one bond ($ka\u2032=ka$ when *n* = 1). For the sake of system comparison, all data will be presented in terms of the net equilibrium constants and the total rates of association. As the number of bonds increases, the region where network chains fully relax prior to undergoing re-association shifts to lower values of *k*_{a}.

#### Case I: Chains are able to relax fully before reassociating

When chains are detached from the network in a system with case I kinetics, they are able to fully relax prior to reattachment. However, in order for the network itself to fully relax, it is necessary for all non-relaxed chains to first form dangling chains. Thus, the total network relaxation should consist of two relaxation processes: the fast relaxation of dangling chains, followed by the slower time scale for chain end dissociation and subsequent relaxation,

where $\varphi i$ represents the fraction of species *i* and *G*(*t*) is a function of both the terminal Rouse relaxation time of the dangling chains (here set to be 1) and 〈*τ*_{d}〉_{Case I}, which is the average time between dissociation events from bridged to dangling chains. It is assumed that looped chains do not contribute to network relaxation except to decrease the fraction of bridged chains at equilibrium, which was also the case in previous work.^{34} This decrease in stress due to the presence of looped chains has been shown to be qualitatively accurate via controlled experimental techniques and theory.^{31,53,54}

For a single bond, the average time between dissociation events in the second exponential in Eq. (22) can be estimated using the imposed kinetic constants to be related to the intermediate rate of dissociation,

Intuitively, increasing the rate of association decreases the terminal relaxation time of the material at constant *K*_{eq}. Similarly, increasing the equilibrium constant to favor the formation of dangling chains decreases the impact of 〈*τ*_{d}〉_{Case I} on the total relaxation of the network by decreasing the value of 〈*τ*_{d}〉_{Case I} towards zero such that the terminal relaxation time approaches 1, as evidenced in Fig. 1.

Graphs of the resultant relaxation spectra from Brownian dynamics simulations show the presence of two relaxation modes (Fig. 1) for values of the kinetic constants corresponding to case I. As the rate of dissociation is increased relative to association, the terminal relaxation time approaches that of a system of dangling chains. In contrast, 〈*τ*_{d}〉_{Case I} $\u226b$ *τ*_{R} results in a relaxation spectrum where the shorter relaxation time is barely visible. This occurs in part because of the disparity between *τ*_{R} and 〈*τ*_{d}〉_{Case I}, where the average time for chain end dissociation is significantly longer than the dangling chain relaxation time. Additionally, decreasing the value of *K*_{eq}, which favors association of dangling chains back into the network, results in a lower total fraction of dangling chains and decreases the magnitude of the first exponential term in Eq. (22). This decrease is evident at low values of dimensionless time in Fig. 1.

It is possible to perform double exponential fits on these relaxation data. However, as the equilibrium constant decreases and 〈*τ*_{d}〉_{Case I} increases, fits to the faster relaxation time become harder for the fitting algorithm to resolve. This behavior is evident in Fig. 1, where the smaller relaxation is barely visible. In order to address this issue, the data were instead fit to a double exponential where one of the terms was fixed with a relaxation time of 1. The resultant time constants from these fits to the relaxation data can be compared to a modified form of Eq. (23), resulting in a master curve for the average dissociation time expressed in terms of the reduced dissociation time (*τ*_{d}*k*_{a}) [Fig. 2(a)],

Equation (24) on its own overestimates the relaxation time. This overestimation occurs due to an oversimplification of 〈*τ*_{d}〉 given that the rate constant for dissociation of bridged chains is not solely dependent upon $kd\u2032$. Instead, the dissociation is also a function of the force of retraction acting on the bridged chain through Bell’s law. Thus, the reduced dissociation time can better be approximated as

where *B*_{avg} represents the average value for Bell’s law term, discussed in detail in Sec. 5 of the supplementary material. This results in a much better fit to the simulation results [Fig. 2(a)]. The non-reduced relaxation time should approach 1 with increasing *K*_{eq}. However, the reduced dissociation time deviates from the 1/*K*_{eq} scaling at higher values of the equilibrium constant.

When more bonds are added, it is necessary to modify the equation to account for the inclusion of multiple bridged states. Instead of dissociating to form dangling chains, bridged chains may either associate to have *k* + 1 bonds or dissociate to have *k* − 1 bonds. If *k* − 1 = 0, the bridged chain completely detaches from the network to form a dangling chain. For the simplest multi-bond case of *n* = 2, bridged chain ends can take on two states: *B*1 and *B*2, where *B*2 corresponds to a fully associated chain end (it cannot associate any further) composed of 2 bonds, and *B*1 is a partially associated state with only 1 bond. *B*1 is the only bridged species that can become a dangling chain. Generalizing to *n* bonds, states 1 through *n* – 1 are only partially associated and state *n* is fully associated.

The addition of these bonds requires that the three possibilities of bridged chain conformation—*B*1, *B*2:*n* − 1, and Bn—are treated separately. In order to detach, a chain may follow many possible kinetic pathways. This results in a tri-diagonal system of equations that can be used to solve for the relaxation time of chains in all states (*k* = 1 to *n*). A detailed derivation of these equations is included in Sec. 5 of the supplementary material. From these equations, it is possible to compute the average time required to dissociate the chain from state *k* to a dangling chain based on the probability that a chain will be present with *k* bonds (*p*_{k}),

When Eq. (26) is solved for multiple values of *n*, the relaxation time increases with increasing bonds to *n* = 3 and then decreases with subsequent bond addition [Fig. 2(d)]. Comparison with data from fits to simulations for n = 2 and n = 5 shows that this equation captures the trend in reduced dissociation time [Figs. 2(b) and 2(c)]. However, the fits do not align as well with multiple stickers as they do with single stickers, likely due to the approximation that *B*_{avg} can be expressed as a constant value that was determined using a Gaussian probability distribution (see Appendix A).

Case I should be challenging to observe experimentally because bond reattachment must be extremely slow to enable full chain relaxation while still having a sufficient equilibrium constant to form a gel. Consequently, this behavior is rarely observed experimentally. The presence of a second relaxation time in physical gels has been seen in diblock copolymer systems, with the slow process being attributed to the Stokes-Einstein diffusion of the micelles and the fast process corresponding to the relaxation of coronal chains.^{55,56} However, telechelic systems typically only exhibit multiple visible relaxation times if there are two or more different association time scales^{30,57} or if the relaxation is measured following a high-strain deformation.^{58} Peters and Lodge used time-resolved small-angle neutron scattering to measure the chain end pullout time of poly[styrene-b-(ethylene-*alt*-propylene)-b-styrene].^{50} The resultant data (*τ*_{R} = 0.14 s, *K*_{eq} ≈ *αχ* = 0.026, and *τ*_{pullout} ≈ *τ*_{d} = 4054 s) could be used to estimate a dimensionless association rate on the order of 10^{−3}. However, Peters and Lodge attributed the presence of additional relaxation times to either the formation of secondary structural ordering or the congestion of a disordered micellar system.^{50}

#### Case II: Chains are not able to relax fully before re-associating

When dangling chains are not able to fully relax before re-attaching to the network, the time scale for relaxation does not separate into the two clear time scales as it does in case I. For every time period *τ*_{d} + *τ*_{a} that passes, the dangling chains are only able to relax during time they spend detached (the time it takes to reattach, *τ*_{a}). Thus, the dangling chain relaxation time *τ*_{R} must be re-normalized by *τ*_{a}/(*τ*_{d} + *τ*_{a}),

This results in a single exponential where the relaxation time is a function of the ratio in Eq. (27), and the coefficient should depend on both the fraction of dangling chains ($\varphi D$) and the fraction of bridged chains ($\varphi B$),

Graphs of the relaxation modulus versus time verify the presence of only a single exponential for values of the kinetic constants corresponding to case II (Fig. 3). The decrease in *G(0)* with decreasing *K*_{eq} is a function of the increased presence of looped chains as *K*_{eq} decreases.

For a single bond, the relaxation time of this single exponential can be approximated using the relaxation time in Eq. (28) as a function of the kinetic constants,

where *f*_{R} is the average value for the Gaussian probability term for bridged chain formation (*f*_{R} = 〈1 − *f*(*r*)〉).

The derivation of the average Gaussian probability term is discussed in detail in Appendix A. While this equation indicates that the average relaxation time when *n* = 1 should not be dependent upon the value of *k*_{a}, the terminal time actually decreases to an asymptotic value with increasing rate of association (Fig. S9-1 of the supplementary material). The prediction overestimates the value from simulation because the prediction assumes that the time to attach from *Bk* to *Bk* + 1 or detach from *Bk* to *Bk* − 1 is constant. However, in reality both the time to attach and the time to detach are dependent upon the location of the chain end.

When Eq. (29) is modified to account for multiple bonds, it is instead necessary to calculate the time it takes for a chain to detach following its initial attachment. The precise calculations are discussed in greater detail in Appendix B, and the resultant equation for the relaxation time becomes

Unlike in case I where increasing the number of bonds resulted in a non-monotonic relationship with the average dissociation time, increasing the number of bonds in case II results in an increase in the terminal time (Figs. 4 and 14). This prediction is surprising given that the overall bond energy is conserved. Experiments often show an increase in the relaxation time with increasing numbers of bonds, but each bond added increases the total *K*_{eq} and decreases the overall *k*_{d}.^{22,50,59} As with a single bond, the predicted terminal time for multiple bonds predicts general trends but does not quantitatively match the data [Fig. 4(b)].

In order to compare to literature studies, it is first necessary to provide context for the specific sorts of materials that exhibit case II kinetics. Associating domains with a single bond state that have equivalent kinetics to case II have been studied in detail by Rapp *et al.* and Craig *et al.* Rapp *et al.* used fluorescence recovery following photobleaching to extract kinetic constants for a variety of telechelic coiled-coil proteins with varying modifications to leucine zipper hydrophobicity. Here, the dissociation rate was assumed to be on the same time scale as the dangling chain relaxation time (k_{d} ∼ 1) with equilibria ranging from 0.04 to 0.5.^{60} Craig *et al.* used small molecule chemistry combined with NMR exchange spectroscopy to determine the association equilibria for specific metal-ligand coordination crosslinkers.^{61,62}

Systems with the sorts of associating domains discussed above exhibit relaxation times that are on the order of those predicted for case II. Feldman *et al.* used quadruple hydrogen-bonding side-chains to study the influence of multiple bond-formers within a single endblock associating group. The presence of these additional bonds at the end of telechelic polymer chains resulted in enhanced mechanical properties and increased network relaxation times.^{63} This idea of multiple bonds enhancing mechanical properties has been applied to supramolecular polymers.^{64} In systems where associations occur via only a single bond state, Shen *et al.* demonstrated a correlation between decreasing pH and increasing the relaxation time for telechelic gels with associating leucine zippers. Decreasing the pH stabilized the leucine zipper aggregates and is synonymous with decreasing *K*_{eq}.^{37} Yesilyurt *et al.* used covalent complexation to illustrate the differences between fast and slow association processes and exhibited dimensionless relaxation times of order 1.^{65} In both cases I and II, the resultant relaxation times—both from the predictive equations and from equilibrium simulations—are within ranges that are seen in experimental studies of triblock copolymers,^{50,66,67} telechelic coiled-coil protein gels,^{37} and metal-coordinate gels.^{57,68}

### Effect of multiple bond states on start-up of flow in shear

The addition of flow adds yet another relative time scale for network rearrangement. For purposes of simplicity, start-up of shear and extension was performed on only a small subset of kinetic constants corresponding to cases I and II at equilibrium—for *n* = 1, 2, and 5 bond states. This combination of kinetic constants illustrates the differences between case I, where dangling chains have time to relax prior to reassociation, and case II, where dangling chains can only relax a fraction of their stress prior to reassociation. Previous work has shown that in systems with a single bond, non-monotonicities develop in the steady-state flow curves of systems with slow association kinetics where only a small fraction of dangling chains are present at equilibrium.^{34} It is important to note that the definition of *Wi* used here is relative to the relaxation time for dangling chains because it is the only time scale that stays constant across all systems. Because of this definition, the crossover from linear to nonlinear behavior is a function of both *Wi* and the terminal relaxation time of the full system for both shear and extensional flows (Fig. S8-6 of the supplementary material).

For the start-up of steady shear, increasing shear rates result in the development of a stress overshoot followed by a decrease to the steady-state value of the stress (Fig. 5) for both case I and case II kinetics. Following the convention in the field, the stress was calculated using only the contribution from the internal energy due to chain extension.^{29,35,69,70} Increasing the number of bond states increases the steady-state stress at shear rates of *Wi* ≤ 1, while for increasing shear rate beyond *Wi* = 1, the networks asymptotically approach a bond number-independent limit (Fig. 6). Thus, at high enough shear rates, these multi-bond systems behave more like a system with a single bond. This is consistent with the intuitive notion that when the deformation rate is fast relative to the bond dynamics, the multiple bond does not dynamically exchange between bonding states in a way that yields different properties; it simply breaks as would a simple bond. Creton *et al.*^{71} and Craig *et al.*^{61,62,72} have seen similar rate-dependent responses in experimental systems.

Intuitively, the application of shear should result in an increased fraction of dangling chains, while the fractions of both bridged and looped chains decrease. However, differences in both the kinetic constants and shear rates cause deviations from this behavior. In some instances, the values of the kinetic constants favor increases in bridged and looped chain fractions.

At low shear rates (*Wi* < 1), increasing the number of bonds does not qualitatively change the behavior of the individual chain fractions (Fig. S7-3 of the supplementary material). However, when the dangling chains have time to fully relax before undergoing an attachment event (case I), the fraction of looped chains in the system increases at higher strains [Fig. 7(a)]. This occurs because when bridged chains undergo their extension-dependent detachment at strains greater than 1, they form dangling chains that are able to fully relax and form looped chains. Indications of this increase in looped chain formation have been seen experimentally for a system that could correspond to case I behavior.^{73} Yekta *et al.* posited that the lack of signal change in fluorescence-probe experiments on aqueous systems of hydrophobically modified ethoxylate-urethane associating polymers under shear is a function of the formation of looped chains from bridged chains.^{73}

At high shear rates (*Wi* > 1), dangling chains are pulled faster than they are able to relax [Fig. 7(b)]. Looped and bridged chains dissociate at strains greater than 1 to become dangling chains. The fraction of looped chains decreases for the duration of the strain, while the fraction of dangling chains slowly increases. However, following the initial decay, the fraction of bridged chains slowly increases to steady state as dangling chains re-form bridged chains.

For a system with faster kinetics corresponding to case II, low shear rates (*Wi* < 1) result in decreases in the looped chain fraction to a steady-state value that is lower than the equilibrium value in the absence of shear [Fig. 7(c)]. In this case, dangling and bridged chains rapidly exchange at steady state with minimal to no decrease in stress as a result of dangling chain relaxation. The increase in average chain extension promotes a net drop in the number of loops. At high shear rates (*Wi* > 1), dangling chains are pulled faster than they are able to relax, and dangling chains rapidly form as bridged chains are broken [Fig. 7(d)]. This results in an overshoot in dangling chain fraction shortly after a decrease in bridged chain fraction. This overshoot corresponds to the overshoot in shear stress. Following the initial overshoot, dangling chains and bridged chains interconvert, causing the number of bridged chains to increase above its initial value. This occurs because looped chains dissociate, and the imposed shear decreases the chances that the loops will reform.

While increasing the number of bonds in the system does not qualitatively change the behavior of chain fractions at low and high shear rates, an increased number of associating bonds result in an increase in both the yield stress and the steady-state value of the stress (Figs. 5 and 6). This increase agrees well with experiments on triblock copolymers, where in addition to increasing the relaxation time, increasing the length of associating groups is known to result in higher shear moduli.^{50,67} However, in these experiments, the effective level of *K*_{eq} is not kept constant. Therefore, these simulation results show something more striking: the fact of having multiple bonds increases viscosity even if the bond strength is not increased. This increase is due to a higher average degree of chain extension in the bridged and dangling chains as the number of bond states increases for both cases (Fig. 8) for smaller values of *n* such that the chain end is able to fully dissociate. When dangling chains are able to relax prior to reassociation into the network, the value of $r2/31/2$ for dangling chains stays near 1 until *Wi* ≥ 1 [Fig. 8(a)]. This results in lower steady-state extensions for both dangling and bridged chains relative to the system with faster association kinetics [Fig. 8(c)]. When dangling chains are not able to fully relax prior to association, the value of $r2/31/2$ for dangling chains is higher than bridged chains. In previous work, this was attributed to the ability of chains to tumble in higher shear flows.^{34} So, while increasing the number of bonds decreases the fraction of elastically effective chains, it results in higher extensions that ultimately increase shear stress.

These differences in chain extension with increasing *n* correspond to differences between the internal energy of the network and the work required to deform it. For a system of permanently linked bridged chains, the work and internal energy increase together such that all work done on the system goes directly toward increasing the internal energy. However, for a system with reversible cross-links, the work and internal energy do not directly overlap. At high strains, the work put into the system to deform the chains is greater than the internal energy stored in the chain ends. The energy dissipated by the system is a function of the difference between the work and internal energy, as well as the energy dissipated during bond breaking. The use of additional bonds to increase the amount of energy dissipated and increase material toughness has been frequently used in experimental work.^{24,27,71} For systems with sacrificial bonds that do not reform, the energy stored in the bond translates to energy dissipation.^{27} However, when the bonds can reform over the time scale of the experiment, the energy dissipated by the breakage of additional bonds becomes dependent on both the rate of imposed flow and the association kinetics.^{71} When an individual bond in state *k* breaks, it breaks due to thermal dissociation under an applied tension. Therefore, the energy dissipated by the bond itself is only the mechanical energy which was applied, not the total bond energy or activation energy which was overcome by a combination of thermal and mechanical energies. Therefore, the energy dissipated is modeled as the energy stored in the bond per Bell’s law term for a system containing *n* bonds where the bond is in state *k*,

If the bond reforms before the chain is able to relax, the internal energy of the bond is effectively conserved, and the reformed bond is under the same tension as previously. Therefore, no net energy is dissipated. However, if the chain is able to relax a sufficient fraction of its stress before reassociating, the force on the chain and thus the internal energy stored in the bond decrease, and the difference contributes to dissipation. Alternatively, if the bond does not reform and the chain remains in the dangling state for the duration of the flow, the entirety of the mechanical energy stored in bond deformation is dissipated.

The total dissipated energy during steady-state shear flow increases with increasing number of bonds for both cases over the range of *Wi* studied (0.0033 ≤ *Wi* ≤ 3.3) (Figs. 9 and 10). In the linear elastic region, when the amount of work put into the system goes directly toward increasing the internal energy of the network, the total dissipated energy is near zero. However, around 100% strain the work put into the system all does not translate to internal energy as the force on the extended chains is high enough to increase the rate of dissociation. During the transition from linear elastic to steady-state behavior, the systems with only a single bond state dissipate more energy. However, increasing the number of bond states from *n* = 1 to *n* = 5 results in higher dissipation at steady state. The degree of difference from *n* = 1 to *n* = 5 is highest at lower shear rates and for systems with slower reaction kinetics.

At shear rates greater than the dangling chain relaxation time, the amount of total dissipation becomes less dependent upon both the relative kinetics and the number of bond states. When the kinetics are slow enough to allow for dangling chains to fully relax as in case I, increasing the shear rate to *Wi* ≥ 1 decreases the differences in total dissipation with increasing bond number (Fig. 10). For fast kinetics, the difference in total dissipation with increasing bond number stays relatively constant regardless of the applied shear rate as dangling chains are only able to relax a small fraction of the stress in the system. In all cases, the total energy dissipated by bond breaking is small enough relative to the total dissipated energy to make the stress calculations used to determine the data in Figs. 7 and 8 a reasonably accurate representation of the total stress in the system. Using both the internal energy from chain extension and the internal energy from bond breaking would not affect the shape of the curves.

### Effect of multiple bond states on start-up of flow in extension

The factors that result in differences between shear stress during start-up for the two cases also affect the resultant network behavior in extensional flows, where increasing numbers of bonds correlate to increased extensional viscosities. The extensional viscosities were calculated following literature conventions, only taking into account the internal energy of the chains.^{29} However, the application of strain rates greater than 0.5 results in breakdown of the network regardless of the system kinetics and number of bonds shortly after a strain of 1 is reached [Figs. 11(c) and 11(d) and Fig. S8-1 of the supplementary material]. This breakdown occurs as a function of the strain on the system such that the extensional viscosity diverges at a strain of 100% for all dimensionless strain rates greater than *Wi* = 0.5; it can be expected that the network will break up prior to the divergence in viscosity. Thus, pulling these telechelic systems with strain rates greater than their characteristic relaxation times will result in much the same behavior regardless of the properties of the associating groups.

At shear rates below *Wi* = 0.5, the system that falls into case I at equilibrium undergoes an overshoot in the extensional viscosity before reaching a steady-state value, while the case II system does not undergo an overshoot. When the kinetics are slow, the behavior is similar to steady shear. The presence of a stress overshoot is partially due to the finite extensibility of the chains, as well as the enhanced dissociation due to Bell’s law term.^{74} In case I, dangling chains are able to fully relax before reassociating, which relaxes the stress [Fig. 11(a)]. However, in the system corresponding to case II, the kinetics of chain association are too fast to allow for the dangling chain to fully relax the stress prior to reattachment [Fig. 11(b)].

NOMENCLATURE | ||

τ_{R} | Rouse time | Dangling chain relaxation time |

$ka\u2032$ | Bond association rate | Rate of association for individual bond states non-dimensionalized by τ_{R} |

$kd\u2032$ | Bond dissociation rate | Rate of dissociation for individual bond states non-dimensionalized by τ_{R} |

k_{a} | Total association rate | Rate of association from dangling to bridged or dangling to looped conformations non-dimensionalized by τ_{R} |

k_{d} | Total dissociation rate | Rate of dissociation from bridged to dangling or looped to dangling conformations non-dimensionalized by τ_{R} |

n | Number of bond states | |

K_{eq} | Equilibrium constant | Equilibrium between dissociation and association (k_{d}/k_{a}) |

$Keq\u2032$ | Partial equilibrium constant | Equilibrium between dissociation and association of one bond state ($kd\u2032/ka\u2032$) |

Wi | Weissenberg number | Dimensionless flow rate for shear ($\gamma \u0307\tau R$) and extension ($\epsilon \u0307\tau R$) |

τ_{d} | Dissociation time | Time between dissociation events represents the time a chain spends attached relative to τ_{R} |

τ_{a} | Association time | Time between association events represents the time a chain spends detached relative to τ_{R} |

〈τ_{d}〉_{Case I} | Average dissociation time (case I) | Average time for a chain to dissociate from an aggregate to form a dangling chain (average time for all associated states to become dangling chains) relative to $\tau R$ |

〈τ_{d}〉_{Case II} | Average dissociation time (case II) | Average time between dissociation events relative to τ_{R} |

〈τ〉_{Case II} | System relaxation time | Terminal time for G(t) when data are fitted to a single exponential decay |

R_{x} | Reaction rate constant | Reaction rate (x denotes association or dissociation and the relevant species) |

β | Bond length | Dimensionless scale for binding in units of $Nb2$ |

N | Kuhn segment number | Number of Kuhn segments |

b | Kuhn segment length | Length of each Kuhn segment |

F_{FENE} | FENE force | Finitely extensible nonlinear elastic force of retraction |

α | Interaction distance | Dimensionless range of interaction between two chain ends in units of $Nb2$ |

f(r) | Probability for loop formation | Gaussian probability with width 3α that a dangling chain will form a loop |

a_{j} | Propensity function | |

$\varphi i$ | Chain fraction | Fraction of species of conformation i |

B_{avg} | Bell’s law average | Bell’s law term averaged |

f_{R} | Average probability of forming a bridged chain | |

H | Spring constant | Dimensional spring constant in three dimensions = 3k_{B}T/Nb^{2} |

ζ | Friction coefficient | |

L | Contour length | Here equivalent to Nb |

$\u2207$u | Velocity gradient | Dimensionless |

t | Time | Dimensionless time |

G(t) | Relaxation modulus | |

r_{jk} | End-to-end vector | Vector location of chain end for species j (= L, D, and B) with k stickers attached |

ΔW or dW | Three-dimensional Wiener process | Random variables drawn from a normal distribution with zero mean and a variance dt |

□* | Dimensional variable | Dimensional form of generic variable □ |

k_{B} | Boltzmann’s constant | |

T | Absolute temperature | |

$\gamma \u0307$ | Shear rate | |

$\epsilon \u0307$ | Extensional rate | |

N_{adj} | Adjusted number of Kuhn segments of the chain | The adjusted number is a function of the number of bond states attached to an association |

p_{k} | Probability that a bridged chain has k stickers attached | |

U_{bond,k} | Internal energy stored in the bonds in a system with k bonds attached |

NOMENCLATURE | ||

τ_{R} | Rouse time | Dangling chain relaxation time |

$ka\u2032$ | Bond association rate | Rate of association for individual bond states non-dimensionalized by τ_{R} |

$kd\u2032$ | Bond dissociation rate | Rate of dissociation for individual bond states non-dimensionalized by τ_{R} |

k_{a} | Total association rate | Rate of association from dangling to bridged or dangling to looped conformations non-dimensionalized by τ_{R} |

k_{d} | Total dissociation rate | Rate of dissociation from bridged to dangling or looped to dangling conformations non-dimensionalized by τ_{R} |

n | Number of bond states | |

K_{eq} | Equilibrium constant | Equilibrium between dissociation and association (k_{d}/k_{a}) |

$Keq\u2032$ | Partial equilibrium constant | Equilibrium between dissociation and association of one bond state ($kd\u2032/ka\u2032$) |

Wi | Weissenberg number | Dimensionless flow rate for shear ($\gamma \u0307\tau R$) and extension ($\epsilon \u0307\tau R$) |

τ_{d} | Dissociation time | Time between dissociation events represents the time a chain spends attached relative to τ_{R} |

τ_{a} | Association time | Time between association events represents the time a chain spends detached relative to τ_{R} |

〈τ_{d}〉_{Case I} | Average dissociation time (case I) | Average time for a chain to dissociate from an aggregate to form a dangling chain (average time for all associated states to become dangling chains) relative to $\tau R$ |

〈τ_{d}〉_{Case II} | Average dissociation time (case II) | Average time between dissociation events relative to τ_{R} |

〈τ〉_{Case II} | System relaxation time | Terminal time for G(t) when data are fitted to a single exponential decay |

R_{x} | Reaction rate constant | Reaction rate (x denotes association or dissociation and the relevant species) |

β | Bond length | Dimensionless scale for binding in units of $Nb2$ |

N | Kuhn segment number | Number of Kuhn segments |

b | Kuhn segment length | Length of each Kuhn segment |

F_{FENE} | FENE force | Finitely extensible nonlinear elastic force of retraction |

α | Interaction distance | Dimensionless range of interaction between two chain ends in units of $Nb2$ |

f(r) | Probability for loop formation | Gaussian probability with width 3α that a dangling chain will form a loop |

a_{j} | Propensity function | |

$\varphi i$ | Chain fraction | Fraction of species of conformation i |

B_{avg} | Bell’s law average | Bell’s law term averaged |

f_{R} | Average probability of forming a bridged chain | |

H | Spring constant | Dimensional spring constant in three dimensions = 3k_{B}T/Nb^{2} |

ζ | Friction coefficient | |

L | Contour length | Here equivalent to Nb |

$\u2207$u | Velocity gradient | Dimensionless |

t | Time | Dimensionless time |

G(t) | Relaxation modulus | |

r_{jk} | End-to-end vector | Vector location of chain end for species j (= L, D, and B) with k stickers attached |

ΔW or dW | Three-dimensional Wiener process | Random variables drawn from a normal distribution with zero mean and a variance dt |

□* | Dimensional variable | Dimensional form of generic variable □ |

k_{B} | Boltzmann’s constant | |

T | Absolute temperature | |

$\gamma \u0307$ | Shear rate | |

$\epsilon \u0307$ | Extensional rate | |

N_{adj} | Adjusted number of Kuhn segments of the chain | The adjusted number is a function of the number of bond states attached to an association |

p_{k} | Probability that a bridged chain has k stickers attached | |

U_{bond,k} | Internal energy stored in the bonds in a system with k bonds attached |

Instead, the extensional viscosity for a system corresponding to case II kinetics slowly approaches steady state at higher strains. This retardation in the increase in extensional viscosity is a function of the generation of unextended dangling chains as looped chains dissociate to form dangling chains that are slowly pulled away from their equilibrium conformation where they can no longer form looped chains. Instead, the dangling and bridged chains interchange—stretching when in the bridged state and relaxing a fraction of the stress in the dangling state. This alternation while looped chains dissociate results in a slow increase in both dangling and bridged chains with increasing strain (Fig. S8-3 of the supplementary material). When the dissipation is calculated for these systems in extension, the maximum dissipation for both fast and slow kinetics occurs at strains just below *Wi* = 0.5. Increasing the number of bonds increases the amount of dissipation in both cases. However, much like in shear, the presence of multiple bonds has a greater impact on dissipation when the system kinetics are slower.

At steady state, the two cases exhibit very different behavior. When the kinetics are slow, the steady-state extensional viscosity continuously decreases with increasing strain rate for *n* = 1 and 2 but increases for *n* = 5 [Fig. 12(a)]. When the kinetics are fast, the steady-state extensional viscosity non-monotonically increases by orders of magnitude and then decreases with increasing strain rate [Fig. 12(b)].

These data agree with other theoretical predictions for extensional flows,^{35,75,76} including a paper by Marrucci *et al.* that compares the results of a theoretical construction involving finitely extensible dumbbells to polystyrene in both the semi-dilute solution and the melt state. In the melt state, the extensional viscosity decreases with increasing strain rate. In semi-dilute solutions, the extensional viscosity decreases to *Wi* = 1 before increasing with increasing strain rate.^{77} While this polymer system is not a network composed of a telechelic associating system, the causes for this behavior—slower kinetics in the melt state and the importance of finite extensibility in the semi-dilute state—could be used to help describe the behaviors seen from this theoretical prediction.

There is not a large amount of tensile data for purely telechelic molecules in extensional flows. However, a number of papers focusing on commercially available associative thickeners have discussed mechanical behavior in both shear and extension. Increasing strain rates on these materials results in increasing steady extensional viscosity values, agreeing qualitatively with the increase seen at steady state for systems with fast kinetics.^{35,78,79} Hotta *et al.* showed extensional data at low strain rates for triblock poly(styrene-b-isoprene-b-styrene) that match qualitatively with equivalently low *Wi* (Fig. S7-5 of the supplementary material).^{80}

## CONCLUSIONS

Using a simple reactive form of the Smoluchowski equation for convection and diffusion, it is possible to provide insight into the relaxation and flow behavior of telechelic associating systems using Brownian dynamics simulations. The ability to add in more than a single association event has enabled the expansion of this theory to describe a wider variety of systems where association and dissociation might occur in stages. The range of parameter space can be divided into two different dynamic regimes. In case I, dissociated chain ends can fully relax prior to re-association into the network. This ability results in two characteristic time scales: one corresponding to dangling chain relaxation and the other corresponding to the total time it takes for a chain end to dissociate from its junction. In case II, dissociated chain ends can only relax a fraction of the way toward equilibrium before re-associating. This results in only a single characteristic time scale that is a function of the fraction of time the chain spends detached relative to the total time the chain spends attaching and detaching. The addition of multiple bonds for case I where dangling chains are able to fully relax before reassociating resulted in a non-monotonic decrease with increasing number of bonds such that the relaxation time was a maximum for *n* = 3 bonds. However, for case II where dangling chains are able to relax a fraction of the way before reassociating, the addition of multiple bond states resulted in longer relaxation times with increasing *n*.

The addition of multiple bond states resulted in a different behavior in the start-up of shear and extension, especially at low shear and strain rates when the kinetics of association were slow. At low *Wi*, the presence of multiple bonds results in increased dissipation that is maximized when the kinetics of association are slower such that dangling chains are able to relax prior to reassociation. Shear and strain rates greater than the dangling chain relaxation time prevent dangling chains from relaxing, and thus the addition of multiple bond states makes less of an impact on network dissipation during flow. These results indicate that for applications where toughness in extensional flows is important, it is best if the material is designed to have a characteristic relaxation time that will result in low shear and strain rates relative to the characteristic relaxation time of the material itself. If this can be done, the addition of multiple bonds in the associating endblock can act as an additional way to dissipate stress.

## SUPPLEMENTARY MATERIAL

See supplementary material for detailed explanation of the methods used in this publication, as well as additional data.

## ACKNOWLEDGMENTS

This project was supported by the U.S. Army Research Office through the Institute of Soldier Nanotechnologies. J.R. was supported by Real Colegio Complutense at Harvard with additional financial support from Project No. FIS2016-78847-P of the MEC. The authors would like to thank Gareth H. McKinley for valuable discussions at the onset of this project.

### APPENDIX A: AVERAGE RATES FOR ASSOCIATION AND DISSOCIATION

While the simple expressions for predicting the terminal time can be fit to simulation data using a single parameter, it is possible to further refine the equations using average reaction rates. The rate of dissociation of bridging chains is mediated by Bell’s law term, while the rate of association from a dangling chain to a bridged chain is mediated by the Gaussian probability of forming a looped chain. It is possible to estimate the average values of these two terms using the probability distribution functions at equilibrium. Assuming that at equilibrium the chains are approximated by a Gaussian end-to-end distance distribution even though the chains are finitely extensible (N = 1000 is sufficiently large), the probability distribution can be written in dimensionless units as

The rate constant for the association from dangling to bridged chains is

where *α* = 0.25. The average value of 1 − *f*(*r*) can be calculated using the equilibrium probability distribution in Eq. (A1),

Thus, the average rate constant for association from a dangling to a bridged chain can be written as

The rate constant for dissociation from bridged to dangling chains or from a bridged chain with *k* bonds attached to a bridged chain with *k* − 1 bonds attached can similarly be solved as

where *β* = 0.35 and *B*_{avg} represents the average of Bell’s law term.

### APPENDIX B: TERMINAL TIME PREDICTION

The average time it takes for the chain to detach from state *Bk* to form a dangling chain can be calculated for the three possible distinct states: *B*1, *B*2:*n* − 1, and *Bn*. When a bridged chain is only attached by one bond (*B*1), it is possible for it to dissociate to form a dangling chain or further associate to a bridged chain with two bonds. The probabilities that these two events could occur can be expressed by

When a bridged chain is attached by 2 to *n* − 1 bonds, it can either dissociate to form a bridged chain consisting of one less bond or associate to form a more bridged chain attached by an additional bond. Thus, the probabilities for association and dissociation remain the same as they were for one bond

However, when a bridged chain is in state *Bn*, the chain is only able to dissociate to form a bridged chain attached by 1 less bond. Equations (B3) and (B4) become

These probabilities can then be used to weight the total time a bridged chain takes to attach or detach. When a chain is in the *B*2:*n* − 1 state, it can associate from *Bk* to *Bk* + 1. The time scale of this is the time it takes to associate $1/ka\u2032$ plus the time it spends in state *Bk* + 1 (*τ*_{k+1}). Similarly, it can dissociate from *Bk* to *Bk* − 1, the time scale for which is a combination of the time between dissociation events $1/Bavgkd\u2032$ and the time it spends in state *Bk* − 1 (*τ*_{k−1}). Combining these with Eqs. (B3) and (B4) results in the expected relaxation time for a bridged chain with 2 to *n* − 1 bonds,

The time scales for a chain in state *B*1 or *Bn* are very similar to Eq. (B7). However, when a chain dissociates from *B*1, it becomes a dangling chain,

Finally, a chain in state *Bn* can only dissociate from *Bn* to *Bn* − 1. So substituting in Eqs. (B5) and (B6) for the probabilities in Eq. (B7) results in the equation for the total time to dissociate from a bridged chain with *n* bonds to a dangling chain,

Equations (B7)–(B9) can be rewritten as a tri-diagonal system of equations by plugging in the probabilities for association and dissociation,

#### Case I: Dangling chain able to relax prior to reassociation

For this particular case where the relaxation spectrum should result in two relaxation times—one at short times that is a function of the chain relaxation time and the other at long times that is a function of the time it takes for chains to detach from the bridged to the dangling state—it is possible to express the probability that a bridged chain has *k* stickers attached (where the stickers are both distinguishable and only allowed to associate ordinally from one end of the bond to the other),

Given this probability, the average time it takes to detach a chain can be written as

The terminal time for the networks in case I can thus be estimated as the sum of the time it takes to detach a chain to form a dangling chain and the relaxation time of that dangling chain (*τ*_{R} = 1) (Fig. 13).

#### Case II: Dangling chain not able to relax fully before reassociating

The tridiagonal systems of equations used to calculate the average time a bond in each state *k* takes to relax can be used to calculate the total relaxation time for case II. However, unlike in case I where it was necessary to calculate the average time it takes for all of the chains to dissociate from state *k* to a dangling chain, it is instead necessary to calculate the time the chain spends between when it first attaches to form a bridged chain until it fully detaches to re-form a dangling chain. Thus, the time between chain dissociations corresponds to the time for a bridged chain with 1 bond to dissociate (*τ*_{d,1}). From Eq. (B10), this corresponds to

This time can then be substituted to get the multi-sticker analog for predicting case II behavior (Fig. 14),

It should be noted that when the average rates are taken into account, it is also necessary to modify the boundary between cases I and II such that for case I the condition is

This modification shifts the boundaries between cases I and II to even lower values of *k*_{a}.