Valence can control the nonexponential viscoelastic relaxation of multivalent reversible gels

Gels made of telechelic polymers connected by reversible cross-linkers are a versatile design platform for biocompatible viscoelastic materials. Their linear response to a step strain displays a fast, near-exponential relaxation when using low-valence cross-linkers, while larger supramolecular cross-linkers bring about much slower dynamics involving a wide distribution of timescales whose physical origin is still debated. Here, we propose a model where the relaxation of polymer gels in the dilute regime originates from elementary events in which the bonds connecting two neighboring cross-linkers all disconnect. Larger cross-linkers allow for a greater average number of bonds connecting them but also generate more heterogeneity. We characterize the resulting distribution of relaxation timescales analytically and accurately reproduce stress relaxation measurements on metal-coordinated hydrogels with a variety of cross-linker sizes including ions, metal-organic cages, and nanoparticles. Our approach is simple enough to be extended to any cross-linker size and could thus be harnessed for the rational design of complex viscoelastic materials.


INTRODUCTION
Soft hydrogels are ubiquitous in biology and dictate the mechanics of cells and tissues [1].Due to their biocompatibility, synthetic hydrogels are thus prime candidates to serve as robust soft tissue implants, although fine control of their viscoelastic properties is crucial for their success in this role [2,3].In simple viscoelastic materials, stress relaxes according to a single exponential with a single relaxation time.This is not however the case for most biological materials such as cells [4], tissues [5], mucus [6], and biofilms [7].Instead, their relaxation is characterized by a broad distribution of relaxation times [8].Such relaxation is often heuristically described by a stretched exponential: where smaller values of the exponent α ∈]0; 1[ denote broader distributions of relaxation time scales [9].Other similarly phenomenological fitting functions include power-law dependences of σ on t [10][11][12] and lognormal distributions of the relaxation times [13,14].
Associative gels, which relax by a succession of binding and re-binding events [15], offer a promising route to design controllable viscoelastic materials.It is thus possible to tune their relaxation time by adjusting the chemical binding energy of their crosslinkers [16,17].Although this chemistry-based approach allows tuning of the overall stress relaxation time as illustrated in Fig. 1(a), less is known about the different approaches to tune the shape of the stress relaxation curve of reversible hydrogels.Accordingly, most existing models for the relaxation of multivalent gels focus on regimes dominated by a single relaxation time scale [18], leading to exponential relaxation [19].Control over the distribution of relaxation time scales could however be achieved in synthetic hydrogels connected with multivalent crosslinkers such as nanoparticles [20], metal-organic cages [21], clay [22], and latex beads [23], which are known to exhibit nonexponential viscoelastic relaxation.Here, we aim to elucidate this emergence of a wide distribution of time scales in materials with high valence crosslinkers to enable the rational design of complex gels [Fig.1(b)].Here we use the term "valence" to designate the number of polymer strands that a crosslinker can bind, a property sometimes also sometimes referred to as their "functionality" [24].
We propose that the emergence of a broad distribution of relaxation time scales arises from microscopic events consisting of the severing of the physical connection between two crosslinkers.We first propose a model where FIG. 1. High-valence crosslinkers yield a slow, potentially complex unbinding dynamics (a) Hydrogels held together by small crosslinkers relax over the time scale associated with the unbinding of a single polymer strand.There, stress decays exponentially in response to a step strain: σ(t)/σ(0) = exp(−t/τ ).Since t/τ = exp(ln t − ln τ ), changes in the time scale τ shift the σ vs. ln t relaxation curve horizontally, but do not alter its shape.(b) In contrast, relaxation events in the presence of high-valence crosslinkers require the simultaneous unbinding of many polymer stands.The associated time scale is long and highly variable depending on the number of strands involved in the "superbond" (grey shade).As a result the stress relaxation of such gels is no longer exponential and the precise shape of the relaxation curve strongly depends on the valence of the crosslinkers.this connection, hereafter termed "superbond", breaks if all its constitutive crosslinkers are detached at the same time [Fig.1(b)].We show that the breaking time of a superbond increases exponentially with the number of strands involved, consistent with previous observation [25].As a result of this strong dependence, small spatial heterogeneities in the polymer concentration may result in widely different relaxation times from one superbond to the next.Such exponential amplification of relaxation times originating from small structural differences forms the basis of models previously used to describe the relaxation of soft glasses [26][27][28].In contrast with these studies, our approach explicitly models the microscopic basis of this amplification.That allows it to not only recover relaxation curves virtually indistinguishable from those discussed in previous studies, but to also predict the influence of temperature and crosslinker valence on the macroscopic stress relaxation observed in the resulting gel.The details of the polymer strand morphology are not central to this influence, and we thus enclose them in a few effective parameters which could be derived from first principles in specialized models related to specific implementations of our basic mechanism.
To confirm these predictions, we conduct experiments on hydrogels with four distinct crosslinker types of different sizes, and find that our model quantitatively reproduces multiple relaxation curves using this small set of microscopic parameters.Finally, we show that several phenomenological fitting functions used in the literature can be recovered as asymptotic regimes of our analytical model.

Model of a single superbond
We first model a single superbond in the simple, experimentally relevant [21] case of strong interactions combined with short polymers, which implies negligible entanglements.In the limit of very large and rigid crosslinkers, the polymer layer around a crosslinker is locally planar, its structure is not affected by small fluctuations of polymer concentration and its thickness fixes the distance between crosslinkers [21,29].We thus use the simplifying assumption that all individual bonds participating in a superbond are identical and non-interacting, an approximation whose validity we ultimately assess through comparisons with experiments.
We model the attachment and detachment of a single polymer strand from a pair of crosslinkers as shown in Fig. 2(a).When both its ends are bound, the strand may or may not connect two different crosslinkers.The corresponding "bridging" and "looping" states have the same energy since we assume the polymer strand to be completely flexible, and we denote by ∆S the entropy difference between them.To transition between these two states, the strand must disconnect one of its ends and form a "dangling" state.The disconnection of the strand in this state implies an energy barrier ∆E that is much larger than the thermal energy k B T = β −1 .This implies that the dangling state is short-lived, and thus need not be explicitly included in our modeling.Our approximation scheme implies that the overall rate ω + and ω − to go from the looping to the bridging state and back are constant.They read where the typical time scale τ 0 takes into account the entropy difference between the looping and dangling state.
At equilibrium, we denote the probability for a single polymer strand to create a bridge as p on = 1 − p off = 1/(1 + e ∆S ).We now consider the dynamics of a superbond involving N polymer strands.Within our approximation of independent attachment and detachment of the individual polymer strands, the superbond undergoes the Markov process illustrated in Fig. 2(b) and the probability P n (t) for n strands to create bridges between the two crosslinkers at time t satisfies the master equation FIG. 2. We model superbond breaking as the disconnection of many independent polymer strands.(a) Disconnecting a single polymer strand requires going through a high-energy, short-lived dangling state (larger arrows indicate faster transitions).The looping and bridging states both have two polymer-crosslinker bonds, and therefore have the same energy.(b) Individual strands in a superbond attach and detach independently, resulting in a one-dimensional random walk in the number n of attached strands [Eq.( 3)].Here we only draw the bridging strands, not the looping strands.
which ensures that the number of bridging polymers can never be greater than N .To determine the rate at which a superbond breaks, we set an absorbing boundary condition P 0 (t) = 0 and define its survival probability as S(t) = N n=1 P n (t).In the limit N ≫ 1 where a large number of strands are involved in the superbond, we show in Sec.S1 of the Supplementary Information that the detachment of the two beads is analogous to a Kramers escape problem.We thus prove that the survival probability decays as a single exponential S(t) = exp(−t/τ N ) [30] with an average detachment time The breaking of the superbond can thus be assimilated to a Poisson process with rate 1/τ N regardless of the initial condition P n (0).The strong, exponential dependence of τ N on N implies that any dispersity in the number of strands involved in a superbond may result in a wide distribution of time scales.

Model of the relaxation of a gel
Two factors influence the dispersity of N .First, its value is constrained by the available space at the surface of each crosslinker, which we model by setting an upper bound N sat on the number of polymer strands (in a loop or a bridge) participating in any superbond.Second, depending on the local density of polymer in the vicinity of the superbond, the actual number of strands present may be lower than N sat .In the regime where the polymer solution surrounding the crosslinkers is dilute, polymer strands are independently distributed throughout the system.As a result, the distribution of local strand concentrations within a small volume surrounding a superbond follows a Poisson distribution.We thus assume that N is also described by a Poisson distribution up to its saturation at N sat : where N would be the average number of strands in a superbond in the absence of saturation and thus depends on the ratio of polymer to crosslinker concentration.Note that the specific form of the distribution used in Eq. ( 5) does not significantly modify our results, as discussed later.
In response to a step strain, we assume that each superbond is stretched by an equal amount and resists the deformation with an equal force prior to breaking.Superbonds may subsequently reform, but the newly formed bonds are not preferentially stretched in the direction of the step strain and therefore do not contribute to the macroscopic stress on average.Denoting by t = 0 the time at which the step strain is applied and by σ(t) the resulting time-dependent shear stress, the progressive breaking of the initial superbonds results in the following stress response function: While the breaking times τ N are unaffected by the applied stress in the linear response regime, nonlinearities could easily be included in our formalism by making ∆S stress-dependent and thus favor strand detachment.The relaxation described in Eq. ( 6) occurs in two stages.At long times t ≫ τ Nsat , few short-lived superbonds remain.Saturated superbonds (N = N sat ) dominate the response and Eq. ( 6) is dominated by the last term of its sum.As a result, the stress relaxes exponentially over time, as seen from the linearity of the log-lin curves of Fig. 3(a) for large values of t.Systems with smaller values of N sat manifest this regime at earlier times; in the most extreme case, the relaxation of a system where superbonds involve at most a single polymer strand (N sat = 1) is fully exponential and extremely fast as compared to systems with higher N sat .Over short times (t ≪ τ Nsat ), stress relaxation involves multiple time scales.This nonexponential regime is apparent on the left of Fig. 3(a).These two regimes have already been reported in several Relationship between the stretch exponent α quantifying the nonexponential character of the relaxation and the microscopic parameter Nsat/ N .Here p off = 0.2.A low Nsat/ N gives an exponential relaxation (α ≃ 1), while a larger Nsat/ N leads to a more complex behavior (α < 1).While α appears to converge to a finite value for large Nsat/ N for the largest values of N , this behavior is contingent on our choice of fitting interval.This issue does not affect the rest of the curves.Large stars correspond to the curves represented in (a).Inset: illustration of the quality of the fits between the heuristic stretched exponential [dashed orange line, Eq. ( 1)] and our prediction [solid blue line, Eq. ( 6)].
While Eq. ( 6) is not identical to the stretched exponential of Eq. ( 1), the inset of Fig. 3(b) shows that they are remarkably close in practice.We thus relate the stretch exponent α to the saturation number N sat by fitting a stretched exponential to our predicted stress response function over the time interval required to relax 90% of the initial stress [Fig.3(b)].The fits are very close matches, and consistently give correlation factors r 2 > 0.98 (see detailed plots in Supplementary Fig. S2).If N sat ≲ 0.5 N then α ≃ 1, indicating a nearlyexponential relaxation.Indeed, in that case superbond saturation occurs well before the peak of the Poisson distribution of N .Physically, this implies that the local polymer concentration surrounding most superbonds is sufficient to saturate them.As almost all superbonds are saturated, they decay over the same time scale τ Nsat .As a result, the material as a whole displays an exponential relaxation.For larger values of N sat , the Poisson distribution is less affected by the saturation and the dynamics is set by the successive decay of superbonds involving an increasing number of strands, implying lower values of α.The larger the value of N , the sharper the crossover between these two regimes.

Experiments
To validate our model of the effect of crosslinker valence on hydrogel relaxation, we perform step-strain experiments of poly(ethylene glycol) (PEG)-based gels involving a range of crosslink valences and compare Eq. ( 6) to the resulting relaxation curves.We implement strand-crosslink bonds based on two sets of metalcooordination chemistry: nitrocatechol-Fe 3+ coordination and bispyridine-Pd 2+ coordination.The first set of gels is made with nitrocatechol-functionalized PEG crosslinked by single Fe 3+ ions with an estimated valence of 3, and by iron oxide nanoparticles.The nanoparticles have a mean diameter of 7 nm with a surface area that allows a valence of ∼ 100 ligands [20].The second set of gels are made with bispyridine-functionalized PEG, wherein bis-meta-pyridine ligands induce self-assembly of gels that are crosslinked by Pd 2 L 4 nanocages with a valence of 4, and bis-para-pyridine ligands induce self-assembly of gels that are crosslinked by Pd 12 L 24 nanocages with a valence of 24 [31].As shown in Fig. 4, these four distinct gel designs result in a broad range of relaxation behaviors.Overall, large-valence gels and lower temperatures result in longer relaxation times, consistent with the illustration of Fig. 1(a).The relaxation curves associated to high-valence crosslinkers are also less steep, consistent with the involvement of a broader distributions of relaxation times and the schematic of Fig. 1(b).
At a more detailed level, our assumption that the dynamics of single polymer strand proceeds independently of its environment implies the existence of a single energy scale ∆E.As a result, we predict that all time scales involved in the relaxation are proportional to exp(−β∆E).We confirm this through a time-temperature collapse shown in the insets of Fig. 4 (see Sec. S4 of the Supplementary Information for details).This collapse provides us with the value of ∆E for each of our four systems, which we report in Table I.The binding energy value ∆E for the Pd 2 L 4 and Pd 12 L 24 gels match, as expected from the fact that they originate from the same composition.On the other hand, the ∆E of the nanoparticles gel cannot be directly compared to the Fe 3+ ion gel because of their dramatic physico-chemical properties as well as the acidic pH used to facilitate mono-coordination in nanoparticles.
To compare the temperature-collapsed curves to our prediction of Eq. ( 6), we fit the parameters p off , τ 0 , N and N sat across multiple temperatures.The resulting fits, shown in Fig. 4, display a good agreement between the theory and experiments across up to 4 orders of magnitude in time scales.The fitted values of N sat are moreover consistent with the estimated valence of the crosslinkers (Table I), which confirms our interpretation of the physical origin of N sat .A more direct relation between the valence and N sat would require additional information about the structure of each category of gels, and in particular the number of nearest neighbors of an individ-

FIG. 4. Stress relaxation function for four experimental systems with increasing crosslinker valences (see Table I for values).
Here we use a log-lin scale [unlike in Fig. 3  TABLE I.Estimated and fitted parameters involved in the comparison between experiment and theory in Fig. 4. The energies are given in units of kBT for T = 300 K. Instead of displaying the parameter τ0, we present the more easily interpreted unbinding time of a single polymer strand at 300K, namely τ1 = τ0e β 300 ∆E /p off .
ual crosslink, which remains debated for these types of gels [32][33][34].The possible clustering of crosslinkers may also influence this relation, as nanocage systems similar to ours [35] may form higher-order nanocage structures.While such complexities as well as possible imperfections in our fitting procedure complicate the literal interpretation of the fitted values of N sat , these nonetheless clearly increase with increasing crosslinker valence.The fit also supports the notion that the mean number of strands per superbond N accounts for the distribution of relaxation time scales in our gels.The Fe 3+ gel thus displays an exponential relaxation consistent with N = N sat = 1.The higher-valence Pd 2 L 4 and Pd 12 L 24 systems have a complex relaxation at early times followed by an exponential behavior, as expected for N ≃ N sat > 1.As expected from our model, the crossover time τ Nsat separating the two regimes is larger in the higher-valence Pd 12 L 24 gel.Finally, the high-valence nanoparticle system shows an extended complex relaxation associated with N < N sat , thus confirming that all the qualitative relaxation regimes discussed in the previous sections are experimentally relevant.

Distribution of relaxation time scales
To further visualize the differences between the responses of our gels, we plot the distributions of relaxation times p(τ ) for our fitted model in Fig. 5.The Fe 3+ gels, which relax according to a single exponential and whose p(τ ) are therefore delta functions, are not represented there.In the Pd 2 L 4 and Pd 12 L 24 systems, a distribution characterized by an initially decreasing distribution of time scales is interrupted by a valence-dependent maximum relaxation time τ Nsat .That time is comprised within the range of time scales observed in Fig. 4, accounting for the crossover to an exponential relaxation within this range.In nanoparticle systems, by contrast, the crossover occurs much later and thus cannot be directly observed in experiments.In all cases, the precise form of the distribution of time scales used in the domain τ < τ Nsat does not critically affect the predicted relaxation curves.Indeed, we show in Sec.S6 of the Supplementary Information that replacing the Poisson distribution of Eq. ( 5) with other distributions with the same mean and variance lead to essentially indistinguishable predictions over experimentally observable time scales.FIG. 5. Distribution of relaxation times corresponding to the theoretical plots of Fig. 4. The distribution associated with the Fe 3+ gel is a delta function for τ /τ1 = 1 and is thus not represented on this graph.The time distributions are given by a Poisson distribution cut for τ = τN sat (vertical lines), as described in Eq. ( 6).The dotted lines represent what these distributions would be in the absence of this saturation.
This emphasizes the robustness of our predictions to the details of that choice of distribution.They are instead primarily determined by the mean and maximum superbond sizes, N and N sat .
In the limit of large N and even larger N sat , the complex relaxation phase of our model characterized by t < τ Nsat may display analytical behaviors identical to some widely used rheological fitting functions.In this regime, the Poisson distribution p(N ) of Eq. ( 5) goes to a Gaussian.Since according to Eq. ( 4) the variable N is essentially the logarithm of the relaxation time τ for N ≫ 1, this results in a log-normal distribution of relaxation time scales: This result adds additional insights to this widely used fitting functional form, as it allows to relate the mean and variance of the distribution to the underlying crosslinkerscale parameters [13,14].It moreover offers a potential molecular-level justification for its use in describing the complex relaxation of systems with multivalent crosslinks.In the alternative case where p(N ) is a decaying exponential, our model results in power-law distributed relaxation time scales and the stress response function takes the form This result may also be presented in terms of the dependence of the storage and loss moduli on the frequency ω in an oscillatory rheology experiment.We thus predict that for γ < 1 The results larger values of γ and detailed derivations of Eqs.(7)(8)(9) are shown in the Supplementary Information.Again, this result has the potential to account for the power-law relaxation observed in many rheological systems [10][11][12], in addition to providing a link to their microscopic constituents.Overall, these results suggest a possible control of the system's rheology through the characteristics of p(N ), which could in turn be modulated through the spatial distribution of the polymer strands and the dispersity of the crosslinkers.

DISCUSSION
Our simple model recapitulates a wide range of rheological behaviors in multivalent systems based on two key superbond parameters: the mean size N and the maximum size N sat .These respectively control the amplitude of the fluctuations in superbond size and the longest superbond relaxation time scale.Prior to the longest relaxation time, the system displays an increasingly nonexponential response for increasing N [Fig.3(b)].Beyond it, it crosses over into exponential relaxation.In contrast with widely used phenomenological fitting parameters, our two variables yield reliable insights into the underlying microscopic dynamics, as demonstrated by the agreement of their fitting values with our a priori knownledge of four experimental systems covering a wide range of values of N and N sat .
Our model bears a mathematical similarity with standard random energy trap models [36].There, a longtailed relaxation emerges from a short-tailed distribution of trap depths due to the exponential dependence of the relaxation times on the trap depths.Similarly, here a nonexponential relaxation emerges from a shorttailed distribution of superbond sizes N [Eq.( 5)] thanks to the exponential dependence of τ N on N [Eq.( 4)].In contrast with trap models, however, our model does not predict a glass transition upon a lowering of temperature.It instead displays a simple Arrhenius time-temperature relation, consistent with the experimental collapses in the insets of Fig. (4).Here again, an additional benefit of our approach is the direct connection between the predicted relaxation and experimentally accessible parameters such as the crosslinker surface (through N sat ).
Our model's focus on the collective aspects of superbond breaking and the characteristics of the crosslinkers implies that it encloses most of the physics of the polymer strands within a few mesoscopic parameters, mainly τ 0 and ∆S.Within our approach, the morphology of the polymer thus does not affect the form of our relaxation, although it may lead to a rescaling of the relaxation times of Eq. ( 4).This formulation remains valid as long as the length and concentration of the polymer strands is low enough that the polymer strands do not become significantly entangled, which could spoil the Poissonian attachment/detachment process of Eq. ( 2).Even in this case however this equation may not be strictly valid, as the polymer layer in a superbond with many bound crosslinkers tends to be more compressed than in one with few.This effect should leads to a smooth (likely power law) dependence of ω ± on N , which would preserve the dominance of the much more abrupt exponential dependence of τ N on N .As a result, while such polymer brush effects could induce corrections in our estimations of the model parameters, the basic mechanism outlined here should still hold in their presence.
Our model reproduces several qualitative characteristics of the rheology of multivalent gels, such as the strong influence of the crosslinker valence, Arrhenius temperature dependence and the transition between a nonexponential and an exponential regime at long times.Due to its simple, widely applicable microscopic assumptions, we believe that it could help shed light on and assist the design of a wide range of multivalent systems.Beyond composite gels, it could thus apply to RNA-protein biocondensates where multivalent interactions between proteins are mediated by RNA strands [37], as well as cytoskeletal systems where filaments linked to many other filaments display a slow relaxation reminiscent of that of our multivalent crosslinkers [38].

Synthesis of Fe3O4 nanoparticles (NPs)
Bare Fe 3 O 4 NPs are synthesized following previously reported methods (Li et al. 2016).100 mg as-synthesized NPs are re-dispersed in 80 mL of 1:1 (v/v) solution of CHCl 3 and DMF, and 100 mg 1-arm PEG-C is added.The mixture is homogenized and equilibrated by pulsed sonication (pulse: 10 s on + 4 s off; power: 125 W) for 1 hour.Then the mixture is centrifuged at 10000 rpm for 10 min to remove any aggregates, and rotary evaporated at 50 • C, 30 mbar to remove CHCl 3 .Then the NP solution is precipitated in 150 mL cold Et 2 O (-20 • C).The precipitate is re-dispersed in H 2 O, and freezedried.The resulting NPs are 7 nm in diameter.

Preparation of the Fe3 + -NC gels
Preparation procedure is similar to a previously reported protocol [39], except that the gel is made in DMSO instead of H 2 O. 50 µL of 200 mg/mL 4-arm PEG-NC solution in DMSO is mixed with 16.7 µL of 80 mM FeCl 3 solution in DMSO (ligand: Fe 3+ molar ratio of 3:1).Then 33.3 µL DMSO and 13.8 µL TEA is added to facilitate deprotonation, and a gel is formed.

Preparation of the Pd2L4 gels
The synthesis of polymer and gel preparation procedures for P2L4 is the same as a reported protocol [31] with minor modifications.The annealing of the Pd 2 L 4 polyMOC gel was done at 60°C for 1 hr instead of 80°C for 4 hr, and 1.05 equivalent of Pd(NO 3 ) 2 . 2 H 2 O (relative to bifunctional polymer ligand) was used instead of 1 equivalent.

Preparation of the Pd12L24 gels
The synthesis of polymer and gel preparation procedures for polyMOC is the same as a reported protocol [31] Preparation of the NP gels Preparation procedure is the same as the reported protocol [40].Briefly, 20 mg PEGylated Fe 3 O 4 NPs (equivalent to 20 mg Fe 3 O 4 core) and 20 mg 4-arm PEG-NC are mixed in a 0.2 M HCl aqueous solution.The solution mixture (pH = 2) is transferred into a mold and sealed, and a solid gel is obtained after curing in a 50 • C oven for 24 hours.

Rheology
Stress relaxation measurements are done on an Anton Paar rheometer with parallel plate geometry (10 mm diameter flat probe for NP gels and polyMOC gels, and 25 mm diameter cone probe for Fe 3+ gels).All tests are done immediately after transferring the gel sample onto the sample stage.A Peltier hood is used for all experiments to control the measurement temperature and prevent solvent evaporation.H 2 O based samples are furthermore sealed with mineral oil before experimentation to reduce the evaporation rate.Relaxation tests were performed by applying a γ = 0.005.step strain for the NP gel, and γ = 0.02 step strain for the other three systems.
Here we present the mathematical proofs of the main results of the manuscript as well as details on the methodology used to analyze the experimental data.

S1. DISTRIBUTION OF SUPERBOND BREAKING TIME AND DERIVATION OF τN
Here we show that the survival probability for the detachment of a superbond (illustrated in main text Fig. 2) containing many polymer strands (N → ∞) asymptotically goes to S(t) = e −t/τ N , where τ N is given by Eq. ( 4) of the main text.We first consider a general one-step process and derive the basic recursion equation used throughout the proof in Sec.S1.1.We solve the recursion in Sec.S1.2 and express the generating function of S(t) as a double sum.In Sec.S1.3, we apply the resulting formula to our particular problem and take the continuum limit of the second sum.Finally, we compute both sums in the N → ∞ limit in Sec.S1.4.Our derivation is adapted from the calculation presented in the appendix of Ref. [1].

S1.1. Backward Kolmogorov equation for the generating function of S(t)
We consider a one-step process, i.e., a stochastic process consisting of transitions between consecutive discrete states on a line, with transition rates r n and g n illustrated in Fig. S1(a).We denote the probability for the particle to be in state k at time t after starting in state n at time 0 by P (k, t|n).We assume an absorbing boundary condition in 0 and a reflecting boundary condition in N , i.e., Inserting these definitions into Eq.(S2) yields which we endeavor to solve for h n (α) in the following.

S1.2. Sum equation for the generating function
We define a rescaled current between sites n − 1 and n This allows us to turn the two-step recursion of Eq. (S4) into one with only one step: which can easily be summed as We now invert Eq. (S5) and use Eq.(S7) to express the finite difference (h n − h n−1 ).We further use the property that h m = h 0 + m n=1 (h n − h n−1 ) and recognize that h 0 = 0 due to Eq. (S1) to obtain Illustration of the similarity of our modeled stress response function with a stretched exponential.We plot the relaxation modulus computed using Eq. ( 6) of the main text for N ∈ [3,5,10,15].For each value of N , we plot four values of Nsat, namely Nsat = 0.1 N , 0.5 N , N and 1.5 N , p off = 0.2.Each plot also mentions the value of the fitted stretch exponent α and the correlation coefficient r 2 .

S2. LINK BETWEEN α AND Nsat/ N
Here establish the connection between the stretch exponent α and the values of N sat / N shown in Fig. 4 of the main text.To mimic the observation of an experimental step strain over a finite time window, we focus our attention on the time interval between t = 0 and t = τ 90 , where τ 90 is the time required to relax 90% of the stress, i.e., σ(τ 90 ) = 0.1 × σ(0).We plot the relaxation curve given by Eq. ( 6) of the main text over this time window, then perform a least-squares fit using a stretched exponential [Eq.(1) of the main text] with α and τ as fitting parameters.As shown in Fig. S2, the agreement is excellent for a large majority of the parameters used.The corresponding value of the fitting parameters (τ and α) for a broader variety of N and N sat is also provided in Fig. S3.This suggests that experimental curves that are well fitted by a stretched exponential could be equally well described by our model.

S3. TIME-TEMPERATURE COLLAPSE
Here we describe the procedure used to determine the binding energy ∆E in the experimental systems discussed in the main text.Equation (4) of the main text implies that the temperature dependence of the stress response function can be eliminated by expressing it as a function of the rescaled time t = te β∆E .This should cause the relaxation curves of a given system at different temperatures to collapse.
For each type of ligand, we have 5 datasets showing the stress relaxation function as a function of time at each different temperature {T To enable the comparison between timerescaled datasets, we first define an interpolating function for the stress relaxation function at each temperature used.We thus compute the set of interpolating coefficients p by perform a least-square fit of the following rational function to the datapoints t . We furthermore define the interval of definition of P (α) (x) as the range over which data is available, i.e., I P (α) = 0, max i t (α) i .
We then perform the collapse of the {T (α) } α∈ [1,4] interpolated curves onto the T (0) curve.To this effect we define the set of rescaling coefficients {a (α) } α∈ [1,4] and performs a separate time rescaling for each temperature: t = te a (α) .For each α ∈ [1,4], we optimise the semidistance between the functions t → P (0) (t) and t → P (α) (te a (α) ) with respect to a (α) .The resulting collapsed curves are shown in Fig. S4 (a,b,c).The optimal rescaling coefficients are plotted as a function of the inverse temperature 1/k B T in Fig. S4 (d,e,f).Consistent with the time-temperature collapse hypothesis, this dependence is affine, and we use the slope of the best fitting line as our value of the binding energy ∆E.

S4. FIT OF THE STRESS RELAXATION FUNCTION TO OUR THEORETICAL PREDICTION
In the main text, we fit the experimental curves with the stress relaxation function predicted by our model.We then represent them on a log-lin scale to allow the simultaneous visualization of short and long time scales.To demonstrate the robustness of our fits, in Fig. S5 we replot these curves in a lin-lin-scale, as well as a lin-log scale that emphasizes intervals of exponential relaxation as straight lines.

FIG. S4.
Collapse of the relaxation modulus: (a-c) Collapsed relaxation modulus of Fe 3+ ions, Pd12L24 nanocages and nanoparticles respectively after the rescaling of the time for an optimized collapse .The curves are represented on a log-lin coordinate system, but the collapsing procedure is performed on a lin-lin scale .(d-f) corresponding rescaling parameters as a function of 1/(kBT ) the inverse temperature.The slope of the line is −∆E and the legend gives the value of ∆E in kBT unit at 300K.

S5. RATIONALIZATION OF THE POISSON DISTRIBUTION OF THE SUPERBOND SIZE p(N )
The polymers used in our experiments are 4-arms polyethylene glycol (PEG).At the end of each arm is a nitrocatechol ligand that allows crosslinker binding.In our model, we assume that the ends of a polymer are always attached to a ligand.For this reason, the diffusion of such a polymer over a distance comparable to the polymer size occurs on a time scale comparable to the time required to rearrange the bonds between crosslinkers, which corresponds to the time required for the relaxation of the stress in the system.Let us consider that the 4-arm PEG are able to diffuse over a volume v during the time of the experiment.We model the spreading of the polymers in the system by discretizing the system into small boxes of volume v between which no polymer exchange occurs over the duration of the experiment.As a result the distribution of the polymers over the boxes is due to the initial preparation of the system.We assume that this processes places each polymer in a random box with equal probability.As a result, the probability that a specific box contains n polymers is given by a Poisson distribution: where ρ PEG is the average concentration of PEG in the system, and vρ PEG is the mean (over the system) number of P EG in a box of volume v. Equation (S21) is the basis for Eq. ( 5) of the main text.

S6. EXPERIMENTAL FIT USING ALTERNATIVE DISTRIBUTIONS OF N
To demonstrate the robustness of our model, we perform the fit of the experimental data using probability distribution different from the Poisson distribution of Eq. ( 5) of the main text.For this, we first define three new distributions with the same mean (and when possible the same variance) than the one of Eq. ( 5) of the main text and keep the Here N = 10 and p off = 0.18 ⇒ γ ≃ 0.0583.(c) Plot of the three distinct asymptotic regimes for a higher value of γ ( N = 10 and p off = 0.935 ⇒ γ = 1.49).The marker at ωτ1 = τ1/τN sat denotes the expected position of the low-frequency crossover, while the high-frequency crossover is expected for ωτ1 ≈ 1.

FIG. 3 .
FIG. 3. (a) Disperse, high-valence superbonds initially display a nonexponential mechanical relaxation, then cross over to an exponential regime when only the saturated superbonds remain.Curves plotted from Eq. (6) with p off = 0.2, N = 10 and different values of Nsat as indicated on each curve.(b) Relationship between the stretch exponent α quantifying the nonexponential character of the relaxation and the microscopic parameter Nsat/ N .Here p off = 0.2.A low Nsat/N gives an exponential relaxation (α ≃ 1), while a larger Nsat/ N leads to a more complex behavior (α < 1).While α appears to converge to a finite value for large Nsat/ N for the largest values of N , this behavior is contingent on our choice of fitting interval.This issue does not affect the rest of the curves.Large stars correspond to the curves represented in (a).Inset: illustration of the quality of the fits between the heuristic stretched exponential [dashed orange line, Eq. (1)] and our prediction [solid blue line, Eq. (6)].
FIG. 4. Stress relaxation function for four experimental systems with increasing crosslinker valences (see Table I for values).Here we use a log-lin scale [unlike in Fig.3(a)] to facilitate the visualization a large range of time scales.Alternate representations are available as Fig. S5.Symbols are experimental datapoints, and the lines are the associated fitting curves.Insets: timetemperature collapsed data obtained by a rescaling t → te β∆E .

P
FIG. S1.Superbond detachment as a Kramers-like barrier-crossing problem.(a) Definition of the rates of the one-step process.(b) Profile of the pseudo-free energy defined in Eq. (S13).Superbond detachment requires the system to fluctuate out of the free energy well to the x = 0 absorbing state, with 1/N playing the role of a temperature.

3 N = 5 N = 10 N = 15 N = 20 FIG
FIG. S3.Best fit values of the stretched exponential parameters τ and α for a range of values of N and Nsat.The right-handside panel is identical to Fig. 4 of the main text.

FIG. S8 .
FIG. S8.Comparison between the storage and loss moduli computed from the exact expression Eq. (S35) (solid lines) and the asymptotic expressions of Eqs.(S36-S38) (dashed lines).(a) Plots in the large Nsat limit (here Nsat = 100), showing a good agreement with the power law regime of Eq. (S37) for two values of N and for constant p off = 0.18 corresponding to γ ≃ 0.116 and γ ≃ 0.0583.(b) Plots for a smaller value of Nsat (Nsat = 30) showing the three distinct asymptotic regimes.Here N = 10 and p off = 0.18 ⇒ γ ≃ 0.0583.(c) Plot of the three distinct asymptotic regimes for a higher value of γ ( N = 10 and p off = 0.935 ⇒ γ = 1.49).The marker at ωτ1 = τ1/τN sat denotes the expected position of the low-frequency crossover, while the high-frequency crossover is expected for ωτ1 ≈ 1.