\section{Introduction}

We may not be able to play soccer like Lionel Messi, but perhaps we can train deep networks to hallucinate imagery suggesting that we can.  Consider the images in Fig.~\ref{fig-intro}.  In this paper we describe research toward synthesizing realistic images or sequences that forecast the appearance of people performing desired actions.  The model can take a pose and a reference image of a person, and generate a novel view of the person in the given pose, while capturing fine appearance details.

Image or video generation has emerged as an important problem in many domains. Within robotics, work has explored predicting the consequences after interactions between an agent and its environment~\cite{finn2016unsupervisedPhysicalInteraction}. In natural language processing, approaches~\cite{mansimov2015generating,reed2016generative} have been proposed to tackle tasks such as text to image or image to text synthesis. Generative models of video sequences~\cite{walker2017poseKnow, Villegas2017LearningToGenerate} and fashion images~\cite{ZhaoWCLF17,MaDPIG17} are a core part of visual understanding and have received renewed attention from the vision community.

In this paper, we focus on generating realistic images of a particular person striking a desired pose. Human body pose is a natural intermediate representation for this generation, and hence utilized in many previous methods for synthesizing human motion and video~\cite{BrandH00,EfrosBMM03,walker2017poseKnow}. We follow in this paradigm, using body poses to generate person images and videos. Simple networks~\cite{walker2017poseKnow} may generate blurry and distorted images.  Stylistic methods~\cite{isola2016pix2pix} have shown great success in generating realistic images, but lack control over the appearance of the generated images. Our task requires the model to be able to generate images of a person with a specific appearance. Inspired by~\cite{chen2017stylebank}, we propose a novel appearance rendering network which encodes appearance into convolutional filters. These filters are operationalized using a fully convolutional network, and utilized in an image-to-image translation structure that transfers the desired appearance to the generated image.

To sum up, we contribute a new state of the art generative model that performs adaptive appearance rendering to create accurate depictions of human figures in these human poses. We demonstrate our model on two applications: fashion image synthesis and multi-person video synthesis involving complex human activities. Comprehensive quantitative results are provided to facilitate the analysis of each component of our model, and qualitative visualizations are shown in complementary to quantitative results.

\section{Related Work}

{\bf Image generation:} Image-to-image translation has achieved great success since the emergence of GANs~\cite{goodfellow2014GAN}. Recent work produces promising results using GAN-based models~\cite{isola2016pix2pix, CycleGAN2017}. Stylized images can be generated by using feed-forward networks~\cite{gatys2015neuralStyleTransfer} with the help of perceptual loss~\cite{johnson2016perceptualLoss}. The recent work of \cite{chen2017stylebank} proposes a structure to disentangle style and content for style transfer. Styles are encoded using a stylebank (set of convolution filters). Visual analogy making~\cite{NIPS2015_5845_analogy, sadeghi2015visalogy} generates or searches for an new image analagous to an input one, based on other previously given example pairs.

{\bf Generation from intermediate representations:} Generation directly in low-level pixel space is difficult and these types of approaches tend to generate blurry or distorted future frames. To tackle this problem, hierarchical models~\cite{walker2017poseKnow, Villegas2017LearningToGenerate, ZhuPrada17, KaracanAEE16, ma2017pose} adopt intermediate representations. This type of approach can alleviate image blur, however the quality of generation largely depends on the the image generation network. Simple generation networks can still produce blurry images as shown in~\cite{walker2017poseKnow}. It's worth noting that Ma et al. propose a two-stage approach to solve a problem similar to our work~\cite{ma2017pose}. However, the main contribution of this work is to propose a simple yet well-performing base model to the problem of generating images conditioned on pose and appearance.

{\bf Video generation:} Data-driven video generation has seen a renaissance in recent years. One major branch of methods uses RNN-based models such as encoder-decoder LSTMs for direct pixel-level video prediction~\cite{ranzato2014videoGenerationBaseline, srivastava2015unsupervisedVideoRepresentationLSTMs, oh2015actionConditionedVideoPrediction, mathieu2015deepBeyondMSE}. These methods successfully synthesized low-resolution videos with relatively simple semantics, such as moving MNIST digits or human action videos with very regular, smooth motion. Subsequent work has attempted to expand the quality of predicted video in terms of resolution and diversity in human activity.  Earlier efforts were focused on optical flow-timescale prediction, further work pushed past into more complex motions (e.g.~\cite{liu2017video}).

In summary, our approach builds on the substantial body of related work in pose analysis and style/analogy-based image generation. We contribute a novel, simple and effective method for adaptive appearance rendering model for image/video generation from human poses.

\begin{figure*}[!tph]
\center
\includegraphics[width=0.85\textwidth]{figures_for_main_paper/sys-flow.pdf}
\caption{Overview of the adaptive appearance rendering network. Given input posemap image and appearance reference image, the appearance is encoded into convolutional filters, then these filters are placed in the image-to-image translation network to transfer the appearance to images of each person striking the desired pose.}
\label{fig:overview}
\end{figure*}
\section{Adaptive Rendering Network}
\label{sec:adaptive_render}
Given pose estimations as intermediate representations, the goal of our model is to synthesize realistic image of person with desired appearance and pose, where the pose of every person is a \textit{posemap} image in which ``white" body joint points are drawn on a black background canvas. Great success has been shown in generating realistic images given a sketch as inputs~\cite{isola2016pix2pix}, which is similar to our task since the posemap image can be treated as a sketch. However, for our task we cannot generate random appearance and we need control over appearance of generated images to make sure that the generated person is wearing the desired clothing.

To accomplish this goal, we propose an adaptive rendering structure where the appearance filters are adaptively computed from an input reference image using a fully convolutional neural network (FCN) and supervised using style loss to enforce similar color distributions as the input reference image. By incorporating this FCN into an image-to-image translation network a realistic image of a target consistent with the desired action and appearance can be generated.

\subsection{Network Structure}

Fig.~\ref{fig:overview} shows our adaptive rendering network (Ada-R Network) architecture, which consists of two branches: an encoder-decoder branch, and an adaptive rendering branch. The network requires two input images: a posemap image, and a reference image which provides the appearance of the same person in a different pose. The goal of the network is to generate a realistic image of a person consistent with the posemap and having appearance consistent with the reference image.

\textbf{Encoder-Decoder}: Instead of training an encoder-decoder network which can reconstruct input images, our encoder-decoder branch shown in Fig.~\ref{fig:overview} is a sketch $\rightarrow$ image model.

We use the same input size and encoder-decoder structure as in~\cite{isola2016pix2pix}: both generator and discriminator use modules of the form convolution-BatchNorm-Relu~\cite{ioffe2015batch}, the encoder consists of convolutional layers with stride 2 and symmetrically the decoder consists of convolutional layers with fractional stride $\frac{1}{2}$.

\textbf{Adaptive Rendering}: The encoder-decoder network takes binary posemap images as inputs which do not contain any information about the appearance of the person.  Hence, we propose to use another network to learn appearance information. By combining these two networks together we are able to generate realistic images of a person wearing the desired clothing. Here we introduce our Ada-R network.

To transfer the desired appearance to the encoder-decoder branch, we replace the last convolutional filter in the encoder-decoder branch with our adaptive appearance transfer filter. The adaptive appearance filter $\bm{K}_{ada \mhyphen app}$ encoding appearance information of a person is derived from an input appearance reference image $\bm{I}_{app}$ using a fully-convolutional network
\begin{equation}
    \bm{K}_{ada \mhyphen app} = FCN(\bm{I}_{app})
\end{equation}
For image generation, the inputs are appearance-posemap image pairs. For video generation, the rendering of one person's posemap sequence only requires one reference image, and the frames could be obtained by performing adaptive appearance rendering frame by frame. The filter is applied to the rendering procedure by
\begin{eqnarray}
    \bm{F} &=& \mathcal{E} (\bm{I}_{pose}) \\
    \bar{\bm{F}} &=& \bm{F} \ast \bm{K}_{ada \mhyphen app}
\end{eqnarray}
where $\mathcal{E}$ is the encoder network, $\bm{I}_{pose}$ is the posemap image and $\ast$ is a regular 2D convolution operator. $\bm{F}$ is the feature map of size $w\times h\times C_{in}$ generated by the encoder network. We set the size of $\bm{K}_{ada \mhyphen app}$ to $k\times k\times C_{in}\times C_{out}$ to be compatible with $\bm{F}$, where $k$, $C_{in}$ and $C_{out}$ denote the kernel size, input channel number and output channel number of this adaptive convolution operation. $\bar{\bm{F}}$ is the feature map after applying the adaptive appearance filter to the feature map $\bm{F}$. The person $\bm{I}_{gen}$ with desired appearance is finally produced with
\begin{equation}
     \bm{I}_{gen} = \mathcal{D} (\bar{\bm{F}})
\end{equation}
where $\mathcal{D}$ is the decoder network.

We release a reference implementation\footnote{\url{https://github.com/wisdomdeng/AdaptiveRendering}} to provide the full details on the model architecture and parameter settings.

\subsection{Loss Function}

Our network is trained in an adversarial setting, where the Ada-R network is the generator $G$, and a discriminator $D$ is introduced to discriminate between real and generated images. The overall architecture is a conditional generative adversarial network conditioned on the reference image. Let $\bm{I}_{goal}$ be the target image that we try to produce, and $\bm{I}_{gen}$ be the image that the Ada-R network generated. The loss of our Ada-R network is defined as
\begin{equation}
\mathcal{L_{CGAN}}(G, D) + \mathcal{L_{T}}
\end{equation}
where $\mathcal{L_{CGAN}}(G, D)$ is the standard adversarial loss.
$\mathcal{L_{T}}$ is our appearance transfer loss, defined below.
\begin{equation}
\mathcal{L_{T}} = \alpha \mathcal{L}_{1}(\bm{I}_{gen}, \bm{I}_{goal}) + \beta \mathcal{L_{C}}(\bm{I}_{gen}, \bm{I}_{goal}) + \gamma \mathcal{L_{S}}(\bm{I}_{gen}, \bm{I}_{app}).
\end{equation}
$\mathcal{L}_{1}$ is an $L_1$ loss that encourages pixel-level agreement between the generated image and the target image, defined as:
\begin{equation}
\mathcal{L}_{1}(\bm{I}_{gen}, \bm{I}_{goal}) = ||\bm{I}_{gen}-\bm{I}_{goal}||.
\end{equation}
$\mathcal{L_{C}}$ and $\mathcal{L_{S}}$ are the content and style loss, defined the same as Gatys et al.~\cite{gatys2015neuralStyleTransfer}.  They
aim to preserve image structure and colour distributions respectively:
\begin{eqnarray}
\mathcal{L_{C}}(\bm{I}_{gen}, \bm{I}_{goal}) = \sum_{l \in {l_{c}}}||F_{l}(\bm{I}_{gen})-F_{l}(\bm{I}_{goal})||^2 \\
\mathcal{L_{S}}(\bm{I}_{gen}, \bm{I}_{app}) = \sum_{l \in {l_{s}}}||G_{l}(\bm{I}_{gen})-G_{l}(\bm{I}_{app})||^2
\end{eqnarray}
where $F_{l}$ is the feature map from layer $l$ of a pretrained VGG-19 network~\cite{Simonyan15VGG}. $l_{c}$ are layers of VGG-19 used to compute the content loss. $G_{l}(\cdot)$ is the Gram matrix which learns the correlations of color distribution given two input images. $l_{s}$ are layers of VGG-19 used to compute the style loss.

The final objective is defined as
\begin{equation}
G^{\star} = \arg \underset{G}{\min}\ \underset{D}{\max}\ \mathcal{L_{CGAN}}(G, D) + \mathcal{L_{T}}
\end{equation}

\section{Experiments}
\label{sec:exp}
We demonstrate our model on the DeepFashion dataset~\cite{liuLQWTcvpr16DeepFashion} and Volleyball dataset ~\cite{Ibrahim_2016_CVPR_volleyball}. For the DeepFashion dataset, the goal is to render a given appearance to different poses of same person. Due to the completeness of poses, we perform comprehensive quantitative evaluations and an ablation study on the DeepFashion dataset. We also demonstrate a novel application of our model on the Volleyball dataset where the goal is to synthesize short videos (5 frames) containing groups of people given the people in the $1^{st}$ frame as appearance reference image.

To train our Ada-R network, we compute content loss at layer \textit{relu4-2} and style loss at layer \textit{relu1-2}, \textit{relu2-2}, \textit{relu3-2}, \textit{relu4-2} and \textit{relu5-2} of the pre-trained VGG-19 network. We set the learning rate to \textit{1e-3}, loss weights are set to bring the $L1$, content and style losses to a similar scale. For fashion dataset $\alpha=100$, $\beta=1e-4$, $\sigma=1e-14$. For Volleyball dataset, $\alpha=100$, $\beta=0.1$, $\sigma=1e-12$. To make the training stable, for DeepFashion dataset, in each iteration the generator is updated three times and the discriminator is updated one time; and for Volleyball dataset, in each iteration the generator is updated twice and the discriminator is updated one time.

We adopt Mean Square Error (MSE), Peak Signal-to-noise Ratio (PSNR), and Structural Similarity (SSIM)~\cite{wang2004image} as evaluation metrics. MSE and PSNR are pixel-wise measurements for quality of reconstructed or generated images. SSIM measures the quality of a generated image by considering the image's structural information instead of only pixel-wise errors. For both datasets, the OpenPose detector~\cite{cao2017realtime} is used to obtain corresponding poses for each target in a given image.

\subsection{Experiments on DeepFashion}
\label{eval:fashion}
DeepFashion contains fashion images of different person IDs, and most IDs contain 4 views: front, side, back and full body. We filtered out IDs without upper body view or with less than 3 views, resulting in 3418 person IDs. To train our model, we use 2395 person IDs, resulting in 14370 pose-appearance pairs. To test our model, we use 1023 person IDs, resulting in 6138 pose-appearance pairs. Original images are resized to $128 \times 128$ for both training and testing.

To analyze the strength of our model as well as the importance of each component in both the network structure and loss, we compare our Ada-R approach with the following approaches:

{\bf Ada-R w/o style loss}: In this baseline, we remove the style loss from the loss function.

{\bf Ada-R w/o content loss}: In this baseline, we remove the content loss from the loss function.

{\bf Ada-R w/o $L_1$ loss}: In this baseline, we remove the $L_1$ loss from the loss function.

{\bf Pose-GAN$^\star$ (PG$^\star$)}: We adopt the generation structure used by Walker et al.~\cite{walker2017poseKnow}, namely the posemap image and appearance image are concatenated as input to the image-to-image translation network and the FCN is removed. The same loss is used as our Ada-R approach.

{\bf Pose-GAN$^\star$ (PG$^\star$) w/o style loss}: We adopt the generation structure of PG and remove the style loss from the loss function.

{\bf Visual Analogy Making$^\star$ (VAM$^\star$)}: We adopt the analogy based generation structure used in~\cite{Villegas2017LearningToGenerate,NIPS2015_5845_analogy}. Same loss is used as our Ada-R approach.

For evaluation, posemaps and images are normalized to $[-1,1]$. To measure the quality of the poses of the generated images, we propose a new \textit{perceptual pose score} which uses a state-of-the-art pose estimator~\cite{cao2017realtime} to extract pose from generated images and compares each pose joint with the corresponding ground truth pose joint using Euclidean distance. Quantitative results showing comparisons among our approach with all baselines are shown in Tab.~\ref{tbl:fashion_measure} and visualizations of all approaches are shown in Fig.~\ref{fig-fashion1} and Fig.~\ref{fig-fashion2}.

\begin{figure}[h]
\center
\includegraphics[width=0.95\textwidth]{figures_for_main_paper/fashion-1-new.pdf}
\caption{Visualization of generation results of Ada-R using different loss terms on DeepFashion.}
\label{fig-fashion1}
\end{figure}

\begin{figure}[h]
\center
\includegraphics[width=0.95\textwidth]{figures_for_main_paper/fashion-2-new.pdf}
\caption{Visualization of generation results of different state-of-the-art approaches on DeepFashion.}
\label{fig-fashion2}
\end{figure}

The quantitative evaluation results show that our \textit{Ada-R} models, except for the ablated one without $L_1$ loss, uniformly outperform \textit{Pose-GAN} and \textit{Visual Analogy Making} models by significant margins across all measures. The ablation study results for \textit{Ada-R} models on different components of the loss function also show that each component contributes positively to the model's performance. Since the target image directly supervises the model training through the content loss and $L_1$ loss, removing one of them, especially the $L_1$ loss, causes the largest degradation of model performance. This degradation in performance is also reflected in visualization results presented in Fig.~\ref{fig-fashion1}. As we can see from Fig.~\ref{fig-fashion1}, removing $L_1$ or content loss results in: (1) the pose of the target in the generated images being inferior (missing arms or hands); (2) images generated being more blurry; (3) the appearance of the generated images showing more errors in details (color transfer to arms, wrong colors, uneven colors in clothes). As seen from Tab.~\ref{tbl:fashion_measure}, removing the style loss from the training objective causes relatively smaller impacts on the model's performance and the model can generate person images striking the correct pose. However it is observed that many examples generated by \textit{Ada-R w/o style loss} do not reflect the colors of the style reference image correctly and the colors tend to be flat and unvaried compared with our full model.

When compared with \textit{PG$^\star$} models that learn the representation for the pose and the style together, \textit{Ada-R} which adopts a two-stream approach and encodes the pose and style separately also shows clear advantages. As shown in Fig.~\ref{fig-fashion2}, \textit{PG$^\star$} tends to generate images that show both the pose of the input posemap and the pose in the style reference image so as to match the color distribution of the style image, resulting in two overlapping targets in the generated images. While \textit{PG$^\star$ w/o style loss} tends to generate blurry images with incorrect poses or appearance as shown by the pose score and examples in Fig.~\ref{fig-fashion2}. The comparison demonstrates the advantage of disentangling the representations of pose and style and targeting each representation with a specific loss. \textit{Visual Analogy Making$^\star$ (VAM$^\star$)}~\cite{NIPS2015_5845_analogy}, also uses disentangled representations for pose and style and we use the same loss function as our model for training. However, we find that it cannot generate images with correct poses most of the time. Failure examples of \textit{Visual Analogy Making$^\star$} are presented in Fig. ~\ref{fig-fashion2}.
\input{fashion_measure.tex}

\subsection{Experiments on Volleyball}
\label{eval:volleyball}
The Volleyball dataset contains sequences of volleyball games. For this dataset, our model is trained to observe players in the $1^{st}$ input frames and predict their future appearances from frame 6 to frame 10. To get training sequences of each player, we run person detection~\cite{ren2015faster} and tracking~\cite{dlib09} to get tracklets of each player in each clip. We follow the data split of original dataset and preprocessing is conducted to filter out instances with less than 10 joints and clips containing less than 10 targets. We get 1262 clips for training and 790 clips for testing. Images of players are cropped and then resized to $256 \times 256$ pixels for both training and testing. For evaluations, we compare our model with state-of-the-art approaches including Pose-GAN and VAM as described in Sec.~\ref{eval:fashion}. Comparisons among our approach and two state-of-the-art approaches are shown in Tab.~\ref{tbl:volleyball_measure} and visualizations of all approaches are shown in Fig.~\ref{fig:volley_visual}.

\input{volleybal_measure.tex}

\begin{figure*}[!tph]
\center
\includegraphics[width=0.9\textwidth]{figures_for_main_paper/volleyball-1-new.pdf}
\caption{Visualizations on Volleyball dataset.}
\label{fig:volley_visual}
\end{figure*}

As shown in both quantitative and qualitative results, our model performs very well on this application. Our model can generate clean and sharp person image sequences with correct appearance generated and consistent poses as in the posemap input, despite the complex scenario in terms of delicate pose changes, motion blur, and appearance variation.

We could also generate frames in original resolution $1280 \times 720$ by replacing players in a frame by generated ones. Examples of the frame synthesis results are shown in Fig.~\ref{fig:volley_whole_frame}, where the generated players are highlighted using blue bounding boxes.

\begin{figure*}[!tph]
\center
\includegraphics[width=0.9\textwidth]{figures_for_main_paper/volleyball-whole-frame-new.pdf}
\caption{Visualizations on whole frame generation on Volleyball dataset. The generated players are highlighted using \textcolor{blue}{blue} bounding boxes in each frame.}
\label{fig:volley_whole_frame}
\end{figure*}
\section{Conclusion}
We proposed a novel approach for synthesizing realistic images or sequences that generate the appearance of people taking the desired poses. Our approach encodes the appearance into convolutional filters. These filters are learned using a fully convolutional network, and utilized in an image-to-image translation structure that transfers the desired appearance to the desired pose. Both quantitative and qualitative results show that our model is superior to state-of-the-art approaches and can generate better images involving complex appearance and better videos involving complex human activities. The success of our model demonstrates that the combination of appearance filters and style loss can render the desired target appearance while adapting to the predicted pose.


\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}