A foundational framework for the mesoscale modeling of dynamic elastomers and
gels

Robert J. Wagner1***Correspondence to<EMAIL_ADDRESS>&
Meredith N. Silberstein2
1Department of Mechanical Engineering, State University of New York at
Binghamton, Binghamton, NY, USA
2Sibley School of Mechanical & Aerosspace Engineering, Cornell University,
Ithaca, NY, USA

Discrete mesoscale network models, in which explicitly modeled polymer chains
are replaced by implicit pairwise potentials, are capable of predicting the
macroscale mechanical response of polymeric materials such as elastomers and
gels, while offering greater insight into microstructural phenomena than
constitutive theory or macroscale experiments alone. However, whether such
mesoscale models accurately represent the molecular structures of polymer
networks requires investigation during their development, particularly in the
case of dynamic polymers that restructure in time. We here introduce and
compare the topological and mechanical predictions of an idealized, reduced-
order mesoscale approach in which only tethered dynamic bonding sites and
crosslinks in a polymer’s backbone are explicitly modeled, to those of
molecular theory and a Kremer-Grest, coarse-grained molecular dynamics
approach. We find that for short chain networks ($\sim$12 Kuhn lengths per
chain segment) at intermediate polymer packing fractions, undergoing
relatively slow loading rates (compared to the monomer diffusion rate), the
mesoscale approach reasonably reproduces the chain conformations, bond kinetic
rates, and ensemble stress responses predicted by molecular theory and the
bead-spring model. Further, it does so with a 90% reduction in computational
cost. These savings grant the mesoscale model access to larger spatiotemporal
domains than conventional molecular dynamics, enabling simulation of large
deformations as well as durations approaching experimental timescales (e.g.,
those utilized in dynamic mechanical analysis). While the model investigated
is for monodisperse polymer networks in theta-solvent, without entanglement,
charge interactions, long-range dynamic bond interactions, or other
confounding physical effects, this work highlights the utility of these models
and lays a foundational groundwork for the incorporation of such phenomena
moving forward.

##  1 Introduction

Synthetic polymers are one of the most diversely applied classes of materials
available to engineers. Due to their extensive sets of design parameters
(e.g., molecular weight; polydispersity; composition; crosslink
concentration/type; targeted incorporation of supramolecular interactions;
swelling with solvent; etc.), polymers exhibit a broad set of emergent
mechanical structures and properties appropriate for diverse applications
including in structural composites (Chung, 2019), renewable energy systems
(Rydz et al., 2021), packaging (Vallejos et al., 2022, Ibrahim et al., 2022),
adhesives (Heinzmann et al., 2016), and biomaterials (Chen et al., 2022). A
particular trait of interest in polymeric design is the inclusion of dynamic
bonds. Polymers with dynamics bonds include covalently adaptable polymers
(Kloxin et al., 2010) such as vitrimers (Yue et al., 2020, Shen et al., 2021,
Hubbard et al., 2021, 2022), and supramolecular polymers (Brunsveld et al.,
2001, Yount et al., 2005) such as physically bonded elastomers and gels
(Vidavsky et al., 2020, Mordvinkin et al., 2021, Xu et al., 2022). The
dissociation of such bonds from highly stressed states, followed by their re-
attachment into lower energy configurations is a dissipative mechanism that
can grant dynamic polymers exceptional toughness (Haque et al., 2011, Gong et
al., 2016, Bai et al., 2018, Li and Gong, 2022), extensibility (Tuncaboylu et
al., 2011, Jeon et al., 2016, Zhang et al., 2019, Cai et al., 2022), self-
healing capabilities (Kersey et al., 2006, Brochu et al., 2011, Li et al.,
2016, 2021), and even mechanosensitive response (Wojtecki et al., 2011, Eom et
al., 2021, Doolan et al., 2023). Furthermore, the additional design parameters
introduced by the inclusion of dynamic bonds (e.g., tether length for
telechelic bonds (Ge and Rubinstein, 2015, Mordvinkin et al., 2021),
coordination number of charged species in metallopolymers (Yount et al., 2005,
Zhang et al., 2020, Vidavsky et al., 2020), etc.) provide polymers with
enhanced tunability. For instance, simply changing the molar concentration or
type of dynamically bonding species in a polymer may yield elastic moduli
spanning multiple orders of magnitude (Xu et al., 2022, Huang et al., 2022).
Therefore, predictive models that can map the highly varied emergent
mechanical properties of such polymers from their ab initio compositions
(prior to synthesis and experimental testing) may greatly benefit researchers
developing new materials requiring certain traits. However, the same rich sets
of constituents and compositional design choices responsible for dynamic
polymers’ highly varied properties also imbue them with hierarchical length
scales and relaxation timescales spanning from $10^{-10}$ to $10^{-3}$ m and
$10^{-15}$ to $10^{6}$ s, respectively, making this mapping process
challenging.

At the macroscale – well above the thermodynamic limit and length scales of
heterogeneities – continuum approaches (James and Guth, 1943, Flory, 1985,
Tanaka and Edwards, 1992, Arruda and Boyce, 1993, Wu and Van Der Giessen,
1993, Miehe et al., 2004, Vernerey et al., 2017) are suitable for mechanical
property predictions, provided proper implicit representation of the pertinent
first order physics such as single-chain mechanics (James and Guth, 1943,
Miehe et al., 2004, Saleh, 2015), swelling (Flory, 1942, Hong et al., 2009,
Bouklas and Huang, 2012), bond dynamics (Leibler et al., 1991, Stukalin et
al., 2013), macroscopic damage (Bouklas et al., 2015, Shen and Vernerey, 2020,
Lee et al., 2023), etc. Although such approaches have been used extensively to
successfully predict polymeric properties and indirectly deduce
microstructural origins of observed traits with a wide breadth of physical
effects (Xu et al., 2022, Bosnjak et al., 2022), they generally rely on some
combination of homogenization assumptions, affine deformations, and mean-field
approximations (C. Picu, 2011). As a result, they are limited in their ability
to directly map composition to microstructure, and then microstructure to
global mechanical properties, as is needed for the predictive design of newly
developed materials.

With the advent and increasing prevalence of powerful computational resources,
as well as the broadening accessibility of open-source software, researchers
have been able to circumvent the need for extensive homogenization by using
high-fidelity, discrete modeling approaches such as molecular dynamics (MD)
(Bergström and Boyce, 2001, Somasi et al., 2002, Doyle and Underhill, 2005) in
frameworks such as LAMMPS (Thompson et al., 2022). Such methods can easily
accommodate dynamic bond kinetics via either deterministic activation of bonds
due to heuristic rules (Goodrich et al., 2018) or, more commonly,
incorporation of Monte Carlo methods that randomly sample the formation or
dissociation of dynamic bonds based on transition state theory (Evans and
Ritchie, 1997, Hoy and Fredrickson, 2009, Stukalin et al., 2013, Amin et al.,
2016, Perego and Khabaz, 2020, Zhao et al., 2022, Wagner et al., 2024). MD
simulations have proven valuable for linking microstructure to bulk properties
of dynamic polymers, and have recently been used to study features such as
enhanced transport (Goodrich et al., 2018, Huang et al., 2023, Taylor et al.,
2024), microstructural relaxation (Yang et al., 2015, Amin and Wang, 2020),
self-healing (Zheng et al., 2021), and polymer reprocess-ability (Zhao et al.,
2022). However, the number of discrete particles that must be modeled to
capture representative volume elements (RVE) of dynamic polymers using MD,
combined with the small vibrational timescales of said particles (on the order
of picoseconds) restricts such models to nanometer and nanosecond domains
(Agrawal et al., 2016, Liu et al., 2020, Zhang et al., 2023). These
spatiotemporal scales regularly reach $10^{2}$ nanometers and $10^{2}$
nanoseconds with the use well-established coarse-graining methods such as
bead-spring representation of Kuhn segments in a polymer chain (Kremer and
Grest, 1990, Somasi et al., 2002) or Brownian dynamics (Doyle and Underhill,
2005) to implicitly model solvent. Yet, the computational cost of such
simulations remains high and there still exists several spatiotemporal orders
of magnitude between the molecular scales of chemical structure (i.e.,
Ångstroms to nanometers and femtoseconds to nanoseconds), and the macroscale
at which experiments are mostly conducted and continuum models are commonly
applied (e.g., upwards of millimeters and milliseconds).

To address this spatiotemporal gap, many researchers have begun investigating
polymeric materials at $10^{2}$ to $10^{3}$ nanometer scales using of a class
of “discrete network models” (DNMs). In these models, crosslinks within
polymer networks are explicitly represented as “nodes”, while the constitutive
chains that link them together are captured implicitly using statistically
derived force-extension relations (Hernández Cifre et al., 2003, Sugimura et
al., 2013, Kothari et al., 2018, Wagner et al., 2021, Wyse Jackson et al.,
2022). Since DNMs avoid representing every atom or lumped Kuhn segment in a
polymer chain, they reduce the number of explicitly tracked particles and
therefore computational cost of simulated polymer networks. Vernerey and
coworkers have recently developed one such set of DNMs for dynamic polymers
that reproduce the viscoelastic mechanical predictions of Transient Network
Theory (TNT) (Vernerey et al., 2017, Wagner et al., 2021), as well as the
experimentally measured mechanical stresses and microstructural traits of both
stable (Wagner et al., 2022) and dynamically crosslinked gels (Wagner and
Vernerey, 2023). However, these DNMs neglect the non-equilibrium effects of
frictional drag due to inter-chain and solvent-chain interactions on the basis
of quasi-static loading conditions. Additionally, these DNMs implicitly model
dynamic bonding sites, often referred to as “stickers” (Leibler et al., 1991,
Mordvinkin et al., 2021), so that they cannot investigate the phenomenon of
bond lifetime renormalization put forth by Stukalin et al. (2013) whereby
stickers may break and reform bonds with the same partner multiple times
before undergoing partner exchange. Without capturing these transient effects,
it remains unclear to what extent the predictions of prior dynamic polymer
DNMs agree with those of nanoscale MD and molecular theory, as well as over
what timescales these DNMs are effectively applicable.

To address these shortcomings and explore the applicable conditions for
dynamic polymer DNMs, we here introduce an MD-consistent DNM (Fig. 1A) greatly
expanding on that of Wagner et al. (2021). This DNM, hereafter referred to as
simply the “mesoscale model”, explicitly tracks only two types of “nodes” to
represent polymers (Fig. 1B). These are (i) the anchoring crosslinks
(representing the locations at which chains are grafted to the polymer
backbones of the network), and (ii) the distal stickers that reversibly
associate with neighboring stickers. Stickers may represent either reversible
covalently bonding (Kloxin et al., 2010, Richardson et al., 2019, Yue et al.,
2020, Wagner and Vernerey, 2023) or physically interacting sites (Vidavsky et
al., 2020, Zhang et al., 2020, Xu et al., 2022, Cai et al., 2022). As in prior
DNMs, node-to-node interactions are captured via idealized implicit pairwise
bond potentials derived from statistical mechanics. We use this model to
investigate how input parameters such as polymer chain lengths, dynamic bond
activation energies, polymer concentrations, and externally applied loading
rates mediate distal sticker exploration and binding kinetics, which in turn
govern network-scale topologies and mechanical stress responses. We examine
these mappings as predicted by not only the mesoscale model, but also a
conventional bead-spring approach in which every Kuhn segment comprising a
polymer chain is explicitly modeled as a node attached to adjacent segments
via a finitely extensible, nonlinear elastic potential (Kremer and Grest,
1990, Somasi et al., 2002, Cruz et al., 2012, Sliozberg and Chantawansri,
2013) (Fig. 1C). In doing so, we identify the applicable regimes and
limitations of the mesoscale model, quantify its computational cost savings
over conventional coarse-grained MD, and explore its ability to model the
mechanics of larger material domains.

Figure 1: Discrete network modeling of dynamic polymers. (A) 3D RVE and close-
up 2D schematic of a branched polymeric network in which backbone chains are
depicted as black, while open and bonded branching side chains are red and
blue, respectively. Attachment and detachment events of stickers are
characterized by their respective rates, $k_{a}$ and $k_{d}$. (B) Schematic of
the mesoscale modeling approach. Only the crosslinks that tether side groups
to the network (grey) and stickers (red) are explicitly modeled as particles.
(C) Schematic of a conventional bead-spring model. Every Kuhn segment is
modeled as a discrete particle attached to neighboring segments by nonlinear
springs. Backbone Kuhn segments are black, branched side-chain Kuhn segments
are dark red, stable crosslinks tethering the branched chains to the network
are grey, and open stickers are red.

In the remainder of this work, we introduce the idealized bead-spring and
mesoscale modeling methods (Section 2); examine the single-chain exploration
and bond dynamics predictions, which shed light on statistical pairwise bond
association theory (Section 3); compare and contrast the network-scale
mechanical predictions and computational costs of the bead-spring and
mesoscale models (Section 4); and then demonstrate how the mesoscale approach
may be employed to probe larger spatiotemporal scales than are easily
accessible using conventional MD (Section 5).

##  2 Discrete modeling methods

Here we detail the methods of both the newly introduced mesoscale model and
the analogous bead-spring model used for validation. Both models are
implemented using the MD framework LAMMPS (Thompson et al., 2022) built with
the “Transient Network Theory” package recently introduced by Wagner et al.
(2021, 2024). Network initiation and data post-processing are conducted using
custom scripts within MATLAB 2022a. Below, we first introduce the equation of
motion used to update both models’ particles’ positions in time as a function
of the forces exerted on them (Section 2.1). We then detail the pairwise bond
association/dissociation rules used to update both models’
attachment/detachment kinetics (Section 2.2). For detailed numerical
initiation procedures of both single-chain and network-scale studies, as well
as lists of pertinent unit conversions and parameters, see Appendices A-B.

###  2.1 Equation of motion

In both the bead-spring and mesoscale models, we update the position of each
explicitly modeled particle, $\alpha\in[1,\mathscr{N}]$, in time, $t$, using
conventional Brownian dynamics. The effects of solvent are implicitly captured
by stochastic forces (that represent thermal fluctuations due to solvent
interactions), and solvent-induced drag forces (proportionate to the
particle’s velocity). Based on relatively low particle masses and high
viscosities of the surrounding mediums (i.e., solvent or adjacent polymers),
the position, $\bm{x}^{\alpha}$, of each particle $\alpha$ is updated
according to an overdamped Langevin equation of motion:

$\gamma^{\alpha}\frac{d\bm{x}^{\alpha}}{dt}=\sum_{\beta}\bm{f}^{\alpha\beta}+\sqrt{2\gamma^{\alpha}k_{b}T}\bm{\eta}^{\alpha}.$
(1)

The term on the left-hand side of Eq. (1) constitutes drag force where
$\gamma^{\alpha}$ is a friction coefficient related to node $\alpha$’s
untethered diffusion coefficient, $D^{\alpha}$, through the Einstein relation
(i.e, $D^{\alpha}=k_{b}T/\gamma^{\alpha}$) (Einstein, 1905, Rubinstein and
Colby, 2003). The first term on the right-hand side of Eq. (1) is the
unbalanced force due to all pairwise polymer interactions between node
$\alpha$ and its interaction neighbors, $\beta$. The final term on the right-
hand side of Eq. (1) captures thermal fluctuations due to solvent interactions
(Rubinstein and Colby, 2003), where $k_{b}T$ is the thermal energy,
$k_{b}=1.38\times 10^{-23}$ J/K is the Boltzmann constant, $T$ is absolute
temperature, and $\bm{\eta}$ stochastically prescribes Gaussian-distributed
thermal noise with a variance $dt^{-1}$. This noise, $\bm{\eta}^{\alpha}$,
satisfies the conditions:

$\langle\bm{\eta}^{\alpha}\rangle_{\mathscr{N}}=0,$ (2)
$\langle\bm{\eta}^{\alpha}(t_{1})\cdot\bm{\eta}^{\alpha}(t_{2})\rangle_{\mathscr{N}}=\delta(t^{\prime}-t).$
(3)

where the operator $\langle\Box\rangle_{\mathscr{N}}$ denotes ensemble
averaging over all $\mathscr{N}$ nodes. Eqs. (2) and (3) respectively convey
that the mean of $\bm{\eta}$ is zero (i.e., no net thermal force) and that
thermal fluctuation is a delta-correlated stationary process in time (i.e.,
there is no temporal correlation between the thermal fluctuation applied on a
given node at any arbitrary time $t$ and subsequent time $t^{\prime}$). Eq.
(1) may be used to explicitly estimate each particle’s discrete displacement,
$\delta\bm{x}$, over a discrete time step, $\delta t$. To achieve numerical
stability, $\delta t$ was set to $4\times 10^{-4}\tau_{0}$ for the bead-spring
model, and $4\times 10^{-3}\tau_{0}$ for the mesoscale model, where
$\tau_{0}=b^{2}/D_{0}$ is the prescribed monomer diffusion timescale
(Einstein, 1905) and $b$ is the Kuhn length characterizing the size of a mer.
Note that no equation of motion is provided for rotational degrees of freedom,
as all particles are treated as point masses, consistent with ideal chain
assumptions of no rotational penalty between bonded segments (Rubinstein and
Colby, 2003, Doi, 2013).

While both models update their nodes’ positions through Eq. (1), the pairwise
forces their nodes experience ($\bm{f}^{\alpha\beta}$) due to both bonded and
non-bonded polymer interactions must be captured distinctly. For the mesoscale
model, bonded interactions consist primarily of the entropic tensile forces of
implicit chains, resulting from their reduced conformational degrees of
freedom as they elongate. Entropic tension of the chain connecting node
$\alpha$ to its $\beta^{th}$ neighbor is computed as
$-\partial\psi_{c}/\bm{r}^{\alpha\beta}$, where $\bm{r}^{\alpha\beta}$ is said
chain’s end-to-end vector, and $\psi_{c}$, is its Helmholtz free energy. Free
energy, $\psi_{c}$, is prescribed via the Padé approximation of ideal Langevin
chains as:

$\psi_{c}=k_{b}T\left\\{\frac{(\bm{\lambda}_{c}^{\alpha\beta})^{2}}{2}-N\log\left[N-(\bm{\lambda}_{c}^{\alpha\beta})^{2}\right]\right\\},$
(4)

where $N$ is the number of Kuhn segments in the chain and
$\bm{\lambda}_{c}^{\alpha\beta}=\bm{r}^{\alpha\beta}/(\sqrt{N}b)$ is the
chain’s stretch (Cohen, 1991). To simplify the mesoscale approach, we invoke
the ideal chain assumption such that all non-bonded polymer interactions
(e.g., excluded volume repulsion, depletion forces, and long-range
interactions) may be neglected.

For the bead-spring model, the primary bonded interactions are the forces
cohering the covalent bonds comprising Kuhn segments. These are computed as
$-\partial\psi_{b}/\bm{r}^{\alpha\beta}$ where the free energy of each Kuhn
segment, $\psi_{b}$, is captured using the “nonlinear” bond style of LAMMPS
given by:

$\psi_{b}=\frac{E(\bm{r}^{\alpha\beta}-b)^{2}}{L^{2}-(\bm{r}^{\alpha\beta}-b)^{2}}.$
(5)

Here $E$ is an energy scale that modulates bond stiffness, $L$ is the finite
length by which the bond may be stretched or compressed from equilibrium, $b$
(the Kuhn length) is the finite rest length of the bond, and
$\bm{r}^{\alpha\beta}$ remains the bond’s end-to-end vector (Rector et al.,
1994). We found that setting $E=800k_{b}T$, and $L=b$ provides ample numerical
stability while accurately reproducing the theoretically predicted end-to-end
distributions of entropic chains (Section 3.1), Kuhn segment length
distributions within $\pm 10\%$ of the prescribed rest length ($b$), and
force-extension relations in agreement with
$\bm{f}(\bm{r})=-\partial\psi_{c}/\partial\bm{r}$ (Appendix C). Dynamic bonds,
which form at the length scale $b$, were also modeled using Eq. (5) for both
modeling approaches, but with $E$ lowered to $100k_{b}T$ to bolster numerical
stability. Careful consideration should be given when prescribing $L$ and $E$
for material-specific dynamic bond types in future applications of this
method. To represent non-bonded, excluded volume interactions in the bead-
spring model, we included a Lennard-Jones potential of the form:

$\psi_{LJ}=4k_{b}T\left[\left(\frac{\sigma_{0}}{d^{\alpha\beta}}\right)^{12}-\left(\frac{\sigma_{0}}{d^{\alpha\beta}}\right)^{6}\right],\;(d^{\alpha\beta}\leq
d_{c})$ (6)

where $\sigma_{0}=2^{-1/6}b$ is the equilibrium distance between two
interacting nodes, and $d^{\alpha\beta}$ is the distance between neighboring
nodes $\alpha$ and $\beta$ within cutoff distance $d_{c}=2\sigma_{0}$ of each
other.

In addition to their differences between pairwise polymer interactions, the
bead-spring and mesoscale approaches must also have different damping
coefficients, $\gamma^{\alpha}$, prescribed to each node. The damping
coefficient prescribed to a particle in the bead-spring model must account for
only one Kuhn segment and may therefore be approximated as
$\gamma_{0}=k_{b}T/D_{0}$, where $D_{0}$ is the diffusion coefficient of an
untethered monomer.111Here, $D_{0}$ is taken as a sweeping parameter using
reasonable values based on experimental literature (Shimada et al., 2005,
Kravanja et al., 2018, Shi, 2021) (see Appendix B). Through the prescription
of $D_{0}$ and $b$, the characteristic time, $\tau_{0}=b^{2}/D_{0}$, it takes
a monomer to diffuse distance $b$, is established. This, in turn, dictates the
timescale of the model. In contrast, the coefficient of a node in the
mesoscale model must also account for the friction of adjacent Kuhn segments
in its attached chain(s). Rouse scaling theory predicts that the frictional
coefficient of a tethered chain with $N$ Kuhn segments scales as
$\gamma^{\alpha}\approx N\gamma_{0}$, (Rouse, 1953, Rubinstein and Colby,
2003). We find that in networks of freely diffusing chains, this relationship
holds (Sections 4 and 5). However, in studies in which one end of a tethered
chain is fixed, we instead find that the mesoscale model most closely
approximates the Rouse sub-diffusion of the bead-spring model when
$\gamma^{\alpha}=N^{2/3}\gamma_{0}$ for the mesoscale nodes (Section 3).

###  2.2 Bond kinetics

Dissociation between attached stickers is prescribed according to Eyring’s
theory (Eyring, 1935) at a rate of:

$k_{d}=\tau_{0}^{-1}\exp\left(-\frac{\varepsilon_{d}}{k_{b}T}\right),$ (7)

where $\tau_{0}=b^{2}/D_{0}$ is the characteristic time it takes a monomer
(i.e., detached sticker) to diffuse its Kuhn length, $b$, and
$\varepsilon_{d}$ is the activation energy of bond dissociation (Eyring,
1935). Analogously to dissociation, bond attachment is prescribed with an
intrinsic rate:

$k_{a}=\tau_{0}^{-1}\exp\left(-\frac{\varepsilon_{a}}{k_{b}T}\right),$ (8)

where $\varepsilon_{a}$ is the activation energy of association. However,
unlike dissociation, bond association is intrinsically predicated on the
ability of detached stickers to encounter one another through their stochastic
diffusion, and so will inherently be influenced by traits such as the open
sticker concentration (Stukalin et al., 2013) and tethered chain length as
discussed in Section 3.1. Nonetheless, assuming no long-range interactions
(Rubinstein and Colby, 2003), encounters are defined as occurring when two
stickers diffuse within one Kuhn length, $b$, of each other (Stukalin et al.,
2013). To determine if an attached set of stickers dissociates, or an
“encountering” pair of open stickers associates, a memoryless Poisson process
is assumed so that the probability of reaction, $P$, evolves in time according
to:

$dP_{i}=1-\exp\left(-k_{i}t\right)dt$ (9)

where the index “$i$” denotes either attachment, “$a$”, or detachment, “$d$”
(Wagner et al., 2021, Wagner and Vernerey, 2023). Integrating Eq. (9) over the
discrete time interval $[t,t+\delta t]$ and checking the resulting probability
against a random number in the uniformly distributed range 0 to 1, permits the
capture of stochastic reactions in both models.

Average effective attachment and detachment rates over the total simulated
time (from $t=0$ to $t=t_{f}$) are computed as:

$\bar{k}_{a}=t_{f}^{-1}k_{s}\int_{0}^{t_{f}}N_{a}\left[(c-c_{a})V\right]^{-1}dt,$
(10)

and:

$\bar{k}_{d}=t_{f}^{-1}k_{s}\int_{0}^{t_{f}}N_{d}(c_{a}V)^{-1}dt,$ (11)

respectively. Here, $k_{s}$ is the sampling frequency; $N_{a}$ and $N_{d}$ are
the discrete numbers of attachment and detachment events at time $t$,
respectively; $c$, $c_{a}$, and $c-c_{a}$ are the total, attached, and
detached chain concentrations at time $t$, respectively; and $V$ is the domain
volume. To ensure adequate temporal resolution, $k_{s}$ was set to twenty
times that of the prescribed monomer oscillation frequency (i.e.,
$k_{s}=20\tau_{0}^{-1}$). We found that increasing $k_{s}$ from
$15\tau_{0}^{-1}$ to $20\tau_{0}^{-1}$ influenced neither the diffusive
behavior of tethered chains nor the binding kinetics between them for either
discrete modeling approach (see Appendix D).

##  3 Single-chain and bond kinetics validation studies

Before utilizing the mesoscale model for network-scale mechanical predictions,
we must first verify that the single-chain statistics of the mesoscale model
agree with both the predictions of the bead-spring model, as well as
prevailing statistical theory for ideal chains (Rubinstein and Colby, 2003,
Doi, 2013). Additionally, we must ensure that the bond association,
dissociation, and partner exchange kinetics predicted by the mesoscale model
reasonably agree with those predicted by the bead-spring model. Model
validation results for single chain exploration, pairwise bond kinetics, and
partner exchange kinetics in ensembles of chains are presented in Sections
3.1-3.3 below.

###  3.1 Tethered chains follow Gaussian statistics and approximate Rouse sub-
diffusion

We conducted tethered single-chain studies to corroborate that individual
chains’ end-to-end conformations and tethered diffusion characteristics agree
between models. To explore the effects of chain length, the number of Kuhn
segments was swept over $N=\\{12,18,24,30,36\\}$. The lower limit of $N=12$
was set adequately high to still observe Gaussian statistics. Meanwhile, the
upper limit of $N=36$ was set to three times the lower limit to ensure that a
comparably high chain length was explored without modeling chains
significantly longer than the entanglement length ($N=35$) predicted by Kremer
and Grest (1990). To also explore the effects of modulating the monomer
diffusivity (which governs the timescale of the model per the relation
$\tau_{0}=b^{2}/D_{0}$), the diffusion coefficient was swept over
$D_{0}=\\{0.125,0.25,0.5,1,2,4,8\\}\times 10^{-10}$ m2 s-1. The central value
of $D_{0}=10^{-10}$ m2 s-1 was selected based on a realistic value for
diffusion of poly(ethylene glycol) oligomers with low degrees of
polymerization (on the order of two mers, which corresponds to one Kuhn
segment) at ambient temperatures in good solvent (Shimada et al., 2005).

To achieve adequate statistical sampling, ensembles of $n_{p}=1331$ and
$n_{p}=125$ non-interacting polymer chains were generated for the mesoscale
and bead-spring models, respectively. The significantly larger ensemble of
chains for the mesoscale model was set arbitrarily high and was easily enabled
by the model’s reduced computational cost, while the sample size for the bead-
spring model was set adequately high to observe convergence in predicted
results. Chains were generated by first initializing their fixed tethering
nodes in a 3D grid at the position set $\\{\bm{x}^{0}\\}$.222For single-chain,
tethered diffusion studies, no inter-chain bond kinetics between the chains’
free ends were included so that the initial spacing between the positions
$\\{\bm{x}_{0}\\}$ is arbitrary. However, to leverage the spatial
parallelization of LAMMPS for higher throughput sampling, the tethering sites
were separated on their grid by a distance of $6Nb$. Simulated polymer chains
were then randomly generated from each tethering site using random-walk
generation procedures according to Eqs. (A1) and (A3) for the bead-spring and
mesoscale iterations of the model, respectively. Once all chains were
initiated, the positions of all non-fixed nodes were updated according to Eq.
(1), and then the end-to-end vectors of all $m\in\left[1,n\right]$ chains were
measured over time. End-to-end vectors are defined as:

$\\{\bm{r}_{m}\\}=\sum_{\alpha}^{\Lambda}\\{\bm{x}_{m}^{\alpha+1}-\bm{x}_{m}^{\alpha}\\},$
(12)

where the maximum node number, $\Lambda$, for each molecule is the number of
Kuhn segments ($\Lambda=N$) for the bead-spring model and two ($\Lambda=2$)
for the mesoscale model (see Fig. 2A).

Fig. 2B-C depict the probability distribution functions (PDFs) of finding a
chain at a given end-to-end stretch, $\lambda_{c}=r/(\sqrt{N}b)$, as predicted
by the mesoscale model (grey histogram), the bead-spring model (black
diamonds), and the analytically derived joint PDF of Eq. (A4) (black curve)
for Gaussian chains with $N=12$ (Fig. 2B) and $N=36$ (Fig. 2C) Kuhn segments.
Recall that stretch is defined as the end-to-end length, $r$, of a chain
normalized by its mean expected length, $\sqrt{N}b$, (as predicted by
Gaussian, random-walk statistics) so that the functional form of the joint PDF
of Eq. (A4) with respect to $\lambda_{c}$ is:

$P(\lambda_{c})=4\pi\left(\frac{3}{2\pi
Nb^{2}}\right)^{\frac{3}{2}}\lambda_{c}^{2}\exp\left(-\frac{\lambda_{c}^{2}}{2}\right),$
(13)

whose variance is unity and is thus independent of $N$ or $b$ (Doi, 2013).
Therefore, providing end-to-end distributions in terms of stretch induces
collapse of the PDFs with respect to the two chain lengths for a normalized
comparison. Note that probabilities of Fig. 2B-C are also re-normalized as
$p=P(\lambda_{c})/\int P(\lambda_{c})d\lambda_{c}$ so that the PDF integrates
to unity over the range $\lambda_{c}\in[0,\infty)$ and thus aligns vertically
with the discrete distributions. Results indicate that the mesoscale model’s
predicted end-to-end distributions are in excellent agreement ($R^{2}\geq
0.99$) with the analytically derived Gaussian PDF. Additionally, the bead-
spring model agrees with both the mesoscale and analytical models when the
characteristic energy scale from Eq. (5) is set to $E=800k_{b}T$.

Figure 2: Tethered diffusion characteristics. (A) Illustration of a tethered
chain, anchored to a fixed node (grey). The sticker (red) diffuses over the
time interval $[t_{0},t]$, from positions $\bm{x}^{\Lambda}_{m}(t_{0})$ to
$\bm{x}^{\Lambda}_{m}(t)$. (B,C) PDFs of chain end-to-end stretch,
$\lambda_{c}$, predicted by the mesoscale model (grey histogram, $n_{p}=1331$
chains), bead-spring model (discrete diamonds, $n_{p}=125$ chains), and
analytical Gaussian distribution of Eq. (13) (black curve), when $N=12$ and
$N=36$. (D) Normalized MSD, $\langle\Delta r^{2}\rangle/b^{2}$, versus time
for $N\in\\{12,18,36\\}$ as predicted by the bead-spring (discrete data) and
mesoscale (solid curves) models. Total time, $t\approx 600\tau_{0}$, is above
the Rouse time, $\tau_{r}$ so that steady state is reached for all $N$. (E)
MSD with respect to time at the order of the Rouse timescale, $0<t<5\tau_{r}$,
for $N\in\\{12,18,36\\}$, predicted by the bead-spring (discrete data) and
mesoscale (solid curves) models. Inset depicts bead-spring data (discrete
data) for $0<t<\tau_{r}$ along with fits from Eq. (15) (dashed curves). (F)
Relative error between the bead-spring and mesoscale models (bead-spring model
as reference) with respect to time. Results of (B-F) are for the median value
of $D_{0}=10^{-10}$ m2 s-1 only, as $D_{0}$ has no meaningful effect on end-
to-end conformations or MSD with respect to normalized time (Fig. E1). All
error bars represent standard error of the mean (SE).

While the PDFs of Fig. 2B-C indicate excellent agreement between the
instantaneous ensemble of chain conformations from both discrete models and
statistical mechanics, they reveal nothing of the exploratory behavior of
individual chains as they diffuse through space. To characterize exploration
in time, we also compute the mean-square displacement (MSD) of the chains’
distal ends according to:

$\langle\Delta
r^{2}(t)\rangle=\langle(\bm{x}_{m}^{\Lambda}(t)-\bm{x}_{m}^{\Lambda}(t_{0}))^{2}\rangle$
(14)

where $\Lambda$ denotes the freely diffusing node at the end of the chain,
$t_{0}$ is the reference time from which MSD is measured, and the operator
$\langle\square\rangle$ denotes ensemble averaging over all $m\in[1,n_{p}]$
chains. Fig. 2D compares the MSD measured for stickers in the mesoscale model
(continuous curves) to those of the bead-spring model (discrete data) over a
duration of $t\approx 600\tau_{0}$. Over this relatively longer timescale, the
MSDs plateau at values that agree between models, regardless of diffusion
coefficient (Appendix E, Fig. E1). However, for the eventual purposes of
investigating bond kinetics that are mediated by the frequency at which two
binding sites enter within short distances (i.e., $b=\sqrt{D_{0}\tau_{0}}$) of
each other, it is paramount to observe exploration behavior over shorter
timescales and smaller spatial areas.

Fig. 2E displays the same MSD data from Fig. 2D, but over the much shorter
timescale of $0<t<5\tau_{r}$ where $\tau_{r}$ is the Rouse time, or the time
it takes a polymer chain to diffuse its own characteristic area, $Nb^{2}$. It
is well established that in the regime $\tau_{0}<t<\tau_{r}$, tethered chains
explore 3D space according to:

$\langle\Delta r^{2}(t)\rangle\approx
b^{2}\left(\frac{t}{\tau_{s}}\right)^{1/2}$ (15)

where $\tau_{s}=\tau_{r}N^{-2}$ is the time it takes a distal sticker at the
end of a chain of $N$ Kuhn segments to diffuse its own characteristic size,
$b^{2}$ (Hult et al., 1999, Rubinstein and Colby, 2003, Stukalin et al.,
2013). Observing Fig. 2E, we see that the mesoscale and bead-spring models are
in relatively good agreement at timescales on the order of $1<t<5\tau_{r}$
(with less than 10$\%$ relative error per Fig. 2F). From the inset of Fig. 2E,
we also see that the bead-spring model is in excellent agreement with the
Rouse model ($R^{2}>0.95$) over the time range $0<t<\tau_{r}$ when the Rouse
time, $\tau_{r}$, is taken as the time at which measured MSD first exceeds
$Nb^{2}$ and the sticker diffusion timescale, $\tau_{s}$, is treated as a
fitting parameter (Fig. E2).333We find $\tau_{s}/\tau_{0}\approx 0.07$,
regardless of $N$ or $D_{0}$ (Fig. E1.H). However, below the Rouse timescale
the mesoscale model under-predicts the MSD predicted by the bead-spring and
Rouse models, with a relative error upwards of 50$\%$ (Fig. 2F). This high
magnitude of relative error is partially due to the fact that MSD approaches
zero when $t\rightarrow 0$. Yet, that the relative error is always negative,
indicates that the mesoscale model consistently under-predicts distal sticker
diffusion over short timescales.

Underprediction of the MSD by the mesoscale model, as compared to the bead-
spring model, could be due to a disagreement in radial MSD (i.e., MSD due to
displacement along the end-to-end direction of the chain), tangential MSD
(i.e., MSD due to displacements normal to the end-to-end direction of the
chain), or some combination of both. The radial and tangential components of
MSD are respectively defined as mean-square change in end-to-end length:

$\langle\Delta
r_{r}^{2}(t)\rangle=\langle\left[r_{m}(t)-r_{m}(t_{0})\right]^{2}\rangle,$
(16)

and mean-square circumferential distance swept by the distal end of the chain:

$\langle\Delta
r_{t}^{2}(t)\rangle=\langle\left[r_{m}(t_{0})\theta_{m}\right]^{2}\rangle,$
(17)

over the time interval $t\in[t_{0},t]$. Here,
$\bm{r}_{m}(t)=\bm{x}^{\Lambda}_{m}(t)-\bm{x}^{0}_{m}(t)$ and
$\bm{r}_{m}(t_{0})=\bm{x}^{\Lambda}_{m}(t_{0})-\bm{x}^{0}_{m}(t_{0})$ are the
end-to-end vectors of chain $m$ at times $t$ and $t_{0}$, respectively; and
$\theta_{m}=\cos^{-1}\left[\frac{\bm{r}_{m}(t_{0})\cdot\bm{r}_{m}(t)}{r_{m}(t_{0})r_{m}(t)}\right]$
is the angle between the end-to-end vectors at time $t_{0}$ and $t$ (see Fig.
3A). Examining Fig. 3B-C, we see that the cumulative radial and tangential
MSDs of the mesoscale model closely mirror those of the bead-spring model,
rarely deviating by more than 10$\%$. However, over the time period at which
dynamic bonding is sampled in future studies ($\tau_{0}/20$), the mesoscale
stickers explore approximately $0.95\pm 0.10$ $b^{2}$ and $2.17\pm 0.06$
$b^{2}$ less than their bead-spring counterparts in the radial and tangential
directions, respectively (Fig. 3D-E). This is true regardless of chain length.
The agreement of long-term ($t\geq\tau_{r}$) exploratory behavior between
models supports that the mean paths of the stickers in both models traverse
similar characteristic trajectories. However, the discrepancy in MSD below the
Rouse time ($t<\tau_{r}$) indicates that the stickers of the bead-spring model
do so with a higher vibrational amplitude on the order of $b^{2}$, and this is
true of both their tangential and end-to-end vibration modes.

Notably, the error is more pronounced for the tangential diffusion mode,
meaning that there are greater deviations between the models’ vibrational
sticker amplitudes in directions normal to chain end-to-end vectors than in-
line with them. This is likely because radial movement is constrained so that
magnitudes of radial MSD are smaller than those of tangential MSD (Fig. 3B-C).
However, it may also arise from the fact that tangential diffusion is mediated
entirely by drag and Brownian forces on the stickers. Both of these depend on
the damping coefficient, which must be set higher for the mesoscale model to
achieve similar MSDs past $\tau_{r}$. In contrast, radial diffusion is also
checked by the entropic chain forces, which are independent of the damping
coefficient and in good agreement between models (Fig. C1.B), thus perhaps
reducing error. In any case, to ensure that discrepancies in diffusive
behavior do not affect bond kinetics, or – by extension – topological network
reconfiguration and network mechanics, we next investigate the emergent bond
attachment and detachment rates of these two models as not only a function of
their intrinsically assigned kinetic rates through Eqs. (7) and (8), but also
the exploratory behavior of their chains.

Figure 3: Radial and tangential tethered diffusion. (A) Illustration of a
tethered chain anchored to a fixed node (grey) whose sticker (red) diffuses
from end-to-end vector $\bm{r}_{m}(t_{0})$ (grey arrow) to $\bm{r}_{m}(t)$
(black arrow). Geometric definitions of $\Delta r_{t}(t)$ (dashed arc) and
$\Delta r_{r}(t)$ (dashed arrow) are illustrated. (B) Radial and (C)
tangential components of MSD (top) with respect to time, $t/\tau_{r}$, for the
bead-spring model (discrete data) and mesoscale model (solid curves). Relative
error (bottom) between models (bead-spring model serves as reference). (D)
Radial and (E) tangential components of average square displacement (top) over
the duration $[t,t+\tau_{0}/20]$, with respect to normalized time for the
bead-spring model (discrete data) and mesoscale model (solid curves). Absolute
error (bottom) between models (bead-spring model serves as reference). (B-E)
Error bars represent SE.

###  3.2 Chain attachment mirrors Bell’s model for detachment

While the kinetic association rate, $k_{a}^{ap}$, prescribed a priori between
two stickers within distance $b$ of each other is set according to Section 2.2
(see Appendix F and Fig. F1 for validation), the actual rate of attachment
must also depend on the probability, $P_{e}$, with which said stickers
encounter one another such that the emergent attachment rate is $k_{a}\propto
k_{a}^{ap}P_{e}$. However, an analytical form of this probability is not
immediately available from the literature. Therefore, to investigate how this
probability evolves in each of the two discrete models while also
interrogating agreement between their associative kinetics, we conducted
studies in which $n_{p}=216$ tethered chains with stickers at their distal
ends were positioned near fixed stickers. The distal stickers and fixed
stickers were then allowed to bond/unbond with each other in mutually
exclusive pairs (Fig. 4A-C). To probe the effects of distance, the separation
length, $d_{ts}$, between each chain’s tethering site and its paired fixed
sticker was swept over the range $d_{ts}\in[0.0125,0.5]Nb$. Neighboring pairs
were separated from each other by a distance greater than $2Nb$, so that no
two chains’ stickers could come within bonding distance, $b$, of each other.
Instead, each tethered chain was only within reach of the fixed sticker
belonging to its pair.

Once all fixed nodes were positioned, the chains were initiated as described
in Section A.1 through Eqs. (A1) to (A3). After initiation, all non-fixed
nodes’ positions were updated according to Eq. (1). Stickers that came within
distances less than $b$ of each other were checked for attachment according to
Eqs. (8-9) and the methods of Section 2.2. The time-averaged rates of
attachment and detachment were then computed according to Eqs. (10) and (11),
respectively. Besides sweeping the separation distance, $d_{ts}$, the number
of Kuhn segments, $N$, and associative activation energies, $\varepsilon_{a}$,
were also swept over the ranges $N=\\{12,18,36\\}$ and
$\varepsilon_{a}=\\{0.01,0.1,1\\}k_{b}T$ to elucidate the effects of chain
lengths and intrinsically prescribed binding rates, respectively. The upper
limit of $\varepsilon_{a}=k_{b}T$ was selected because over computationally
viable time domains (on the order of $10^{3}\tau_{0}$ to $10^{5}\tau_{0}$ for
the bead-spring model with longer chains) we found that the number of discrete
attachment events no longer provided adequate statistical sampling sizes when
$\varepsilon_{a}\sim 10k_{b}T$. Meanwhile, the lower limit of
$\varepsilon_{a}=0.01k_{b}T$ was selected because, at this activation energy
scale, the intrinsic attachment rate approaches the monomer diffusion
frequency ($k_{d}^{ap}\rightarrow\tau_{0}^{-1}$) so that sampling lower
activation energies ($\varepsilon_{a}<0.01k_{b}T$) becomes redundant. The bond
kinetics sampling rate, $k_{s}$ and data output frequency were both set to
$20\tau_{0}^{-1}$.

Fig. 4D-E indicates that the ensemble-averaged attachment rates,
$\bar{k}_{a}$, measured from both the bead-spring and mesoscale models are in
reasonable agreement with one another across separation distances (here
characterized by the chain stretch, $\lambda_{c}=d_{ts}/(\sqrt{N}b)$, required
for attachment), chain lengths (through $N$)444Results from three unevenly
spaced values of $N$ are presented based on an observed nonlinear relation
between $k_{a}$ and $N$, which reveals that as $N$ increases the sensitivity
of $k_{a}$ to $N$ decreases., and for two disparate activation energies
($\varepsilon_{a}=\\{0.01,1\\}k_{b}T$). As expected, increasing the bond
activation energy reduces the emergent attachment rate as indicated by the
lower values of $\bar{k}_{a}$ from Fig. 4E when compared to those of Fig. 4D.
Intuitively, increasing chain length (by increasing $N$) diminishes the
attachment rate. This is attributed to the fact that longer chains have larger
available exploration volumes and therefore are statistically less likely to
encounter neighboring sticker sites at any given moment. Note that the models’
predicted values of $\bar{k}_{d}$ are consistently in excellent agreement with
the value set a priori through Eq. (7), which remains constant with respect to
$\lambda_{c}$ since $k_{d}$ is intentionally made independent of chain stretch
in Eq. (7).

Figure 4: Bond kinetics of a single tethered chain. (A-C) Illustration of a
tethered chain whose free end is a sticker (red) that may bind or unbind with
a fixed sticker (fixed red node to right). Distance between the fixed
tethering node (grey) and fixed sticker (red) is denoted $d_{ts}$. This
sticker pair is shown (A) in an arbitrary detached initial state, (B)
immediately after an attachment event, and (C) after the subsequent detachment
event. (D-E) Average normalized attachment rates, $\bar{k}_{a}\tau_{0}$, with
respect to normalized separation distance, $\lambda_{c}=d_{ts}/(\sqrt{N}b)$
when the associative bond energy is (D) $\varepsilon_{a}=0.01k_{b}T$ and (E)
$\varepsilon_{a}=k_{b}T$. Data is provided for both the bead-spring (circles)
and mesoscale (triangles) models. Error bars represent SE. Best fits of Eq.
(19) for the bead-spring (solid curves) and mesoscale (dashed curves) data are
also displayed, treating prefactor, $A$, as a fitting parameter. Data is also
shown for the LAMMPS implementation of Eq. (19) when $A$ is taken as the
average value from the fits to the bead-spring and mesoscale models. (F)
Prefactor, $A$, for all activation energies and chain lengths. Error bars
represent the 95$\%$ confidence interval. (F) Goodness of fit between the
discrete models and scaling theory, characterized by $R^{2}$.

Interestingly, the average attachment rates of Fig. 4D-E follow Gaussian
relations with respect to separation distance, $d_{ts}$. Based on this
observation, we postulate that the encounter probability, $P_{e}$, between the
two stickers is proportionate to the probability, $P(d_{ts})$, of finding a
Gaussian chain at end-to-end length $d_{ts}$ through Eq. (A4). We also
logically postulate that $P_{e}$ scales directly with the characteristic
volume, $b^{3}$, within which a sticker “encounters” and may therefore bind to
a neighbor. Thus, the encounter rate scales as:

$P_{e}\propto b^{3}\left(\frac{2}{3}\pi
Nb^{2}\right)^{-\frac{3}{2}}\exp\left(-\frac{3}{2}\frac{d_{ts}^{2}}{Nb^{2}}\right).$
(18)

Substituting Eq. (18), along with the definition of $k_{a}^{ap}$ through Eq.
(8), into the postulated relation $k_{a}\propto k_{a}^{ap}P_{e}$, writing the
relation in terms of chain stretch, $\lambda_{c}=d_{ts}/(\sqrt{N}b)$, and then
simplifying predicts that the emergent attachment rate scales as:

$k_{a}\propto\tau_{0}^{-1}\left(\frac{2}{3}\pi
N\right)^{-\frac{3}{2}}\exp\left[-\frac{\varepsilon_{a}+\psi_{a}(\lambda_{c})}{k_{b}T}\right],$
(19)

where $\psi_{a}(\lambda_{c})\approx\frac{3}{2}{k_{b}T}\lambda_{c}^{2}$ is the
Helmholtz free energy of an ideal chain at relatively low end-to-end
stretches, $\lambda_{c}\sim 1$ (i.e., within the Gaussian regime).

Significantly, Eq. (19) predicts that the attachment rate scales with the
inverse exponent of not only the intrinsic bond activation energy,
$\varepsilon_{a}$, but also the Helmholtz free energy, $\psi_{a}$, of the
entropic chain that would exist if attachment occurred. While we here
prescribe force-independent bond dissociation for simplicity, Eq. (19) also
mirrors Bell’s model for force-dependent slip-bond dissociation (Bell, 1978)
which states that the rate of bond detachment is given by:

$k_{d}^{Bell}=K\exp\left[-\frac{\varepsilon_{d}-\psi_{a}(\lambda_{c})}{k_{b}T}\right]$
(20)

where $\psi_{a}$ is the Helmholtz free energy of the already attached chain
and the prefactor $K$ is an attempt frequency generally taken as
$\tau_{0}^{-1}$ (Evans and Ritchie, 1997, Song et al., 2021, Wagner and
Vernerey, 2023). In essence, Eq. (19) predicts that the stretch-dependent
Helmholtz free energy of a polymer chain additively increases the effective
energy barrier of attachment for a binding site located at its distal end.
Meanwhile, Eq. (20) states that the same stretch-dependent Helmholtz free
energy subtractively decreases the effective energy barrier for bond
dissociation. These relations between kinetic rates and single-chain Helmholtz
free energies emerge from statistical mechanics when forward and reverse
reaction rates are functionally derived using the end-to-end state and
transition state partition functions of polymer chains (Buche and Silberstein,
2021). Furthermore, they are consistent with the experimental and theoretical
findings of investigators such as Guo et al. (2009) or Bell and Terentjev
(2017) who studied the association of polymer-tethered ligands to fixed
receptor sites.

The theoretical scaling predictions from Eq. (19) are also in good agreement
with both models investigated here if we multiply $k_{a}$ from Eq. (19) by a
dimensionless prefactor, $A$, which is treated as a fitting parameter (see
Fig. 4F). Values of $A$ fitted to the mesoscale model are in agreement with
those fitted to the bead-spring model for every combination of activation
energy and chain length investigated (Fig. 4F) further illustrating agreement
between the two discrete approaches. While values of $A$ range from $0.53\pm
0.08$ to $1.45\pm 0.23$, they are generally less than unity suggesting that
Eq. (19) tends to over-predict measured attachment rates from the discrete
approaches. Furthermore, there are consistent trends whereby the largest chain
length and the highest activation energy correspond to greater values of $A$.
The former trend may indicate that the scaling proportionality $k_{a}\propto
N^{-3/2}$, which states that $k_{a}$ is proportionate to the cubic inverse of
the maximum allowable stretch, $\lambda_{c}^{max}=N^{1/2}$, is incomplete in
Eq. (19). Alternatively, it may suggest that there exists some unknown source
of error (common to both discrete models), which is more pronounced for
shorter chain lengths due to their higher sensitivity of $k_{a}$ to $N$.
Meanwhile, the latter trend (that $A$ increases to above unity as the
$\varepsilon_{a}$ increases) indicates that for systems with higher activation
energy, the scaling theory trends towards under-predicting the attachment
rate. We hypothesized that this trend results from bond-history dependent
sticker exploration. Namely, that bonds that recently detached are
statistically more likely to undergo repeat attachment with the same neighbors
over shorter timescales, a notion put forth by Stukalin et al. (2013) and
investigated as it pertains to partner exchange rates in Section 3.3. However,
independently plotting the rate of first-time attachment rates,
$\bar{k}_{a,1}$, versus overall attachment rates, $\bar{k}_{a}$, (including
repeat events) reveals little change in the observed kinetics for the
simulation times here (Fig. G1).

While the exact physics necessitating the inclusion of prefactor, $A$, remain
unclear, no fitting is required to match the variance of Eq. (19) (see also
Eq. (A4)). Indeed, the high $R^{2}$ values of Fig. 4G (all $>0.9$) confirm
that the scaling theory can be fit appropriately across all activation
energies and chain lengths simply by modulating $A$. This is true irrespective
of chain length or activation energy (Fig. 4D-E). This also remains true when
a set of two tethered chains capable of binding at their distal ends are
modeled within various tether-to-tether distances, $d_{tt}$, of one another
(Fig. G2) and allowed to bond/unbond. In this case, the number of Kuhn
segments in the would-be chain after attachment is simply $2N$ so that the
Helmholtz free energy in terms of end-to-end (i.e., tether-to-tether) distance
becomes $\psi_{a}=\frac{3}{4}\frac{k_{b}T}{Nb^{2}}d_{tt}^{2}$, but remains
$\psi_{a}=\frac{3}{2}k_{b}T\lambda_{c}^{2}$ in terms of chain stretch per Eq.
(19). Substituting this free energy into Eq. (19) and then fitting only $A$
nicely reproduces the average two-chain, pairwise attachment rates measured
from both discrete models (Fig. G2). Importantly, this lack of needed
adjustment to the variance confirms the key interpretation of Eq. (19) that
the rate of attachment is penalized by the Helmholtz free energy that a chain
must attain in order to stretch sufficiently for bond attachment.

Evidently, Eq. (19) provides additional means for coarse-graining whereby only
the backbone sites (at which side chains are grafted into a polymer network)
are explicitly modeled. Then, the bonds formed between tether sites by the
telechelic association of distal stickers can be implicitly captured as
energetic potentials via Eq. (4) that act with probabilities set through Eqs.
(9) and (19). Furthermore, the maximum bond interaction length may be set to
the maximum stretch, $\lambda_{c}\approx 1.5$, at which chains observably
associate (observe Fig. 4D-E and Fig. G2.D for one- and two-chain systems,
respectively). We gauged this prospective coarse-graining method by
implementing Eq. (19) into LAMMPS and conducting an analogous study in which
only the fixed nodes are modeled and the attachment between them is checked
through Eq. (19). This study reproduces the Gaussian relation between $k_{a}$
and $\lambda_{c}$ observed for the bead-spring and mesoscale model predictions
(Fig. 4D-E, discrete diamonds), with considerable reduction in computational
cost (Table G1). Since it circumvents the need to track stickers’ oscillations
within and outside of distance $b$ from one another, this implicit approach
allows for dynamic bonding check frequencies, $k_{s}$ on the order of chain
oscillation rates (i.e., $k_{s}\sim\tau_{r}^{-1}\sim N^{-2}\tau_{0}^{-1}$)
instead of $k_{s}\sim 20\tau_{0}^{-1}$. In this study, $k_{s}$ is arbitrary
since the tethers are fixed, but here we set $k_{s}=\tau_{0}^{-1}$ yielding
two to three orders of magnitude reduction in computational run time as
compared to the bead-spring model, and one to two orders of magnitude
reduction as compared to the mesoscale model (Table G1), without observed
deviation in emergent kinetic rates on a pairwise basis (Fig. 4D-E).

Despite the significant cost savings of this scaling theory-based approach, it
is limited in a number of ways. First, it requires calibration of $A$ using
either the bead-spring or mesoscale model. Second, $k_{s}$ is still limited by
the Rouse diffusion rate, $\tau_{r}^{-1}$, it takes for chains to migrate
within or outside of $\sim 1.5\sqrt{N}b$ of each other. Finally, this implicit
method loses information about the exact position of open stickers. Therefore,
it is fundamentally unable to reproduce predictions from the bead-spring and
mesoscale models about partner exchange or the phenomena of renormalized bond
lifetime as discussed in Section 3.3. Implementation of Eq. (19) should
therefore be relegated to modeling systems with dilute dynamic bond
concentrations for which partner exchange is rare; stickers should otherwise
be explicitly modeled.

###  3.3 Bond dynamics and partner exchange kinetics agree between discrete
models

Having confirmed good agreement between the bead-spring and mesoscale models’
predictions of pairwise bond kinetics for single and two-chain systems, we
next investigate whether the complex bond kinetics emerging within ensembles
of chains agree between models. Once again, we investigate the rates of bond
attachment and detachment. However, attachment events for an ensemble of
chains (e.g., a network) at equilibrium can be further sub-divided into two
types. The first type is repeat attachment whereby two stickers dissociate and
then reattach (without finding new partners during the interim) after some
time $\tau_{rpt}$ (Fig. 5, from $t_{1}$ to $t_{2}$). The second type is
partner exchange whereby stickers bond to new partners after detaching from
old ones after some comparatively longer detached lifetimes $\bar{\tau}_{exc}$
(Fig. 5, from $t_{3}$ to $t_{4}$). Based on this distinction, Stukalin et al.
(2013) introduced the renormalized bond lifetime, $\bar{\tau}_{rnm}$, or the
time from when a pair of stickers newly bonds (Fig. 5, $t_{0}$) to when one or
both of the stickers in the partnership attaches to a new neighbor (Fig. 5,
$t_{4}$). Evidence suggests that it is the renormalized bond lifetime or
inverse exchange rate (not just the detachment rate) that dictates the
reconfigurational stress relaxation time in dynamic polymers (Stukalin et al.,
2013).

Figure 5: Illustration of bond lifetimes. A schematic of four tethered chains
(labeled “i”-“l”), whose anchoring nodes (grey) are configured in a fixed grid
(with nearest neighbor separation distance, $d$), displays detachment, repeat
attachment, and partner exchange events. The definitions of attached bond
lifetime, $\tau_{a}$, detached bond lifetime prior to repeat attachment,
$\tau_{rpt}$, detached bond lifetime prior to partner exchange, $\tau_{exc}$,
and renormalized bond lifetime, $\tau_{rnm}$ are also illustrated.

To probe partner exchange kinetics and investigate agreement between the
models, we simulated $n_{p}=343$ tethered polymer chains with stickers at
their free ends. Their fixed ends were positioned in a cubic lattice within a
periodic RVE. Prior work has demonstrated that sticker concentrations, $c$,
significantly influences exchange kinetics (Stukalin et al., 2013). To
investigate this effect, the lattice constant (i.e., separation distance,
$d=c^{-1/3}$, between tether nodes) was swept, which is equivalent to the
average separation distance between neighboring stickers, $\bar{d}$, in an
ensemble of isotropically oriented chains. To ensure realistic values of $d$
for a given chain length, $d$ was set based on the prescribed chain length,
via $N$, and polymer packing fraction, $\phi$, according to $d=b[\pi
N/(6\phi)]^{1/3}$, assuming that each Kuhn segment occupies an approximate
volume of $\pi b^{3}/6$. To evaluate whether chain length influences exchange
kinetics, the number of Kuhn segments per chain was swept over
$N=\\{12,18,36\\}$ for consistency with Section 3. Meanwhile, $\phi$ was swept
over the range $\phi\in[0.01,0.52]$ to capture a breadth of values likely
encompassing those of real gels and elastomers (Rubinstein and Colby, 2003,
Marzocca et al., 2013, Wagner et al., 2022). The upper limit of $\phi$ was set
based on empirical estimates of polymer packing in systems at ambient
conditions (see Appendix H for details), as well as our present models’
limiting ideal chain assumption, which invokes that chains do not interact via
volume exclusion or cohesive forces (e.g. Van der Waals). In reality, as
polymer free volume is decreased, a higher resultant concentration of inter-
chain interactions will constrain the conformations available to each chain,
likely decreasing their initial attachment and partner exchange rates due to
greater subdiffusion, while increasing their effective stiffness due to
reduction in entropy. Therefore, we impose that $\phi\leq 0.52$ based on the
understanding that the validity of the ideal chain assumption is improved at
lower packing fractions with fewer inter-chain interactions.

With $d$ computed for each combination of $N$ and $\phi$, the tether nodes
were positioned and then chain initiation was conducted per the methods of
Section A.1. Once initiated, the nodes comprising the tethered chains were
stepped through time according to Eq. (1) and their distal ends were allowed
to bond/unbond according to Eqs. (7-9). Stickers were only allowed attach to
one neighbor at a time, thus enforcing the simple conditions of mutually
exclusive pairwise binding. Attachment activation energy was set to
$\varepsilon_{a}=0.01k_{b}T$ (a low value), to induce rapid associative
kinetics (i.e., that $k_{a}^{ap}\approx\tau_{0}^{-1}$) and therefore
interrogate the agreement between models with high temporal resolution demands
through $k_{s}$. The activation energy for detachment was set to
$\varepsilon_{d}=7k_{b}T$ (a relatively high value so $k_{d}\approx
10^{-3}\tau_{0}^{-1}$) to ensure that bond lifetimes were comparable to the
detached lifetimes, but not so long that detachment (and by extension repeat
attachment/partner exchange) was rarely observed over the simulated duration
of $8\times 10^{3}\tau_{0}$. Data for this study was output with a frequency
of $10\tau_{0}^{-1}$, which we found ensures adequate detection of repeat
attachment events. Adequate sampling of average kinetic rates ($n\geq 100$)
was achieved under these conditions. Reasonable sampling of bound and unbound
sticker lifetimes was achieved at most chain concentrations. However, as few
as $n<10$ fully attached bond lifetimes were observed for very low chain
concentrations. This is because attachment/detachment kinetic rate
measurements require observation of only one state transition, whereas full
bond lifetime measurements require the observation of two (i.e., start and end
state transitions). Additionally, bond lifetimes are power law distributed so
that mean lifetimes are sensitive to outliers. For these reasons we emphasize
interpretation of rates in our discussion below, reserving mean bond lifetimes
for qualitative comparison across event types.

To characterize bond kinetics, average attachment rates, $\bar{k}_{a}$, were
computed according to Eq. (10) (Fig. 6A).555 First-time attachment events are
excluded from the pool of attachment events from which $\bar{k}_{a}$ was
computed, as their inclusion here significantly alters the relationship
between $\bar{k}_{a}$ and $d/(Nb)$. In contrast their exclusion led to
consistent $\bar{k}_{a}$ versus $d/(Nb)$ relations, regardless of whether
sampling was conducted over the entire simulation duration, or only the
simulation after approximately steady state fractions of attached/detached
chains were reached. The attachment rates were further partitioned into
partner exchange rates, $\bar{k}_{exc}$ (Fig. 6B), and repeat attachment
rates, $\bar{k}_{rpt}$ (Fig. 6C). These were also computed using Eq. (10),
except the number of attachment events, $N_{a}$, was simply replaced by the
number of exchange events, $N_{exc}$, or repeat events, $N_{rpt}$,
respectively. Values of average detached lifetime, $\bar{\tau}_{d}$, average
detached lifetime prior to partner exchange, $\bar{\tau}_{exc}$, and average
detached lifetime prior to repeat attachment, $\bar{\tau}_{rpt}$, were
directly measured from both models (insets of Fig. 6). The average detachment
rate, $\bar{k}_{d}$ (computed using Eq. (11)) and average fractions of
attached and detached chains at steady state ($f_{a}$ and $f_{d}$,
respectively) are provided for reference (Fig. I1), with steady state defined
as occurring once $df_{a}/dt=-df_{d}/dt<10^{-4}\tau_{0}^{-1}$. Average values
of attached bond lifetime, $\bar{\tau}_{a}$, and renormalized bond lifetime,
$\bar{\tau}_{rnm}$ are also provided for reference (Fig. I2). Importantly,
when plotted with respect to chain concentration, $c=d^{-3}$, the curves for
all outputs reported in Fig. 6 approximately collapse into single curves
irrespective of chain length (Fig. I3). This is consistent with the
predictions of Stukalin et al. (2013) and provides additional evidence that
bond kinetics within ensembles of chains are governed by sticker
concentration, independent of chain length. Distances, $d$, may be also be
interchanged with chain stretches, $\lambda_{c}=d/(\sqrt{2N}b)$, in Fig. 6 to
facilitate a more direct comparison to the results of Section 3. However, we
here visualize results in terms of $d/(Nb)$, which most clearly elucidates
reasonable agreement between the bead-spring and mesoscale models across
distinct chain lengths.

Figs. 6 and I1-I2 demonstrate that the bead-spring and mesoscale models’
predictions of kinetic rates, steady state bond fractions, and bond lifetimes
are in reasonable agreement across all tether separation distances, $d$, and
chain lengths, $N$. As in prior sections, the dissociation rates universally
match the value set a priori (Fig. I1.A). However, for the purposes of this
study we focus on the associative kinetics, which divulge additional
information about repeat attachment versus partner exchange. The magnitudes of
the associative kinetic rates are generally on the order of
$10^{-3}\tau_{0}^{-1}$, while corresponding detached bond lifetimes are on the
order of $10\tau_{0}$ to $10^{3}\tau_{0}$, confirming that adequate bond
kinetic sampling was established per Appendix D. As expected, attachment rates
decrease with increasing chain separation and chain length. However,
attachment rates no longer exhibit Gaussian scaling with respect to separation
distance as they did on a pairwise attachment basis through Eq. (19) in
Section 3.2. Instead, the attachment rates for ensembles of chains evolve as:

$\bar{k}_{i}=k_{i}^{semi}\left(\frac{d_{max}}{d}-1\right),$ (21)

where $d_{max}$ is the maximum empirically predicted tether-to-tether
separation distance at which attachment occurs, $k^{semi}_{i}$ is the
attachment rate when $d=d_{max}/2$, and $k_{i}$ denotes either overall
attachment rate, $k_{a}$, partner exchange rate, $k_{exc}$, or repeat
attachment rate, $k_{rpt}$. Fit curves are displayed in Fig. 6, and empirical
values of maximum chain stretch, $\lambda_{c}^{max}=d_{max}/(\sqrt{2N}b)$,
minimum chain concentration, $c_{min}b^{3}=(b/d_{max})^{3}$, and minimum
packing fraction, $\phi_{min}=\pi Nb^{3}/(6d_{max}^{3})$, at which associative
kinetics occur, are provided in Fig. I4.A-C. Fig. I4.D-F provides the
characteristic attachment rates $k_{a}^{semi}$, $k_{exc}^{semi}$, and
$k_{rpt}^{semi}$, for each event type.

Figure 6: Associative ensemble bond kinetics. (A) overall bond attachment,
$\bar{k}_{a}$, (B) exchange attachment, $\bar{k}_{exc}$, and (C) repeat
attachment, $\bar{k}_{rpt}$, rates with respect to normalized chain
separation, $\bar{d}/(Nb)$. Insets in (A-C) display average detached
lifetimes, $\bar{\tau}_{d}$, detached lifetimes prior to partner exchange,
$\bar{\tau}_{exc}$, and detached lifetimes prior to repeat attachment,
$\bar{\tau}_{rpt}$, respectively. Error bars represent standard error of the
mean. The empirical model of Eq. (21) is fit to all sets of bead-spring and
mesoscale data (solid and dashed curves, respectively).

Maximum chain stretch at which attachment occurs is generally on the order of
$\lambda_{max}\approx 1.5$, consistent with the findings of Sections 3.2.
Neither minimum chain concentration nor minimum packing fraction at which
attachment kinetics occur appear to vary significantly with chain length
(within the 95% confidence interval), nor do the empirically fit values of
$k_{a}^{semi}$, $k_{exc}^{semi}$, and $k_{rpt}^{semi}$ (Fig. I4). Again, this
supports the notion that sticker concentration dictates bond kinetics
(Stukalin et al., 2013). While any effects of chain length on these empirical
parameters are muted, agreement between models appears strongest for the
shortest chain length, which is visually expressed in the attached and
detached chain fractions (Fig. I1.B-C). Although reasonable agreement between
models is achieved, the steady state fraction of attached chains is
consistently higher for the bead-spring model than mesoscale model at lower
chain concentrations (i.e., higher values of $d$). This reveals a disparity in
the attachment rates at low chain concentrations that is not otherwise obvious
(due to their low magnitudes in Fig. 6A). Specifically, it suggests that the
attachment rates of the bead-spring model outstrip those of the mesoscale
model for long chains at low concentrations, which has meaningful network-
scale consequences discussed further in Section 4.

Aside from facilitating comparison of the two models, the results of Fig. 6
also highlight how Ångstrom and nanometer binding length scales – if not
calibrated on a bond-specific basis – could drive inaccuracies that propagate
upwards in length and time for either of these approaches. For instance, the
exchange rates at lower separation distances (or higher concentrations) in
Fig. 6B are greater than the corresponding repeat attachment rates in Fig. 6C.
This result is unexpected (Stukalin et al., 2013), but has been carefully
confirmed and is an artifact of our models’ assumptions. It arises from the
fact that the maximum sticker-to-sticker attachment distance is set equal to
the equilibrium length of dynamic bonds (both are $b$). Dynamic bond lengths
are normally distributed according to Fig. C1.A. Their length is nominally $b$
plus or minus the standard deviation of this distribution, so that their
dissociation does not necessarily leave their stickers within immediate
reattachment range ($\leq b$) of each other. Since the number of free stickers
available for partner exchange is typically greater than the number of
stickers available for repeat attachment (one), exchange kinetics occur more
frequently and this result is more pronounced at higher chain concentrations.
Thus, when applying these models to evaluate long-term network structures or
stress responses, an initial, high-resolution modeling approach (e.g., density
functional theory) or appropriate leveraging of prior work is strongly
recommended to justify the prescription of associative length scales.
Prescription of a variable, short-ranged binding probability based on the
competition between binding potential and kinetic energy may also be
considered.

##  4 Network-scale Mechanical Response

To evaluate whether the bead-spring and mesoscale models return reasonable
predictions of network-scale mechanical stress response, we simulated polymer
networks undergoing uniaxial extension and stress relaxation while monitoring
the stress response. All simulations described in this section were run on the
the 80-core Ice Lake (ICX) compute nodes of Stampede3, using the built-in MPI
capability of LAMMPS. First, we generated networks of $n_{p}=60$ branched
polymeric chains inside of periodic RVEs. Each polymer has five branched side
groups of length $Nb$ tethered to its backbone at even intervals of $Nb$ (Fig.
7A). Branches also terminate the chains as illustrated in Fig. 7A. For
detailed network initiation procedures see Appendix A.2. Once initiated, each
networks’ nodes were stepped through time according to Section 2.1 and the
distal stickers of their side chains were allowed to dynamically attached to
and detach from one another according to Section 2.2. In order to reasonably
reproduce the MSDs, bond kinetics, and initial network topologies of the bead-
spring model, we find that the damping coefficient must be set to the expected
value of $\gamma=Nz\gamma_{0}/2$, based on Rouse theory (Rouse, 1953,
Rubinstein and Colby, 2003). Here, $z$ is the number of chains attached to
each crosslink ($z=3$ for tethering crosslinks and $z=1$ for stickers) and the
multiple of $1/2$ enforces that each node of the mesoscale model encapsulates
half of the Kuhn segments of each chain it is attached to.

Once the networks were formed and equilibrated for some time, $t_{eq}$, the
RVEs were uniaxially extended to a maximum network stretch of $\lambda=3$, at
a constant true strain rate, $\dot{\varepsilon}$. Stretch was applied by
deforming the boundaries of the RVE according to the volume conserving
deformation gradient $\bm{F}(t)=$
diag$(e^{\dot{\varepsilon}t},e^{-\dot{\varepsilon}t/2},e^{-\dot{\varepsilon}t/2})$,
where $\dot{\varepsilon}$ denotes the true strain rate in the direction of
stretch. Once full extension was reached, the deformation gradient
$\bm{F}(t)=$ diag$(3,1/\sqrt{3},1/\sqrt{3})$ was sustained for
$t_{rlx}=220\tau_{0}$ so that stress relaxation responses could be evaluated.
The loading history in the direction of extension is plotted in Fig. 7B. To
attain adequate sampling, the results of $n=15$ simulations were ensemble
averaged for every parameter combination described below. These samples
consisted of five initial network structures equilibrated over three different
durations, $t_{eq}=\\{220,275,330\\}\tau_{0}$, for each parameter combination.

Figure 7: Network-scale parameter space. (A) Initiated bead-spring (top) and
mesoscale (bottom) model RVEs comprised of $n_{p}=60$ polymers where red nodes
indicate the positions of stickers. (B) Loading history applied to the RVEs
with respect to time. Snapshots of the bead-spring RVE (top) visualize the
loading in time. The orthonormal basis
$\\{\bm{e}_{1},\bm{e}_{2},\bm{e}_{3}\\}$ is denoted with the principle
direction of stretch as $\bm{e}_{1}$. (C-D) Parameter spaces indicating the
ranges of investigated (C) spatially and (D) temporally related inputs.
Snapshots of the network in (C) illustrate simulated RVEs of the bead-spring
model at the extreme chain lengths ($N=\\{12,36\\}$) and packing fractions
($\phi=\\{0.2,0.5\\}$). The grey and yellow regions in (D) indicate where
linear TNT with constant bond kinetic rates predicts steady state stress onset
and elastic-like response during loading, respectively.

To evaluate mechanical response, Cauchy stress was computed throughout
simulations as the potential energy-governed component of the virial stress:

$\bm{\sigma}=V^{-1}\sum_{\alpha}\sum_{\beta}\bm{r}^{\alpha\beta}\otimes\bm{f}^{\alpha\beta}-\pi\bm{I}.$
(22)

where $V$ is the total RVE volume, $\bm{r}^{\alpha\beta}$ denotes the end-to-
end vector spanning from node $\alpha$ to its attached neighbor $\beta$,
$\bm{f}^{\alpha\beta}=-\partial\psi/\partial\bm{r}^{\alpha\beta}$ is the force
between $\alpha$ and $\beta$, and the isotropic pressure term $\pi\bm{I}$
represents the volume conserving hydrostatic stress from excluded volume
interactions. Since a Lenard-Jones potential is included in the bead-spring
model but no such interactions are modeled between the implicit chains of the
mesoscale model, $\pi\bm{I}$ is neither expected nor intended to necessarily
match between approaches. Rather, the focus of this study is on each model’s
prediction of the change in stress due to reorientation, stretch, and
reconfiguration of its chains in response to loading. Thus, after
equilibration, initial stress is taken as $\bm{\sigma}=\bm{0}$ so that
$\pi\bm{I}$ is equal and opposite to the first term from Eq. (22). For both
models, $\alpha\beta$ pairs are treated as the polymer chains comprised of $N$
Kuhn segments that chemically attach tether-tether or tether-sticker pairs so
that virial stress is computed at the network or “mesh-scale”. Thus, the free
energy of chains is represented by Eq. (4), based on the models’ agreement
between force extension relations in Fig. C1.B. It is also possible to
directly measure the virial stress in the bead-spring model by carrying out
the sum in Eq. (22) over all Kuhn segments (Fig. J1). This approach delivers
good agreement with mesh-scale estimates of virial stress for the bead-spring
model (except at very high loading rates when
$\dot{\varepsilon}=0.1\tau_{0}^{-1}$ per Fig. J1.C). However, to reduce noise
and improve clarity for discussion, we present results within this section
using the mesh-scale virial stress.

These network models maintain the assumptions of ideal, monodisperse chains
without long-range interactions. For these simplified conditions, we
investigate the effects of four key parameters that may influence model
agreement. The first two parameters – related to spatial occupancy of the
chains – are chain length, via $N=\\{12,36\\}$, and polymer packing fraction,
$\phi=\\{0.2,0.5\\}$ (Fig. 7C). The limits of $N$ and upper limit of $\phi$
were selected for consistency with Section 3. The lower limit of $\phi=0.2$
was selected because percolated networks were not always observed when
$\phi<0.2$. The other two swept parameters – related to timescales of the
simulations – are the true strain rate,
$\dot{\varepsilon}=\\{0.01,0.1\\}\tau_{0}^{-1}$, and detachment rate,
$k_{d}=\\{0,0.1\\}\tau_{0}^{-1}$ (Fig. 7D). The range of $\dot{\varepsilon}$
was selected because it results in microsecond-scale simulations, but still
permits ensemble sampling using the bead-spring model at reasonable
computational cost. The range of $k_{d}$ was selected because setting
$k_{d}=0$ begets near-permanent network structures after equilibration
(allowing only for attachment of residually dangling stickers), while setting
$k_{d}=0.1\tau_{0}^{-1}$ induces near-complete stress relaxation over
relaxation times comparable to the loading times. The intrinsic attachment
rate was held constant at a high value, $k_{a}\approx\tau_{0}^{-1}$, to ensure
that solid-like networks structures were maintained.

We began by simulating both models at the extremes of these four parameter
ranges. For all cases investigated, the mesoscale model results in at least an
86% reduction in computational runtime (that scales with the number of polymer
chains as $\%\Delta t_{cpu}=1-0.81n_{p}^{-0.44}$), and an approximately 90%
reduction in data storage requirements (Fig. J2). These savings immediately
highlight the utility of the mesoscale model for exploring larger
spatiotemporal domains (discussed further in Section 5). To compare the
models’ predictions, Figs. 8 and J3 display the normal Cauchy stress,
$\sigma_{11}$, in the direction of loading with respect to time for every
combination of $\phi$, $k_{d}$, and $\dot{\varepsilon}$ when $N=12$ and
$N=36$, respectively. Results are presented in SI units of MPa with respect to
$\mu$s, demonstrating that reasonable orders of magnitude stress are predicted
by these models when compared to rubbery elastomers (Qi et al., 2003,
Vatankhah-Varnosfaderani et al., 2017) or small-mesh gels (Shibayama et al.,
2019).666Magnitudes of stress may be renormalized by simply adjusting the Kuhn
length per Appendix B. Error between models is expressed with respect to time
directly beneath each plot as the difference between mesoscale and bead-spring
model stress, normalized by the peak bead-spring model stress. This measure
(used throughout the remainder of this section) avoids overemphasizing error
at small values of stress (e.g., prior to loading or after relaxation), while
still allowing for meaningful error comparison between systems with different
magnitudes of absolute stress (e.g., networks at higher versus lower packing
fractions). Positive error indicates that the mesoscale model predicts higher
values of stress than the bead-spring model. For permanent networks loaded at
$\dot{\varepsilon}=0.01\tau_{0}^{-1}$ and comprised of chains with $N=12$ Kuhn
lengths between crosslinks (Fig. 8A), the models are in reasonable agreement
at both packing fractions. Peak error of $16\%$ between the models’ stress
predictions occurs at the end of the loading phase and the models converge to
similar values of elastic stress with less than $2\%$ error between them. This
suggests similar initial network structures after equilibration.

Even in the absence of bond detachment (Fig. 8A-B), stress relaxation is still
observed in both models due to frictional drag that retards conformational
chain relaxation. Due to drag during loading, both models predict elastic-like
responses of monotonically increasing stress with no reduction in Young’s
modulus (Fig. 8C), even in the regime in which linear TNT with constant bond
kinetics predicts onset of steady state stress (i.e., when the Weissenberg
number, $\dot{\varepsilon}/k_{d}<0.5$, Fig. 7D). This is attributed to the
fact that loading rates are relatively high
($\dot{\varepsilon}\not\ll\tau_{0}^{-1}$) and higher loading rates impart
greater particle velocities that increase drag effects. Indeed, increasing the
loading rate from $\dot{\varepsilon}=0.01\tau_{0}^{-1}$ to
$\dot{\varepsilon}=0.1\tau_{0}^{-1}$ while maintaining $k_{d}=0$ and $N=12$
(Fig. 8B), imparts higher peak mechanical stresses. Importantly, this increase
in loading rate also increases peak error between the approaches from $16\%$
(Fig. 8A) to $68\%$ (Fig. 8B). This implicates drag forces as a major
potential source of error between models and warrants further investigation
into the effects of not only loading rate, but also chain length since damping
coefficient depends on $N$.

Observing Fig. J3, the mesoscale model under-predicts the stress response of
the bead-spring model at all combinations of $\dot{\varepsilon}$, $k_{d}$, and
$\phi$ when $N=36$ (Fig. J3). Notably, very little change in mechanical
response is observed as $k_{d}$ is increased from $0$ to $0.1\tau_{0}^{-1}$,
suggesting that the stress response is not topologically governed, but is
rather almost entirely governed by drag for these longer-chained networks
(Fig. J3). Therefore, this disparity between models is likely attributable to
differences in the localization of drag forces within each approach that are
exacerbated as $N$ is increased. Since drag forces increase with $N$ in Eq.
(1) (i.e., $\gamma^{\alpha}d\bm{x}^{\alpha}/dt\propto N$) but entropic tensile
forces decrease with $N$ (i.e., $\bm{f}^{\alpha\beta}\propto N^{-1}$), any
differences between models’ drag effects may become more pronounced for long-
chained networks. Differences inevitably arise because the bead-spring model
distributes drag forces evenly along the lengths of its chains, whereas the
mesoscale model concentrates these forces at crosslink sites. To better
understand these effects while mapping the regimes in which the mesoscale
model is best suited, the remainder of this section provides deeper
investigation into the effects of chain length (Section 4.1), loading rate
(Section 4.2), and detachment rate (Section 4.3). The polymer packing fraction
is not investigated in further detail since little difference is observed in
the error between models when $\phi$ is increased from $0.2$ to $0.5$ (Fig.
8).

Figure 8: Mechanical response at the parameter extremes when $\bm{N=12}$.
Normal Cauchy stress in the direction of loading for the bead-spring and
mesoscale models and a modified measure of relative error (mesoscale model
stress minus bead-spring model stress as a percentage of peak stress bead-
spring model stress), are plotted with respect to time at both polymer packing
fractions ($\phi=\\{0.2,0.5\\}$), loading rates
($\dot{\varepsilon}=\\{0.01,0.1\\}\tau_{0}^{-1}$), and bond detachment rates
($k_{d}=\\{0,0.1\\}\tau_{0}^{-1}$) when $N=12$. Horizontal dashed lines on
error plots denote zero error.

###  4.1 Short chains impart topological and mechanical agreement

To investigate the effects of chain length, we initiated and deformed
permanent networks ($k_{d}=0$) comprised of chains with
$N=\\{12,18,24,30,36\\}$ Kuhn segments according to the procedure and loading
described above when $\dot{\varepsilon}=0.01\tau_{0}^{-1}$ and $\phi=0.2$.
Cauchy stress and relative error are plotted with respect to time in Fig.
9A-B. Fig. 9B reveals that error is generally lowest (never exceeding 8% and
dropping to less than 5% after stress relaxation) when $N=12$, and magnitudes
of relative error increase with $N$.

Figure 9: Effects of chain length. (A) Cauchy stress in the direction of
loading and (B) error between the mesoscale and bead-spring models (relative
to the peak error of the bead-spring model) with respect to time when
$\phi=0.2$, $k_{d}=0$, $\dot{\varepsilon}=0.01\tau_{0}^{-1}$ and
$N=\\{12,18,36\\}$. (C) Cauchy stress (normalized by peak stress) with respect
to time. Solid and dashed curves represent best fits of Eq. (23) to the bead-
spring and mesoscale data, respectively. (D) Stored stress, $\sigma_{e}$
(solid purple), and dissipated stress, $\sigma_{d}$ (dashed purple), along
with $\alpha$-relaxation rate, $k_{\alpha}$, (solid green) for the mesoscale
model (triangles) and bead-spring model (circle) with respect to $N$. Error
bars represent the 95% confidence interval from nonlinear least-squares
regression analysis.

Since no bond dissociation is allowed here ($k_{d}=0$), error in stress
responses between models must originate from differences in initial topology,
frictional damping, or both. Observing Fig. 9A, the mesoscale model under-
predicts bead-spring model stress when $N=36$ even near the end of stress
relaxation, strongly suggesting differences in network structure. We find that
mesoscale networks of long chains ($N=36$) – although containing the same
relative fractions of attached/detached bonds as bead-spring networks (Fig.
J4A) – express far higher average degrees of non-load transmitting bonds
between branches of the same molecules (Fig. J4.B). This is attributed to much
slower diffusion of the mesoscale (versus bead-spring) nodes during network
initiation and equilibration (see Fig. J5 for MSD data). It appears that the
slower diffusion of the stickers and tethers in the mesoscale model
predisposes the networks comprised of the highest chain lengths to clustering
and intra-molecular bonding, and that this effect is sharply more pronounced
when $N=36$ (Fig. J4.B). Fewer percolated load paths in the mesoscale model
results in reduced strain energy storage (You et al., 2024) and under-
predicted mechanical stress as compared to the bead-spring model. This
emphasizes the importance of replicating associative bond kinetics between
models and affirms network structure as a secondary effect of frictional
damping.

To better understand the effects of chain length on frictional damping, we
isolate the stress relaxation response. Fig. 9C displays the stress relaxation
response normalized by peak stress, $\sigma_{11}^{max}$. We may fit a set of
two parallel Maxwell elements to this response following the form:

$\sigma_{11}=\sigma_{11}^{max}\left[f_{\alpha}\exp(-k_{\alpha}t)+(1-f_{\alpha})\exp(-k_{r}t)\right]$
(23)

where $f_{\alpha}$ represents the relative fraction of stress dissipated due
to chain relaxation or “$\alpha$-relaxation”, $k_{\alpha}$ represents the
$\alpha$-relaxation rate mediated by drag forces (not to be confused with
attachment rate, $k_{a}$), and $k_{r}$ represents the rate of
reconfigurational relaxation driven by $k_{d}$. Since $k_{d}=0$, $k_{r}$ must
also be zero so that Eq. (23) simplifies to a Zener element model where
$\sigma_{d}=\sigma_{11}^{max}f_{\alpha}$ becomes the dissipated stress,
$\sigma_{e}=\sigma_{11}^{max}(1-f_{\alpha})$ becomes the stored or
“equilibrium” stress, and $k_{\alpha}$ is definitively mediated by drag forces
alone. Thus, deviations in $k_{\alpha}$ between models expose differences in
the characteristic chain relaxation rates due to drag. Meanwhile, differences
in $\sigma_{e}$ broadcast differences in the initial, permanent network
structures, whereby fewer load carrying, inter-molecular attachments results
in smaller values of $\sigma_{e}$ (You et al., 2024). Fig. 9D displays fit
values of $\sigma_{d}$, $\sigma_{e}$, and $k_{\alpha}$ with respect to $N$.
The $\alpha$-relaxation rates, $k_{\alpha}$, are consistently higher for the
mesoscale model than bead-spring model. This faster conformational relaxation
indicates that the mesoscale model, which lumps $Nz/2$ Kuhn segments into
every node, under-represents frictional resistance of intermediate bead-spring
chains, despite its slower diffusive exploration of crosslinks and stickers
due to localization of drag forces. However, these effects appear to minimally
effect stress response when $N=12$, for which diffusion kinetics (Fig. J5) and
thus network topologies (see $\sigma_{e}$ in Fig. 9D) are in reasonable
agreement. Therefore, recommended practice is to limit the mesoscale model’s
implicit chains to contour lengths around $12b$.

###  4.2 Slow loading bolsters mechanical agreement

To investigate the effects of loading rate, we initiated and deformed networks
comprised of chains with $N=12$ and $\phi=0.2$, at loading rates of
$\dot{\varepsilon}=\\{0.01,0.02,0.03,0.06,0.1\\}\tau_{0}^{-1}$. To isolate the
effects of loading rate on frictional $\alpha$-relaxation we again disable
dissociation by setting $k_{d}=0$. Cauchy stress and relative error between
models are plotted with respect to time in Fig. 10A-B. As expected based on
prior sections, peak absolute error always occurs at the end of loading and
increases monotonically with respect to loading rate. The predictions of both
models always converge to values of equilibrated stress within 5% error of one
another under these conditions, further substantiating agreement in initial
network topologies between approaches.

To explore clear loading rate-dependence of the stress relaxation behavior, we
isolate and fit the normalized stress relaxation response (Fig. 10C-D) with a
stretched Kohlrausch-Williams-Watts exponential decay function (Richardson et
al., 2019) of the form:

$\sigma_{11}=\sigma_{11}^{max}\left\\{f_{\alpha}\exp\left[-(k_{\alpha}t)^{\chi})\right]+(1-f_{\alpha})\exp(-k_{r}t)\right\\},$
(24)

where $k_{r}$ remains zero, $k_{\alpha}$ is the nominal $\alpha$-relaxation
rate, and the parameter $\chi\in[0,1]$ quantifies the degree of variable
relaxation rate in time. If $\chi=1$ there is no variability, while if
$\chi<1$ the relaxation rate decreases with time. The stretched exponential
fit is chosen because networks loaded at faster rates experience greater chain
stretches that correspond to disproportionately high entropic tension of their
nonlinear chains (Fig. C1.B). These disproportionately high forces drive the
networks’ crosslinks into a relaxed state more quickly, resulting in a faster
initial relaxation response that slows with time.

Figure 10: Effects of loading rate. (A) Normal Cauchy stress in the direction
of loading and (B) percent relative error between the bead-spring and
mesoscale models plotted with respect to time when $\phi=0.2$, $N=12$,
$k_{d}=0$, and $\dot{\varepsilon}=\\{0.01,0.03,0.1\\}\tau_{0}^{-1}$. (C)
Cauchy stress, normalized by peak stress, during stress relaxation with
respect to time. (D) Fit values of $f_{\alpha}$ (purple), $\chi$ (orange), and
$k_{\alpha}$ (green) with respect to $\dot{\varepsilon}$ for the bead-spring
(circles) and mesoscale (triangles) models.

For both models, $\chi$ decreases with increasing loading rates, as expected,
indicating more variable $\alpha$-relaxation rates when $\dot{\varepsilon}$ is
higher. The relative fractions of dissipated stress, $f_{\alpha}$, and nominal
relaxation rates, $k_{\alpha}$, both increase with respect to loading rate,
also as expected. Comparing the modeling approaches, we see that the bead-
spring model consistently exhibits lower relative fractions of dissipated
stress, $f_{\alpha}$, consistent with the observation that peak stresses of
the mesoscale model are greater than those of bead-spring model despite
subsequently reaching close equilibrium stresses. While the nominal relaxation
rates, $k_{\alpha}$, of the models are in good agreement at the slowest
loading rate, $k_{\alpha}$ for the mesoscale model increase more dramatically
with respect to strain rate. This rate dependence is consistent with the
mesoscale model’s under-prediction of frictional resistance to conformational
change. These results suggests that the mesoscale model should be reserved for
loading cases in which $\dot{\varepsilon}\leq 0.01\tau_{0}^{-1}$, which is
both realistic and necessary for most polymers given that molecular
oscillations are typically on the order of $10^{12}$ to $10^{14}$ Hz
(Herzberg, 1955). Since slower loading rates require longer loading times,
they are also more readily achieved by the mesoscale model than bead-spring
approach due to its lower computational cost (see Section 5).

###  4.3 Model agreement is achieved across detachment rates

To investigate the effects of detachment rate, we initiated and deformed
networks comprised of chains with $N=12$ and $\phi=0.2$, while sweeping
$k_{d}=\\{0,0.001,0.003,0.01,0.03,0.1\\}\tau_{0}^{-1}$. The applied loading
rate was held at $\dot{\varepsilon}=0.01\tau_{0}^{-1}$ based on the results of
Section 4.2. Cauchy stress and relative error between models are plotted with
respect to time in Fig. 11A-B. Generally, the models are in reasonable
agreement across detachment rates (Fig. 11A), although the frictional effect
that causes slight over-prediction of mesoscale model stress remains present.
While relative error appears to increase with $k_{d}$ (because magnitudes of
stress also decrease with $k_{d}$), maximum absolute values of error between
models remain around 0.2 MPa regardless of $k_{d}$ (Fig. 11B). These results
suggest that $k_{d}$ does not drive significant deviation between the models.

To quantify the effect of $k_{d}$ on relaxation rates, we once more fit Eq.
(24) to the normalized stress relaxation response (Fig. 11C). However, $k_{r}$
is now treated as a fitting parameter that is expected to scale directly with
$k_{d}$. Fit parameters are plotted with respect to $k_{d}$ in Fig. 11D-E.
Examining Fig. 11D, $0.8<\chi\leq 1$ for all $k_{d}$ of both models,
indicating that chain stretch causes little variability in the relaxation rate
for the prescribed loading rate777Since $k_{d}$ is deliberately made force-
independent here, variability in relaxation rate cannot arise from force-
dependent detachment and can instead only be attributed to stretch-dependent
conformational relaxation. Additionally, $\chi$ is uncorrelated with $k_{d}$,
consistent with findings from Section 3.2 that chains attach to one another
within the linear regime of force-stretch (i.e., $\lambda<1.5$). While the
reconfigurational relaxation rate, $k_{r}$, scales directly with $k_{d}$ (Fig.
11E, $p<0.01$ for both models), the $\alpha$-relaxation rate, $k_{\alpha}$, is
uncorrelated with $k_{d}$ ($p>0.1$ for both models) since it is mediated only
by the damping coefficient. Nonetheless, one might expect the fraction,
$f_{\alpha}$, of stress dissipated due to conformational relaxation to
decrease as $k_{r}$ increases. However, $f_{\alpha}$ is also independent of
$k_{d}$ (Fig. 11D, $p>0.2$ for both models). This underscores the coupling
between configurational and conformational network relaxation that both models
capture, whereby every detachment event results in new conformational degrees
of freedom for chain relaxation (Stukalin et al., 2013, Yu et al., 2014,
Wanasinghe et al., 2022, Wagner and Vernerey, 2023). Despite the complex
coupling between $k_{r}$ and $k_{\alpha}$, Fig. 11 broadly suggests that
modulating $k_{d}$ does little to alter agreement in stress response between
models. In the following section, we leverage this finding, as well as the
reduced computational cost of the mesoscale model (Fig. J2), to simulate large
spatiotemporal domains and explore the regime in which $k_{d}\ll\tau_{0}^{-1}$
as is generally the case in real materials at operational temperatures
(Richardson et al., 2019, Chen et al., 2019).

Figure 11: Effects of detachment rate. (A) Normal Cauchy stress in the
direction of loading and (B) percent relative error between the bead-spring
and mesoscale models are plotted with respect to time when $\phi=0.2$, $N=12$,
$\dot{\varepsilon}=(0.01)\tau_{0}^{-1}$ and
$k_{d}=\\{0,0.01,0.1\\}\tau_{0}^{-1}$. (C) Cauchy stress normalized by peak
stress during stress relaxation with respect to time. (D-E) Fit values of
$f_{\alpha}$ (purple) and $\chi$ (orange), as well as (E) and $k_{\alpha}$
(green) and $k_{r}$ (blue), all with respect to $k_{d}$ for the bead-spring
(circles) and mesoscale (triangles) models.

##  5 Mesoscale access to larger spatiotemporal domains

In Section 4, we interrogated the effects of chain length, loading rate, and
detachment rate at spatiotemporal scales readily accessible to the bead-spring
model ($\sim 10$ nm and $\sim 1$ $\mu$s). Importantly, this restricts loading
rates to fast regimes in which $\dot{\varepsilon}$ is on the order of 1% to
10% of $\tau_{0}^{-1}$ and frictional effects are dominant. By extension, this
restricts worthwhile examination of the effects of detachment rates to regimes
in which $k_{d}$ is roughly 10% of $\tau_{0}^{-1}$, since crossover from
predominantly elastic to viscous behavior tends to occur at loading rates on
the order of $\dot{\varepsilon}\sim k_{d}$. Here we demonstrate how the
mesoscale model’s computational cost savings may be put towards modeling
relatively large networks with dimensions on the order of 100 nm, over time
domains on the order of 100 $\mu$s (i.e., $10^{5}\tau_{0}$), thus lessening
the gap between molecular theory and experimentally relevant scales.

One feature observed in many dynamic polymers is the ability to undergo
maximum global stretch ratios, $\lambda_{max}$, well over $10\times$ without
failure (Zhang et al., 2019, Cai et al., 2022, Xu et al., 2022) at relatively
slow loading rates relative to $k_{d}$ (i.e., $\dot{\varepsilon}/k_{d}<0.5$).
However, to capture this effect in a Lagrangian RVE requires simulating slow
loading rates ($\dot{\varepsilon}<k_{d}\ll\tau_{0}^{-1}$), at loading
timescales up to $t_{load}=\ln(\lambda_{max})/\dot{\varepsilon}$. It also
requires simulating RVEs with large enough initial dimensions, $L_{0}$, so
that incompressible uniaxial extension does not result in domain widths less
than the maximum allowable bond length (i.e., $L_{0}\geq
Nb\sqrt{\lambda_{max}}$). For instance, to model a network undergoing uniaxial
stretch up to $\lambda_{max}=20$, the initial RVE should have dimensions of no
less than $4.6Nb$, which for the networks investigated here ($N=12$),
corresponds to $L_{0}\geq 34$ nm. Furthermore, supposing
$\dot{\varepsilon}/k_{d}=0.5$, and $k_{d}\sim 10^{-3}\tau_{0}^{-1}$, then
deformation would need to be applied for a duration of
$t_{load}=\ln(20)/[0.5\times 10^{-3}\tau_{0}^{-1}]\approx 6\times
10^{5}\tau_{0}$. For the bead-spring model, whose time step is $\delta
t=4\times 10^{-4}\tau_{0}$, this would require on the order of $10^{9}$
discrete iterations. While such simulations become impractical for the bead-
spring approach, especially as $\lambda_{max}$ is increased or
$\dot{\varepsilon}$ is decreased, the mesoscale model’s reduced timestep and
explicitly tracked number of particles renders modeling such deformations
easily achievable.

To demonstrate the capabilities of the mesoscale model, we simulated networks
with initial RVE dimensions of $L_{0}\approx 70$ nm, undergoing
incompressible, uniaxial extension to $\lambda_{max}=20$ per Fig. 12. As in
Section 4, all simulations in this section were run using the ICX compute
nodes of Stampede3 and LAMMPS MPI capability. Networks were initiated and
equilibrated for $t_{eq}=220\tau_{0}$ per the procedure of Appendix A.2. Based
on the results of Section 4 the networks’ polymer chains are comprised of
$N=12$ Kuhn lengths per chain segment, and the packing fraction was set to
$\phi=0.2$ so that $n_{p}=5\times 10^{3}$ polymer chains are needed to attain
the desired RVE size. In order to diminish the effects of frictional damping
and meaningfully interrogate the effects of loading rate, we set
$k_{d}=10^{-3}\tau_{0}^{-1}$ and swept the Weissenberg number over
$\dot{\varepsilon}/k_{d}=\\{0.125,0.25,0.5,1\\}$. Simulations were run until
$\lambda_{max}=20$ was reached, or chains within the network approached
stretches of $\sqrt{N}$ causing divergence in force and numerical
instability.888Numerical instability may be circumvented by enforcing force-
dependent bond dissociation (Puthur and Sebastian, 2002, Shen and Vernerey,
2020, Lamont et al., 2021, Buche and Silberstein, 2021, Wagner et al., 2021,
Song et al., 2021, Buche et al., 2022, Mulderrig et al., 2023). Here we focus
on the rate-dependence originating from constant bond detachment rates and
stopped simulations when chain lengths reached $r\approx Nb$ (at network
stretches of $\lambda=10.5$ and $\lambda=4.5$ when
$\dot{\varepsilon}/k_{d}=0.5$ and $\dot{\varepsilon}/k_{d}=1$, respectively)

Fig. 12B displays the mechanical stress response of the simulated networks
with respect to time for all four Weisenberg numbers. As expected, the
mechanical response is highly rate-dependent. However, unlike the rate-
dependence of Section 4.2, here we observe the onset of steady state stress
for values of $\dot{\varepsilon}/k_{d}<0.5$, in accordance with TNT (Vernerey
et al., 2017). This is representative of many dynamic polymers (Jeon et al.,
2016, Cai et al., 2022) for which the rate of energy dissipation equilibrates
with the rate of work input at sufficiently slow loading speeds. In contrast,
when $\dot{\varepsilon}/k_{d}\geq 0.5$, an elastic-like response of
monotonically increasing stress is observed. The true stress-strain responses
for networks loaded at rates $\dot{\varepsilon}/k_{d}\geq 0.5$ assume the
nonlinear shapes characteristic of rubbery, hyperelastic materials (Treloar,
1943), and this effect is more pronounced when $\dot{\varepsilon}/k_{d}=1$
than $\dot{\varepsilon}/k_{d}=0.5$. Fig. 12C display snapshots of sample
networks at a stretch of $\lambda=4$ when $\dot{\varepsilon}/k_{d}=0.125$ and
$\dot{\varepsilon}/k_{d}=1$ to visibly convey differences in chain stretch for
a network undergoing steady state creep (top) versus monotonically increasing
stress (bottom). This study successfully reproduces a typical rate-dependent
dynamic polymer response by probing the $10^{-7}$ m scale for over 90 $\mu$s
in the longest cases, and does so in the span of just 7 wall-clock hours.

Figure 12: Effects of loading rate during monotonic large deformation at
constant true strain rate. (A) Initiated mesoscale model RVE comprised of
$n_{p}=5\times 10^{3}$ polymers, where red nodes indicate the positions of
stickers. (B) Normal Cauchy stress in the direction of loading with respect to
true strain, $\varepsilon=\ln(\lambda)$, when $\phi=0.2$, $N=12$,
$k_{d}=10^{-3}\tau_{0}^{-1}$, and
$\dot{\varepsilon}/k_{d}=\\{0.125,0.25,0.5,1\\}$. Data represents the average
of $n=3$ samples and shaded regions represent S.E. of the mean. (C) Snapshots
of two RVEs normal to the $\\{\bm{e}_{i}\\}_{i=1,2}$ plane at a stretch of
approximately $\lambda=4$ when $\dot{\varepsilon}/k_{d}=0.125$ (top) and
$\dot{\varepsilon}/k_{d}=1$ (bottom). Bond colors represents chain stretch in
the direction of applied extension ($\bm{e}_{1}$).

Dynamic mechanical analysis (DMA) frequency sweeps are another common approach
for probing strain rate dependence of polymers. We therefore simulated
networks of $n_{p}=100$ polymers with $N=12$ and $\phi=0.2$ undergoing
oscillatory shear. For elastomers and gels, pure shear loading conditions are
typically applied via parallel plate rheology. To replicate pure shear within
the orthonormal RVEs of the mesoscale model, we applied plane strain in
directions $\bm{e}_{1}$ and $\bm{e}_{2}$ following
$\epsilon_{11}=\epsilon_{0}\sin(2\pi\omega t)$ and
$\epsilon_{22}=-\epsilon_{11}$, respectively, where $\epsilon_{0}$ is the
strain amplitude and $\omega$ is the angular frequency. For gels and other
soft materials, $\epsilon_{0}$ is typically on the order of $0.01$ to maintain
linearity; however, we set $\epsilon_{0}=0.05$ to enhance the signal-to-noise
ratio. These loading conditions reasonably replicate pure shear for small
strains so that shear strain and stress may be approximated using in-plane
transformations as
$\epsilon_{12}\approx-(\frac{\epsilon_{11}-\epsilon_{22}}{2})\sin(2\theta)$
and $\sigma_{12}\approx-(\frac{\sigma_{11}-\sigma_{22}}{2})\sin(2\theta)$,
where $\theta\approx\pi/4$ is the orientation of maximum shear with respect to
the the principle basis, $\\{\bm{e}_{1},\bm{e}_{2}\\}$. Snapshots of an RVE
undergoing oscillatory pure shear are shown in Fig. 13A, while a sample of the
corresponding applied strain and steady state stress response are depicted in
Fig. 13B-C. Storage modulus
($G^{\prime}=\sigma_{peak}/\epsilon_{0}\cos(\delta)$) and loss modulus
($G^{\prime\prime}=\sigma_{peak}/\epsilon_{0}\sin(\delta)$) are calculated
from such stress and strain data, where $\sigma_{peak}$ is the peak value of
the measured shear stress and $\delta$ is the phase shift between the applied
strain and resulting stress.

Figure 13: Results of numerical DMA. (A) Snapshots of an RVE of $n_{p}=100$
polymers normal to the $\\{\bm{e}_{i}\\}_{i=1,2}$ plane (1) in the undeformed
state, (2) at $\epsilon_{11}=0.05$, and (3) at $\epsilon_{11}=-0.05$. (B)
Sample loading history ($\epsilon_{12}$ versus time) over $2$ $\mu$s when
$f=(10^{-2})\tau_{0}^{-1}$. (C) Cauchy stress, $\sigma_{12}$, resulting from
the load history of (B). (D) Storage and loss moduli with respect to
normalized frequency, $f\tau_{0}$ for $k_{d}=(10^{-4})\tau_{0}^{-1}$ (top)
$k_{d}=(10^{-3})\tau_{0}^{-1}$ (center), and $k_{d}=(10^{-2})\tau_{0}^{-1}$
(bottom). Loss tangent, $\tan(\delta)$ is provided in faint grey. Vertical
dotted-dashed lines denote $k_{rlx}$ (i.e., the frequency at which
$\tan\delta=1$). Data is ensemble averaged from $n=15$ samples at each
frequency and detachment rate. Shaded region in (C) and error bars in (D)
represent SE.

Storage and loss moduli, as well as corresponding loss tangents
($\tan\delta$), for the ensemble average of $n=15$ networks are plotted with
respect to frequency ($f=2\pi\omega$) in Fig. 13D, for detachment rates of
$k_{d}=\\{10^{-4},10^{-3},10^{-2}\\}\tau_{0}^{-1}$. We choose to sweep
detachment rate because $k_{d}$ is tunable in many material systems via
crosslinker chemistry (e.g., cationic species and valency in metallopolymers)
(Kloxin et al., 2010, Richardson et al., 2019, Vidavsky et al., 2020, Zhang et
al., 2020). All $n=15$ networks were independently equilibrated and simulated
at $17$ distinct frequencies. Cumulatively, each sample network was simulated
for a duration in excess of $400$ $\mu$s over roughly 11 wall-clock hours,
exhibiting this method’s practicable access to large timescales. Evaluating
the results, storage and loss moduli are consistently on the order of $10$ to
$10^{2}$ kPa, which is typical of a gel (Okamoto et al., 2011) or soft
elastomer (Chen et al., 2015, Shabbir et al., 2016, Zhang et al., 2020, Peter
et al., 2021). At all frequencies the loss tangent is relatively high
($\tan(\delta)>0.5$), signifying that the networks generally undergo a large
degree of non-affine chain relaxation regardless of $k_{d}$.

Trends in $G^{\prime}$ and $G^{\prime\prime}$ are exemplary of a linear
viscoelastic material in that high frequencies beget predominantly elastic
response ($G^{\prime\prime}/G^{\prime}=\tan(\delta)<1$), while low frequencies
impart viscous behavior ($\tan(\delta)>1$). The loading frequency at which
this transition occurs (Fig. 13.D, dotted-dashed lines) is interpretable as
the relaxation rate, $k_{rlx}$ (Rubinstein and Colby, 2003). As the detachment
rate increases (Fig. 13.D, dashed lines), so too does $k_{rlx}$, and no
transition is observed for the slowest detachment rate
($k_{d}=10^{-4}\tau_{0}^{-1}$) for which $\tan(\delta)$ is always less than
unity (indicative of rubbery response). The dependence of $k_{rlx}$ on $k_{d}$
confirm that this model is capturing the transition regime wherein relaxation
is mediated by bond detachment, rather than $\alpha$-relaxation (Chen et al.,
2015, Shabbir et al., 2016, Zhang et al., 2020) as in Section 4.2. For the two
detachment rates with observable transitions
($k_{d}=\\{10^{-3},10^{-2}\\}\tau_{0}^{-1}$), $k_{d}/k_{rlx}\approx 6$
affirming that the relaxation rate is slower than $k_{d}$. This is attributed
to bond lifetime renormalization (Stukalin et al., 2013) and is representative
of many experimental dynamic polymers (Rubinstein and Colby, 2003, Chen et
al., 2015, Shabbir et al., 2016, Zhang et al., 2020) whose rheological
responses are best captured by the sticky Rouse model wherein chain relaxation
is retarded by the presence of dynamic binding sites (Leibler et al.,
1991).999To measure partner exchange rates directly requires high frequency
data output ($k_{s}\sim 0.1\tau_{0}^{-1}$) so that data curation becomes
untenable for detailed bond exchange studies at these timescales. The
mesoscale model’s capture of reconfigurational relaxation at small strains and
slow loading rates, as well as its ability to model stretch-dependent
relaxation spectra (Yu et al., 2014) at large strains and fast loading rates
in Section 4, spotlight its capacity for predicting physically realistic
trends in polymers across many decades of time.

##  6 Conclusion

We have formulated an idealized mesoscale model for dynamic polymers at the
single-chain and network scales and conducted detailed investigation on its
validity. Unlike comparable prior studies, which implicitly modeled
dynamically bonding stickers (Wagner et al., 2021, Wagner and Vernerey, 2023),
the mesoscale model introduced here explicitly tracks sticker positions.
Furthermore, whereas prior studies primarily compared mesoscale models to
macroscale theory or experiments, we have taken a bottom-up approach in which
the mesoscale model was compared to a Kremer-Grest coarse-grained MD approach
(Kremer and Grest, 1990), as well as to statistical mechanics-based molecular
theory where possible (Rubinstein and Colby, 2003, Stukalin et al., 2013). Our
method revealed that the effective bond activation energy for telechelic
association of polymer chains is penalized by the Helmholtz free energy of the
end-to-end stretch they must assume for attachment. This relation is predicted
by statistical mechanics (Buche and Silberstein, 2021), has been discovered
experimentally (Guo et al., 2009), and mirrors Bell’s well-established model
for force-dependent bond dissociation (Bell, 1978). While this finding may
offer a route by which to implicitly model stickers under specific conditions
(e.g., low sticker concentration), our bottom-up approach highlights the
importance of explicitly tracking stickers.

We found that seemingly small differences in bead-spring versus mesoscale
sticker MSDs below the Rouse time culminate in distinct partner exchange and
repeat attachment rates (Stukalin et al., 2013), which in turn drive
measurable structural differences for networks of long chains and low chain
concentrations between models. These differences were reduced, but not
eliminated, by setting the damping coefficient for a mesoscale node to
$\gamma^{\alpha}=Nz\gamma_{0}/2$ following Rouse subdiffusion (Rubinstein and
Colby, 2003, Stukalin et al., 2013), or $\gamma^{\alpha}=N^{2/3}\gamma_{0}$
for chains tethered to a fixed point. Based on these findings, we recommend
maintaining implicit chains with on the order of $N=12$ segments when
employing these mesoscale methods, as seen in many polymers with high dynamic
binding site concentrations (Colby et al., 1998, Chen et al., 2013, Zhang et
al., 2020, Xu et al., 2022, Xie et al., 2024). Future work could explore
incorporating intermediate nodes along the lengths of mesoscale chains to
simulate high molecular weight systems. We found that these models predicted
mechanical network responses in excellent agreement with each other for
networks comprised of short-chained ($N=12$), monodisperse polymers loaded at
strain rates below 1% of the monomer diffusion rate. Furthermore, they nicely
predicted the chain length, loading rate, and detachment rate-dependent trends
expected of dynamic elastomers and gels (Jeon et al., 2016, Zhang et al.,
2019, Xu et al., 2022, Cai et al., 2022). The mesoscale model achieved these
predictions with a 90% reduction in computational time and data storage
requirements compared to the bead-spring model. Via these savings, the
mesoscale approach offers access to spatiotemporal scales not readily
available using conventional MD approaches including large Lagrangian
deformations and long, variable-timescale numerical experiments (e.g.,
frequency sweeps).

When applying this mesoscale method over longer durations, we emphasize the
importance of accurately representing a material’s structural evolution. Any
discrepancies in the rates of bond restructuring and conformational changes
are liable to cause error stack-up over timescales significantly longer than
the renormalized bond lifetime. A major focus of future work for this
mesoscale approach that may improve long-term structural accuracy is capturing
crucial inter-chain and polymer-solvent interactions. Real materials such as
gels and elastomers host innate homogenization mechanisms such as osmotic
pressure (Flory, 1942, 1985) or excluded volume interactions between chains
(Zimm et al., 1953). Both of these mechanisms culminate in a material-scale
pressure ($\pi\bm{I}$ in Eq. 22) that resists clustering of chains (Mordvinkin
et al., 2021) and mediates bulk volumetric deformation in real materials and
conventional MD approaches. In the case of gels, a recently introduced coarse-
grained method for capturing osmotic pressure as a function of the
$\chi$-parameter may be incorporated into this framework (Flory, 1942, Doi,
2013, Wagner et al., 2022). However, more work is required to develop
mechanisms that encapsulate the effects of steric interactions between chains,
especially as seen in denser polymers such as melts and elastomers, or highly
entangled polymers with high molecular weights. Indeed, these same inter-chain
interactions are responsible for the dissipative entanglements that partially
cohere and greatly toughen many high molecular weight polymers (Sun and
Faller, 2006, Ge et al., 2013, Schieber and Andreev, 2014, Masubuchi, 2014, Ge
et al., 2018, Kim et al., 2021, Steck et al., 2023, Shi et al., 2023). This
mesoscale model may serve as a foundational framework into which researchers
may incorporate novel reduced-order methods for such mechanisms and provide a
powerful tool for the predictive design of dynamic polymers.

##  Appendix A Model initiation procedures.

Here we describe the numerical initiation procedures for all simulations of
the various studies in the main text. Initiation procedures for the chains in
tethered diffusion, bond kinetics, and bond exchange studies of Section 3 are
described in Section A.1. The procedures used to generate RVEs for network-
scale mechanics studies in Sections 4-5 are described in Section A.2. All
initiation procedures were carried out using custom codes written in
MATLAB2022a.

###  A.1 Single-chain initiation procedures

To validate single-chain characteristics (e.g., end-to-end distributions) and
bond kinetics (e.g., attachment, detachment, and exchange rates), various
tethered chain studies were conducted in Section 3, wherein an end of each
chain was fixed while the rest of it diffused following Eq. (1). To initiate
the chains in these studies, arrays of fixed, tethering nodes were first
generated at the position set $\\{\bm{x}^{1}_{m}\\}$ where the index $m$
denotes the molecule number. The detailed configuration (e.g., number of
nodes, their relative position to each other, etc.) of these initial tethering
nodes are described in detail on a case-by-case basis in Sections 3.1-3.3.
However, in every case, simulated polymer chains were randomly generated from
these seeded tethering sites.

For the bead-spring iterations of these studies, this was achieved using a 3D
random walk approach whereby new beads were subsequently positioned at
$\\{\bm{x}_{m}^{\alpha}\\}$ based on the positions of the previous beads of
their chain, $\\{\bm{x}_{m}^{\alpha-1}\\}$, according to:

$\bm{x}_{m}^{\alpha+1}=\bm{x}_{m}^{\alpha}+b\bm{\hat{r}}^{\alpha},$ (A1)

until the desired number of Kuhn segments, denoted by index $\alpha$, was
obtained (i.e., $\alpha\in[1,N]$). Here, the step size is that of a Kuhn
length, $b$, and the directional vector,
$\bm{\hat{r}}^{\alpha}=\bm{r}^{\alpha}/r^{\alpha}$, was determined according
to:

$\bm{r}^{\alpha}=(\sin\varphi^{\alpha}\cos\theta^{\alpha})\bm{\hat{e}}_{1}+(\sin\varphi^{\alpha}\sin\theta^{\alpha})\bm{\hat{e}}_{2}+(\cos\varphi^{\alpha})\bm{\hat{e}}_{3}$
(A2)

where $\\{\bm{e}\\}=\\{\bm{\hat{e}}_{1},\bm{\hat{e}}_{2},\bm{\hat{e}}_{3}\\}$
denotes the orthonormal basis of the simulation, and the polar, $\theta$, and
azimuthal, $\varphi$, angles were randomly sampled from the uniform
distributions $\theta\in[0,2\pi)$ and $\varphi\in[0,\pi]$, respectively.

For the mesoscale model, only one additional node was appended to each
tethering site, which represents the distal end of the chain. This node was
positioned according to:

$\bm{x}_{m}^{2}=\bm{x}_{m}^{1}+\bm{r}_{m},$ (A3)

where the direction of vector $\bm{r}_{m}$ is again assigned per Eq. (A2), but
the step-size, $|\bm{r}_{m}|$ is randomly sampled from the 3D joint PDF for
finding a Gaussian chain with end-to-end vector, $\bm{r}$, given by
(Rubinstein and Colby, 2003):

$P(\bm{r})=\left(\frac{2}{3}Nb^{2}\pi\right)^{-\frac{3}{2}}\exp\left(-\frac{3\bm{r}^{2}}{2Nb^{2}}\right).$
(A4)

Once initiated, all non-tethered nodes’ positions were updated according to
Eq. (1) and the details of Section 2.1.

###  A.2 Network initiation

To generate stable initial conditions for simulations in Sections 4-5,
particle positions and bond configurations were initiated for both models
using a custom code written in MATLAB2022a. First, a cubic RVE of dimensions
$V=L^{3}$ was sized to achieve the desired polymer packing fraction, $\phi$,
per the relation:

$L=\left\\{\phi^{-1}\left[\frac{\pi}{6}b^{3}n_{p}(2n_{t}-1)N\right]\right\\}^{1/3},$
(A5)

where the quantity $\frac{\pi}{6}b^{3}n_{p}(2n_{t}-1)N$ is the total occupied
volume of polymer assuming that the volume of a single Kuhn segment is
approximately $\pi b^{3}/6$. Additionally, $n_{p}$ (varied) is the number of
polymers in the network, $n_{t}=5$ is the number of tethered side chains for
each polymer, and $N$ is the number of Kuhn segments between crosslinks, so
that $n_{p}(2n_{t}-1)N$ is the total number of Kuhn segments in the network
based on the branched configuration of Fig. 7A. The RVE was positioned with
its center at Cartesian coordinates $\bm{X}_{0}=\left[0,0,0\right]$ and its
initial boundaries spanning in each direction by $\pm L/2$.

Once the RVE was positioned, $n_{p}n_{t}$ particles representing all tethering
sights of side chains on the polymers’ backbones were positioned using a
Poisson growth process. The first particle was positioned at $\bm{X}_{0}$.
Tether particles were then seeded in series at positions
$\bm{X}_{i}=\bm{X}_{i-1}+\bm{u}_{i}$ where $i\in\left[1,n_{p}n_{t}\right]$
denotes the tether index and $\bm{u}_{i}$ is a randomly selected displacement
vector. Direction vectors of $\bm{u}$ were assigned using the spherically
uniform sampling procedure from Eq. (A2). The norms of $\bm{u}$ were randomly
assigned from the uniform distribution
$u\in\left[0.78,1.41\right]c_{t}^{-1/3}$, where $c_{t}=n_{p}n_{t}/V$ is the
tether concentration so that $c_{t}^{-1/3}$ is their nominal separation.
Particles were not permitted to be within a distance of less than
$0.78c_{t}^{-1/3}$ from nearby neighbors. If a particle was seeded outside of
the boundaries of the RVE, a new seeding branch was begun by selecting an
earlier seed particle at random. This process was carried out until the domain
was occupied by $n_{p}n_{t}$ tether particles within the RVE boundaries. Any
particles seeded outside of the RVE were removed. Once all tether crosslink
positions were initiated, their radial distribution was homogenized by
introducing an arbitrary soft, pairwise repulsive force of the form:

$\bm{f}_{r}=E\left(\frac{1}{\sigma_{r}}-\frac{\sigma_{r}^{\alpha}}{d^{\alpha+1}}\right)\hat{\bm{d}},$
(A6)

and then allowing their positions to equilibrate using an overdamped steepest
descent algorithm and periodic boundary conditions according to Wagner et al.
(2021). Here, $E$ is an energy scale to modulate the magnitude of force,
$\sigma_{r}$ is the characteristic length scale over which this force acts,
$\alpha$ is a scaling parameter that controls the stiffness of the force, and
$d$ is the distance between particles. Eq. (A6) is phenomenological and merely
used to erase any process-specific artifacts of the particle initiation
procedure. To homogenize tether positions, we set $E=1.25k_{b}T$,
$\sigma_{r}=1.2c_{t}^{-1/3}$, and $\alpha=2$.

After tether homogenization, we polymerized the backbones of all $n_{p}$
polymers in series. For each polymer, an initial node was selected at random,
and then a chain was added between it and its nearest unbound neighboring
tether while observing RVE periodicity. This was carried out $n_{t}-1$ times
for all $n_{p}$ polymers. To equilibrate polymer segments and mitigate
stochastically initiated chains with lengths in excess of $Nb$, the system was
then equilibrated using the repulsive forces of Eq. (A6) and linear tensile
forces between all attached crosslinks according to:

$\bm{f}_{t}=\frac{3k_{b}T}{Nb^{2}}\bm{r},$ (A7)

where $\bm{r}$ is the end-to-end vector representing each chain segment. While
Eq. (A7) is expressed in physical units, it is phenomenologically applied like
Eq. (A6) simply to erase initiation history. Again, equilibration was carried
out per Wagner et al. (2021) while enforcing periodic boundaries.

Next, side chains were grafted to each of the polymer backbones’ $n_{t}$
tethering sites using the procedure described above via Eqs. (A3) and (A4).
Another step of phenomenological equilibration was then carried out per the
procedure above (i.e., applying the forces of Eqs. (A6) and (A7)) to
homogenize all sticker and tether positions. Finally, $N$ Kuhn segments were
placed at the linearly interpolated positions between the ends of every
attached tether-tether and tether-sticker pair. To ensure stable positioning
of each Kuhn segment with respect to its bonded neighbors (which occurs when
they are approximately a distance of $b$ apart from one another), the soft
pairwise repulsion of Eq. (A6) was again applied but with its characteristic
length scale set to $\sigma_{r}=b$. A final step of equilibration (without any
pairwise tensile forces) was conducted to homogenize the Kuhn segments. This
procedure was executed identically for both the bead-spring and mesoscale
approaches so that initial crosslink distributions were statistically
equivalent between models; however, the mesoscale models’ intermediate Kuhn
segments between stickers and tether crosslinks were retroactively removed and
replaced by direct sticker-to-sticker or sticker-to-tether connections with
implicit pairwise bond potentials through Eq. (4). Although the initiation
procedure was the same between approaches (with the additional step of Kuhn
segment removal for the mesoscale), distinct seeds were used for the mesoscale
and bead-spring models to properly evaluate statistical agreement between
random samples. With the initial chain positions established, the periodic
boundaries were unwrapped in accordance with the requirements of LAMMPS and
input files were automatically generated using MATLAB2022a to enact the
loading criteria described in each section of Sections 4-5.

##  Appendix B Unit conversion and parameter selection

Conversions between the arbitrary units of the discrete model and SI units are
provided in Table B1. Conversions are prescribed directly for the fundamental
units of temperature, time, and length. However, since the model is
overdamped, particle masses were not prescribed and a conversion for the
Boltzmann constant was prescribed instead. These four conversions are used to
derive conversions for other pertinent units such as energy, force, and
stress. While these conversions are provided for reference, results are
generally provided in normalized units throughout the work unless specified
otherwise.

Table B1: Discrete model unit conversions.

Model parameters are listed in Table B2 in both SI and arbitrary model units.

Table B2: Discrete model parameters.

Justifications for the parameter values of Table B2 are as follows:

  * •

Temperature: Temperature was set to 293 K based on typical ambient conditions.

  * •

Kuhn length: The Kuhn length was set to 0.7 nm based on the order of
experimentally estimated values of Kuhn lengths for polymers such as
poly(ethylene glycol) in near-theta solvent (Ahlawat et al., 2021, Liese et
al., 2017, Lee et al., 2008).

  * •

Number of Kuhn segments per polymer chain segment: The number of Kuhn segments
was swept over the range $N\in[12,36]$ (corresponding to chain lengths of
$Nb\in[7.5,22.2]$ nm) per the justification of Section 3.1.

  * •

Molecular weights: Molecular weights, $M_{w}$, were considered for the
branched polymers investigated in Sections 4-5, as related to the selected
values of $N$. In full network-scale studies, in which $n_{t}-1$ chain
segments of $N$ Kuhn segments are linked in series to form a polymer’s
backbone, with $n_{t}$ side-branches grafted to their sides (also of $N$ Kuhn
segments), the total number of Kuhn segments per polymer is
$N_{p}=(2n_{t}-1)N$. Given the estimate that each Kuhn segment is comprised of
two mers (Liese et al., 2017), and the estimated molecular weight of ethylene
glycol is $M_{eg}\approx 44.1$ kg mol-1 (Wagner et al., 2022), then the
molecular weight per polymer chain may be estimated as $M_{w}=2N_{p}M_{eg}$.
Thus, setting $n_{t}=5$ (Fig. 7A), then the molecular weights for polymers
with $N=12$ and $N=36$ Kuhn segments are $M_{w}=9.53$ kDa and $M_{w}=28.6$
kDa, respectively. These are reasonable values as compared to many low
molecular weight polymers and gels.

  * •

Diffusion rate: The diffusion coefficient, $D_{0}$, was used to set the
monomer damping coefficient, $\gamma_{0}=k_{b}T/D_{0}$, and monomer diffusion
timescale, $\tau_{0}=b^{2}/D_{0}$. It was swept over the range
$D_{0}\in\\{0.125,0.25,0.5,1,2,4,8\\}\times 10^{-10}$ m2 s-1, to study its
effects on MSD (Fig. E1). However, it had no meaningful effect on MSD
characteristics and merely renormalized the models’ timescales. Hence,
thereafter it was fixed at $D_{0}=10^{-10}$ m2 s-1 per the justification of
Section 3.1.

  * •

Bond detachment activation energies: Bond detachment activation energy,
$\varepsilon_{d}$, was set as specified within each study of Sections 3-5.
However, it was broadly set as a multiple of $k_{b}T$ in the range
$\varepsilon_{d}\in[2.3,\infty)k_{b}T$ such that
$k_{d}\in[0,0.1]\tau_{0}^{-1}$ through Eq. (7), depending on the need of each
study.

  * •

Bond attachment activation energies: As with $\varepsilon_{d}$, bond
attachment activation energy, $\varepsilon_{a}$, was set as specified within
each study of Sections 3-5 to multiples of $k_{b}T$. Besides when it was swept
over $\varepsilon_{a}\in[0.01,1]k_{b}T$ in Section 3, it was held at
$\varepsilon_{a}=0.01k_{b}T$, begetting a fast intrinsic attachment rate of
$k_{a}^{ap}\approx\tau_{0}^{-1}$ per Eq. (8) and resulting in percolated
network structures.

##  Appendix C Calibrating the bead-spring potential

To verify that the chains of the bead-spring model reproduced the force-
extension relations of the ideal Langevin chains approximated by Eq. (4), we
conducted a simple study in LAMMPS. Chains modeled as $N+1$ beads attached by
springs per the description of Section 2 were modeled at various end-to-end
lengths, $r$, in the range $r\in[1,0.95N]b$, while their average end-to-end
forces were measured. To achieve numerically stable chains with initial end-
to-end separations of approximately $b$, chains were initiated using MATLAB
2022b. One end of each chain was fixed at Cartesian coordinates,
$\bm{x}^{0}=(0,0,0)$, and the other was initially positioned at
$\bm{x}^{N}=(0.95Nb,0,0)$. The positions of $N-1$ beads were then linearly
interpolated between the endpoints, $\bm{x}^{0}$ and $\bm{x}^{N}$, and then
adjacent beads were connected to generate $N$ Kuhn segments. These $N-1$
intermediate beads were then offset in the yz-plane by some amount, $(0,v,w)$,
to break axial symmetry, where $v$ and $w$ are random displacements in the
range $\pm b\sqrt{0.05}$. The positions of the particles were then
equilibrated using arbitrarily soft, pairwise repulsive potentials between all
beads within distance $b$ of each other, and soft harmonic potentials between
all connected beads. The soft pairwise repulsive potential is given by