# Structural features and nonlinear rheology of self-assembled networks of
crosslinked semiflexible polymers

Saamiya Syed  College of Technology, University of Houston, Houston, TX 77204
Center for Theoretical Biological Physics, Rice University, Houston, TX 77005
Fred C. MacKintosh  Center for Theoretical Biological Physics, Rice
University, Houston, TX 77005  Department of Chemical and Biomolecular
Engineering, Rice University, Houston, TX 77005  Department of Chemistry, Rice
University, Houston, TX 77005  Department of Physics & Astronomy, Rice
University, Houston, TX 77005  Jordan L. Shivers  Center for Theoretical
Biological Physics, Rice University, Houston, TX 77005  Department of Chemical
and Biomolecular Engineering, Rice University, Houston, TX 77005

###### Abstract

Random networks of semiflexible filaments are common support structures in
biology. Familiar examples include the fibrous matrices in blood clots,
bacterial biofilms, and essential components of the cells and tissues of
plants, animals, and fungi. Despite the ubiquity of these networks in
biomaterials, we have only a limited understanding of the relationship between
their structural features and highly strain-sensitive mechanical properties.
In this work, we perform simulations of three-dimensional networks produced by
the irreversible formation of crosslinks between linker-decorated semiflexible
filaments. We characterize the structure of networks formed by a simple
diffusion-dependent assembly process and measure their associated steady-state
rheological features at finite temperature over a range of applied prestrains
encompassing the strain-stiffening transition. We quantify the dependence of
network connectivity on crosslinker availability and detail the associated
connectivity dependence of both linear elasticity and nonlinear strain
stiffening behavior, drawing comparisons with prior experimental measurements
of the crosslinker concentration-dependent elasticity of actin gels.

##  I Introduction

The formation of living things involves the energy-intensive assembly of
complex structures from scarce resources. Yet, the demanding physical
processes of life require robust materials that remain intact and functional
under significant and often repetitive applied stresses and strains. Quasi-
one-dimensional or filamentous structures address these challenges efficiently
by supporting significant tensile stresses with minimal material cost.Burla et
al. (2019a) Biological polymers and fibers are generally also
semiflexible,Janmey (1991) meaning that they resist modes of deformation that
induce bending, such as applied compression. In examples spanning a wide range
of length scales, including information-storing DNA and RNA, actin and
intermediate filaments in the cell cytoskeleton, and extracellular collagen
and elastin fibers in tissues, these mechanical features consistently play
crucial roles in biological function.

The cytoskeleton and extracellular matrix are quintessential examples of
disordered networks, a common class of higher-order structure in living
materials. Building upon the qualities of their underlying semiflexible
filaments, these networks act as responsive elastic scaffolds that resist
potentially damaging deformation while leaving ample space for the transport
and storage of functional components, such as interstitial fluids and cells.
Unlike conventional elastic solids, their mechanical properties are scale-
dependentTyznik and Notbohm (2019); Proestaki et al. (2019) and extremely
sensitive to changes in applied stress or strain,Onck et al. (2005); Picu et
al. (2018); Ban et al. (2019) to which they respond with dramatic stiffening,
alignment, and changes in local filament density, enabling essential
biological phenomena such as long-range force transmission by cellsNotbohm et
al. (2015); Wang et al. (2015); Han et al. (2018); Alisafaei et al. (2021);
Grill et al. (2021) and the contraction of muscle fibers.Huxley (1969); Gunst
and Fredberg (2003) Recent work has suggested that the nonlinear viscoelastic
properties of these and related fibrous networks are governed by an underlying
mechanical phase transitionSharma et al. (2016a); Jansen et al. (2018)
associated with the onset of stretching-dominated rigidity under tension-
generating (i.e. shear or extensile) applied strains. Within this framework,
elastic properties are predicted to exhibit a power law dependence on applied
strain in the vicinity of a critical strain, and growing nonaffine
(inhomogeneous) strain fluctuations are expected to drive a significant
slowing of stress relaxation as the transition is approached.Shivers et al.
(2022) The location of this transition (the magnitude of the critical strain)
depends sensitively on key features of the underlying network
architecture.Wyart et al. (2008); Broedersz et al. (2012)

Extensive simulation-based studies have explored the nonlinear rheological
properties of disordered networks of crosslinked stiff or semiflexible
polymers,Onck et al. (2005); Picu (2011); Carrillo et al. (2013); Freedman et
al. (2017) in some cases with realistic three-dimensional geometries produced
either by physical assembly processesHuisman et al. (2007); Kim et al. (2009);
Dobrynin and Carrillo (2011); Müller et al. (2014, 2015); Žagar et al. (2015);
DeBenedictis et al. (2020); Grill et al. (2021) or artificial generation
procedures.Huisman et al. (2008); Lindström et al. (2010); Huisman and
Lubensky (2011); Huisman et al. (2010); Lindström et al. (2013); Amuasi et al.
(2015) However, efforts to specifically connect network structure with strain-
controlled critical behavior have typically focused on simplified random
spring networks.Wyart et al. (2008); Vermeulen et al. (2017); Jansen et al.
(2018); Rens et al. (2018); Merkel et al. (2019); Shivers et al. (2019, 2020);
Burla et al. (2020, 2019b); Dussi et al. (2020); Tauber et al. (2022) In
spring networks, the critical strain coincident with the onset of stretching-
dominated mechanics can be tuned by changing the average connectivity $z$,
defined as the average number of bonds joined at each network junction.Wyart
et al. (2008) For a network of crosslinked filaments, $z$ is controlled by the
typical number of crosslinks formed per filament, approaching an upper bound
of $z\to 4$ in the high crosslink-density limit.Broedersz et al. (2012) In
biopolymer gels, small changes in the concentration of available crosslinkers
can drive dramatic changes in rheological properties,Broedersz and MacKintosh
(2014) including changes in the linear elastic modulus Piechocka et al. (2016)
and shifts in the critical strain corresponding to the onset of stretching-
dominated mechanics.Gardel et al. (2004) These changes are naively consistent
with a tendency of the connectivity to increase with crosslinker
concentration. However, the quantitative relationships between crosslinker
concentration and connectivity, and the structural characteristics of
assembled networks more generally, remain mostly unknown. Improving our
knowledge of how the concentration-dependent microscopic structural details of
self-assembled networks translate into strain-dependent macroscopic
rheological properties is essential for understanding the forces at play in
important biological processes, such as wound healing and cancer metastasis,
and for the design of biomimetic synthetic materials.Romera et al. (2017)

Figure 1:  Model system and rheological approach. (a) Coarse-grained
semiflexible filaments are decorated with randomly assigned sticky sites
(light blue spheres in the image on the right) with coverage fraction $p$.
When two sticky sites meet, they are connected by a permanent crosslink (royal
blue dumbbell). Solutions of these filaments diffusively self-assemble into
percolating disordered networks with macroscopic elasticity. (b) After network
assembly, the structure is fixed and the strain-dependent rheological
properties are determined from the fluctuating shear stress $\sigma(\gamma,t)$
measured under constant macroscopic simple shear strain $\gamma$.

In this study, we consider a system composed of coarse-grained semiflexible
filaments that diffusively self-assemble into a system-spanning network
through the formation of permanent inter-filament crosslinks. We begin with
randomly positioned free filaments with a specified coverage fraction $p$ of
"sticky" or linker-decorated sites. The linker coverage fraction serves as a
proxy for the ratio of crosslinker and filament concentrations in a real
system. Allowing diffusive motion to proceed, we add permanent crosslinks
(short elastic bonds) between sticky sites whose pairwise distance decreases
beneath a designated crosslink formation distance. After the rate of
diffusion-driven formation of new crosslinks becomes sufficiently slow, we
halt the assembly process and analyze the structure of the fixed network,
measuring the dependence of both the connectivity and the average inter-
crosslink contour length on the linker coverage fraction.

We then transition to the rheology stage, in which we observe the network’s
steady state behavior under constant simple shear strain $\gamma$. We obtain
time-series measurements of the thermally fluctuating shear stress
$\sigma(\gamma,t)$ in the mechanically equilibrated state, as shown
schematically in Fig. 1b. Repeating these measurements over a range of strains
for each set of input parameters, we compute relevant elastic quantities such
as the differential shear modulus $K_{\mathrm{eq}}$ and the critical strain
$\gamma_{c}$, at which the macroscopic system transitions between bending-
dominated and stretching-dominated mechanical regimes. Appropriate physical
parameters are chosen to enable comparison with prior experimental
measurements of the elasticity of irreversibly crosslinked networks of
F-actin,Gardel et al. (2004) an essential component of the cytoskeleton of
eukaryotic cells. Using the same time-series measurements, we then
characterize the dynamics of stress relaxation via time correlations in the
stress fluctuations. Building upon recent work,Shivers et al. (2022) we
demonstrate that the excess differential viscosity, a measure of energy
dissipation reflected in a system’s finite-temperature stress fluctuations, is
directly proportional to the same system’s corresponding quasistatic, athermal
differential nonaffinity, a measure of the strain heterogeneity induced by a
small perturbation in the quasistatic, athermal limit. Since, in disordered
networks, the quasistatic nonaffinity is highly strain-dependent and reaches a
maximum at the critical strain, analogous behavior is expected in the excess
differential viscosity and the slowest viscoelastic relaxation time. Our
simulations confirm this expectation, providing crucial insight into
potentially measurable effects of nonaffine fluctuations, which have generally
proven challenging to experimentally quantify.

##  II Model definition and network assembly

We imagine a system that begins as a solution of free semiflexible filaments
covered to some extent with bound linker proteins that are capable of
dimerizing to form permanent elastic crosslinks. In experiments, the linker
coverage could be controlled by varying the relative concentration of
crosslinking protein and filament monomer, as in Ref. 47. If the filament
concentration and linker coverage in such a solution are sufficiently high,
the formation of crosslinks between diffusively migrating filaments eventually
produces a macroscopic network.

To capture this behavior, we turn to a simplified computational model. We
specify the number of filaments $N_{f}$, filament length $\ell_{f}$, and
filament length density $\rho$, which together determine the size of the
periodic simulation box, $L=(N_{f}\ell_{f}/\rho)^{1/3}$. To ensure that the
(initially straight) filaments do not span the entire system, we require
$L>\ell_{f}$, or equivalently $\ell_{f}<\sqrt{N_{f}/\rho}$. We further specify
the number of evenly-spaced nodes per filament, $n$, which determines the
filament bond rest length $\ell_{0}=\ell_{f}/(n-1)$. The total number of nodes
in the system is then $N=N_{f}n$. We then designate a fraction $p$ of the
nodes, chosen randomly, as sticky, or capable of forming a crosslink bond with
another node of the same type. Crosslinks are permanent bonds with a shorter
rest length $\ell_{0,c\ell}$ than the filament bonds. All bonds, including
crosslinks, are treated as harmonic springs with stretching stiffness $\mu$,
and harmonic bending interactions with stiffness $\kappa=k_{B}T\ell_{p}$ act
between adjacent filament bonds. Here, $\ell_{p}$ denotes the filament
persistence length. In terms of the $3N$-dimensional vector of node positions
$\mathbf{x}$, we write the elastic potential energy of this system as

$U(\mathbf{x})=\frac{\mu}{2}\sum_{ij}\frac{\left(\ell_{ij}-\ell_{ij,0}\right)^{2}}{\ell_{ij,0}}+\frac{\kappa}{2}\sum_{ijk}\frac{\left(\theta_{ijk}-\theta_{ijk,0}\right)^{2}}{\ell_{ijk,0}}$
(1)

in which $\ell_{ij}=\left|\mathbf{x}_{j}-\mathbf{x}_{i}\right|$,
$\theta_{ijk}$ is the angle between bonds $ij$ and $jk$,
$\ell_{ijk}=(\ell_{ij}+\ell_{jk})/2$, and the subscript $0$ denotes rest
values. The first sum is taken over all bonds $ij$, including crosslink bonds,
and the second over all pairs of connected bonds $ij$ and $jk$ along the
backbone of each filament. Because we are considering crosslinked networks
with very low filament volume fractions, we deem it acceptable to ignore
steric interactions between filaments. Neglecting inertia, as is appropriate
for the timescales studied here, the system obeys the overdamped Langevin
equation,Allen and Tildesley (2017); Frenkel and Smit (2002) such that the
forces acting on all nodes satisfy

$\mathbf{F}_{N}+\mathbf{F}_{D}+\mathbf{F}_{B}=\mathbf{0},$ (2)

in which the terms on the left represent the network forces, drag forces, and
Brownian forces, respectively. The force due to the network is

$\mathbf{F}_{N}=-\frac{\partial U(\mathbf{x})}{\partial\mathbf{x}}.$ (3)

Nodes are subjected to a Stokes drag force

$\mathbf{F}_{D}=-\zeta\frac{\partial\mathbf{x}}{\partial t}$ (4)

with drag coefficient $\zeta=6\pi\eta_{s}a$, in which $\eta_{s}$ is the
solvent viscosity and $a$ is an effective node radius of half the filament
bond length, $a=\ell_{0}/2$. Finally, the Brownian force is

$\mathbf{F}_{B}=\sqrt{2\zeta k_{B}T}\dot{\mathbf{w}},$ (5)

in which each component of $\dot{\mathbf{w}}(t)$ is a Gaussian random variable
with zero mean and unit variance.Allen and Tildesley (2017) Equation 2 can be
rewritten as

$\frac{\partial\mathbf{x}}{\partial t}=-\frac{1}{\zeta}\frac{\partial
U(\mathbf{x})}{\partial\mathbf{x}}+\sqrt{\frac{2k_{B}T}{\zeta}}\dot{\mathbf{w}}(t).$
(6)

During each timestep in the assembly stage, we check whether the distance
separating any pair of sticky nodes has decreased below a specified crosslink
formation distance $r_{c}$, in which case we connect the two with a crosslink
bond. Each sticky node can form a maximum of one crosslink, and connections
between crosslinks and filaments are treated as freely hinging. In other
words, we do not include a bending potential between adjacent crosslink and
filament bonds. We halt the assembly stage after a total time
$\tau_{\mathrm{a}}$ has elapsed, ensuring that $\tau_{\mathrm{a}}$ is suitably
long to neglect the small fraction of crosslinks forming at longer times.
Here, we use $\tau_{\mathrm{a}}=6\times 10^{7}\Delta t$, corresponding to
approximately one minute in real units. For the parameters specified in Table
S1 (see Supplemental Material), we deem this $\tau_{\mathrm{a}}$ sufficiently
long for the rate of crosslink formation to become negligible, meaning that
the network structure has stabilized. Note that crosslink formation does not
take place after the assembly stage; during the rheology stage, described in
the next section, the network topology remains fixed.

After assembly, we analyze the structural features of each network. To
determine the connectivity $z$, or average number of connections at each
network junction, we consider a reduced version of the simulated network. Each
pair of crosslinked nodes in the original network corresponds to a single node
in the reduced network, and each filament section between two crosslinked
nodes in the original network corresponds to an edge (see sketch in Fig. S1
and further details in Supplemental Material, Section S1). Dangling ends, or
filament sections in the original network connected to only one crosslink, are
therefore neglected. This is a reasonable choice, as dangling ends do not
contribute to the network’s zero-frequency elastic response. We then calculate
the connectivity from the number of edges $n_{\mathrm{edges}}$ and nodes
$n_{\mathrm{nodes}}$ present in the reduced network structure as
$z=2n_{\mathrm{edges}}/n_{\mathrm{nodes}}$.

In practice, we construct systems with filament length
$\ell_{f}=9\;$\mathrm{\SIUnitSymbolMicro m}$$ and filament persistence length
$\ell_{p}=17\;$\mathrm{\SIUnitSymbolMicro m}$$, chosen to approximate F-actin,
with filament length per volume $\rho=2.6\;$\mathrm{\SIUnitSymbolMicro
m}$^{-2}$ (for F-actin, this corresponds to a concentration of
$c_{A}=1.6\;$\mathrm{\SIUnitSymbolMicro}\mathrm{M}$$). Additional parameters
are specified in Table S1 (see Supplemental Material). Unless stated
otherwise, reported measurements throughout this work are averaged over three
randomly generated network samples, and error bars represent $\pm 1$ standard
deviation. All simulations are performed using the open-source molecular
dynamics simulation tool LAMMPS.Thompson et al. (2022)

We first consider the effects of adjusting the crosslinker coverage fraction
$p$ on the assembled network structure. Varying $p$ from $0.4$ to $0.9$, we
find that the connectivity of the fully assembled network structures ranges
from $z\in[2.8,3.6]$, increasing monotonically with $p$ (see Fig. 2a). These
values are similar to those measured for collagen and fibrin networks in prior
experimental work.Lindström et al. (2010); Beroz et al. (2017); Jansen et al.
(2018); Xia et al. (2021) Next, we compute the average inter-crosslink contour
length $\ell_{c}$, which we define as the average length of filament sections
in the original network corresponding to edges in the reduced network. As
shown in Fig. 2b, we find that the average inter-crosslink contour length
$\ell_{c}$ decreases monotonically with $p$.

Figure 2:  Linker coverage-dependence of key structural characteristics of
assembled networks. With increasing linker coverage fraction $p$, we observe
(a) increasing average connectivity $z$ and (b) decreasing average inter-
crosslink contour length $\ell_{c}$.

##  III Nonlinear rheology

After the assembly stage is halted, we proceed to the rheology stage, in which
the strain-dependent steady-state viscoelastic properties of the system are
measured via the time-dependent shear stress as the system fluctuates about
the mechanically equilibrated state under a fixed shear strain. We impose a
constant macroscopic simple shear strain $\gamma$ using Lees-Edwards periodic
boundary conditions Lees and Edwards (1972) and, using the conjugate gradient
method, initially obtain the energy-minimizing configuration of the network,
corresponding to the mechanically equilibrated state at $T=0$. Then, we evolve
the system according to Eq. 6 at temperature $T>0$, specified in Table S1 in
Supplemental Material, for a total run time $\tau_{\mathrm{tot}}=3\times
10^{7}\Delta t$, discarding the first half of the trajectory to avoid
initialization effects. For a given configuration of the system, we compute
the instantaneous virial stress tensor,

$\sigma_{\alpha\beta}=\frac{1}{2V}\sum_{i>j}{f_{ij\alpha}r_{ij\beta}}$ (7)

in which $\mathbf{r}_{ij}=\mathbf{x}_{j}-\mathbf{x}_{i}$ and $\mathbf{f}_{ij}$
is the force acting on node $i$ due to its interaction with node $j$. Because
we will focus on macroscopic simple shear oriented along the $x$-axis with
gradient direction $z$, we now define $\sigma\equiv\sigma_{xz}$ to simplify
notation. Once we have obtained a time series of the shear stress
$\sigma(\gamma,t)$ at strain $\gamma$, we calculate the time-averaged shear
stress $\langle\sigma(\gamma,t)\rangle_{t}$. After repeating this procedure
over many strains in the interval $\gamma\in[0,1]$, we compute the strain-
dependent differential shear modulus,

$K_{\mathrm{eq}}(\gamma)=\frac{\partial\langle\sigma(\gamma,t)\rangle_{t}}{\partial\gamma},$
(8)

which measures the apparent stiffness of the sample under macroscopic strain
$\gamma$ in response to an infinitesimal additional shear strain step. For
sufficiently small strains, this yields the corresponding linear shear
modulus, $G_{0}=\lim_{\gamma\to 0}K_{\mathrm{eq}}(\gamma)$.

Figure 3:  (a) Time-averaged shear stress $\langle\sigma(\gamma,t)\rangle_{t}$
as a function of applied shear strain $\gamma$, for a single network sample
with linker coverage fraction $p=0.9$. (b) Differential shear modulus
$K_{\mathrm{eq}}=\partial\langle\sigma(\gamma,t)\rangle_{t}/\partial\gamma$.
The linear shear modulus $G_{0}$ and critical strain $\gamma_{c}$ are
indicated by dashed lines, and the solid line is a fit to the equation of
state from Ref. 15, described in Supplemental Material, Section S3.

In Fig. 3a and b, we plot the mean stress $\langle\sigma(\gamma,t)\rangle_{t}$
and differential shear modulus $K_{\mathrm{eq}}$ as a function of strain for a
single network sample with the parameters specified in Table S1 (see
Supplemental Material) and linker coverage fraction $p=0.9$. We determine the
critical strain $\gamma_{c}$, which indicates the transition between the
bending-dominated and stretching-dominated mechanical regimes, as the
inflection point of the $\log{K_{\mathrm{eq}}}$ vs. $\log{\gamma}$ curve. To
avoid issues associated with differentiating noisy data, we can alternatively
find $G_{0}$ and $\gamma_{c}$ by fitting the entire $K_{\mathrm{eq}}$ vs.
$\gamma$ curve to an Ising-like equation of state discussed in prior
work,Sharma et al. (2016a) as we describe in further detail in Supplemental
Material, Section S3. Such a fit is shown in Fig. 3b. For the data presented
here, both methods yield effectively equivalent values of $G_{0}$ and
$\gamma_{c}$.

Figure 4:  Effects of linker coverage on features of the linear and nonlinear
stiffness in shear. a) Differential shear modulus $K_{\mathrm{eq}}$ for
networks with varying linker coverage fraction $p$ as a function of strain and
(b) as a function of shear stress. Increasing the linker coverage fraction $p$
drives (c) an increase in the linear modulus $G_{0}=\lim_{\gamma\to
0}K_{\mathrm{0}}$, with $G_{0}\propto p^{3}$, and (d) a decrease in the
critical strain $\gamma_{c}$ as $\gamma_{c}\propto p^{-1}$.

In Fig 4, we report the strain-dependence of the differential shear modulus
for networks with varying linker coverage fraction $p$. As $p$ increases, the
critical strain evidently decreases (stiffening occurs at lower applied
strains) and the linear shear modulus increases. These observations agree
qualitatively with both the observed crosslinker concentration-dependence of
the linear modulus $G_{0}$ and rupture strain $\gamma_{\mathrm{max}}$ (which
we expect to be proportional to $\gamma_{c}$) of F-actin gels reported in Ref.
47, and with the connectivity dependence of the linear modulus and critical
strain observed in spring network simulations.Wyart et al. (2008) The same
differential shear modulus data are plotted in Fig. 4b as a function of
stress. The linear modulus $G_{0}$ and critical strain $\gamma_{c}$ values
extracted from these curves are plotted as functions of the linker coverage
fraction $p$ in Fig. 4c and d, respectively. We find that each exhibits a
power law dependence on $p$, with $G_{0}\propto p^{3}$ and $\gamma_{c}\propto
p^{-1}$. In Ref. 47, in which the analogous quantity $R$ (the ratio between
the concentrations of actin and the crosslinking protein scruin) is varied,
the authors report $G_{0}\propto R^{2}$ and $\gamma_{c}\propto R^{-0.6}$.
These observations may be consistent with ours, provided that the linker
coverage fraction in our system maps to the experimental crosslinker
concentration ratio in Ref. 47 as $p\propto R^{x}$ with $x\sim 0.6$.
Separately, the authors of Ref. 55 reported that, for networks of actin
filaments crosslinked by heavy meromyosin, $G_{0}\propto R^{1.2}$ and
$\gamma_{c}\propto R^{-0.4}$, consistent with our observed dependence of both
quantities on $p$ if $x\sim 0.4$. However, meaningfully relating $p$ and $R$
will require further investigation. Here, for example, no two sticky sites on
the same filament can reside closer than a distance $\ell_{0}$
($1\;$\mathrm{\SIUnitSymbolMicro m}$$ with our parameters) from one another.
This is obviously not the case in real F-actin networks.

Figure 5:  The linear shear modulus $G_{0}$ (a) increases dramatically with
connectivity, shown here for systems with varying linker coverage fraction
$p$, and (b) decreases with increasing average inter-crosslink contour length
$\ell_{c}$. For small $\ell_{c}$ (high $p$), we find
$G_{0}\propto\ell_{c}^{-3}$, consistent with the naive scaling predicted for
affine, isotropic semiflexible filament networks.MacKintosh et al. (1995)
Deviation from this scaling appears at larger $\ell_{c}$. (c) The critical
strain $\gamma_{c}$, in contrast, decreases approximately linearly with $z$. A
linear fit suggests $\gamma_{c}\to 0$ as $z\to z^{*}\sim 4$, roughly
consistent with the upper bound for $z$ in the high crosslink density
limit.Broedersz et al. (2012) (d) The critical strain increases approximately
linearly with increasing average inter-crosslink contour length,
$\gamma_{c}\propto\ell_{c}$, again in apparent agreement with the scaling
expected for the naive affine model.

In Fig. 5, we plot the same extracted linear shear modulus and critical strain
measurements as functions of the $p$-dependent structural quantities discussed
in the previous section, the average network connectivity $z$ and average
inter-crosslink contour length $\ell_{c}$. Notably, we find that the linear
shear modulus $G_{0}$ decreases with increasing inter-crosslink contour length
approximately as $G_{0}\propto\ell_{c}^{-3}$. This exponent is consistent with
the analytically predicted dependence for an affinely deforming isotropic gel
of semiflexible filaments, $G_{0}\propto\rho
k_{B}T\ell_{p}^{2}/\ell_{c}^{3}$.MacKintosh et al. (1995) This agreement may
be coincidental, as, for our range of parameters, the typical inter-crosslink
contour lengths are presumably too small for bending modes to be properly
resolved. This warrants further investigation in simulations with reduced
coarse graining, i.e. smaller $\ell_{0}/\ell_{f}$. In Fig. 5c, we show that
$\gamma_{c}$ approaches 0 as $z$ approaches an intercept value $z^{*}\sim 4$,
in apparent agreement with the upper limit of $z$ for crosslinked networks
with high crosslinking density.Broedersz et al. (2012) However, we expect some
dependence of this limit on the total filament length $\ell_{f}$ and the
bending rigidity $\kappa$, which both necessarily influence the structure of
the assembled network. We also find that the critical strain increases
linearly with the average inter-crosslink contour length $\ell_{c}$ (see Fig.
5d), in apparent agreement with the expected scaling of the critical strain
with contour length for an affinely deforming isotropic gel of semiflexible
filaments, $\gamma_{c}\propto\ell_{c}/\ell_{p}$; as with the dependence of
$G_{0}$ on $\ell_{c}$, this may be a coincidence.

Further information is contained in the fluctuations of the instantaneous
stress about its average value,Wittmer et al. (2013)

$\delta\sigma(\gamma,t)=\sigma(\gamma,t)-\langle\sigma(\gamma,t)\rangle_{t},$
(9)

from which we compute the stress fluctuation autocorrelation function,

$C(\gamma,\tau)=\frac{V}{k_{B}T}\langle\delta\sigma(\gamma,t)\delta\sigma(\gamma,t+\tau)\rangle_{t},$
(10)

which reveals useful details about the time-dependence of energy
dissipation.Wittmer et al. (2015a, b)

Figure 6:  The relaxation of stress fluctuations slows dramatically near the
critical strain as rearrangements become increasingly nonaffine. (a) The shear
stress fluctuation autocorrelation function $C(\gamma,\tau)$ measured at
representative strains below (purple), near (red), and above (orange) the
critical strain $\gamma_{c}$ (see the corresponding $K_{\mathrm{eq}}$ data in
the inset of panel (b)) reveals this slowing, which is captured quantitatively
by (b) a peak in the apparent slowest relaxation time $\tau_{R}$ (Eq. 11).
These data correspond to the same sample as in Fig. 3 (see inset), and the
labeled strains are indicated by the same colors in all panels. (c) At long
lag times, the time-dependent apparent excess viscosity
$\delta\eta(\gamma,\tau)$ computed with Eq. 12 grows much larger near
$\gamma_{c}$ than elsewhere. This behavior is especially clear when (d) the
apparent excess zero-shear viscosity $\delta\eta(\tau=\tau_{\mathrm{max}})$ is
plotted as a function of $\gamma$, revealing a peak at the $p$\- (and thus
$z$-) dependent critical strain. These observations are supported by
independent measurements of the strain-dependent quasistatic, athermal
differential nonaffinity $\delta\Gamma_{\infty}$, which is quantitatively
related to $\delta\eta$ by Eq. 14 (solid blue line).

In Fig. 6a, we show representative $C(\gamma,\tau)$ data for a single sample
(the same sample as in Fig. 3) under applied strains below, near, and above
the critical strain $\gamma_{c}$. We find that the stress fluctuations decay
slowly when the macroscopic applied strain is near $\gamma_{c}$, in contrast
to a much faster decay for strains below or above the critical regime. We can
quantify this strain-dependent change in relaxation dynamics by integrating
the shear stress autocorrelation function $C(\gamma,\tau)$ over a sufficiently
long range of lag times $\tau$. In practice, the range of lag times for which
we can reliably estimate $C(\gamma,\tau)$ is limited by the simulation run
time. Assuming $C(\gamma,\tau)$ is known for lag times below a maximum
$\tau_{\mathrm{max}}$, we can estimate the system’s slowest relaxation time as

$\tau_{\mathrm{R}}(\gamma,\tau_{\mathrm{max}})=\frac{\int_{0}^{\tau_{\mathrm{max}}}{\tau
C(\gamma,\tau)d\tau}}{\int_{0}^{\tau_{\mathrm{max}}}{C(\gamma,\tau)d\tau}}$
(11)

In the limit $\tau_{\mathrm{max}}\to\infty$, this converges to the true
slowest relaxation time $\tau_{\mathrm{R,0}}(\gamma)$. In Fig. 6b, we plot the
apparent slowest relaxation time $\tau_{\mathrm{R}}$ for the same system, with
marker colors indicating the strains plotted in Fig. 6a. Note that the
corresponding stiffening curve from Fig. 3 is also shown in the inset. We
observe that $\tau_{\mathrm{R}}$ grows substantially as the critical strain
$\gamma_{c}$ is approached from either side, reaching a maximum at
$\gamma_{c}$. In fact, we see consistent behavior – that is, growth of the
slowest relaxation time by an order of magnitude or more with a peak at
$p$-dependent critical strain $\gamma_{c}$ – in networks over the entire range
of $p$ considered. This is shown in Fig. 7d, in which the applied strain is
normalized by the $p$-dependent critical strain $\gamma_{c}$, revealing maxima
in $\tau_{\mathrm{R}}$ at $\gamma/\gamma_{c}=1$ for all $p$.

Figure 7:  Shared features of stiffening and slowing down at the structure-
dependent critical strain. (a) For networks with varying crosslinker coverage
fraction $p$, the stiffening regime of the normalized differential shear
modulus $K/G_{0}$ collapses onto a single curve when plotted vs.
$\gamma/\gamma_{c}$. Lines correspond to fits to the equation of state
described in Ref. 15 (see Supplemental Material, Section S3), and the
associated fit parameters are plotted as a function of $p$ in Fig. S2. (b) The
mean-squared stress fluctuations
$\langle\delta\sigma(\gamma,t)^{2}\rangle_{t}$ are largest at the $p$\- (and
thus $z$-) dependent critical strain. Solid lines correspond to the difference
between the measured strain-dependent affine and equilibrium differential
shear moduli, $K_{\mathrm{aff}}$ and $K_{\mathrm{eq}}$, according to Eq. 16.
Both the (a) apparent excess zero-shear viscosity $\delta\eta$ and (b)
apparent slowest relaxation time $\tau_{\mathrm{R}}$ likewise exhibit large
peaks at the $p$-dependent critical strain $\gamma_{c}$.

The stress fluctuation autocorrelation function is likewise related to the
system’s lag-dependent excess differential viscosity Wittmer et al. (2015a)
$\delta\eta(\gamma,\tau)\equiv\eta(\gamma,\tau)-\eta_{s}$ measured at strain
$\gamma$,

$\delta\eta(\gamma,\tau)=\int_{0}^{\tau}C(\gamma,\tau^{\prime})d\tau^{\prime},$
(12)

which we plot in Fig. 6c for the same labeled strain values. It is clear that
$\delta\eta(\gamma,\tau)$ grows far more dramatically near the critical strain
than elsewhere. In recent work,Shivers et al. (2022) it was suggested that the
low-frequency or “zero-shear” excess differential viscosity
$\delta\eta_{0}(\gamma)=\lim_{\tau\to\infty}\delta\eta(\gamma,\tau)$ of
prestrained disordered networks is controlled by nonaffinity. Specifically,
$\delta\eta_{0}$ was shown to be related to the quasistatic, athermal
differential nonaffinity $\delta\Gamma_{\infty}(\gamma)$, defined as

$\delta\Gamma_{\infty}(\gamma)=\lim_{\delta\gamma\to
0}\frac{1}{N\ell_{0}^{2}\delta\gamma^{2}}\sum_{i}\left|\delta\mathbf{x}_{i}^{\mathrm{NA}}\right|^{2}$
(13)

in which the sum is taken over all network nodes and the vector
$\delta\mathbf{x}_{i}^{\mathrm{NA}}$ represents the nonaffine component of the
displacement vector of node $i$ under a macroscopic, quasistatically applied
strain step $\delta\gamma$. The viscosity-nonaffinity relationship in Ref. 17
can be stated concisely as

$\delta\eta_{0}(\gamma)=\frac{N}{V}\zeta\ell_{0}^{2}\delta\Gamma_{\infty}(\gamma).$
(14)

The left side of Eq. 14 reflects the dynamics of the stress fluctuations at
finite temperature, while the right side reflects the heterogeneous nature of
the strictly quasistatic deformation field, which can be obtained by comparing
energy-minimized, $T=0$ system configurations under varying applied $\gamma$.

Simulations have suggested that, for disordered filament networks, the
differential nonaffinity generically reaches a maximum at the structure-
dependent critical strain,Rens (2014); Sharma et al. (2016b) at which the
system macroscopically transitions between bending-dominated (or floppy, for
$\kappa=0$) and stretching-dominated regimes. Thus, according to Eq. 14, we
should generically see a proportional peak in the excess differential
viscosity at $\gamma_{c}$. Testing this prediction requires measuring the
quasistatic, athermal differential nonaffinity. To do so, for a given network
at strain $\gamma$, we first obtain the energy-minimizing configuration at
$T=0$ using the conjugate gradient method. To this configuration, we then
apply a small additional strain step $d\gamma=0.01$, after which we repeat the
energy minimization procedure. Comparing the positions in the energy-
minimizing configurations at $\gamma$ and $\gamma+d\gamma$, we compute the
quasistatic differential nonaffinity $\delta\Gamma_{\infty}$ using Eq. 13. In
Fig. 6d, we plot $\delta\Gamma_{\infty}$ along with the apparent excess
differential zero-shear viscosity for the same system, demonstrating excellent
agreement with Eq. 14. To demonstrate that this behavior is preserved as the
structure is varied, we show in Fig. 7c that a clear peak in $\delta\eta$
occurs at $\gamma/\gamma_{c}=1$ over the entire range of $p$ considered.

We can also use the stress fluctuations to determine the differential
relaxation modulus,Wittmer et al. (2015c)

$K(\gamma,\tau)=K_{\mathrm{eq}}(\gamma)+C(\gamma,\tau),$ (15)

which quantifies the time-dependent apparent stiffness of the system at
prestrain $\gamma$ in response to an instantaneous additional strain
step.Rubinstein and Colby (2006); Doi and Edwards (1988) For sufficiently
short times, the differential relaxation modulus approaches an upper limit of
$\lim_{\tau\to 0}K(\gamma,\tau)=K_{\mathrm{aff}}(\gamma)$, corresponding to
the apparent stiffness of the energy-minimized, athermal equivalent of the
system at strain $\gamma$ under an instantaneous, homogeneous infinitesimal
shear strain step. Consequently, the equilibrium stiffness of the system under
applied strain $\gamma$ is simply $K_{\mathrm{aff}}(\gamma)$ reduced by
$C(\gamma,0)$,Squire et al. (1969); Lutsko (1989); Wittmer et al. (2013) or
equivalently

$K_{\mathrm{eq}}(\gamma)=K_{\mathrm{aff}}(\gamma)-\frac{V}{k_{B}T}\langle\left(\delta\sigma(\gamma,t)\right)^{2}\rangle_{t}.$
(16)

We find that this relationship provides a useful estimate of the mean squared
stress fluctuations, as shown in Fig. 7b, and we observe that the mean squared
stress fluctuations robustly reach a maximum at $\gamma/\gamma_{c}=1$.

##  IV Conclusions

We have investigated the assembly and mechanical testing of disordered
networks of crosslinked semiflexible polymers via Brownian dynamics
simulations. Such networks serve as essential mechanical constituents of a
wide variety of biological materials spanning many length scales. We explored
the structural and rheological consequences of varying the crosslinker
coverage fraction $p$, an analog of the experimental ratio between the
concentrations of crosslinker and filament proteins. Using physical parameters
intended to mimic the cytoskeletal polymer F-actin, we measured the effects of
increasing linker coverage fraction $p$ on the average connectivity $z$ of
self-assembled networks, observing a corresponding increase in connectivity
between $z\in[2.8,3.6]$.

We then investigated the relationship between these changes in connectivity
and associated changes in the strain-dependent rheological properties of
assembled networks, obtaining extended time series measurements of the
fluctuating shear stress $\sigma(\gamma,t)$ for systems held under fixed shear
strain $\gamma$. Analyzing many such trajectories gathered over a range of
applied strains $\gamma$, we computed the strain-dependent differential shear
modulus $K_{\mathrm{eq}}$, from which we extracted both the linear shear
modulus $G_{0}$ and the critical strain $\gamma_{c}$. We described the
dependence of these quantities on the crosslinker coverage fraction $p$,
demonstrating qualitative agreement with the experimentally observed
dependence of the same quantities on crosslinker concentration ratio $R$
reported in Ref. 47. Specifically, we found that increasing $p$ produces
inherently stiffer networks (having an increased linear shear modulus $G_{0}$)
that simultaneously exhibit an earlier tendency to strain-stiffen (having a
lower critical strain $\gamma_{c}$). Drawing comparisons between the scaling
of both $G_{0}$ and $\gamma_{c}$ with $p$ and the analogous experimentally-
observed scaling of the same quantities with $R$, we suggested a simple power
law relationship between $p$ and $R$. We next described the apparent
dependence of these elastic features on the $p$-dependent structural
quantities $z$ and $\ell_{c}$, notably showing that the critical strain
decreases linearly with $z$. A linear fit suggests that $\gamma_{c}$
approaches $0$ as the connectivity $z$ approaches a limiting value $z^{*}\to
z_{c}$ near $4$, the theoretical upper bound for $z$ at high crosslinking
density.

We then extended our observations beyond strictly elastic properties by
analyzing the stress fluctuation autocorrelation function $C(\gamma,\tau)$,
which revealed the development of extremely slow dynamics in systems subjected
to applied shear strains near the critical strain $\gamma_{c}$. From
$C(\gamma,\tau)$, we obtained estimates of the slowest relaxation time
$\tau_{\mathrm{R},0}$ and the excess differential zero-shear viscosity
$\delta\eta_{0}$, both of which we showed are consistently maximized at the
$p$-dependent critical strain, irrespective of the details of the underlying
network. Building upon results from Ref. 17, we demonstrated that the excess
differential viscosity in these finite temperature systems is controlled
quantitatively by the athermal, quasistatic differential nonaffinity
$\delta\Gamma_{\infty}$, a measure of the inherent tendency of the strained
network to deform heterogeneously. Importantly, our results suggest that one
should expect measurable dynamical signatures of transition-associated
nonaffine fluctuations to appear in semiflexible polymer networks with
biologically relevant elastic properties, e.g. those of the F-actin
cytoskeleton, under physiologically relevant applied strains. In other words,
one should generically expect to observe slowing stress relaxation in
biopolymer networks under prestrain levels near the critical strain
$\gamma_{c}$ marking the macroscopic transition between bending-dominated and
stretching-dominated elasticity. Since $\gamma_{c}$ is controlled by the
average connectivity $z$, which is controlled in turn by the availability of
crosslinking sites (here, $p$), our results suggest a route to experimentally
control the strain-dependence of major features of both network elasticity and
the dynamics of stress relaxation.

###### Acknowledgements.

This study was supported in part by the National Science Foundation Division
of Materials Research (Grant No. DMR-1826623) and the National Science
Foundation Center for Theoretical Biological Physics (Grant No. PHY-2019745).
This work began as a summer research experience for undergraduates (REU)
project supported by the Frontiers in Science (FIS) program of the Center for
Theoretical Biological Physics.

## References

  * Burla et al. (2019a) F. Burla, Y. Mulla, B. E. Vos, A. Aufderhorst-Roberts, and G. H. Koenderink, Nature Reviews Physics 1, 249 (2019a).
  * Janmey (1991) P. A. Janmey, Current Opinion in Cell Biology 3, 4 (1991).
  * Tyznik and Notbohm (2019) S. Tyznik and J. Notbohm, Mechanics of Materials 138, 103155 (2019).
  * Proestaki et al. (2019) M. Proestaki, A. Ogren, B. Burkel, and J. Notbohm, Experimental Mechanics 59, 1323 (2019).
  * Onck et al. (2005) P. R. Onck, T. Koeman, T. van Dillen, and E. van der Giessen, Physical Review Letters 95, 178102 (2005).
  * Picu et al. (2018) R. C. Picu, S. Deogekar, and M. R. Islam, Journal of Biomechanical Engineering 140, 021002 (2018).
  * Ban et al. (2019) E. Ban, H. Wang, J. M. Franklin, J. T. Liphardt, P. A. Janmey, and V. B. Shenoy, Proceedings of the National Academy of Sciences 116, 6790 (2019).
  * Notbohm et al. (2015) J. Notbohm, A. Lesman, P. Rosakis, D. A. Tirrell, and G. Ravichandran, Journal of The Royal Society Interface 12, 20150320 (2015).
  * Wang et al. (2015) H. Wang, A. S. Abhilash, C. S. Chen, R. G. Wells, and V. B. Shenoy, Biophysical Journal 107, 2592 (2015).
  * Han et al. (2018) Y. L. Han, P. Ronceray, G. Xu, A. Malandrino, R. D. Kamm, M. Lenz, C. P. Broedersz, and M. Guo, Proceedings of the National Academy of Sciences 115, 4075 (2018).
  * Alisafaei et al. (2021) F. Alisafaei, X. Chen, T. Leahy, P. A. Janmey, and V. B. Shenoy, Soft Matter 17, 241 (2021).
  * Grill et al. (2021) M. J. Grill, J. Kernes, V. M. Slepukhin, W. A. Wall, and A. J. Levine, Soft Matter 17, 10223 (2021).
  * Huxley (1969) H. E. Huxley, Science 164, 1356 (1969).
  * Gunst and Fredberg (2003) S. J. Gunst and J. J. Fredberg, Journal of Applied Physiology 95, 413 (2003).
  * Sharma et al. (2016a) A. Sharma, A. J. Licup, K. A. Jansen, R. Rens, M. Sheinman, G. H. Koenderink, and F. C. MacKintosh, Nature Physics 12, 584 (2016a).
  * Jansen et al. (2018) K. A. Jansen, A. J. Licup, A. Sharma, R. Rens, F. C. MacKintosh, and G. H. Koenderink, Biophysical Journal 114, 2665 (2018).
  * Shivers et al. (2022) J. L. Shivers, A. Sharma, and F. C. MacKintosh, arXiv:2203.04891 (2022).
  * Wyart et al. (2008) M. Wyart, H. Liang, A. Kabla, and L. Mahadevan, Physical Review Letters 101, 215501 (2008).
  * Broedersz et al. (2012) C. P. Broedersz, M. Sheinman, and F. C. MacKintosh, Physical Review Letters 108, 078102 (2012).
  * Picu (2011) R. C. Picu, Soft Matter 7, 6768 (2011).
  * Carrillo et al. (2013) J.-M. Y. Carrillo, F. C. MacKintosh, and A. V. Dobrynin, Macromolecules 46, 3679 (2013).
  * Freedman et al. (2017) S. L. Freedman, S. Banerjee, G. M. Hocky, and A. R. Dinner, Biophysical Journal 113, 448 (2017).
  * Huisman et al. (2007) E. M. Huisman, T. van Dillen, P. R. Onck, and E. Van der Giessen, Physical Review Letters 99, 208103 (2007).
  * Kim et al. (2009) T. Kim, W. Hwang, H. Lee, and R. D. Kamm, PLOS Computational Biology 5, e1000439 (2009).
  * Dobrynin and Carrillo (2011) A. V. Dobrynin and J. M. Y. Carrillo, Macromolecules 44, 140 (2011).
  * Müller et al. (2014) K. W. Müller, R. F. Bruinsma, O. Lieleg, A. R. Bausch, W. A. Wall, and A. J. Levine, Physical Review Letters 112, 238102 (2014).
  * Müller et al. (2015) K. W. Müller, C. Meier, and W. A. Wall, Journal of Computational Physics 303, 185 (2015).
  * Žagar et al. (2015) G. Žagar, P. R. Onck, and E. Van Der Giessen, Biophysical Journal 108, 1470 (2015).
  * DeBenedictis et al. (2020) E. P. DeBenedictis, Y. Zhang, and S. Keten, Macromolecules 53, 6123 (2020).
  * Huisman et al. (2008) E. M. Huisman, C. Storm, and G. T. Barkema, Physical Review E 78, 051801 (2008).
  * Lindström et al. (2010) S. B. Lindström, D. A. Vader, A. Kulachenko, and D. A. Weitz, Physical Review E 82, 051905 (2010).
  * Huisman and Lubensky (2011) E. M. Huisman and T. C. Lubensky, Physical Review Letters 106, 088301 (2011).
  * Huisman et al. (2010) E. M. Huisman, C. Storm, and G. T. Barkema, Physical Review E 82, 061902 (2010).
  * Lindström et al. (2013) S. B. Lindström, A. Kulachenko, L. M. Jawerth, and D. A. Vader, Soft Matter 9, 7302 (2013).
  * Amuasi et al. (2015) H. E. Amuasi, C. Heussinger, R. L. Vink, and A. Zippelius, New Journal of Physics 17, 083035 (2015).
  * Vermeulen et al. (2017) M. F. Vermeulen, A. Bose, C. Storm, and W. G. Ellenbroek, Physical Review E 96, 053003 (2017).
  * Rens et al. (2018) R. Rens, C. Villarroel, G. Düring, and E. Lerner, Physical Review E 98, 62411 (2018).
  * Merkel et al. (2019) M. Merkel, K. Baumgarten, B. P. Tighe, and M. L. Manning, Proceedings of the National Academy of Sciences 116, 6560 (2019).
  * Shivers et al. (2019) J. L. Shivers, J. Feng, A. Sharma, and F. C. MacKintosh, Soft Matter 15, 1666 (2019).
  * Shivers et al. (2020) J. L. Shivers, S. Arzash, and F. C. MacKintosh, Physical Review Letters 124, 038002 (2020).
  * Burla et al. (2020) F. Burla, S. Dussi, C. Martinez-Torres, J. Tauber, J. van der Gucht, and G. H. Koenderink, Proceedings of the National Academy of Sciences of the United States of America 117, 8326 (2020).
  * Burla et al. (2019b) F. Burla, J. Tauber, S. Dussi, J. van der Gucht, and G. H. Koenderink, Nature Physics 15, 549 (2019b).
  * Dussi et al. (2020) S. Dussi, J. Tauber, and J. van der Gucht, Physical Review Letters 124, 018002 (2020).
  * Tauber et al. (2022) J. Tauber, J. van der Gucht, and S. Dussi, The Journal of Chemical Physics 156, 160901 (2022).
  * Broedersz and MacKintosh (2014) C. P. Broedersz and F. C. MacKintosh, Reviews of Modern Physics 86, 995 (2014).
  * Piechocka et al. (2016) I. K. Piechocka, K. A. Jansen, C. P. Broedersz, N. A. Kurniawan, F. C. MacKintosh, and G. H. Koenderink, Soft Matter 12, 2145 (2016).
  * Gardel et al. (2004) M. L. Gardel, J. H. Shin, F. C. MacKintosh, L. Mahadevan, P. Matsudaira, and D. A. Weitz, Science 304, 1301 (2004).
  * Romera et al. (2017) M. F.-C. Romera, R. P. M. Lafleur, C. Guibert, I. K. Voets, C. Storm, and R. P. Sijbesma, Angewandte Chemie International Edition 56, 8771 (2017).
  * Allen and Tildesley (2017) M. P. Allen and D. J. Tildesley, _Computer Simulation of Liquids_ (Oxford University Press, 2017).
  * Frenkel and Smit (2002) D. Frenkel and B. Smit, _Understanding Molecular Simulation: From Algorithms to Applications_ (Academic Press, 2002).
  * Thompson et al. (2022) A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, et al., Computer Physics Communications 271, 108171 (2022).
  * Beroz et al. (2017) F. Beroz, L. M. Jawerth, S. Münster, D. A. Weitz, C. P. Broedersz, and N. S. Wingreen, Nature Communications 8, 16096 (2017).
  * Xia et al. (2021) J. Xia, L.-H. Cai, H. Wu, F. C. MacKintosh, and D. A. Weitz, Proceedings of the National Academy of Sciences 118, e2020541118 (2021).
  * Lees and Edwards (1972) A. W. Lees and S. F. Edwards, Journal of Physics C: Solid State Physics 5, 1921 (1972).
  * Tharmann et al. (2007) R. Tharmann, M. M. A. E. Claessens, and A. R. Bausch, Physical Review Letters 98, 88103 (2007).
  * MacKintosh et al. (1995) F. C. MacKintosh, J. Kas, and P. A. Janmey, Physical Review Letters 75, 4425 (1995).
  * Wittmer et al. (2013) J. P. Wittmer, H. Xu, P. Polińska, F. Weysser, and J. Baschnagel, The Journal of Chemical Physics 138, 12A533 (2013).
  * Wittmer et al. (2015a) J. P. Wittmer, H. Xu, O. Benzerara, and J. Baschnagel, Molecular Physics 113, 2881 (2015a).
  * Wittmer et al. (2015b) J. P. Wittmer, H. Xu, and J. Baschnagel, Physical Review E 91, 022107 (2015b).
  * Rens (2014) R. Rens, PhD Thesis, University of Amsterdam (2014).
  * Sharma et al. (2016b) A. Sharma, A. J. Licup, R. Rens, M. Vahabi, K. A. Jansen, G. H. Koenderink, and F. C. MacKintosh, Physical Review E 94, 042407 (2016b).
  * Wittmer et al. (2015c) J. P. Wittmer, I. Kriuchevskyi, J. Baschnagel, and H. Xu, The European Physical Journal B 88, 242 (2015c).
  * Rubinstein and Colby (2006) M. Rubinstein and R. Colby, _Polymer Physics_ (Oxford University Press, 2006).
  * Doi and Edwards (1988) M. Doi and S. F. Edwards, _The Theory of Polymer Dynamics_ (Oxford University Press, 1988).
  * Squire et al. (1969) D. R. Squire, A. C. Holt, and W. G. Hoover, Physica 42, 388 (1969).
  * Lutsko (1989) J. F. Lutsko, Journal of Applied Physics 65, 2991 (1989).