\section{Introduction}

Living cells are known to exhibit unusual mechanical properties including an internal nonlinear stiffening under external stress, in which their stiffness can increase by orders of magnitudes \cite{Trepat2007592}, reversible softening under compression~\cite{Chaudhuri2007}, viscoplasticity~\cite{Bonakdar2016} and poroelasticity~\cite{Kimpton2015}. A long-standing puzzle is related to the surprisingly slow stress relaxation that has been measured in living cells\SC{~\cite{Navarro1996, Fabry2001,Stamenovic2004,Desprat20052224,Bursac2005557,Balland2006,Stamenovic2007}}. 
This stress relaxation reveals more than just a long relaxation time, but also a broad spectrum of relaxation times, with
a dynamic modulus that varies with frequency as a power-law with exponent $\beta$ in the range of $\sim0.1- 0.3$~\cite{Fabry2001}.
It has been argued that this may be related to the soft glassy rheology (SGR) model \cite{Sollich1998,Fabry2001,Bursac2005557}, although the relevance and validity of this in cell mechanics remains unclear. 

Most of the mechanical properties of living cells originate in the cytoskeleton, a dynamic network composed of crosslinked biopolymers, which gives the cell its shape and rigidity~\cite{Hall1998}. One reason for its dynamic nature is that many of the crosslinking proteins, for example $\alpha$-actinin, form transient bonds with both finite binding and unbinding rates~\cite{Schiffhauer20161473}. Such crosslinking proteins, denoted as transient crosslinkers, introduce a distinct type of stress relaxation in semiflexible polymer networks, since the unbinding of crosslinkers allows the networks to flow at long time~\cite{Ward20084915,Lieleg2008,Lieleg20094725,Broedersz2010,Yao2013,Muller2014}. Previous theory and experiments involving reconstituted biopolymer networks have revealed a characteristic scaling exponent of $\beta=1/2$ for the frequency-dependent linear shear modulus~\cite{Broedersz2010,Yao2013,Muller2014}. Unlike the Rouse model, in which the same exponent appears in the high-frequency regime~\cite{RouseJr.19531272,gennes_1990,doi}, transient networks only exhibit the $1/2$ exponent in the low-frequency regime, indicating a different mechanism of stress relaxation from the Rouse model. Moreover, the Rouse model applies to flexible polymers, while the analogous high-frequency regime for semiflexible polymers such as actin is known theoretically and experimentally to exhibit a $3/4$ exponent~\cite{Gittes1997,Gittes1998,Morse1998, Morse19987030,Koenderink2006,Broedersz2014, Pritchard2014}. 

Within the cytoskeleton there are also molecular motors that generate internal stresses~\cite{Alberts,Bray2001,Markovich2019}, which may alter the rheological properties of living cells~\cite{Wang2002, Stamenovic2004, FERNANDEZ2006}.  Recently, it has been shown experimentally that the apparent scaling exponent of the linear shear modulus in reconstituted actin networks with transient crosslinks can be further reduced and systematically varied over the range of $0.1\lesssim\beta\le0.5$ by an applied external stress~\cite{Mulla2019}. 
\SC{It is well known that applying external or internal stress on a permanent biopolymer networks can cause nonlinear stiffening~\cite{Gardel1301,Storm2005191,Mizuno2007,Koenderink2009} and a reduction in the high-frequency exponent from 3/4 to 1/2~\cite{ Granek1997,Caspi1998,Rosenblatt2006, Mizuno2007,Majumdar2008}}. When applied to transient networks, aside from the stress-stiffening response, external stress  solidifies the network~\cite{Yao2013}, suppressing both stress relaxation and reducing the apparent exponent of the frequency-dependent nonlinear shear modulus~\cite{Mulla2019}. This provides a possible explanation for the weak scaling exponent observed in living cells, as living cells are intrinsically under internal stress generated by molecular motors or other active processes~\cite{Alberts,Julicher20073,Roh-Johnson20121232, Chen2020}, although a microscopic model for how stress qualitatively changes the dynamics is still lacking.

In this paper, we develop such microscopic theory for the mechanical response of transient-crosslinked semiflexible networks under applied stress. Based on equilibrium thermodynamics, we show that the dynamics of transiently crosslinked semiflexible polymer networks can be described by a mode-dependent mobility. We analytically derive the form of this mobility in transient networks, and calculate the nonlinear modulus under external stress. Our theoretical prediction of the nonlinear modulus quantitatively agrees with experiment data of reconstituted cytoskeletal networks. We also show that external stress naturally leads to a weak frequency dependence by suppressing fluctuations of bending modes. 

\begin{figure}[h]
	\centering
	\includegraphics[scale=0.35]{illustration.pdf}
	\caption{(a) Schematic figure of a transient network, \SC{in which polymers (lines) are connected by transient crosslinkers (circles). The unbinding and rebinding of a single crosslinker (purple circle) relaxes its adjacent segments.} (b) The relaxation of a long-wavelength mode requires successive unbinding events (\SC{green circle -> blue circle -> brown circle in sequence }). }
	\label{Fig.1}
\end{figure}


\section{Overview}
\SC{Transient crosslinkers introduce a distinct type of stress relaxation: when a crosslinker unbinds, it relaxes the stress on adjacent polymer segments (the crosslinker later bind on a different position), see Fig.~\ref{Fig.1}(a). In order for a long polymer to relax, multiple successive unbinding events are required (Fig.~\ref{Fig.1}(b)), resulting in a relaxation time much longer than the timescale of a single unbinding event~\cite{Broedersz2010,Broedersz2014}. } 


\SC{We propose a microscopic model that accounts for the effect of unbinding/rebinding of transient crosslinkers.} In our model, that is described in detail below, we show that the relaxation of semiflexible polymer networks can be decomposed into the relaxation of the independent bending modes (see Sec.~\ref{section:3}). The relaxation of each mode follows a mode-dependent mobility, $M_{qq}$, \SC{($q$ being the wave number)}  which leads to non-trivial dependences on the frequency of the linear and nonlinear viscoelastic moduli. 
In Fig.~\ref{Fig.rgm} (a) we show a schematic diagram of the various regimes  for a transient-crosslinked biopolymer network in the ($\omega,\sigma$) phase-space, where $\omega$ is the frequency and $\sigma$ is the applied shear stress. With only transient crosslinkers in the system, reptation of finite molecular weight polymers is expected to lead to liquid-like behavior on the longest timescales. In the present model we have focused on stress relaxation that is entirely due to transient-crosslinking, and we therefore consider the limit of high molecular weight and timescales less than the reptation time. On these timescales, we find that the stress relaxation can be devided into two regimes: in the low-frequency regime the stress relaxation is governed by the transient behavior of the crosslinkers, while in the high-frequency regime the network behaves as if the crosslinkers were permanent~\cite{Broedersz2010}. The two regimes are separated by a characteristic frequency, $\omega_c(\sigma)$, which depends on $\sigma$ (see Sec.~\ref{section:5}, Eq.~(\ref{e37})).  These two different regimes appear as a result of the different $q$-dependence of the mode-dependent mobility $M_{qq}$, for $q>\pi/\ell_c$ and $q<\pi/\ell_c$, where $\ell_c$ is the average spacing between crosslinkers in the network. The relaxation of bending modes with $q>\pi/\ell_c$ is dominated by the solvent viscosity, as in permanent networks~\cite{Gittes1998, Morse1998}. However, the relaxation of bending modes with $q<\pi/\ell_c$ is limited by the transient nature of the  crosslinkers, leading to $M_{qq}\sim q^2$ (see Sec.~\ref{section:4}). This quadratic dependence in $M_{qq}$ results in a linear modulus $G(\omega) \sim \omega^{1/2}$. Moreover,  when  external stress is applied, the network may stiffen nonlinearly, where the characteristic stress, $\sigma_c(\omega)\sim \omega^{1/3}$, governs the transition from the linear to the nonlinear stiffening regimes (see Sec.~\ref{section:5}, Eq.~(\ref{e22})). When the network stiffens nonlinearly, \SC{the differential shear modulus $K=K'+iK''=\diff \sigma/\diff \gamma $ is used to characterize the viscoelastic behavior, where $\gamma$ is the shear deformation. }~\SC{As we show below (see Sec.~\ref{section:6}) and sketched Fig.~\ref{Fig.rgm}(b), the apparent exponent with which $G$ depends on $\omega$ is reduced with applied stress and can become arbitrarily small at high stress. This is consistent with recent experiments~\cite{Mulla2019} and may provide an explanation for the week dependence seen in living cells~\cite{Navarro1996, Fabry2001,Stamenovic2004,Desprat20052224,Bursac2005557,Balland2006,Stamenovic2007}.} 
\begin{figure}[t]
	\centering
	\includegraphics[scale=0.4]{regime_nn.pdf}
	\caption{(a) Schematic regime diagram of a transient network as function of frequency $\omega$ and prestress $\sigma$. The extremely low-frequency regime is dominated by reptation, which is not considered in our model. The rest of the diagram consists of a transient and permanent regimes, separated by $\omega_c(\sigma)$, which is the characteristic frequency, see Eq.~(\ref{e37}).  $\sigma_c$ and $\sigma_p$ are the characteristic stresses for the nonlinear-stiffening in the transient and permanent regimes, respectively. (b) Schematic curves of \SC{the differential elastic modulus $K'$} versus $\omega$, for small prestress (linear transient and linear permanent), intermediate prestress (non-linear transient and linear permanent) and large prestress (nonlinear permanent).  The corresponding power-law regimes are indicated. With stress, the low-frequency regime is not predicted to be a strict power law, although we show that it may appear to be so. }
	\label{Fig.rgm}
\end{figure}
\section{Model}
\label{section:3}
To predict the mechanical response of a biopolymer network under stress, we begin by describing a single semiflexible polymer under tension. For simplicity we first discuss a polymer moving on a 2D plain where the transverse deformation is limited in one direction, and then extend the result to polymers in 3D. The Hamiltonian for such polymers with length $\ell$ under tension $F$ is~\cite{Gittes1998,Morse1998,Broedersz2014}:

\begin{equation}
\begin{aligned}
H=\frac \kappa 2 \int dx \,\left(\frac {\partial^2u}{\partial x^2}\right)^2 + \frac F 2 \int \diff x \left(\frac{\partial u}{\partial x}\right)^2 .
\end{aligned}
\label{e1}
\end{equation}
Here $\kappa$ is bending rigidity and $u(x)$ is the transverse displacement at position $x$. The first term in Eq.~(\ref{e1}) is the bending energy of the chain, while the second term is the work done by the external force, where $\Delta \ell = \int \diff x (\partial u/\partial x)^2/2$ is the contraction of the end-to-end distance due to bending fluctuations. This transverse displacement can be decomposed using Fourier series to a series of bending modes, $\{u_q\}$, with $q$ being the wave number:
\begin{equation}
\begin{aligned}
u(x)=\sqrt{\frac{ 2}{ \ell}}\sum_q u_q {\sin}(qx) \qquad \left(q=\frac{n\pi}\ell\right) ,
\end{aligned}
\label{e2}
\end{equation}
where the dimension of $u(x)$ is set to be $[L]^{1/ 2}$ to simplify the expression of the Hamiltonian. Using Eq.~(\ref{e2}), the Hamiltonian in Eq.~(\ref{e1}) is diagonalized and assumes a simple quadratic form: 
\begin{equation}
\begin{aligned}
H= \frac{1}{2}\sum_q\left(\kappa q^4+F q^2\right)u_q^2.
\end{aligned}
\label{e3}
\end{equation}
The dynamics of the amplitudes of all bending modes follow a standard Langevin equation with $\{u_q\}$ as variables (this model is usually referred to as Model A) \cite{Hohenberg1977435}, 

\begin{align}
&\frac{du_q}{dt}=-\sum_p M_{qp}(\{u_s\})\frac{\partial H}{\partial u_p}+\eta_q\notag
\\&=-\sum_p M_{qp}(\{u_s\}) \left(\kappa p^4+F p^2\right)u_p+\eta_q,\label{e4}
\end{align}
where $M_{qp}(\{u_s\})$ is the generalized mobility matrix, $\{u_s\}$ denotes the set of amplitudes of all bending modes, \{$u_{\pi/\ell}$, $u_{2\pi/\ell}$, ...\}, and $\eta_q$ is a thermal Gaussian white noise with zero mean and variance  $\langle \eta_p(t)\eta_q(t')\rangle=2k_BT M_{pq}\delta(t-t')$, where $k_B$ is  Boltzmann constant and $T$ is the temperature.

In general, $M_{pq}$ is a function of all $\{u_s\}$, and Eq.~(\ref{e4}) is nonlinear. However, the transverse displacement of semiflexible polymers is small ($u(x) \ll \ell$), and as we shall see below, when applying an external stress these displacements are even smaller. Therefore, within our framework, Eq.~(\ref{e4}) should be linearized, leading to a constant mobility matrix $M_{pq}$. The evolution of the correlation function of mode $q$, $\langle u_q(0)u_q(t)\rangle$, is then


\begin{align}
&\frac{d}{dt}\langle u_q(0)u_q(t)\rangle=\left\langle u_q(0)\frac{du_q(t)}{dt}\right\rangle\notag\\&=-\sum_pM_{pq} \left(\kappa p^4+F p^2\right)\langle u_q(0)u_p(t)\rangle\notag\\&=-M_{qq} \left(\kappa q^4+F q^2\right)\langle u_q(0)u_q(t)\rangle,\label{e5}
\end{align}
where $\langle ... \rangle$ denotes average over thermal noise realizations. In the last equation we use the fact that $H(\{u_q\})$ is diagonal, implying that there are no correlations between different bending modes, i.e. $\langle u_p(0)u_q(t) \rangle \sim \delta_{pq}$. Thus, the correlation function follows a simple exponential decay, where the variance at $t=0$ is obtained from Eq.~(\ref{e3}) using the equipartition theorem (assuming the system is in equilibrium at $t=0$)~\cite{Gittes1998}, such that :
\begin{equation}
\begin{aligned}
\langle u_q(0)u_q(t)\rangle=\frac{k_B T}{\kappa q^4+F q^2}e^{-M_{qq}({\kappa q^4+F q^2})t}.
\end{aligned}
\label{e6}
\end{equation}
For every bending mode $q$, there is only one parameter associated with its relaxation process, $M_{qq}$, which is the mobility for mode $q$. This mode-dependent mobility naturally emerges  from the linearization of the Langevin equation and can thus be generally applied to any  semiflexible polymer networks. Once the mode-dependent mobility is known, the dynamics of the network is determined. 

Notice that Eq.~(\ref{e6}) shows that the variance of $u_q$ decreases with increasing $F$, thus tension reduces bending fluctuations. This is consistent with our assumption that linearizing Eq.~(\ref{e4}) is always valid for semiflexible polymers. For flexible polymers, although their Hamiltonian can also be diagonalized to a quadratic form (Rouse Model~\cite{RouseJr.19531272}), the amplitude of the transverse fluctuations is large and the linearization of Eq.~(\ref{e4}) is not generally valid. 

Next, we use the correlation function of Eq.~(\ref{e6}) to calculate the correlation function of the end-to-end distance under tension $F$, $C_F(t)\equiv\langle \delta \ell(0)\delta \ell(t)\rangle$. Here $\delta \ell = \langle \Delta \ell \rangle - \Delta \ell$ is the projected end-to-end extension of the polymer with respect to its rest length  for small $u$ and $\Delta \ell = \int \diff x (\partial u/\partial x)^2/2$. This relation leads to a simple formula for $C_F(t)$~\cite{Gittes1998}:

\begin{align}
&C_F(t)= \sum_q q^4\langle u_q(0)u_q(t)\rangle^2\notag\\
&=\sum_q \frac{(k_B T)^2 \,q^4}{\left(\kappa q^4+F \,q^2\right)^2} \exp{\left[-2M_{qq} {\left(\kappa q^4+ F q^2\right)}t\right]}. \label{e8}
\end{align}
The Fourier-Transform of the correlation function is then used to compute the end-to-end power spectrum, $P_F(\omega)\equiv |\delta \ell(\omega)|^2$,

\begin{equation}
\begin{aligned}
P_F(\omega)=\sum_q \frac{(k_B T)^2 \,q^4}{\left(\kappa q^4+F \,q^2\right)} \frac{4M_{qq} }{\omega^2+4M_{qq}^2 {\left(\kappa q^4+ F q^2\right)}^2}. 
\end{aligned}
\label{e27}
\end{equation}
Using the fluctuation-dissipation theorem, we can relate this power spectrum to the response function $\chi(\omega;F)$: $\ell \chi''(\omega;F)=\omega P_F(\omega)/2k_B T$, which together with the  Kramers-Kronig relations gives
\begin{equation}
\begin{aligned}
\!\!\!\!\chi(\omega;F) &\!= \! \sum_q\!  \frac{2k_B TM_{qq}\left({\kappa q^4+F q^2}\right)  }  {\ell\left(\kappa q^2+F\right)^2\left[2M_{qq}{\left(\kappa q^4+F q^2\right)} -i\omega\right] }.\!
\end{aligned}
\label{e18}
\end{equation}
This response function describes the mechanical response of the chain, given that the tension is slightly perturbed around $F$. For polymers in 3D, one should add a factor of two to the right-hand-side of Eq.~(\ref{e18}) in order to account for the transverse displacement in two directions. In the rest of the paper we will use the 3D result. As we will show later, the macroscopic modulus of the entire network can be derived using the response function of a single polymer, Eq.~(\ref{e18}). 

\section{Mode-dependent Mobility }
\label{section:4}
\subsection{Mode-dependent Mobility for Transient-crosslinked Networks}

In the previous section we have shown that the dynamics of a single polymer in any biopolymer network is well described by a mode-dependent mobility $M_{qq}$, while its specific form depends on  the network structure. In this section we first derive $M_{qq}$ for transient networks in the small-$q$ limit in a simple and intuitive way, and then detail a general method that can be used to derive $M_{qq}$ for any $q$.  

Lets us consider a transient netwrok in the hydrodynamic limit (i.e. long-wavelength limit). For small $q$, we can Taylor expand $M_{qq}$,
\begin{equation}
\begin{aligned}
M_{qq} = a_0+a_2 q^2 +a_4 q^4...,
\end{aligned}
\label{e24}
\end{equation}
where $a_n$ denotes the coefficient of the $n$-degree term. The form of $M_{qq}$ is constrained by polymer symmetries. Therefore, since the polymer does not have a preferred transverse direction (it can have a preferred longitudinal direction though due to the polarity~\cite{Alberts}) $M_{qq}$ must be an even function of $q$. The term $a_0$ stands for the mobility for the $q=0$ mode (infinite wavelength). For transient networks, the value of $a_0$ must be $0$. The reason for this can be seen from Fig.~\ref{Fig.1}(b), where we sketch the relaxation process of a long-wavelength mode. In order to relax such a mode, multiple successive unbinding events are required, indicating that transient crosslinkers impose stronger limitations to long-wavelength modes, leading to smaller mobilities. In order to relax the infinite-wavelength mode there should be infinite successive binding events, each of them takes finite time, therefore leading to zero mobility. Hence, the leading term in Eq.~(\ref{e24}) is the quadratic term, and for small $q$ we have $M_{qq}=a_2 q^2$. 

So far we have shown that the transient nature of the network results in a quadratic dependence of $M_{qq}$. However, this quadratic dependence should only be valid for $q\ll q_c$, where $q_c = \pi/\ell_c$ is the characteristic wave number with wavelength $\ell_c$. 
Bending modes with $q\gg q_c$ are not limited by transient crosslinking, and their relaxation is determined by the substrate mode-independent viscosity (rather than the networks itself), which corresponds to a constant mobility $M_0$. Together, we can approximate $M_{qq}$ for all wavelengths:
\begin{equation}
\begin{aligned}
M_{qq}=
\left\{ 
\begin{aligned}
&{a_2 q^2} \qquad&(q \leq q_c )\\
&M_0  \qquad&(q >q_c)
\end{aligned}
\right.. 
\end{aligned}
\label{e16}
\end{equation}
This mobility is discontinuous at $q=q_c$, because we simply separate the bending modes into a crosslink-limited and viscous-dominated parts. In fact,
for bending modes with $q\sim q_c$, both the transient crosslinkers and the substrate viscosity contribute to the stress relaxation, and we anticipate a transition in $M_{qq} $ from the quadratic dependence to the constant mobility. We determine the dependence of $M_{qq}$ for $q\sim q_c$  in Sec.~\ref{section:4}.~B.

Although this form of the mobility is sufficient for predicting the macroscopic modulus (see Sec.~\ref{section:5}), a microscopic understanding of the parameter $a_2$ is important for understanding the physics of the model we use. In general, the mobility (and thus also $a_2$) can depend on $F$, for example due to catch/slip bond. For simplicity we neglect this effect hereafter (see Sec.~\ref{section:6} for further discussion). We can then consider the dynamics of a polymer for $F=0$ and calculate $a_2$. The relaxation time of bending mode $q$, $\tau_q$, can be read from Eq.~(\ref{e6}). When $F=0$, it is $\tau_q = 1/(\kappa q^4 M_{qq})$. For small $q$, we have $\tau_q = 1/(a_2 \,\kappa q^6)$. On the other hand, the relaxation process for small $q$ is limited by binding/unbinding of crosslinkers, which is characterized by a single timescale $\tau_{\rm off}$ (the unbinding time, as appropriate for strong crosslinkers that spend most of their time in the bound state~\cite{Broedersz2010}), therefore, we  also have $\tau_q \sim \tau_{\rm off}$. From dimensional analysis, the coefficient $a_2$ can be written in terms of microscopic parameters:
\begin{equation}
\begin{aligned}
a_2 = c\, \frac{\ell_c^6}{\kappa \tau_{\rm off}},
\end{aligned}
\label{e31}
\end{equation}
where $c$ is a dimensionless constant. Here $\ell_c$ is the average spacing between crosslinkers which appears as it is the only lengthscale associated with transient crosslinkers. 

In order to find the value of $c$, we use the (mean-field) result of Ref.~{\cite{Broedersz2010}} for the response function $\chi$ for $F=0$ (see Appendix A for a complete derivation). In this case, the relaxation rate of each bending mode, $\omega_r(q)$, can be written as \cite{Gittes1998}
\begin{equation}
\begin{aligned}
\omega_r(q) = \frac{2\kappa q^4}{\xi(q)}.
\end{aligned}
\label{e9}
\end{equation}
where $\xi(q)=1/M_{qq}$ is the mode-dependent friction. Since a mode with longer wavelength must have longer relaxation time, $\omega_r(q)$ must increase monotonically, suggesting the existence of the inverse function $q(\omega_r)$. The slowest relaxation rate is $\omega_r(q=\pi/\ell)$ ($\ell$ being the polymer length), corresponding to the longest wavelength mode, which must vanish as $\ell\to\infty$. Therefore, in the long-chain limit we always have $\omega \gg \omega_r(q=\pi/\ell)$. This allows us to approximate the summation in Eq.~(\ref{e18}) with an integral (as will be done in the rest of the paper), leading to:

\begin{eqnarray}
\chi(\omega;F=0)&\simeq& \frac{2k_B T}{\pi\kappa}\int_0^{\infty} \diff q \, \frac{2a_2{ q^2}  }  { \omega_r(q) -i\omega} \nonumber\\&=&\frac{a_2^{1/ 2}k_B T}{3\kappa^{3/ 2	}}\omega^{-1/2}(1+i),\label{e26}
\end{eqnarray}
where the $-1/2$ exponent is consistent with Ref.~\cite{Broedersz2010}. In Ref.~\cite{Broedersz2010}, the linear response function $\chi(\omega;F=0)$ is derived using a mean-field theory (see appendix for details-- note the factor of two for a polymer in 3D):
\begin{align}
\chi(\omega;F=0)&=0.0036\frac{2k_BT\ell_c^3}{\pi \kappa^2}\int_{-\infty}^{\infty} \,\frac{\diff q}{q^2-2i\omega \tau_{\rm off}}\notag\\&=0.0036\frac{k_BT\ell_c^3}{\tau_{\rm off}^{ 1/ 2} \kappa^2}\omega^{- 1/ 2}(1+i). \label{e28}
\end{align}
Comparing this response function with our result (Eq.~(\ref{e26})) gives the analytical expression of the coefficient $a_2$
\begin{equation}
\begin{aligned}
a_2 =  0.00012 \frac{\ell_c^6}{\kappa \tau_{\rm off}}. 
\end{aligned}
\label{e30}
\end{equation}

Although $\ell_c$ and $\kappa$ can be measured experimentally, the unbinding time $\tau_{\rm off}$ of the crosslinker is usually unknown. Therefore, replacing $a_2$ with Eq.~(\ref{e30}) does not reduce the number of fitting parameters. However, Eq.~(\ref{e30}) gives a microscopic understanding of the coefficient $a_2$, and can be further used to calculate $\tau_{\rm off}$ (see Sec.~\ref{section:5}).

\subsection{Mode-dependent Mobility for Generic Networks}
Heretofore we have derived the mode-dependent mobility for transient networks. We find that a quadratic dependence in the small-$q$-limit naturally emerges as a result of the transient nature of the network. However, the mobility for $q\sim q_c$ remains unclear, and the same derivation does not apply to other networks, despite the fact that the mode-dependent mobility can be applied to any semiflexible polymer networks. In this subsection we provide a general method for deriving the mode-dependent mobility for any network, provided that the correlation function of the end-to-end distance is known. The latter can be found analytically, using numerical simulations, or measured experimentally. We then use this method to derive the mode-dependent mobility of transient networks for $q\sim q_c$.

\begin{figure}[t]
	\centering
	\includegraphics[scale=0.5]{fig2.pdf}
	\caption{Mode-dependent mobility, $M_{qq}$, as a function of $q$, calculated from Eq.~(\ref{e34}) with $M_0=100a_2\kappa q_c^4$. The form of the mobility (including the cusp) is further discussed in Appendix B. }
	\label{Fig.mdm}
\end{figure}
The correlation function, $C_F(t)$, is related to the mode-dependent-mobility, $M_{qq}$, via Eq.~(\ref{e8}). Once $C_F(t)$ of a given network is obtained, $M_{qq}$ can be calculated by inverting Eq.~(\ref{e8}). For simplicity, and as the mobility is assumed to be independent of $F$, we only consider $F=0$. Following the same reasoning as in the previous section (see paragraph before Eq.~(\ref{e26})) we approximate the summation in Eq.~(\ref{e8}) with an integral:
\begin{align}
C_0(t)\equiv C_{F=0}(t)&= \frac{2\ell (k_B T)^2}{\pi \kappa^2}\int \frac{\diff q}{q^4}  \,\exp{\left[-\omega_r(q)t\right]}\notag
\\&\!\!\!\!\!\!\!\!=\frac{2\ell (k_B T)^2}{\pi \kappa^2}\int \frac{\diff \omega_r}{q^4}  \,\frac{\diff q}{\diff \omega r}\exp{\left[-\omega_r t\right]}. \label{e10}
\end{align}
Surprisingly, we find that the correlation function at vanishing force is proportional to the Laplace transform of $q^{-4}\diff q/\diff \omega_r$. The inverse Laplace transform of Eq.~(\ref{e10}) gives a differential equation for $q(\omega_r)$,
\begin{equation}
\begin{aligned}
\mathcal{L}^{-1} \{C_0\}(\omega_r)
=\frac{2\ell (k_B T)^2}{\pi \kappa^2q^4}\,\frac{\diff q}{\diff \omega r},
\end{aligned}
\label{e11}
\end{equation}
where 
\begin{equation}
\begin{aligned}
\mathcal{L}^{-1} \{C_0\}(\omega_r)
=\frac{1}{2\pi i}\lim_{T\to \infty}\int_{\gamma -i T}^{\gamma+i T}e^{s\omega_r}C_0(s) \diff s. 
\end{aligned}
\label{e12}
\end{equation}
Here $\gamma $ is any real number greater than the real part of all singularities of $C_0(s)$. Using Eq.~(\ref{e11}) and the correlation function $C_0(t)$, one can derive the relaxation rate $\omega_r(q)$. Then the mode-dependent friction is easily found using Eq.~(\ref{e9}). 


We now use this method to derive the mode-dependent mobility of transient networks for all $q$'s. As discussed above, the stress relaxation of transient networks is governed by two different mechanisms: a slow relaxation relying on the transient crosslinkers and a fast relaxation dominated by the substrate viscosity. We assume the timescales for the two processes to be separated, i.e. $\tau_{\rm off}\gg \tau_{\rm per}$ or $a_2 \ll2 M_0 q_c^4$, where $\tau_{\rm per} = 1/(2\kappa M_0 q_c^4)$ is the longest relaxation time governed by the substrate viscosity.  Under this assumption, the contribution to the correlation function from the two relaxation processes is additive \cite{Broedersz2010}:
\begin{align}
C_0(t)&=&\frac{2(k_BT)^2\ell}{\pi \kappa^2}\int_{0}^{\infty}\frac{\diff q}{q^4} \,\exp{\left(-2a_2 {\kappa q^6}t\right)}\notag\\&\,\,+&\frac{2(k_BT)^2\ell}{\pi \kappa^2}\int_{q_c}^{\infty}\frac{\diff q}{q^4} \,\exp{\left(-2M_0 {\kappa q^4}t\right)}, \label{e13}
\end{align}
where the first term is the contribution due to the transient nature of the crosslinkers. The second term is the classic result of the correlation function of a semiflexible polymer with mobility $M_0$ where the average crosslinking distance is $\ell_c$~\cite{Gittes1998}. 
The inverse Laplace transform of Eq.~(\ref{e13}) is:
\begin{align}
& \mathcal{L}^{-1}\{C_0\}(\omega_r)=\frac{2(k_BT)^2\ell}{\pi \kappa^2}\int_{0}^{\infty}\frac{\diff q}{q^4} \,\delta{\left(\omega_r-2a_2 {\kappa q^6}\right)}\notag\\ &+\frac{2(k_BT)^2\ell}{\pi \kappa^2}\int_{q_c}^{\infty}\frac{\diff q}{q^4} \,\delta{\left(\omega_r-2M_0 {\kappa q^4}\right)}\notag
\\& =\frac{(2a_2)^{ 1/ 2}(k_BT)^2\ell}{3\pi \kappa^{ 3/ 2}}\omega_r^{-3/ 2}\notag\\&+\frac{(2M_0)^{ 3/ 4}(k_BT)^2\ell}{2\pi \kappa^{ 5/ 4}}\omega_r^{-7/ 4}\Theta\left(\omega_r-2M_0 \kappa q_c^4\right), \label{e14}
\end{align}
where $\Theta(x)$ is the heaviside function.
Substituting Eq.~(\ref{e14}) into Eq.~(\ref{e11}) and integrating from $q$ to infinity, with the boundary condition $\omega_r(q\to \infty)=\infty$ gives:
\begin{align}
&\!\!\!\!\!(2a_2\kappa)^{1/2}\omega_r^{-1/2}+{(2M_0\kappa)^{ 3/ 4}}\omega_r^{-3/ 4}\Theta\left(\omega_r-2M_0 \kappa q_c^4\right)\notag\\&\!\!\!\!\!=q^{-3}-q_c^{-3}\Theta\left(-\omega_r+2M_0 \kappa q_c^4\right).\label{e15}
\end{align}

Because we assume $a_2\ll(2M_0 \kappa q_c^4)$, the solution of  Eq.~({\ref{e15}) is
\begin{equation}
\begin{aligned}
\omega_r(q)=
\left\{ 
\begin{aligned}
&{2 a_2\kappa}q^{6}{\left[1-\left(\frac{q}{q_c}\right)^{3}\right]^{-2}} \!\!&(q \leq q_c-\Delta q )\\
&2M_0\kappa q^4  \!\!&(q >q_c-\Delta q ),
\end{aligned}
\right.
\end{aligned}
\label{e32}
\end{equation}
where $\Delta q = {{a_2}^{ 1/ 2}q_c^2}/({3M_0^{1/ 2}})$. 
The mobility $M_{qq}$ can then be found using Eq.~(\ref{e9}),

\begin{align}
M_{qq}=
\left\{ 
\begin{aligned}
&{a_2 q^2}{\left[1-\left(\frac{q}{q_c}\right)^{3}\right]^{-2}} &(q \leq q_c-\Delta q )\\
&M_0  &(q >q_c-\Delta q).
\end{aligned}
\label{e34}
\right.
\end{align}
In Fig.~\ref{Fig.mdm} we plot the mobility $M_{qq}$ as a function of $q$. For $q\ll q_c$, $M_{qq} $ shows a quadratic dependence on $q$, in agreement with our previous analysis for the small $q$ limit. When $q$ approaches $q_c$, $M_{qq}$ increases dramatically until it reaches $M_0$. Although this cusp in $M_{qq}$ near $q=q_c$ appears unphysical, it is essential within our model with a single lengthscale $\ell_c$ for the appearance of the plateau in the linear modulus~\cite{Broedersz2010} (see Appendix B for a mathematical proof). 



\begin{figure}[t]
	\centering
	\includegraphics[scale=0.41]{Fig_2nn.pdf}
	\caption{(a) A plot of $K'$ as function of the external stress $\sigma$ for different frequencies using experimental data of reconstituted actin networks (taken from Ref.~\cite{Mulla2019}). In (b) we plot the rescaled $K'$, vs the rescaled prestress for the various frequencies. The purple curve is the theoretical prediction of Eq.~(\ref{e20}). The best fit values are $\gamma_1 = 1.71 \,{\rm Pa \, s^{1 /3}}$, $\gamma_2 = 35.9\,{\rm Pa \, s^{ 1/ 2}}$. (fit was done using only the $K'$ data.)}
	\label{Fig.2}
\end{figure}

\section{Dynamic Modulus} 
\label{section:5}
In this section we derive the dynamic modulus of the transient network using the mode-dependent mobility calculated above. We focus on the transient nature of the network and neglect the detailed network structure by assuming all the polymers are aligned in the same direction, i.e.,~an effective  1D network. A more complete discussion of the 3D case is deferred to future publication~\cite{unpublished}. When an external prestress $\sigma$ is applied on an aligned network, all the polymers feel the same tension, $F=\sigma/\rho$, where $\rho$ is the polymer length per unit volume~\cite{Broedersz2014}. When  $\sigma$ is perturbed by $\diff\sigma$, the tension is also perturbed, $\diff F =\diff\sigma/\rho $, leading to an extension, $\diff \ell = \ell \chi \diff F$. Therefore, the perturbation in stress causes a perturbation in strain, $\diff \gamma = \diff \ell/\ell= \diff \sigma\chi /\rho$, and the differential modulus is: 
\begin{equation}
\begin{aligned}
K(\omega;\sigma)\equiv \frac{\diff \sigma}{\diff \gamma}=\frac{\rho}{\chi(\omega;F=\frac{\sigma}{\rho})}.
\end{aligned}
\label{e17}
\end{equation}




We start by calculating $K$ in the low-frequency regime, which is governed by the transient crosslinkers, \\$\omega\ll\omega_c(\sigma)$, where $\omega_c(\sigma)$ is the characteristic frequency separating the transient regime and the permanent regime (see Fig.~\ref{Fig.rgm}). The mathematical definition of $\omega_c(\sigma)$ will be given later in this section. In this regime the contribution from long-wavelength modes dominates, allowing us to only use the long-wavelength part of $M_{qq}$, i.e.,~$M_{qq}=a_2q^{2}$. Substituting this into Eq.~(\ref{e18}), we have: 
\begin{align}
\chi(\omega;F) &=\frac{2k_B T(2a_2)^{ 1 /2}}{\pi\kappa^{ 3/ 2}\omega^{ 1/ 2}} \notag\\&\times \bigintssss_0^{\infty} \frac{r^4 \diff r} {(r^2+\frac{F}{F_c(\omega)})\left[r^6+\frac{F}{F_c(\omega)}r^4-i\right] },\label{e19}
\end{align}
where $r = (2\kappa a_2/\omega)^{ 1/ 6}q$ and $F_c(\omega)= (\kappa^2\omega/2a_2)^{1/3}$ is a frequency-dependent characteristic tension. Using the expression in Eq.~(\ref{e19}), the differential modulus is calculated using Eq.~(\ref{e17}). We separate the dynamic modulus into the storage modulus, $K'$, and the loss modulus $K''$, using $K=K'+ i K''$. The dynamic modulus can be expressed in a simple form:


\begin{align}
&K'(\omega;\sigma) = K_c(\omega){\rm Re}\left[g\Big({\sigma}/{\sigma_c(\omega)}\Big)\right]\notag\\
&K''(\omega;\sigma) = K_c(\omega){\rm Im}\left[g\Big({\sigma}/{\sigma_c(\omega)}\Big)\right],\label{e20}
\end{align}

where 

\begin{equation}
\begin{aligned}
\left[g(x)\right]^{-1}= {\bigintssss \dfrac{(q^6+ {x}q^4)\diff q}{(q^2+ {x})^2[(q^6+xq^4)-i]} },
\end{aligned}
\label{e21}
\end{equation}
is a function describing the stress-stiffening behavior of semiflexible polymers, and  
\begin{equation}
\begin{aligned}
\sigma_{c}(\omega)= \gamma_1 \omega^{\frac 1 3}\quad;\quad
K_c(\omega) = \gamma_2 \omega^{\frac 1 2},
\label{e22}
\end{aligned}
\end{equation}
are the characteristic stress and modulus, with $\gamma_1=\rho(\kappa^2/2a_2)^{1/3}$ and $\gamma_2 =\pi \rho\kappa^{ 3/ 2}/(2^{ 3/ 2}a_2^{ 1/ 2}k_BT)$. These characteristic stress and  modulus agree with what has been observed in reconstituted actin networks crosslinked by $\alpha$-actinin-4, and was explained using a phenomenological theory~\cite{Mulla2019}. By rescaling the experimental data using $K_c$ and $\sigma_c$, we collapse the  data of the storage modulus onto a single curve, which is fitted with our theoretical prediction of Eq.~(\ref{e20}), see  Fig.~\ref{Fig.2}. Here we fit the data using only two fitting parameters, $\gamma_1$ and $\gamma_2$, for all curves~(i.e. after fitting one line we predict all the rest)~\footnote{The collapse of the experimental data was previously done in Ref.~{\cite{Mulla2019}} using a phenomenological theory.}. 

We continue with calculating the differential modulus in the entire frequency regime. To do so, we substitute  Eq.~(\ref{e34}) into Eq.~(\ref{e18}) and Eq.~(\ref{e17}), leading to the differential modulus, which is the central result of this paper:
\begin{equation}
\begin{aligned}
&K(\omega;\sigma)=[\tilde\chi_t(\omega;\sigma)+\tilde\chi_p(\omega;\sigma)]^{-1},
\end{aligned}
\label{e37}
\end{equation}
where $\tilde \chi_t$ and $\tilde \chi_p$ are contributions from the transient and permanent modes, respectively. The two terms  are
\begin{align}
\tilde\chi_t(\omega,\sigma)&= D_t\int_0^{1-\Delta r}\diff r \Bigg[\frac{r^4[1-r^3]^{-2}}{\omega_{t}r^2+\omega_{t\sigma}}\notag\\&\times\frac{ 1 }  {r^4[1-r^3]^{-2}({\omega_t r^2+\omega_{t\sigma} }) 
	-i\omega}\Bigg]\notag\\
\tilde \chi_p(\omega,\sigma)&=D_p\int_{1-\Delta r}^{\infty}\diff r \Bigg[ \frac{r^2 }  {\omega_{p} r^2+\omega_{p\sigma}}\notag\\&\times\frac{1 }  {\omega_{p} r^2+\omega_{p\sigma} -i\omega}\Bigg],\label{e36}
\end{align}
where $r = q/q_c$, $\Delta r = \Delta q /q_c$, $D_t =0.11\gamma_2^{-1}\tau_{\rm off}^{-3/2}$, $\omega_t=0.23/\tau_{\rm off}$, $\omega_{t\sigma}=0.38\gamma_1^{-1}\tau_{\rm off}^{-2/3}\sigma$, $D_p=2.09 \gamma_2^{-1}\tau_{\rm off}^{1/2}\tau_{\rm per}^{-2}$, $\omega_p=1/\tau_{\rm per}$ and $\omega_{p\sigma}=1.63 \gamma_1^{-1}\tau_{\rm per}^{-1}\tau_{\rm off}^{1/3}\sigma$.  Here $\Delta q = {{a_2}^{ 1/ 2}q_c^2}/({3M_0^{1/ 2}})$ as defined after Eq.~(\ref{e32}). The nonlinear modulus for any given $\omega$ and $\sigma$ is then determined by four parameters: $\gamma_1$, $\gamma_2$, $\tau_{\rm off}$ and $\tau_{\rm per}$. 

The transient and permanent regimes are determined by the dominant term in  Eq.~(\ref{e37}) (when $\tilde\chi_t$ dominates the network is in the transient regime and when $\tilde\chi_p$ dominates it is in the permanent regime). We then define the characteristic frequency $\omega_c$ to be the one satisfying $|\tilde \chi_t(\omega_c,\sigma)|=|\tilde\chi_p(\omega_c,\sigma)|$. For small stress, numerical analysis suggests $\omega_c(\sigma=0)=1.26\omega_t=0.29/\tau_{\rm off}$, indicating that the unbinding time is the timescale separating the transient regime \SC{ from the permanent regime}. For larger stress $\omega_c$ is decreasing and it vanishes as the stress exceeds a threshold (see Fig.~\ref{Fig.rgm}). This suggests that for very large stress the transient behavior vanishes, and the permanent-regime plateau expands and covers the entire low-frequency regime. However, such a shift is only obvious when $\sigma>\sigma_p$, where $\sigma_p$ is the characteristic stress for the permanent regime, and $\omega_c$ can be regarded as a constant in the linear regime (see Fig.~\ref{Fig.rgm}). 

In the transient regime, i.e. $\omega\ll\omega_c(\sigma)$, $\tilde\chi_t$  dominates Eq.~(\ref{e37}), and the dynamic moduli are well described by Eq.~(\ref{e20}). For small prestress we have $K\sim \omega^{1/2}$, same as the linear modulus. When the prestress increases the term $r^4(\omega_t r^2+\omega_{t\sigma})$ is comparable to $\omega$, and the network  stiffens nonlinearly. The characteristic stress for this stiffening is $\sigma_c$, see Eq.~(\ref{e22}). Since $\sigma_c\sim \omega^{1/3}$, the characteristic stress is larger for higher frequency, and the frequency dependence of $K$ is also weakened by the prestress. In fact, the apparent scaling exponent of $K$ can be written as a function of $\sigma/\sigma_c$, see Sec.~\ref{section:6} for further discussion. 

When $\omega\gg\omega_c(\sigma)$, $\tilde\chi_p$  dominates  Eq.~(\ref{e36})  and the network is in its permanent regime. In this case the network behaves as the well-studied permanent network~\cite{Gittes1998,Morse1998}. There is a noteworthy frequency in Eq.~(\ref{e37}) associated with the permanent regime, $\omega_p+\omega_{p\sigma}$, which is the relaxation rate for the bending modes with $q$ slightly larger than $q_c$. As we assume $\tau_{\rm off}\gg \tau_{\rm per}$, we have $\omega_c(\sigma)\ll\omega_p+\omega_{p\sigma}$, and we expect a plateau in $K'$ between the two frequencies (see Ref.~\cite{Gittes1998} and Appendix B). In fact, when $\omega_c(\sigma)\ll\omega\ll\omega_p+\omega_{p\sigma}$ in Eq.~(\ref{e37}), we have $K\approx C_1 + (C_2/\omega)i$, where $C_1$ and $C_2$ do not depend on $\omega$.  The different scaling for $K'$ and $K''$ in this plateau regime results in different behavior in the transition between the transient and permanent regimes, see discussion on Fig.~\ref{Fig.3} below. For $\omega\gg\omega_p+\omega_{p\sigma}$, we have $K\sim \omega^{3/4}$, in agreement with analytical results for permanent networks~\cite{Gittes1998, Morse1998}. The network in the permanent regime nonlinearly stiffens when the prestress reaches the characteristic stress, $\sigma_p$, which is defined as the prestress satisfying $\omega_p=\omega_{p\sigma}$. 







 

 We then determine the four parameters ($\gamma_1, \gamma_2, \tau_{\rm off}, \tau_{\rm per} $) by fitting the experiment data of Ref.~\cite{Mulla2019}. We find $\gamma_1$ and $\gamma_2$ by fitting the $K'$ data with Eq.~(\ref{e20}) (Fig.~\ref{Fig.2} (b)), then fix these values and fit the $K''$ data to find $\tau_{\rm off}$ and $\tau_{\rm per}$ (Fig.~\ref{Fig.3} (b)), where we minimize the sum of the squared deviation of the theoretical prediction (Eq.~\ref{e37}) from the experimental data for all curves simultaneously (each curve with different $\sigma$). The same values of $\gamma_1$, $\gamma_2$, $\tau_{\rm off}$ and $\tau_{\rm per}$ are used for the entire family of curves. Using the best-fit values of the parameters, we find that $\tau_{\rm off} =0.0085\,{\rm s}$, in agreement with previous result~\cite{Mulla2019}. The value of $\tau_{\rm per}$ may be inaccurate, since the experimental data is only in the transient and intermediate regimes. Changing the value of $\tau_{\rm per}$ only slightly affects the predicted moduli, as long as $\tau_{\rm per}\ll \tau_{\rm off}$. 

As shown in Fig.~\ref{Fig.3}(a), the fitting for the $K'$ data is excellent. Our theory predicts that increasing the prestress results in an increase of $K'$ together with a weakening of its frequency dependence, as is also observed experimentally. For $K''$ the theory does not agree well with the experimental data. Although our theory shows the same qualitative features as the experiments, including the strengthened loss modulus and decreased scaling exponent for $K''(\omega)$, our predicted  $K''$ is always smaller than experimentally observed in the low-frequency regime (see Fig.~\ref{Fig.3}(b)). We believe that one reason for this is the faster increase of $K'$ with prestress (compared to $K''$). Then, for high prestress, the ratio $K''/K'$ can be close to $0.1$, making it experimentally hard to get accurate  measurement of $K''$ (see inset of Fig.~\ref{Fig.3}(b)).  This is a known issue that was discussed in Ref.~\cite{Velankar2007}. The other reason for the disagreement is that the mobility we used for $q\sim q_c$ can be inaccurate: the mobility for $q\gg q_c$ and $q\ll q_c$ is well understood, but the transition between the two parts  was not explored in this work. The prediction for $K''$ is more affected by this transition because the two regimes are well separated for $K''$: $K''$ increases with frequency in the transient regime but decreases $\sim\omega^{-1}$  in the permanent regime.  On the contrary, $K'$ undergoes an insignificant change because after increasing with frequency in the transient regime it reaches a plateau in the permanent regime, and is thus less affected by the transition itself. 
\begin{figure}[t]
	\centering
	\includegraphics[scale=0.93]{Fig3_n.pdf}
	\caption{Differential storage modulus $K'$ (a) and differential loss modulus $K''$ (b) as measured in the experiments of Ref.~\cite{Mulla2019} (x symbols), together with the  theoretical fitted curves of Eq.~(\ref{e37}). The best-fit values of the parameters are: $\gamma_1 = 1.71 \,{\rm Pa \, s^{\ 1/ 3}}$, $\gamma_2 = 35.9\,{\rm Pa \, s^{ 1/ 2}}$, $\tau_{\rm off}=0.0085\, {\rm s}$ and $\tau_{\rm per} = 2.0\times 10^{-14} \,{\rm s\,}$. In the inset of (b) we plot the ratio $K''/K'$ measured experimentally as function of $\sigma$. }
	\label{Fig.3}
\end{figure}






\section{Discussion and Conclusion}
\label{section:6}

We have proposed a general theory for the nonlinear modulus of biopolymer networks, which describes the dynamics of the networks using a mode-dependent mobility, $M_{qq}$. In transient-crosslinked networks, we find that $M_{qq}\sim q^2$ for the long-wavelength modes, indicating that the relaxation of  these modes is slowed down by transient crosslinkers. This explains both the $1/2$ scaling exponent of the frequency-dependent linear modulus reported in Ref.~\cite{Broedersz2010}, and the glassy-like relaxation in the presence of prestress~\cite{Mulla2019}. Our theory suggests that the transient nature of crosslinkers in  biopolymer networks is the cause of the apparent weak exponent of the frequency-dependent modulus found in living cells~\cite{Fabry2001}.

\begin{figure}[t]
	\centering
	\includegraphics[scale=0.6]{fig5.pdf}
	\caption{Apparent scaling exponent $\beta$ as function of stress measured in experiment, together with our theoretical prediction of Eq.~(\ref{e23}). Experimental data is taken from Ref.~\cite{Mulla2019}, with $\omega=0.31 {\rm Hz}$.}
	\label{Fig.5}
\end{figure}

To further analyse how the prestress affects the frequency dependence of the nonlinear modulus, we calculate the apparent scaling exponent $\beta(\omega,\sigma)$ around a given frequency $\omega$ when the network is subjected to a prestress $\sigma$:

\begin{equation}
\begin{aligned}
\beta(\omega;\sigma)=\frac{\partial \ln K'(\omega;\sigma)}{\partial \ln \omega}\,.
\end{aligned}
\label{e23}
\end{equation}
In the transient regime where Eq.~(\ref{e20}) is valid, $\beta$ is only a function of $\sigma/\sigma_c(\omega)$.  In Fig.~\ref{Fig.5} we plot   $\beta$ for varying $\sigma/\sigma_c(\omega)$ and fixed $\omega=0.31\,{\rm Hz}$, and find that the apparent scaling exponent $\beta$ decreases with increasing prestress. Further, $\beta$ can be  tuned to any value between $0$ and $1/2$, depending on the prestress strength. This might explain the weak scaling exponent found in living cells. The cytoskeleton, which gives the cell its rigidity, is formed by many transient crosslinking proteins, and is effectively a prestressed transient biopolymer network, where molecular motors are responsible for the prestress \cite{Schiffhauer20161473}. Moreover, if the crosslinkers binding/unbinding process is out of  equilibrium, an internal stress can be created even in the absence of molecular motors~\cite{Chen2020}. Experiments on living cells have also observed reduced scaling exponent for increasing internal motor stress~\cite{Stamenovic2004,Kollmannsberger20113127}. 


The weak scaling exponent of the frequency-dependent shear modulus in living cells has been discussed repeatedly during the past two decades ~\cite{Fabry2001,Stamenovic2004,Desprat20052224,Bursac2005557,Balland2006}. Most of the previous works try to explain this weak exponent using the soft glassy rheology (SGR) model, which describes the dynamics of soft materials that have structural disorder and metastability ~\cite{Fabry2001,Bursac2005557}. These materials are out of equilibrium, as thermal energy is insufficient to drive the systems across the energy barrier to reach equilibrium~\cite{Navarro1996, Sollich1998,Cloitre2003}. Although living cells are also out of equilibrium,  there is no direct evidence that the origin of the reported weak exponent in living cells is the same as that of soft glassy materials. \SC{Also, the scaling exponent in the SGR model does not change with prestress, which does not agree with most experiments on living cells as well as the more recent experiments on reconstituted networks ~\cite{Stamenovic2004,Kollmannsberger20113127,Mulla2019}. }


In this paper, we describe the dynamics of transient-crosslinked biopolymer networks and provide the microscopic understanding of the nonlinear transient regime. Although our work treats a system close to thermal equilibrium, we find that it exhibits the same weak-scaling phenomena. Here the weak exponent is a result of the coupling between multiple relaxation times, which comes from the relaxation of infinite bending modes that are slowed down by the transient crosslinkers. Our theory suggests that a distribution of long relaxation times can exist in an equilibrium system with only short obvious timescales, where the long relaxation times come from the collective behavior of multiple crosslinks.  Unlike glassy systems with metastability, our system does have a single longest relaxation time. This relaxation time is not the system intrinsic time-scale $\tau_{\rm off}$, which is related to the binding/unbinding process, rather it is the relaxation time for the longest-wavelength mode, $q = \pi/\ell$ (Eq.~\ref{e6}): $\tau_r \sim \tau_{\rm off}\ell^6/\ell_c^6$. When the polymers have sufficiently long contour length $\ell$, $\tau_r$ can become very long such that the system is in a nearly glassy state. 


In this work we assumed that the mode-dependent mobility and $a_2$ for transient networks are independent of the external stress. This simple assumption is sufficient to explain the experimental data of Ref.~\cite{Mulla2019}, although in general, $a_2$ may be a function of the polymer tension. For example, the lifetimes of the bound states may depend on the mechanical force exerted on the crosslinkers, which can lead to a change in $a_2$. For most protein complexes this unbinding time decreases under external force, as they are effectively dragged off. Such an effect is called slip bond. The opposite behavior, catch bond, also exists in some crosslinking proteins~\cite{Schiffhauer20161473}, including $\alpha$-actinin-4 that is used in Ref.~\cite{Mulla2019}. To account for this change in the unbinding time, we can replace the constant $a_2$ with $a_2(\sigma)$, where the rest of the derivation remains unaffected. For slip bond, we expect $a_2$ to increase with $\sigma$ and the exponent to be smaller than in Fig.~\ref{Fig.5}, while for catch bond we expect the opposite. 






Our results are focused on transient networks, but the mode-dependent mobility theory we present is applicable for  other biopolymer networks, such as entangled networks or networks with multiple types of crosslinkers.  Different networks correspond to different mode-dependent mobility, $M_{qq}$. Our theory suggests a one-to-one relation between $M_{qq}$ and the correlation function of a single filament. Since  the linear modulus of the network is usually related to the single-filament correlation function~\cite{Gittes1998}, our theory also implies an underlying relation between the linear and the nonlinear modulus in biopolymer networks. 



\begin{acknowledgements}
	\emph{Acknowledgements}:
	This work 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). The authors
	acknowledge helpful discussions with G.~Koenderink and Y.~Mulla.
\end{acknowledgements} 






