Theory of Semiflexible Filaments and Networks

Abstract

We briefly review the recent developments in the theory of individual semiflexible filaments, and of a crosslinked network of such filaments, both permanent and transient. Starting from the free energy of an individual semiflexible chain, models on its force-extension relation and other mechanical properties such as Euler buckling are discussed. For a permanently crosslinked network of filaments, theories on how the network responds to deformation are provided, with a focus on continuum approaches. Characteristic features of filament networks, such as nonlinear stress-strain relation, negative normal stress, tensegrity, and marginal stability are discussed. In the new area of transient filament network, where the crosslinks can be dynamically broken and re-formed, we show some recent attempts for understanding the dynamics of the crosslinks, and the related rheological properties, such as stress relaxation, yield stress and plasticity.

1. Introduction

What is today routinely called “semiflexible filaments”, 30 years ago was more frequently referred to as polymer chains with high bending modulus. The classical theory of polymers starts from the ideal chain models (such as beads and springes, or freely-jointed rods) and the Edwards Hamiltonian [1] that accounts for the connectivity of monomers:

$$H=\frac{3k_{B}T}{2b^2}\sum _{j=1}^N{({\mathbf{r}}_j-{\mathbf{r}}_{j-1})}^2=\frac{3k_{B}T}{2b}{\int }_0^{L}ds{\left|\frac{d\mathbf{r}\left(s\right)}{ds}\right|}^2,$$

(1)

where b is the Kuhn length (essentially the monomer size) and N the number of monomers, each positioned at $\mathbf{r}\left(s\right)$ along the chain. The contour length of this chain is $L=Nb$. This Gaussian Hamiltonian gives the average end to end distance of the chain $\langle {[\mathbf{r}\left(N\right)-\mathbf{r}\left(0\right)]}^2\rangle =b^{2}N$ in full analogy with free diffusion laws. It is also the basis for the dynamic theories such as Rouse model and its versions. However, very soon people realized that adding the bending rigidity to this ideal chain would not only lead to many interesting theoretical results and analogies, but also would better describe many actual polymer chains and filaments. We review this in detail in the next section.

As with ordinary, fully flexible polymer chains, the main avenue of further research went towards the dense systems, when the chains begin to interact with each other. The interactions range from the plain excluded volume for two monomers approaching contact, starting with a pair-interaction term $\sim \delta [\mathbf{r}\left(s\right)-\mathbf{r}\left(s^{\prime }\right)]$ added to the Edwards Hamiltonian, to various long-range interactions leading to additional physics and phase transformations. In the case of chains with bending rigidity, at higher density one often encounters the nematic liquid crystal phases [2,3,4,5], where chains fold into tight hairpins with long straight segments aligned [6,7,8]. For very stiff chains, when backfolding (hairpins) is not allowed, the analogy with vortex lines in superconductors has led to another successful stream of theoretical research using the transverse chain excursions away from the mean centerline ${\mathbf{r}}_{\perp }\left(s\right)$ as the field variable [9,10]. The issues of semiflexible chains confined in limited and shaped volume were recently reviewed [11].

Here we will not touch these very interesting topics. Our aim is to examine the properties of an isolated filament (subject to various external forces and extensions), and then proceed to the crosslinked networks of such filaments, in all cases remaining in the regime of low total volume fraction. This is closely corresponding to many biological scaffolding systems, which underpins the current active interest in this class of materials.

2. A Single Semiflexible Filament

Semiflexible filaments can be categorized into natural, mainly biological ones such as DNA, actin, collagen, microtubules, intermediate filaments, etc., and the synthetic filaments such as carbon nanotubes or stiff polymer chains. The main physical feature of this class of polymer chains is that there is a finite bending rigidity, such that their mechanical properties are determined by the competition between the elastic potential energy of bending and the entropic free energy arising from thermal fluctuations of their conformations. There is a possibility to have elastic energy of extension added to a model, but our first aim is to examine filaments that are inextensible.

The wormlike chain (WLC) model, proposed by Kratky and Porod in 1949 [12], is widely used for describing the elasticity of a single semiflexible filament, since the work of Fixman and Kovac in 1973 [13]. As sketched in Figure 1, a filament with its contour length as L is coordinated by $\mathbf{r}\left(s\right)$, where $0\le s\le L$ is the arc-length coordinate along the chain. By assuming a finite bending modulus of the filament, its Hamiltonian can be expressed as

$$\begin{array}{ccc}\hfill H\left[\mathbf{r}\right(s\left)\right]& =& \frac{\kappa }{2}{\int }_0^{L}ds{\left(\frac{d^2\mathbf{r}}{d{s}^2}\right)}^2,\hfill \end{array}$$

(2)

where $d^2\mathbf{r}/d{s}^2$ denotes for local curvature of a curve, and κ is the bending modulus. In classical elasticity, the bending modulus of a filament is given by the Young modulus of its material and the shape of its cross-section: $\kappa =YI$, where the second moment of filament cross-section $I=\int y^{2}dxdy$; for an elastic tube of outer radius a and inner radius b (the wall thickness $a-b$) the bending modulus will be $\kappa =(\pi /4)Y(a^4-b^4)$. Persistence length, defined as $l_p=\kappa /k_{B}T$, is a key quantity for describing whether a chain is regarded as flexible ($L\gg l_p$), semiflexible ($L\sim l_p$), or rigid ($L\ll l_p$). As the persistence length is a function of the temperature, a chain can become more flexible under high temperature, and vice versa. As noticed from how the conformation is parameterized in WLC model by the coordinate $\mathbf{r}\left(s\right)$, the semiflexible chain is naturally assumed as locally inextensible, i.e., ${(d\mathbf{r}/ds)}^2=1$, for all s.

As mentioned at the beginning, except for the bending energy defined above, the other important energy source of a semiflexible chain is thermal fluctuation of its conformations, or simply the conformational entropy. This can be more clearly understood by formulating the free energy of a semiflexible chain, $F\left(\mathbf{R}\right)$, as a function of its end-to-end separation, $\mathbf{R}=\mathbf{r}\left(L\right)-\mathbf{r}\left(0\right)$, which is useful in analyzing its mechanic properties, such as the force-extension relation. Alternatively, it can be re-expressed as a function of the non-dimensional ’end-to-end factor’ $\mathbf{x}=\mathbf{R}/L$. The free energy of a semiflexible chain is

$$\begin{array}{c}\hfill F=-k_{B}Tlnp\left(\mathbf{x}\right),\end{array}$$

(3)

where the probability $p\left(\mathbf{x}\right)$ is the partition function at a constrained end-to-end distance, expressed by the standard ’Edwards path integral’:

$$\begin{array}{c}\hfill p\left(\mathbf{x}\right)\propto \int \mathcal{D}\mathbf{r}\left(s\right)e^{-H\left[r\right]/k_{B}T}\delta \left[{\left(\frac{d\mathbf{r}}{ds}\right)}^2-1\right],\mathrm{with}\mathbf{x}=\frac{1}{L}{\int }_0^{L}ds\left(\frac{d\mathbf{r}}{ds}\right),\end{array}$$

(4)

where the local inextensibility is incorporated as a delta-function constraint in the path integral. However, it is difficult to obtain an analytic free energy expression of an inextensible chain defined above: although the bending Hamiltonian is quadratic in the field $\mathbf{r}\left(s\right)$, the inextensibility constraint makes this classical problem highly non-linear. One of the ways for resolving this difficulty is by using the concept of mean extensibility (often called the ’mean field approximation’ in this context), where the total length of the chain L is conserved, while local fluctuations are allowed. Although mathematically this is an approximation – physically this limit is perhaps a more reasonable model of a semiflexible chain or filament. The idea of mean extensibility is simple: the Fourier transformation of δ-function in Equation (4) is

$$\begin{array}{ccc}\hfill \delta \left[{\left(\frac{d\mathbf{r}}{ds}\right)}^2-1\right]& =& \int \mathcal{D}\varphi \left(s\right)e^{i\varphi \left(s\right)\left[{\left(d\mathbf{r}/ds\right)}^2-1\right]},\hfill \end{array}$$

(5)

where $\varphi \left(s\right)$ is introduced as an auxiliary field, imposing the original constraint at each point s along the filament. If we replace $\varphi \left(s\right)$ with a constant value ϕ (the ’mean field’), then the chain can be extensible locally, but with the average inextensibility will be maintained: $\langle {(d\mathbf{r}/ds)}^2\rangle =1$. The validity of this approximation can be checked by calculating the fluctuations of the local length. As discussed in Ref. [14], if a chain is divided into N segments with equal length as $L/N$, then the fluctuation in the local arc-length of an individual segment, Δ, can be expressed as ${\langle {\Delta }^2\rangle }^{1/2}=L{\left[(N-1)(N+1)\right]}^{1/2}/N$, which is equal to $L/N$ for $N\to \infty$; the fluctuations are well-bounded. With this approximation, the probability function (4) can be reduced to a compact integral in infinite limits:

$$\begin{array}{c}\hfill p\left(\gamma \right)=\int d\varphi e^{i\gamma \varphi }{\left(\frac{\sqrt{i\varphi }}{sin\sqrt{i\varphi }}\right)}^{3/2},\end{array}$$

(6)

where the single non-dimensional parameter $\gamma =\kappa (1-x^2)/2k_{B}TL=(l_p/2L)(1-x^2)$ measures the deviation of the chain from its full extension, and at the same time – encodes the effective stiffness of the filament [14]. It turns out that a simple analytical expression

$$p\left(\gamma \right)=exp\left[-{\pi }^2\gamma -\frac{1}{\pi \gamma }\right]$$

(7)

captures the full range of Equation (6). A numerical calculation of the probability function is shown in Figure 2a. This numerical integration is quite non-trivial; it uses the fact that the imaginary part of the integrand in Equation (6) is antisymmetric in ϕ, while the (symmetric) real part has several points on the ϕ-axis where an artefact of discontinuous sign-reversal occurs in numerical evaluation. So the integral has to be evaluated piece-wise, correcting for these artefacts. The analytical interpolation Equation (7) is also plotted in Figure 2a, shown to follow the numerical result with high precision. The close match across the whole range of γ suggests that this is, in fact, an exact analytical result of integration (which we simply didn’t know how to calculate). This allows an immediate evaluation of the free energy of semiflexible filament $F_{\mathrm{c}}\left(x\right)$:

$$F_{\mathrm{c}}\left(x\right)=k_{B}T{\pi }^{2}c(1-x^2)+\frac{k_{B}T}{\pi c(1-x^2)},$$

(8)

where the non-dimensional stiffness parameter $c=l_p/2L$ is the single factor that determines the properties of an inextensible filament. Note that this expression has no approximations (apart from the ’mean-field’ approximation for an auxiliary field ϕ), and so it is equally valid in the limit of flexible chains ($c\ll 1,x\ll 1$), giving the classical entropic-spring result, in the limit of highly stretched flexible chain ($c\ll 1,x\to 1$), and also in the limit of rigid athermal elastic rod ($c\gg 1,x\to 1$).

2.1. Models of the Force-Extension Relation

Over the years, several famous analytical forms of the force-extension relation of a semiflexible chain as a function of its end-to-end distance have been developed (comparison shown in Figure 3). We discuss these force-extension expressions here, while the corresponding free energy in each case could be obtained by an integral of the force over the extension. In particular, when we compare with the numerical result of Equation (6) in Figure 2a, the constrained free energy expressions are obtained by integration, and then raised to exponent to get a non-normalized probability plotted.

Marko-Siggia model [10]. In 1995, Marko and Siggia proposed an approximate force-expansion relation as:

$$\begin{array}{c}\hfill f\left(x\right)=\frac{k_{B}T}{l_p}\left(x+\frac{1}{4{(1-x)}^2}-\frac{1}{4}\right).\end{array}$$

(9)

It originates from the high-extension limit of the WLC model, which was interpolated such that it also matches the low-force limit of entropic spring. This expression works very well for chains near full extension, and has a long and successful history of application in many experimental situations [15]. For instance, most values of the persistence length of various biological filaments in the literature are obtained by fitting force-extension curves with the Marko-Siggia expression (9). When we examine this expression, it turns out that it is valid for a long, flexible chain in the highly stretched limit; in other words, at $\gamma \to 0$ in Equations (7) and (8).

Ha-Thirumalai model [16]. In 1997, Ha and Thirumalai proposed another approximate force-expansion relation by numerically solving:

$$\begin{array}{ccc}\hfill f\left(x\right)& =& 2\lambda x,\mathrm{where}\lambda \left(f\right)\mathrm{satisfies}:\frac{f^2}{4{\lambda }^2}+\frac{3}{4}\sqrt{\frac{2k_{B}T}{l_p\lambda }}coth\left(\frac{\sqrt{\lambda }L}{\sqrt{2k_{B}T{l}_p}}\right)=1.\hfill \end{array}$$

(10)

This expression, as well as the Marko-Siggia model, and also the Odijk model [17], all predict the limiting behavior of $f\to (k_{B}T/l_p){[1-x]}^{-2}$ at $x\to 1$, with a slightly different numerical factor each. This is an inherent feature of the inextensible WLC model in the limit of high extension when the transverse excursions of the chain are small (certainly no ’overhangs’ are possible in the chain configuration). However, slightly different values of $l_p$ would result by equally successfully fitting the force-extension curves of chains and filaments by these expressions. In terminology of Equations (7) and (8), these limits correspond to the regime $\gamma \ll 1$, that is: a flexible chain at high extension. The validity of Ha-Thirumalai model is expanded to a wider range of γ compared with that by Marko and Siggia, as shown in Figure 2a.

MacKintosh-Käs-Janmey model [18]. Another way of evaluating the equilibrium elastic response of filaments to tensile force f, starting from the same WLC model, is by applying equipartition to each bending mode in Fourier space, which is nicely reviewed in [19] by Broedersz and Mackintosh, and in [20] by Palmer and Boyce. Again, using the mean-field constraint on chain extensibility, but a different way of performing the summation of Fourier modes, the end-to-end factor is expressed as a function of applied force:

$$\begin{array}{c}\hfill x=1-\frac{L}{l_p{\pi }^2}\sum _n\frac{1}{n^2+{\xi }^2}=1-\frac{L}{l_p{\pi }^2}\left(\frac{\pi \xi coth\left(\pi \xi \right)-1}{2{\xi }^2}\right),\end{array}$$

(11)

where the shorthand $\xi =(L/\pi )\sqrt{f/k_{B}T{l}_p}$. The force-extension relation can be obtained by numerically inverting the relation between x and f (Palmer and Boyce [20] offer a rather accurate analytical interpolation using Padè approximation, valid in the tensile deformation regime). In the limit of high extensions, the divergent tensile force follows the same characteristic limit as all other models: $f\to (k_{B}T/4l_p){[1-x]}^{-2}$. A very important factor, missing in earlier models, is the natural length of the filament: the force-free average end-to-end extension we labelled as $x_0$ in Figure 2b. Palmer and Boyce give it as $x_0=1-L/6l_p$ for this model (e.g., at $L=l_p$ we have $x_0\approx 0.83$). In the linear regime (at small deformations around $x_0$) the filament responds with the spring constant 90 $k_{B}T{l}_p^2/L^4$.

Blundell-Terentjev model [14]. Equations (7) and (8) constitute the model of Blundell and Terentjev, which has an advantage of being fully analytical and capturing the right physics in both the low and the high temperature limits, as well as in tension and compression/buckling regimes up to the highly bent elastica limit [21]. The force-extension relation is written as

$$\begin{array}{c}\hfill f\left(x\right)=-\frac{{\pi }^2{k}_{B}T{l}_p}{L^2}x+\frac{4k_{B}T}{\pi l_p}\frac{x}{{(1-x^2)}^2},\end{array}$$

(12)

where an internal elastic energy of a bent filament is obtained in the limit $T\to 0$ for a fixed bending modulus $\kappa =k_{B}T{l}_p$. As always for theories based on WLC, the finite extension limit gives the divergent force scaling: $f\to (k_{B}T/\pi l_p){[1-x]}^{-2}$. For a very short or stiff filament ($L\ll l_p$), the internal potential dominates over the entropic effects, and the filament can be treated as a rigid athermal rod. The zero-force average natural extension is $x_0^2=1-2L/{\pi }^{3/2}{l}_p$ (e.g., at $L=l_p$ we have $x_0\approx 0.80$). As shown in Figure 2b, if the ratio $c=l_p/2L$ is below the critical value of ${\pi }^{-3/2}$, the minimum of the free energy will be located at $x=0$, i.e., the filament can be regarded as a flexible chain and reproduce a Gaussian limit upon first-order expansion (with the entropic spring constant $4k_{B}T/\pi l_{p}L$).

2.2. Mechanical Response of a Semiflexible Filament

The characteristic response of a chain or filament to a tensile force has been widely explored in the literature and remains the basis of determining the persistence length $l_p$, usually by fitting the M-S model. For the filaments that are stiff enough to have a non-zero average natural extension $x_0$, the response to a compression force also becomes relevant. When a rigid elastic rod is compressed, it will buckle if the compression force exceeds what is known as the Euler buckling threshold force $f_c={\pi }^2\kappa /L^2$ [22]. Similar response also appears in the case of a semiflexible filament, although in this case the thermal fluctuations result in a finite range of buckling force, and the associated hysteresis. The free energy $F_c\left(x\right)$ given in Equation (8) has the minimum located at $x_0$. The value of the free energy at $x=0$ decreases with the increase of the compression force until it becomes a metastable state; in other words, the Euler buckling of a semiflexible filament is a typical first order phase transition [19,21]. The equilibrium threshold force for buckling is obtained as

$$\begin{array}{c}\hfill f_c^{*}=\frac{{\pi }^2{k}_{B}T{l}_p}{L^2}\left[1-\alpha {\left(\frac{L}{l_p}\right)}^{\nu }\right],\end{array}$$

(13)

where the bending modulus κ is replaced by the persistence length $k_{B}T{l}_p$, and $\alpha =1.11$ and $\nu =0.56$ are fitting (interpolation) parameters [21]. The threshold force reduces to the classical Euler form for an athermal elastic rod with $l_p\to \infty$.

So far we only discussed the physics of semiflexible filaments with a fixed contour length L, which was the basis of the WLC statistical model. However, especially at large tension forces, the real filament may be elastically stretched by extending the bonds holding its monomers in a chain. In actin filaments, these bonds are the physical interactions (hydrophobic and hydrogen bonds) holding the G-actin units together, with a characteristic energy of ∼25–30 $k_{B}T$; a similar bonding energy is found to hold the β-sheet peptides in an amyloid filament [23,24]. A stretched DNA double helix would elongate by modifying hydrogen-bond interactions between base pairs. In all of these cases, a finite stretching modulus was indeed detected by experiment.

When a semiflexible filament is highly stretched, its end-to-end factor will be close to unity, and an approximation can be implemented in Equation (8): $1-x^2\simeq 2(1-x)$, making the force-extension relation at $x>x_0$:

$$\begin{array}{ccc}\hfill f\left(x\right)& \approx & -\frac{{\pi }^2{k}_{B}T{l}_p}{L^2}+\frac{k_{B}T}{\pi l_p{(1-x)}^2}\equiv -f_c+\frac{k_{B}T}{\pi l_p{(1-x)}^2}.\hfill \end{array}$$

(14)

An effective spring constant in charge of the elasticity of a filament can be defined by $K_e=df/Ldx$, producing

$$\begin{array}{c}\hfill K_e=\frac{2k_{B}T}{\pi L{l}_p{(1-x)}^3}\equiv k_e\sqrt{{(1+f/f_c)}^3},\end{array}$$

(15)

where the second version of the expression is obtained implementing the Equation (14), with the shothand notation for the constant $k_e=2{\pi }^{7/2}{k}_{B}T{l}_p^2/L^4$, and the Euler critical force $f_c$. If in addition, the material making up the filament has its own intrinsic Young’s modulus Y, an effective mechanical spring constant has to be added: $k_m\sim Y{a}^2/L$, with a as radius of the filament. For actin filament of length $10\mathsf{\mu }$m, this spring constant is about 5 pN/nm, while a microtubule of the same length has a spring constant around 15 pN/nm. Considering the possibility of filament extension, the full length now becomes $L_{\mathrm{eff}}=L+f/k_m$ under a large tension f, and force-extension relation (with the approximation that x is close to 1) can be expressed as

$$\begin{array}{ccc}\hfill x& =& \left(1+\frac{f}{k_{m}L}\right)\left[1-\sqrt{\frac{k_{B}T}{\pi l_p(f+f_c)}}\right].\hfill \end{array}$$

(16)

In other words, the mechanic response becomes softer than entropic one if the filament is experiencing a very large stretching force, while the entropic mode is softer if a filament is stretched by a smaller force. Similar expressions were obtained by Odijk [17], Lipowsky et al. [25], and Holzapfel and Ogden [26].

There is a crossover region bridging the entropic stretching regime with the mechanical stretching regime, see Figure 4. One can make an analogy of the extensible semiflexible filament to two elastic springs connected in series, giving the effective spring constant $K\left(f\right)$ as:

$$\begin{array}{c}\hfill \frac{1}{K\left(f\right)}=\frac{1}{k_e\sqrt{{(1+f/f_c)}^3}}+\frac{1}{k_m}.\end{array}$$

(17)

If the tensile stretching force is small, i.e., $f<f_c$, the spring constant of the filament is approximately a constant $K_0=1/(k_e^{-1}+k_m^{-1})$, which defines linear entropic response regime. If the tensile force is very large, $f>f_c{(k_m/k_e)}^{2/3}$, the filament extension is dominated by the mechanical (enthalpic) stretching, while the entropic elasticity essentially diverges. If the width of ’nonlinear-entropic’ regime is sufficiently large, the entropic part in Equation (17) scales as $f^{3/2}$; however, if the filament material is more easily stretched, the crossover could be narrow and much less distinct.

3. Network of Crosslinked Filaments

Semiflexible filaments can be crosslinked together to form filament bundles and filament networks as illustrated in Figure 5. In biological context, the crosslinks are frequently various biological motors, for example, myosin and kinesin. Though physical properties of an individual semiflexible filament are relatively well studied, understanding how a semiflexible filament bundle or a network responds to an applied deformation is not a trivial question. For bundles, there is a wormlike bundle (WLB) model, proposed and developed by Heussinger et al. [27,28], where the main idea is that the free energy of a semiflexible bundle not only includes bending and stretch energy of the filaments, but also shear energy due to sliding of neighboring filaments with respect to each other. This model works well with the filament bundles which possess a hierarchical structure, such as microtubules. Filament networks, with a more complex structure than that of an individual chain or a hierarchical bundle, show various emergent phenomena, such as negative normal stress of a sheared network or marginal stability of a filament network with reduced number of crosslinks. Compared with the widely used WLB model for semiflexible bundles, there are different approaches for a filament network, which are constructed in various ways. In this section we will mainly focus on reviewing models on filament networks.

3.1. Models of the Filament Network Elasticity

Storm et al. model [29]. Storm et al. proposed a theoretical model of equilibrium semiflexible network, where all the filaments are assumed to be deformed affinely. Suppose that $\mathbf{E}$ is the deformation tensor, and $\rho \sim 1/{\xi }^2$ is the filament length per unit volume (which is a measure of network density). After deformation, the length density of the filaments per unit volume crossing a plane perpendicular to j axis changes to $\rho E_{jk}{n}_k/\det \mathbf{E}$, where $\mathbf{n}$ is the orientation of a filament end-to-end vector in an initially undeformed network, and $\det \mathbf{E}$ measures the relative volume change of the deformed network. The tension acting on the filament bridging two neighboring crosslinks can be denoted by $f\left(\right|\mathbf{E}·\mathbf{n}|-1)$, where $|\mathbf{E}·\mathbf{n}|-1$ represents for the axial strain of the filament, and the specific form of the tension force can be obtained from the preferred force-extension relationship discussed in previous section. The ith component of the tension can be calculated as $f\left(\right|\mathbf{E}·\mathbf{n}|-1)E_{il}{n}_l/|\mathbf{E}·\mathbf{n}|$, and the ij component of the stress tensor can be obtained by following the Doi-Edwards construction of adding all contributions of the filaments weighed by the filament length crossing the j plane [1]:

$$\begin{array}{ccc}\hfill {\sigma }_{ij}& =& \frac{\rho }{\det \mathbf{E}}〈f\left(\right|\mathbf{E}·\mathbf{n}|-1)\frac{E_{il}{n}_l{E}_{jk}{n}_k}{|\mathbf{E}·\mathbf{n}|}〉.\hfill \end{array}$$

(18)

A simple application of the theory is for the simple shear, where the sample experiences a deformation:

$$\begin{array}{c}\hfill \mathbf{E}=\left(\begin{array}{ccc}1& \gamma & 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right),\end{array}$$

(19)

where γ is shear strain. In the case of small shear strain (linear elastic response), the stress defined in Equation (18) will be simplified as [30]:

$$\begin{array}{ccc}\hfill {\sigma }_{ij}& =& \rho 〈f\left(\gamma n_x{n}_y\right)n_i{n}_j〉+O\left({\gamma }^3\right),\hfill \end{array}$$

(20)

where $\det \mathbf{E}=1$ for the simple shear. Mackintosh et al. [18] applied the linear regime of the force-extension relation provided in Equation (11), which makes the $xy$ component of shear stress:

$$\begin{array}{c}\hfill {\sigma }_{xy}=\rho k_{B}T\frac{90l_p^2}{L^3}\gamma 〈n_x{n}_y{n}_x{n}_y〉.\end{array}$$

(21)

Both here and in Equation (18) the averaging is over the statistical ensemble of filaments, having random orientations in the network. The shear modulus of a filament network consisting of isotropically distributed filaments can be obtained as $G_0=6\rho k_{B}T{l}_p^2/L^3$, as a linear function of crosslink density. For a given concentration of monomer in the system $c\approx 1/a{\xi }^2$, the distance to find one chain to be intersected with another, or the distance between two neighboring entanglements, is calculated to be about $L\sim l_p^{1/5}{\left(ac\right)}^{-2/5}$, which gives the scaling of the shear modulus $G_0\sim k_{B}T{l}_p^{7/5}{\left(ac\right)}^{11/5}$ [18]. There is usually a linear regime under small deformation and a nonlinear regime under large deformation in the elasticity of a filament network. By defining the crossover strain and stress from the linear regime to the nonlinear one as ${\gamma }_t=L/6l_p$ and ${\sigma }_t=\rho k_{B}T{l}_p/L^2$, respectively, one can obtain the relationship between the linear shear modulus and the crossover stress as $c^{1/2}{G}_0\sim {\sigma }_t^{3/2}$, recalling that $\rho \sim c$ [31].

Heussinger-Frey model [32,33]. In contrast to perfect crystalline solids or the isotropic homogeneous rubbery networks, in filament networks the non-affine deformations cannot strictly be avoided. That is, network local structure would relax on the scale of a single mesh unit to lower the local energy and achieve global equilibrium. In 2006, Heussinger and Frey proposed a theoretical model for semiflexible filament network where the filaments experience non-affine deformations, which is also named as floppy mode model (sketched in Figure 6).

Assume a 2D filament network with area A, which consists of N filaments with uniformly distributed orientations and positions. The isotropic distribution of angles between two intersecting filaments gives $P\left(\theta \right)=sin\theta /2$, and the random distribution of the mesh size (taken as approximately equal to the contour length for a stiff filament connecting two neighboring crosslinks) can be written as $P\left(\xi \right)=exp[-\xi /\langle \xi \rangle ]/\xi$. Here the average mesh size can be related to the line density as $\langle \xi \rangle =\pi /2\rho =\pi A/2L_{f}N$, where $L_f$ is the full length of a filament, rather than the distance between two neighboring crosslinks. For subfilaments with radius a, which connect two neighboring crosslinks, they have a high stretching stiffness $k_{\parallel }=4\kappa /\xi a^2$, and a relatively softer bending stiffness $k_{\perp }=3\kappa /{\xi }^3$, using the measure of native filament bending modulus $\kappa =k_{B}T{l}_p$. From numerical simulations by Wilhelm and Frey [34], it was found that the effective shear modulus of such a network has two distinct regimes, characterized by a length scale $l_c=L_f{\left(\delta \rho L_f\right)}^{-\nu }$ with $\nu \approx 2.84$. For filaments with radius $a\gg l_c$, the network responds affinely and the elastic response is mainly arises from the stretching deformation, and the shear modulus $G_{aff}\sim k_{\parallel }$ independently of the total filament length; for filaments with radius $a\ll l_c$, the non-affine bending deformations dominates in the mechanic response of the network, and the shear modulus depends on the total fiber length as $G_{na}\sim k_{\perp }\left(\langle \xi \rangle \right){(L_f/\langle \xi \rangle )}^{\mu -3}$ with $\mu \approx 6.67$.

In Heussinger-Frey model, the filament is relatively long compared with the mesh size, $L_f>\xi >a$, i.e., a filament can be shared in several meshes. A subfilament of length ξ connecting two crosslinks typically stores the energy $w_b\left(\xi \right)\simeq \kappa {\delta }_{na}^2/{\xi }^3\simeq \kappa L_f^2/{\xi }^3$. Then the soft bending response is incorporated, in charge of the local relaxation of one subfilament with respect to its neighbors, and the average bending energy stored in a whole filament consisting n subfilaments (line density defined as $\rho =n/L_f$) can be expressed as

$$\begin{array}{c}\hfill \langle W_b\rangle \simeq n{\int }_{l_{min}}^{\infty }d\xi P\left(\xi \right)\frac{\kappa {\delta }_{na}^2}{{\xi }^3},\end{array}$$

(22)

where ${\delta }_{na}$ is the non-affine displacement which is assumed to scale as the total length of the filament (the affine displacement is on the length scale of the mesh size), and $l_{min}$ is introduced as the cutoff as the mesh with small size does not contribute to the bending energy obviously, due to their high bending stiffness $k_{\perp }\sim {\xi }^{-3}$. In other words, the filaments with length smaller than $l_{min}$ will relax from their floppy mode deformation ${\delta }_{na}$, which reduces their bending energy from $w_b\left(l_{min}\right)$ to zero. However, this type of relaxation leads to the movement of the whole filament, which results in a floppy mode of the whole filament. By equalizing $w_b\left(l_{min}\right)$ to the average bending energy $\langle W_b\rangle$, $l_{min}\simeq 1/{\rho }^2{L}_f$ and $\langle W_b\rangle \simeq \kappa /L_f{\left(\rho L_f\right)}^6$ are obtained. Moreover, the shear modulus $G_{na}\simeq \rho /L_f\langle W_b\rangle \sim {\rho }^7$, which closely coincides with the simulation result as ${\rho }^{6.67}$. The crossover component $\nu =3$ can be obtained by equalizing bending and stretching energy $\langle W_b\rangle =\langle W_s\rangle \simeq \kappa L_f/r^2$.

3.2. Continuum Models of Semiflexible Filament Network

For rubber elasticity [35,36], there is a widely used class of theories formulated by constructing an effective continuum composed by periodically repeating blocks. In these models the free energy of the network can usually be explicitly expressed as a function of the invariants of the deformation tensor $\mathbf{E}$: $I_1={\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2$, $I_2={\lambda }_1^2{\lambda }_2^2+{\lambda }_2^2{\lambda }_3^2+{\lambda }_3^2{\lambda }_1^2$ and $I_3={\lambda }_1^2{\lambda }_2^2{\lambda }_3^2$, where ${\lambda }_{1,2,3}$ are the eigenvalues of $\mathbf{E}$ (see Figure 7 for illustration). Such continuum models are often named unit-cell models, where the cells, or the repeating blocks, are imaginary spheres for 1-chain model, cubes for 3- and 8-chain models, or tetrahedra for the 4-chain model. The macroscopic elastic free energy of the network is obtained by multiplying the free energy of an individual cell by the number of cells, whereas any non-affinity is assumed to be confined within a cell. Similarly, based on the energy of an individual filament as discussed above, several continuum models of semiflexible filament network are proposed.

8-chain model. In 2008, Palmer and Boyce [20] applied the 8-chain model for a filament network, where the body-centered cubic cell is constructed with eight filaments connected from the center to all eight lattice points in the corners. The lengths of the edges of the cell without any deformation are equal to ξ, as we define the mesh size, while the connecting chains are shorter by a factor $\sqrt{3}/2$. The important feature of the 8-chain model is that on deformation all eight chains in the cubic cell are stretched in exactly the same way, with the length changing from $\sqrt{3}\xi /2$ to simply $\sqrt{I_1}\xi /2$, since the cell is aligned along the principal axes of deformation tensor $\mathbf{E}$ by construction. Then the free energy density of the network is given directly by summing all equivalent single-filament expressions:

$$\begin{array}{ccc}\hfill F_{8\mathrm{c}}\left(I_1\right)& =& n{F}_{\mathrm{chain}}\left(\sqrt{I_1/3}\xi \right).\hfill \end{array}$$

(23)

where n is the number of crosslinked filament segments per unit volume, and the explicit energy form of an individual filament, $F_{\mathrm{chain}}$, can be chosen from the models introduced in the previous section. Though 8-chain model succeeds in fitting stress-strain relationship in sheared filament networks, as illustrated in next section, it fails to explain negative normal stress due to the intrinsic isotropy in the 8-chain model. The negative normal stress, as well as other important effects such as the marginal stability, are significant physical properties that need to be addressed.

1-chain model. In 2013, Unterberger et al. [37] applied the 1-chain model to describe the equilibrium elasticity of the filament network. In this model, one end of a filament is fixed at the center of an imaginary spherical cell, while the other end is on the sphere surface with an arbitrary, isotropic orientation $(\theta ,\phi )$. When deformed (with the deformation tensor $\mathbf{E}$), the sphere will be deformed into an ellipsoid, where the lengths of the three semi-axes are $\xi {\lambda }_1$, $\xi {\lambda }_2$ and $\xi {\lambda }_3$, respectively. The elastic energy density of the 1-chain network model can be expressed as an orientational average of a deformed filament,

$$\begin{array}{c}\hfill F_1\mathrm{c}=n\int \frac{sin\theta d\theta d\phi }{4\pi }{F}_{\mathrm{chain}}\left[\tilde{\lambda }(\theta ,\phi )\xi \right],\end{array}$$

(24)

where $\tilde{\lambda }(\theta ,\phi )=\sqrt{{sin}^2\theta ({cos}^2\phi {\lambda }_1^2+{sin}^2\phi {\lambda }_2^2)+{cos}^2\theta {\lambda }_3^2}$ is the deformation ratio at orientation $(\theta ,\phi )$. If taking the ’entropic spring’ limit of the Gaussian chain, the average $F_1\mathrm{c}$ reduces to the neo-Hookean model for rubber, expressed as $n{k}_{\mathrm{B}}T(2{\xi }^2/3\pi l_{p}L)[{\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2]$, with the prefactor as the shear modulus, which is the basis of the classical rubber elasticity. However, with the more complicated filament deformation energy, the lack of an explicit form of $F_1\mathrm{c}$ as a function of the strain invariants limits the applications of 1-chain model.

3-chain model. In 2016, Meng and Terentjev [38] applied the 3-chain model, where an imaginary cubic cell is constructed with the lattice points representing the crosslinking sites, and the edges aligned along the principle directions of the deformation tensor $\mathbf{E}$. Three orthogonal chains at each lattice point in one cell are linked with their end-to-end vectors along the edges, and the equilibrium mesh size is ξ. On deformation, the lengths of three perpendicular edges at one lattice point become ${\lambda }_1\xi$, ${\lambda }_2\xi$ and ${\lambda }_3\xi$, respectively. Then the free energy density of such a network can be explicitly expressed as a function of the strain invariants by applying the free energy form of an individual filament defined in Equation (8):

$$\begin{array}{c}\hfill F_{3\mathrm{c}}\left(\left\{{\lambda }_{i=1,2,3}\right\}\right)=\frac{n}{3}\sum _{i=1,2,3}{F}_{\mathrm{chain}}\left({\lambda }_i\xi \right)=\frac{n{k}_{\mathrm{B}}T}{3}\left[{\pi }^{2}c\left(3-x^2{I}_1\right)+\frac{3-2I_1{x}^2+I_2{x}^4}{\pi c\left(1-I_1{x}^2+I_2{x}^4-I_3{x}^6\right)}\right],\end{array}$$

(25)

where $c=l_p/2L$ and $x=\xi /L$, as used in Equation (8). When the elastic energy is expressed as a function of strain invariants $I_{1,2,3}$, the stress tensor of an incompressible material (with $I_3=1$) can be obtained as [39,40]:

$$\begin{array}{ccc}\hfill {\sigma }_{ij}& =& 2\left[\left(\frac{\partial F}{\partial I_1}+I_1\frac{\partial F}{\partial I_2}\right)C_{ij}-\left(I_1\frac{\partial F}{\partial I_1}+2I_2\frac{\partial F}{\partial I_2}\right)\frac{{\delta }_{ij}}{3}-\frac{\partial F}{\partial I_2}{C}_{ik}{C}_{kj}\right]-P{\delta }_{ij},\hfill \end{array}$$

(26)

where $\mathbf{C}=\mathbf{E}{\mathbf{E}}^{\mathrm{T}}$ is the left Cauchy strain tensor, and P the Lagrangian multiplier for incompressibility, the value of which determined by the boundary conditions. This is also known as constitutive relationship, or stress-strain relationship. Given the compact analytical form of the elastic energy Equation (25), the equilibrium stress-strain relation can now be derived for any applied deformation.

3.3. Applications: Simple Shear, Network Stability and Negative Normal Stress

Consider the simple shear defined in Equation (19), i.e., a deformation when a point in Cartesian coordinates $(x,y,z)$ in an original material will change to $(x+\gamma y,y,z)$ after being sheared with shear strain γ. The incompressibility is satisfied automatically in this deformation, with $I_3=1$ and the other two strain invariants $I_1=I_2=3+{\gamma }^2$. The shear stress obtained from the general constitutive relation (26) of the 3-chain and the 8-chain model are, respectively [38]:

$$\begin{array}{ccc}\hfill {\left[{\sigma }_{xy}\right]}_{3\mathrm{c}}& =& \frac{2}{3}n{k}_{\mathrm{B}}T\gamma x^2\left[\frac{\left(1-x^4\right)}{c\pi {\left[1-\left(2+{\gamma }^2\right)x^2+x^4\right]}^2}-c{\pi }^2\right],\hfill \end{array}$$

(27)

$$\begin{array}{ccc}\hfill {\left[{\sigma }_{xz}\right]}_{8\mathrm{c}}& =& \frac{2}{3}n{k}_{\mathrm{B}}T\gamma x^2\left[\frac{9}{c\pi {\left[3-\left(3+{\gamma }^2\right)x^2\right]}^2}-c{\pi }^2\right].\hfill \end{array}$$

(28)

Both 3-chain and 8-chain models fit the shear-experiment data from Refs. [29,41] equally well, as shown in Figure 8a. We will need more stringent tests to distinguish which of these two models (both offering compact analytical expressions) is working better in filament networks. The stiffness, c, and the initial end-to-end factor, $x=\xi /L$, obtained by fitting with 3-chain model are a bit smaller than those in 8-chain model, but both in the semiflexible regime. Note that in rubber elasticity, there are two types of shear modulus, the nominal shear modulus $G\left(\gamma \right)={\sigma }_{xy}\left(\gamma \right)/\gamma$, and the differential shear modulus $K\left(\gamma \right)=\partial {\sigma }_{xy}/\partial \gamma$ [19,29]. The linear shear modulus of the network, $G_0$, is equal to the differential modulus K for $\gamma \to 0$.

Intrinsic heterogeneity in the 3-chain model, reflected in the fact that the filament lying along the maximum principle stretch direction dominates the response of the whole cubic in semiflexible networks (due to the high nonlinearity of the filament tensile elasticity), is supposed to be closer to the realistic physics of filament network. In Figure 8b, experimental data for a wide variety of semiflexible filaments is fit by Equation (27). Fitted parameters $(c,x)$ in

Table 1. Fitting parameters $(c,x)$ for collagen, actin, vimentin, fibrin and neurofilament data obtained from experiment [29]. Also shown are the linear modulus $G_0$ extracted from the original data and used for scaling in Figure 8, the literature values of $l_p$ with their source, and the calculated mesh size $\xi =l_p(x/2c)$.

	${\mathit{G}}_{\mathbf{0}}\mathbf{\left(}\mathbf{Pa}\mathbf{\right)}$	c	x	${\mathit{l}}_{\mathit{p}}\mathbf{\left(}\mathsf{\mu }\mathbf{m}\mathbf{\right)}$	$\mathit{\xi }\mathbf{\left(}\mathsf{\mu }\mathbf{m}\mathbf{\right)}$
collagen	13.8	1.44	0.85	20.0 [42]	5.9
actin	95.2	1.36	0.85	17.7 [43]	5.5
vimentin	3.82	0.34	0.57	1.0 [44]	0.83
fibrin	18.9	0.25	0.40	0.50 [29]	0.40
neurofilament	2.83	0.14	0.15	0.45 [45]	0.24

, along with other parameters for different filament network coincides with those reported in experiments. An interesting point is that by calculating $x_0=\sqrt{1-c^{*}/c}$ for different filaments listed in Table 1, it is found that the fitted x is close to, but slightly smaller than $x_0$, indicating all of the filaments in Figure 8 are very slightly pre-compressed, rather than pre-stretched in the equilibrium state of the network.

As shown in Table 1, the ratio of fitted parameters $c/x=l_p/2\xi$ remains close to 1, for these semiflexible filaments with various effective stiffness. Moreover, it will become clear that all the networks we examined here lie very near the stability boundary, when discussing about the network stability (Table 1). It remains as a problem whether the fact that the mesh size is close to the filament persistence length is an unintended result of different crosslinking density in experiments [29], or is a relevant and universal biological feature.

An important scaling law in sheared filament network is about its differential modulus and the shear stress, given by $K\sim {\sigma }^{3/2}$, which was found in a number of experiments and theories [31,41,46,47,48,49]. In 3-chain model, the differential shear modulus $K\left(\gamma \right)$ can be expanded in powers of shear strain γ:

$$\begin{array}{c}\hfill K\left(\gamma \right)=\frac{\partial {\sigma }_{xy}}{\partial \gamma }\approx \frac{2}{3}n{k}_{\mathrm{B}}T{x}^2\left[\frac{1+x^2}{\pi c{(1-x^2)}^3}-{\pi }^{2}c\right]+4n{k}_{\mathrm{B}}T{x}^4\frac{1+x^2}{c\pi {(1-x^2)}^5}{\gamma }^2.\end{array}$$

where the first term represents the linear shear modulus $G_0$ in charge of the linear response, and the second term represents the first strain-dependent correction, or the shear-hardening part of the modulus. The crossover between the linear and nonlinear response is obtained when these two terms are comparable with each other. In the limit of a stiff network with $x\to 1$, the crossover, or transition shear strain can be approximated as ${\gamma }_t\sim (1-x^2)/\sqrt{6}$, and the corresponding stress is ${\sigma }_t=G_0{\gamma }_t\sim 2\sqrt{2}n{k}_{\mathrm{B}}T/\left[3\sqrt{3}c\pi {(1-x^2)}^2\right]$, see Figure 9a. This stress ${\sigma }_t$ increases with the temperature as ${\left(k_{\mathrm{B}}T\right)}^2/\kappa$, which coincides with the observations in [49].

When the shear strain approaches the divergence point, i.e., $\gamma \to 1/x-x$ in Equation (27), the stress and the differential modulus can be approximated as,

$$\begin{array}{ccc}\hfill \sigma \left(\gamma \right)& \simeq & n{k}_{\mathrm{B}}T\frac{2\gamma x^2(1-x^4)}{3\pi c{\left[1-(2+{\gamma }^2)x^2+x^4\right]}^2},\hfill \end{array}$$

(29)

$$\begin{array}{ccc}\hfill K\left(\gamma \right)& \simeq & n{k}_{\mathrm{B}}T\frac{2(1-x^4)\left[x^2-(2-3{\gamma }^2)x^4+x^6\right]}{3\pi c{\left[1-(2+{\gamma }^2)x^2+x^4\right]}^3}.\hfill \end{array}$$

(30)

Both expression diverge due to the same vanishing denominator, while maintaining the obvious scaling relation $K\sim {\sigma }^{3/2}$. We see this limit exposed clearly in Figure 9b for several very different actin networks. In this plot, the data points are the same as in Figure 9a, and the predicted curves of $K\left(\sigma \right)$ is obtained by differentiating Equation (27), while using parameters $G_0,c$ and x fitted in Figure 9a, for each dataset. Note that this 3/2 scaling at large deformation, acting as an intrinsic property of filament network, arises from mechanic properties of an individual filament, as

$$\begin{array}{c}\hfill \frac{d\sigma }{d\gamma }\sim \frac{df}{dx}\sim \frac{1}{{\left(1-x^2\right)}^3}\sim f^{3/2}\mathrm{when}x\to 1.\end{array}$$

(31)

One can easily check it with the force-extension relations given in the previous section. Note that in many cases, such as in Unterberger et al. [37], the filament network retains transient activity, which results in stress-softening and plasticity at higher stress due to the breakage of the crosslinks (this is also the case in experiments [50] before the actin filaments formed bundles, or with F-actin crosslinked by filamin). This type of filament network with transient crosslinks will be discussed in the next section. The actual data for the actin network of Storm et al. [29] also deviates from the 3/2 scaling. In fact, the authors of [29] develop their own theory invoking various additional factors (e.g., filament extensibility) to account for this data. However, the basic 3-chain model, assuming permanent crosslinks, bulk incompressibility and inextensible filaments, evidently fits both $\sigma \left(\gamma \right)$ and $K\left(\sigma \right)$ data very well. This fact might indicate that, for these values of c and x, the final crossover to this 3/2 characteristic scaling regime would occur at an even higher stress (at which point the network would probably not survive in practice).

Network stability and non-affine effects. Although the filament network can be pre-compressed, force-free or pre-stretched at its reference state, i.e., $x<,=,\mathrm{or}>x_0$, there is a condition for having a mechanically stable network, given by the linear shear modulus $G_0\ge 0$. In the 3-chain model, using the expression for $G_0$ in Equation (29), this condition translates into:

$$\begin{array}{c}\hfill c\le \frac{1}{{\pi }^{3/2}}\frac{\sqrt{1+x^2}}{{(1-x^2)}^{3/2}}.\end{array}$$

(32)

Note that the right-hand side in Equation (32) is always greater than $c^{*}=1/{\pi }^{3/2}$, as defined in section 2 and illustrated in Figure 2b. Therefore, for sufficiently flexible chains with $c<c^{*}$, the network is always stable; indeed we always expect $G_0>0$ in the rubbery network. On the other hand, for filament networks composed of the more stiff filaments with $c\ge c^{*}$, the network strands have to be crosslinked with $x=\xi /L$ exceeding the pre-tension threshold given by Equation (32), in order for the network to have a finite equilibrium shear modulus. In other words, there usually are forces acting on the crosslinked filaments in the network in order to be mechanically stable with a non-zero shear modulus. This notion is familiar from the “tensegrity” concept in biology and engineering [51]. For stiff athermal filament network, the window of pre-tension between the linear modulus $G_0=0$ at $x\approx 1-1/\pi {\left(2c\right)}^{2/3}$, and $G_0\to \infty$ at $x\to 1$ (assuming inextensible chains) is very narrow. The stability boundary of the network is plotted in Figure 10. The condition for $G_0=0$ gives the equilibrium case labelled as $\gamma =0$ in the phase diagram. The full analysis [38] shows that the magnitude of the shear strain modifies the stability condition, as represented by the additional marginal stability curves for $\gamma =0.2$ and $0.5$. This shift means that the region of mechanical stability of filament network expands on increasing deformation, matching what the recent experimental data of Sharma et al. states [52].

As mentioned above, non-affine deformation in the local context of individual filament junctions and mesh cells, plays an important role when the mesh size is very large compared with the radius of the filament, or when the crosslink density is low. Note that in continuum models, such as 1-, 3-, or 8-chain model discussed here, the “mesh cells” are all deformed in the same way, but differently from the material, as these cells are constructed in orientation aligned with the principal stretching direction of the material, which is sometimes also referred to as “non-affinity”. But here we refer to the concepts of microscopic “affinity” or “non-affinity” for local deformation of individual filaments or mesh cells, rather than the global one between the orientation of a mesh cell and the macroscopic deformation of the material.

When the filament is compressed or stretched with a small deformation strain, the elasticity of the filament can be regarded as approximately linear. In other words, the filament in this case can be approximately as a linear spring with optimal end-to-end length as ${\xi }_0=L\sqrt{1-2L/{\pi }^{3/2}{l}_p}$ and spring constant $\alpha =2{\pi }^{7/2}{k}_{B}T{c}^2$ from Equation (12), i.e., the free energy of the filament becomes approximately $\frac{1}{2}\alpha {(\xi -{\xi }_0)}^2$, where ξ is the end-to-end length of the filament. Note that in a crosslinked network, the length of the filament strand is not necessarily equal to the zero-force (optimal) one, which means the filament can be either pre-stretched, pre-compressed or remain force-free. In an isotropic material experiencing a small shear deformation, the affine part of the linear shear modulus can be calculated by Born-Huang theory as [53,54]:

$$G_0=\frac{1}{V}\sum _{ij}(\alpha {\xi }_{e,ij}-f_{ij}){\xi }_{e,ij}{n}_{x,ij}{n}_{x,ij}{n}_{y,ij}{n}_{y,ij},$$

(33)

where (ij) denotes for a bonded pair, ${\xi }_{e,ij}$ is the equilibrium end-to-end length of the filament bridging points (ij) in a network, $f_{ij}$ is the pre-tension force on the filament bridging the neighboring points, and $n_{ij}$ are unit vectors along the $\left(ij\right)$ line. For such a filament, the force-extension relation gives $f_{ij}=\alpha ({\xi }_{e,ij}-{\xi }_{0,ij})$. By substituting this linear force-extension relation back into Equation (33), the shear modulus takes the form $G_0=(z{N}_c/2V)\alpha {\xi }_e{\xi }_0·\langle n_x{n}_x{n}_y{n}_y\rangle$, where $N_c$ is the number of crosslinks, and z is the number of nearest crosslinks connected with a chosen crosslink, or the average network connectivity. The averaging is over the statistical ensemble of orientations of bond vectors ${\mathbf{n}}_{ij}$.

The non-affine part of the linear shear modulus reduces the modulus by the local force relaxation. Under the same assumptions of harmonic bonds, the non-balanced local relaxation force acting on a crosslink $\left(i\right)$ can be expressed as [54] ${\mathsf{\Xi }}_i=-{\sum }_j(\alpha {\xi }_{e,ij}-f_{ij})n_{ij}{n}_{x,ij}{n}_{y,ij}$. In this case, the cancellation in the bracket leaves no pre-tension contribution and the local force only depends on the force-free distance ${\xi }_0$: ${\mathsf{\Xi }}_i=-{\sum }_j\alpha {\xi }_{0,ij}{n}_{ij}{n}_{x,ij}{n}_{y,ij}$. The detailed derivation [55,56] gives the resulting non-affine part of the modulus: $G_{\mathrm{na}}=(3N_c/V)\alpha {\xi }_0^2·\langle n_x{n}_x{n}_y{n}_y\rangle$. Combining these two parts, we obtain the estimate of effective linear shear modulus:

$$G=G_0-G_{\mathrm{na}}=\frac{N_c}{V}\alpha (\frac{z}{2}{\xi }_e{\xi }_0-3{\xi }_0^2)·\langle n_x{n}_x{n}_y{n}_y\rangle .$$

(34)

If chains are distributed isotropically in the network, then $\langle n_x{n}_x{n}_y{n}_y\rangle =1/15$ [55,56]. The important conclusion of this analysis is that the shear modulus is proportional to $G_0(1-6{\xi }_0/z{\xi }_e)$. That is, a connectivity of $z\ge 6$ is required in the network crosslinked with force-free filaments to be mechanically stable (a classical Maxwell criterion for the assumer central forces). Yet when the filaments are linked with pre-tension (i.e., ${\xi }_e>{\xi }_0$), the point of marginal stability becomes increasingly lower: $z\ge 6{\xi }_0/{\xi }_e$ (again, the “tensegrity” effect at play).

Note that in this schematic analysis of the role of non-affinity in filament network we used a random orientation of filaments as a basic assumption: strictly, $\langle n_x{n}_x{n}_y{n}_y\rangle =0$ for a 3-chain cubic cell, indicating that the 3-chain model cannot be used with nonaffine-effects in a simple way.

Negative normal stress. When being sheared, filament networks are usually found to show negative normal stress [57,58,59]. This is highly unusual and characteristic of the so-called auxetic materials, while the ordinary polymer networks (rubbers and gels) all show the positive normal stress, often visualized as the so-called Weissenberg effect. Most of the experiments until now reporting normal stress measurements in filament networks were conducted in the rotating cylindrical geometry of a standard rheometer, as shown in Figure 11a. Suppose the height and the radius of the undeformed cylinder are $h_0$ and $R_0$, respectively. When being deformed, they change to $h={\lambda }_h{h}_0$ and $R={\lambda }_R{R}_0$, with ${\lambda }_{h,R}$ as the stretch ratio along the longitudinal and the radial directions, respectively, with assumed incompressibility ${\lambda }_h{\lambda }_R^2=1$. In the cylindrical coordinate system, the point at $(r,\theta ,z)$ has its position change to $(r/\sqrt{{\lambda }_h},\theta +\mathsf{\Theta }z/h_0,{\lambda }_{h}z)$, after rotating the top plate with angle Θ (pure shear, rather than simple shear, is used for discussions about the normal stress, because of the uncertain boundary conditions in simple shear case [60]), and the deformation tensor is written as,

$$\begin{array}{c}\hfill \mathbf{E}=\left(\begin{array}{ccc}1/\sqrt{{\lambda }_h}& 0& 0\\ 0& 1/\sqrt{{\lambda }_h}& \gamma \left(r\right)/\sqrt{{\lambda }_h}\\ 0& 0& {\lambda }_h\end{array}\right),\end{array}$$

(35)

where $\gamma \left(r\right)=r\mathsf{\Theta }/h_0$ is the shear strain, depending on the radial position in this parallel-plate geometry. In this case, the first two strain invariants are, $I_1={\lambda }_h^2+[2+{\gamma }^2]/{\lambda }_h$ and $I_2=2{\lambda }_h+[1+{\gamma }^2]/{\lambda }_h^2$. Given the shear strain at the outermost surface, ${\gamma }_0=\gamma \left(R_0\right)$, the total free energy can be expressed as a function of the stretching ratio ${\lambda }_h$ along the z axis, after integration over radius,

$$\begin{array}{c}\hfill F({\lambda }_h;{\gamma }_0)={\int }_0^{R_0}dr{F}_{3\mathrm{c}}({\lambda }_h,r;{\gamma }_0).\end{array}$$

(36)

Then upon minimization of the free energy, the equilibrium ${\lambda }_h$ can be obtained [38]. When the cylinder radius becomes smaller (or larger) upon shear, with its height becoming correspondingly larger (or smaller) – this phenomenon is called the positive (or negative) Poynting effect [61,62]; this geometric effect with stress-free top plate corresponds to the positive (negative) normal stress required to maintain the fixed plate separation in a more common rheometry experiment [57,58].

Figure 11b shows the results for two model materials with different filament stiffness and pre-tension [38]. A flexible network ($c=0.1,x=0.1$) has a positive Poynting effect, or positive normal stress is required to counter the expansion along the height (the Weissenberg effect), where the mesh size is much smaller than the contour length of the subfilaments connecting the neighboring crosslinks ($x_0\ll 1$) thus the entropic energy is important. With the chains being stretched along the principal extension direction in the shear geometry of Figure 11a, the chains in other two principle directions are being compressed, leading to the material contraction along the radial and circumferential directions. Conversely, if the mesh size is comparable to the contour length of the subfilaments, $x_0\le 1$, the material can behave in a negative Poynting effect manner, with contracted height and expanded radius, or exerting the negative normal stress is required to maintain the fixed height. This is because the force acting on the chains directed along the principal extension in Figure 11a is close to a divergence if $x_0\to 1$, and it causes less energy when the material is compressed in the longitudinal direction, rather than stretched. Similar reason works also for a network of more stiff filaments (see red curve for $c=0.3>c^{*},x=0.6$). The magnitude of normal stress in a fixed height geometry can be estimated from the linear relationship ${\sigma }_h\sim 3G_0({\lambda }_h-1)$, which works for the case where the normal strain is small and the Young modulus is $3G_0$. Taking the fibrin in Table 1 as an example, the normal stress is about 11.5 Pa under a shear strain of ${\gamma }_0=0.5$. These are very close to the observations in experiments [57,58]. The magnitude of ${\sigma }_h\simeq 20$ Pa for the same shear also matches the actin results obtained in [37].

Figure 11c gives the full ’phase diagram’ of the stiffness-pretension parameter space ($c,x$), with phase boundaries separating stable/unstable regions and positive/negative normal stress regions. There is a not shown weak dependence of this boundary on the magnitude of the applied shear strain, due to the inherent non-linearity of stress-strain response. Figure 11c shows that a sheared flexible network with $c<c^{*}$ usually shows a positive normal stress, while a semiflexible filament network has a negative negative normal stress, especially when the filaments are crosslinked with increasing pre-tension.

The above model explains the negative normal stress in equilibrium. de Cagny et al. [63] recently observed that the normal stress could be positive at short times after being sheared, and then becomes negative after a sufficient long relaxation time, which is dependent on the mesh size of the network. Clearly this question is yet far from being fully resolved.

4. Transient Filament Network with Breakable Crosslinks

In contrast to a permanently crosslinked network such as classical rubber, biological crosslinks in a filament network can be dynamic. The classical collagen network is held together by hydrogen-bonded triple-helical segments that are known to dynamically respond to an applied tension. The F-actin network is often held together by molecular motors (such as myosins) or physically bonded proteins (such as filamins). Self-assembled biological filaments themselves can be broken from or re-crosslinked to the network. In synthetic networks, one frequently finds dynamic crosslinks formed by physical bonds (such as hydrophobic interactions or hydrogen bonds, Figure 12). Actually, the concept of a transient crosslink here is very similar to the “trap” defined in soft glassy systems [64,65], where semiflexible chains are assumed to be trapped in a glassy environment, and also to the dynamic “entanglement” in concentrated polymer solutions [1].

Due to the transient feature of such crosslinks, the resulting network shows characteristic rheological properties, such as a high stress relaxation at given constant strain, a yield stress upon ramp deformation, or the elastic-plastic transition. As transient filament network is a better analogy to the filament network in many biological systems, there are several relevant theories developed in the recent years from different viewpoints.

In 2007, Kroy et al. proposed the glassy wormlike chain model [66,67], and then used it to study the viscoelastic properties of “sticky filament network” [68]. One assumes there is a certain crosslink breaking rate β, and also a rate of crosslink re-connection, ρ. In adiabatic approximation, the rates can be expressed by the activation laws (Figure 13):

$$\begin{array}{c}\hfill \beta {\tau }_0=e^{[-E+(x_t-x_b)f]/k_{B}T},\rho {\tau }_0=e^{[-E+U-(x_u-x_t)f]/k_{B}T},\end{array}$$

(37)

where ${\tau }_0$ is an intrinsic time scale of attempts due to the thermal fluctuations, E is the energy barrier for breaking the crosslink, U is the energy penalty for the broken crosslink, f is the pulling force, $x_b$ is the coordinate of crosslinked state, $x_t$ is the coordinate of the barrier, and $x_u$ is the coordinate of the uncrosslinked state (crosslink broken). The fraction of the complete crosslinks is interpreted via the effective change in filament length: $n\equiv L_0/L\left(t\right)$ (the fraction of broken crosslinks is $1-n$), where $L_0$ is the minimum contour length of the filament between two crosslinks, and $L\left(t\right)$ is the contour length of the filament between crosslinks at time t, which can change with the breakage and the formation of the crosslinks. Then the evolution equation of the fraction of the complete crosslinks is

$$\begin{array}{c}\hfill \dot{n}{\tau }_0{e}^{E/k_{B}T}=-\left(e^{[U-(x_u-x_t)f]/k_{B}T}+e^{(x_t-x_b)f/k_{B}T}\right)n+e^{[U-(x_u-x_t)f]/k_{B}T}.\end{array}$$

(38)

Assuming that the linear shear modulus is proportional to the density of active crosslinks, the main difference in the complex modulus between the ordinary filament and the glassy semiflexible one in transient network is in the relaxation time ${\tau }_n$, which is a constant ${\tau }_{n0}$ in the first case, while for glassy semiflexible filament it depends on time via the $L\left(t\right)$ relaxation:

$$\begin{array}{ccc}\hfill {\tau }_n& =& \left\{\begin{array}{c}{\tau }_{n0},{\lambda }_n<L\left(t\right)\hfill \\ {\tau }_{n0}{e}^{E({\lambda }_n/L\left(t\right)-1)/k_{B}T},{\lambda }_n\ge L\left(t\right)\hfill \end{array}\right.,\hfill \end{array}$$

(39)

with ${\lambda }_n$ as wavelength of the nth mode. Microscopic complex modulus for transverse point excitations of a generic filament element was defined in this model as

$$\begin{array}{c}\hfill g\left(\omega \right)=g^{\prime }\left(\omega \right)+i{g}^{″}\left(\omega \right)=1/\left[\frac{L^3}{l_p{\pi }^4}\sum _{n=1}^{\infty }\frac{1}{(n^4+n^{2}f/f_c)(1+i\omega {\tau }_n)}\right],\end{array}$$

(40)

where f is filament tension, and $f_c$ the Euler buckling force for rigid rod as defined previously.

In 2008, Lieleg et al. [69,70] performed a rheological experiment on actin network crosslinked by rigor heavy meromyosins (HMM), which is yet another kind of transient crosslinks. They reported local stress relaxation and energy dissipation in an intermediate elasticity-dominated frequency regime. In order to explain this, they developed a semi-phenomenological model by introducing a breakage rate of the crosslinks, ${\beta }_0\approx 0.09{\mathrm{s}}^{-1}$, which is a simpler version of the glassy filament network model above. By assuming a breakage rate of crosslinks, the number of the crosslinks decays exponentially with time as $N\left(t\right)\sim N\left(0\right)exp(-{\beta }_{0}t)$, which can be translated into the frequency space by Fourier transformation, $\widehat{N}\left(\omega \right)$, where ω is the oscillation frequency in the experiment. With the rubber-elasticity relation $G\propto N$, the storage and the loss modulus are written as functions of the oscillation frequency:

$$\begin{array}{c}\hfill G^{\prime }\left(\omega \right)=G_0-a\frac{N{\beta }_0}{{\beta }_0^2/4{\pi }^2+{\omega }^2}+b{\left(\frac{\omega }{{\omega }_0}\right)}^{3/4}\mathrm{and}{G}^{″}\left(\omega \right)=c\frac{N\omega }{{\beta }_0^2/4{\pi }^2+{\omega }^2}+d{\left(\frac{\omega }{{\omega }_0}\right)}^{3/4},\end{array}$$

(41)

where the terms with fitting parameters a and c arise from the breakage of the crosslinks, via $\widehat{N}\left(\omega \right)$, the ${\omega }^{3/4}$ terms with parameters b and d arise from the fluctuation of single filaments [30,71], and the time scale of the relaxation mode is set by ${\omega }_0$ which is dependent on the solvent viscosity. There is a characteristic frequency in experiments, ${\omega }_{\min }$, at which the dissipation has a minimum. When the frequency ω is larger than ${\omega }_{\min }$, the storage modulus is almost a constant, while decreasing at lower frequencies due to the breakage of crosslinks at longer time scales. This semi-phenomenological model provides a basic understanding in the viscoelastic response of a transient network at different time scales.

In 2010, Broedersz et al. [72] proposed a model of transient filament networks based on crosslink-governed dynamics, where Monte Carlo simulations and an analytic approach are combined for explaining the rheological properties of the network. There are three time scales in the model: the thermal relaxation time for chains after breaking off a crosslink, ${\tau }_{\mathrm{eq}}$, the re-crosslink time, ${\tau }_{\mathrm{r}}$, and the crosslink breakage time, ${\tau }_{\mathrm{b}}$. In this model they are assumed to obey a relation: ${\tau }_{\mathrm{eq}}\ll {\tau }_{\mathrm{r}}\ll {\tau }_{\mathrm{b}}$. The elastic energy of the network can be expressed as a sum of the bending energy and the stretching energy of all the subfilaments connecting the pairs of crosslinks:

$$\begin{array}{c}\hfill F=\frac{1}{L}\sum _n\left[\frac{\kappa }{2}{|\Delta {\mathbf{t}}_n|}^2+\frac{{\mu }_{\mathrm{th}}}{2}{\left(|\Delta {\mathbf{r}}_n|-L\right)}^2\right],\end{array}$$

(42)

where L is the length of a strand between the two crosslinks, which is assumed to be the same in the whole network, the sum is taken for all the crosslink positions ${\mathbf{r}}_n$ (with the “mesh size” $\Delta {\mathbf{r}}_n={\mathbf{r}}_{n+1}-{\mathbf{r}}_n$), and ${\mathbf{t}}_n$ is unit tangent vector of nth crosslink (sketched in Figure 14. In this case, the nth segment curvature is $|\Delta {\mathbf{t}}_n|/L$, defining the first term (the bending energy). The second term in Equation (42) denotes for the elastic energy of the spring, where the stretching modulus ${\mu }_{\mathrm{th}}$ is directly related to the discussion below Equation (11) in Section 2. In order to calculate the polymer displacement due to the unbinding and rebinding, they separate the local equilibration step into a move to the mechanical equilibrium position, with an added stochastic contribution around this equilibrium. Mechanical relaxation step from the initial position ${\mathbf{r}}_n^i$ to the local equilibrium position ${\mathbf{r}}_n^{eq}$ is determined by

$$\begin{array}{c}\hfill {\left.\frac{\partial F}{\partial {\mathbf{r}}_n}\right|}_{{\mathbf{r}}_n={\mathbf{r}}_n^{eq}}=0.\end{array}$$

(43)

Actually, this numerical method can describe non-affine deformation in filament networks. By solving Equation (43) discretely and taking the continuum long-wavelength limit, the leading order evolution equations with thermal fluctuations are

$$\begin{array}{c}\hfill {\tau }_b\frac{\partial r_{\parallel }}{\partial t}=\frac{l_c^2}{2}\frac{{\partial }^2{r}_{\parallel }}{\partial x^2}+{\mathbf{e}}_x·{\mathit{\xi }}_{\perp },{\tau }_b\frac{\partial {\mathbf{r}}_{\perp }}{\partial t}=\frac{l_c^2}{2}\frac{{\partial }^2{\mathbf{r}}_{\perp }}{\partial x^2}+{\mathit{\xi }}_{\perp }.\end{array}$$

(44)

where x denotes for the position of the cross-linker along the average direction ${\mathbf{e}}_x$ of the relaxed filament, and $r_{\parallel }$ and ${\mathbf{r}}_{\perp }$ are the longitudinal and transverse deflections of the filament with respect to its average direction. The noise term ${\mathit{\xi }}_{\perp }$ possesses both thermal effects and also the local buckling contributions from thermal compression. Note that the noise actually depends on the local state of the polymer and couples the longitudinal and the transverse motion of a filament as described Equation (44). In the limit of small stretch modulus, $\mu \ll {\mu }_{\mathrm{th}}$, Equation (44) decouple and the resulting transverse contribution to the power spectrum ($C\left(\omega \right)=\langle |\delta l\left(\omega \right){|}^2\rangle$ with $\delta l$ as extension of length) approaches $C_{\perp }\sim {\omega }^{-7/4}$, in agreement with simulation. In this limit, the transverse bending dynamics are effectively those of a stiff filament fluctuating in a viscous solvent, which gives ${\langle \left|\delta l\left(t\right)\right|}^2\rangle \sim t^{3/4}$ [30,71]. Recently, this Langevin-type model was applied to active network systems [73].

In 2013, Unterberger et al. [74] applied a generalized Maxwell model to the transient filament network for describing its viscoelastic properties, as shown in Figure 15. In addition to the pure elastic part of the material, the model adds another m parallel Maxwell elements, each of which consists of an elastic spring and a viscous dashpots. In studies of viscoelasticity and relaxation, this is a common method of introducing the spectrum of relaxation times that is often necessary to describe complex relaxation processes.

Rather than using Cauchy stress $\mathit{\sigma }$, as defined in Equation (26), they applied the second Piola-Kirchhoff stress, $\mathbf{S}$, which is related with Cauchy stress as: $\mathbf{S}=\sqrt{\det \mathbf{C}}·{\mathbf{E}}^{-1}·\mathit{\sigma }·{\mathbf{E}}^{-T}$. There are three distinct contributions to the total stress of the material, which is expressed as

$$\begin{array}{c}\hfill \mathbf{S}={\mathbf{S}}_{vol}^{\infty }+{\mathbf{S}}_{iso}^{\infty }+\sum _{n=1}^m{\mathbf{Q}}_n.\end{array}$$

(45)

where the first term denotes the equilibrium volumetric part, the second term denotes the equilibrium isochoric-elastic part, and the third term reflects the isochoric relaxation (non-equilibrium) part. The first two contributions are typical elastic ones. Here we will only focus on the third part, which is closely related to the “transient” characteristics of the filament network. The volume-preserving non-equilibrium stress tensors ${\mathbf{Q}}_n$, are defined as a function of internal variables, ${\mathsf{\Gamma }}^n$, through the dissipative potential, ${\mathsf{{\rm Y}}}_{iso}$,

$$\begin{array}{ccc}\hfill {\mathbf{Q}}_n& =& {\left(\det \mathbf{E}\right)}^{-2/3}(\mathbf{1}-{\mathbf{B}}^{-1}\otimes \mathbf{B}):\frac{2\partial {\mathsf{{\rm Y}}}_{iso}^n(\overline{\mathbf{B}},\overline{\mathsf{\Gamma }})}{\partial \overline{\mathbf{B}}},\hfill \end{array}$$

(46)

where $\mathbf{B}$ is the right Cauchy tensor $\mathbf{B}={\mathbf{E}}^T\mathbf{E}$, $\overline{\mathbf{B}}=\det {\mathbf{E}}^{-2/3}\mathbf{B}$ is the isochoric part of the right Cauchy tensor, and $\overline{\mathsf{\Gamma }}$ is the isochoric part of the Γ. A set of differential equations and initial conditions for the non-equilibrium part of the nth Maxwell element, ${\mathbf{Q}}_n$ is introduced to describe the transient property of the material, as

$$\begin{array}{c}\hfill {\dot{\mathbf{Q}}}_n+\frac{{\mathbf{Q}}_n}{{\tau }_n}={\theta }_n{\dot{\mathbf{S}}}_{iso}^{\infty },\mathrm{with}{\dot{\mathbf{Q}}}_n(t=0)={\left[{\left(\det \mathbf{E}\right)}^{-2/3}(\mathbf{1}-{\mathbf{B}}^{-1}\otimes \mathbf{B}):\frac{2\partial {\mathsf{{\rm Y}}}_{iso}^n(\overline{\mathbf{B}},\overline{\mathsf{\Gamma }})}{\partial \overline{\mathbf{B}}}\right]}_{t=0},\end{array}$$

(47)

where ${\theta }_n$ are non-dimensional free energy parameters, related to the relaxation times ${\tau }_n$. Convolution integrals provided a solution for Equation (47) of the nth Maxwell element as

$$\begin{array}{c}\hfill {\mathbf{Q}}_n={\left.{\mathbf{Q}}_n\right|}_{t=0}{e}^{-t/{\tau }_n}+{\int }_0^{t}d{t}^{\prime }{e}^{-(t-t^{\prime })/{\tau }_n}{\theta }_n{\dot{\mathbf{S}}}_{iso}^{\infty }\left(t\right),\end{array}$$

(48)

which basically describes how the energy of the material dissipates with time due to the transient crosslinks, and the time scale for the nth Maxwell element is given by its the relaxation time, ${\tau }_n$.

One could instead apply a dynamic theory of transient rubbery networks [75,76], which was based on classical Gaussian (neo-Hookean) chain response, to the case of filament network. This strategy is similar to the theory of equilibrium filament network, Equations (25) and (27), where the specific form of a non-linear elastic response of the filament strand was implemented in a continuum model. By combining the free energy form of a permanent filament network with the transient characteristics of its crosslinks, we could work out the time-dependent response to various modes of dynamic deformation.

We shall use the crosslink dynamics very similar with those in glassy wormlike filament network model, Equation (37). The rate of a chain to break from a crosslink depends on the force f acting on the chain as: $\beta ={\beta }_0{e}^{fa/k_{\mathrm{B}}T}$, where ${\beta }_0$ is the breakage rate for a force-free chain, and a is the monomer size and taken as a characteristic length of a bond. As in Equation (37), ${\beta }_0{\tau }_0=e^{-E/k_{B}T}$. The dangling chains are able to re-crosslinked into the network with a rate ρ, determined by the barrier $E-U$ (cf. Figure 13), although here we will assume that the free filament would be able to relax and re-crosslink in the force-free state (which is the basis of the stress relaxation in the material).

Suppose the total number of filaments in the network is $N_{\mathrm{tot}}$, with $N_c\left(0\right)$ of them crosslinked at the time $t=0$; $N_b\left(0\right)=N_{tot}-N_c\left(0\right)$ chains are un-crosslinked. The breakage rate of the crosslinks is a function of the local force acting on the given chain, rather than as a constant,and thus can be only defined locally and dynamically: it depends on both the current state at time t and the crosslinking state at time $t^{\prime }$ that defines the reference state for the given chain. We denote the breakage rate of the crosslinks at time t formed at time $t^{\prime }$ as ${\beta }_{t;t^{\prime }}$. After an infinitesimal time interval, $\Delta t$, the number of the originally crosslinked chains, will decrease to $N_c\left(0\right)[1-{\beta }_{0;0}\Delta t]\simeq N_c\left(0\right)exp[-{\beta }_{0;0}\Delta t]$, where ${\beta }_{0;0}$ denotes the breakage rate at time $t=0$. Meanwhile, there will be $N_b\left(0\right)\rho \Delta t$ new filaments re-crosslinked into the network in the same time. By combining these two parts, the total number of the crosslinked chains at time $t=\Delta t$ then becomes: $N_c\left(0\right)exp[-{\beta }_{0;0}\Delta t]+N_b\left(0\right)\rho \Delta t$. After another time interval, i.e., at $t=2\Delta t$, the number of the crosslinked chains surviving from the beginning decreases further to $N_c\left(0\right)exp[-{\beta }_{0;0}\Delta t-{\beta }_{\Delta t;0}\Delta t]$; the number of chains which were crosslinked during the first time interval decreases to $N_b\left(0\right)\rho \Delta texp[-{\beta }_{\Delta t;\Delta t}\Delta t]$; finally, the number of newly crosslinked chain during the second time interval is $N_b(\Delta t)\rho \Delta t$. The total number of filaments that are crosslinked (i.e., a part of the elastic network) at time $t=2\Delta t$ is $N_c\left(0\right)exp[-{\beta }_{0;0}\Delta t-{\beta }_{\Delta t;0}\Delta t]+N_b\left(0\right)\rho \Delta texp[-{\beta }_{\Delta t;\Delta t}\Delta t]+N_b(\Delta t)\rho \Delta t$. By repeating this process, the total number of the crosslinked chains after m time intervals will become

$$\begin{array}{c}\hfill N_c(m\Delta t)=N_c\left(0\right)exp[-\sum _{i=0}^{m-1}{\beta }_{t=i\Delta t;0}\Delta t]+\sum _{j=0}^{m-1}{N}_{\mathrm{b}}(j\Delta t)\rho \Delta t{e}^{-{\sum }_{k=j+2}^m\beta (k\Delta t;[j+1]\Delta t)\Delta t},\end{array}$$

(49)

which can be expressed in continuous form as

$$\begin{array}{c}\hfill N_c\left(t\right)=N_c\left(0\right)e^{-{\int }_0^t\beta (t^{\prime };0)d{t}^{\prime }}+{\int }_0^t{N}_{\mathrm{b}}\left(t^{\prime }\right)e^{-{\int }_{t^{\prime }}^t\beta (t^{\prime \prime };t^{\prime })d{t}^{\prime \prime }}\rho d{t}^{\prime },\end{array}$$

(50)

where the first term on the right hand side represents the number of the originally crosslinked chains surviving from the beginning, and the second term represents those chains that were crosslinked between $t=0$ and t, and still remain crosslinked [76]. For the first term, the elastic reference state remains the original undeformed configuration of the network at $t=0$. In contrast, the tension force that determines the rate $\beta (t_1,t_2)$ in the second term depends on the reference state at a time $t_2$ when the particular filament was re-crosslinked into the network at its force-free conformation.

As the reference states of the filaments are defined according to when they are crosslinked, the elastic energy density of the system is

$$\begin{array}{ccc}\hfill F_{\mathrm{t}.\mathrm{n}}\left(t\right)& =& e^{-{\int }_0^t\beta (t^{\prime };0)d{t}^{\prime }}{F}_{\mathrm{p}.\mathrm{n}}(t;0)+{\int }_0^t\frac{N_{\mathrm{b}}\left(t^{\prime }\right)}{N_c\left(0\right)}\rho e^{-{\int }_{t^{\prime }}^t\beta (t^{\prime \prime };t^{\prime })d{t}^{\prime \prime }}{F}_{\mathrm{p}.\mathrm{n}}(t;t^{\prime })d{t}^{\prime },\hfill \end{array}$$

(51)

where the free energy form of the permanent network, $F_{\mathrm{p}.\mathrm{n}}$, can be chosen from the continuum models showed in the previous section, e.g., Equation (25). As with the rate of breaking, the expression for $F_{\mathrm{p}.\mathrm{n}}(t_1;t_2)$ depends on the deformation at time $t_1$, measured with respect to the reference state (at $t_2$) when the fraction of the chains that is counted towards this contribution was re-crosslinked. The corresponding deformation tensor, measured at $t_1$ with respect to the reference state at $t_2$ is: $\mathbf{E}(t_1,t_2)=\mathbf{E}(t_1,0){\mathbf{E}}^{-1}(t_2,0)$. With the formal expression for the dynamic free energy of the material, one can obtain the stress-strain relationship and other rheological characteristics of the transient network, once the regime of imposed deformation is specified. Two simple applications of the theory are explored below:

Stress relaxation. If the transient network is subjected to a constant uniaxial extension with stretch ratio λ, imposed instantaneously at $t=0$ (stretch ratios along the other two orthogonal directions are $1/\sqrt{\lambda }$ from the incompressibility condition), then the tensile stress along the stretching direction decays with time purely due to the breakage of crosslinks: any crosslinks that are re-established during this process remain in the force-free state and the second term in Equation (51) does not contribute to the response at all. The rate of this exponential decay depends on λ but remains constant over the course of relaxation:

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{ela}}={\sigma }_0{e}^{-\beta t},\end{array}$$

(52)

where ${\sigma }_0$ is the instant stress at $t=0$, which is equal to $d{F}_{\mathrm{p}.\mathrm{n}}/d\lambda$. The breakage rate of the crosslink in a primitive cubic cell of a 3-chain model, ${\beta }_0{e}^{fa/k_{\mathrm{B}}T}$, can be calculated by taken the force in the equation as $max\left[\right|f\left(\lambda x\right)|,|f(x/\sqrt{\lambda })\left|\right]$, which is the maximum force acting on the three orthogonal chains of a crosslink. The explicit form of the force can be calculated by Equation (12). Assuming relatively stiff filaments with $c>c^{*}$, originally crosslinked near their force-free point (as most biological filaments studied in Figure 10 appear to be), the rate of relaxation takes the form

$$\begin{array}{c}\hfill \beta \left(\lambda \right)\approx \frac{1}{{\tau }_0}{e}^{-E/k_{B}T}·exp\left[\frac{16{\pi }^2{l}_p}{5L}{(\lambda -1)}^2\right],\end{array}$$

(53)

where the additional exponential factor arising from the tensile force on a filament is approximated for small deformations $(\lambda -1)\ll 1$ to allow easy linearization. This characteristic exponential stress relaxation is a feature of all transient networks with a single value of the bare rate of crosslink breaking (a single value of the breaking barrier E). This is often the case in biological networks, as well as in the new vitrimer materials [77], while in physically crosslinked synthetic polymer networks one often finds a broad distribution of barriers, and the breaking rates. In such systems, the observed stress relaxation is usually the stretched exponential [76].

Ramp deformation. This deformation regime occurs when the material undergoes a continuously increasing uniaxial stretch, i.e., the stretch ratio changes with time as $\lambda =1+\dot{\gamma }t$, with a constant rate $\dot{\gamma }$. Contributions from the re-crosslinks can be usually ignored, since this rate remains constant while the rate of breaking continuously increases with time [78]. In this case the tensile engineering stress can be obtained as

$$\begin{array}{ccc}\hfill \sigma & =& e^{-{\int }_0^t\beta (t^{\prime };0)d{t}^{\prime }}\frac{\partial F_{\mathrm{p}.\mathrm{n}}}{\partial \lambda }.\hfill \end{array}$$

(54)

The 3-chain model is chosen for $F_{\mathrm{p}.\mathrm{n}}\left(\lambda \right)$ as an example. The network can respond in three distinct regimes: increase stress as an elastic solid, plastically flow as a fluid, or exhibit a necking instability when the two regimes coexist along the length of the sample, as illustrated in Figure 16a. The transitions between these regimes depend on how the imposed strain rate $\dot{\gamma }$ compares with the bare breaking rate ${\beta }_0$ (at zero tensile force on a crosslink). These regimes arise from the competition between the intrinsic non-linear elasticity of the material and the plasticity due to the breakage of the crosslinks.

Figure 16b illustrates a typical predicted response. When the reduced strain rate is large, the material just responses as pure elastic one as the chains have no time to break from the crosslinks. Note that the stress rapidly falls to zero when the strain approaches its critical value where the free energy of the network diverges, as the filaments experience infinitely large forces which break the filaments from the crosslinks. When the reduced strain rate is small, the stress of the material first increases due to elasticity of the network with its crosslinks still largely intact, and then reaches the yield point, following which the stress decreases (the network fluidises) due to an increasing number of broken crosslinks. When the reduced strain rate is moderate, the stress of the material first increases then decreases after a yield poin—however, the stress of the material will increase again due to a non-negligible number of the crosslinks surviving from the beginning and the increasingly high force some of these filaments exert. In this intermediate case the material shows the classical necking instability, where a strongly stretched region (plastic zone) and a weakly stretched region (elastic zone) coexist at the same level of stress across the sample.

5. Conclusions

This is a brief overview of the current state of theoretical understanding of physical properties of semiflexible filaments and their crosslinked networks, permanent or transient. This is a broad field with well over 30 years of active research and seminal contributions, and we could not possibly present all the references that perhaps ought to be here; our aim was to illustrate the key concepts as we see them presently, and hopefully assign these concepts to the right original work. In all cases we tried to leave the reader with the “final result”—a compact analytical expression for the relevant physical system, which can be manipulated and fitted to experimental data. The three main sections of this review start from the filaments themselves, examining their mechanical response (essentially force-extension, although the compression becomes a relevant matter for stiffer and athermal filaments). Here the main features are the non-linearity and divergence of the force-extension, the average natural extension of the stiffer filaments, and their metastable buckling on compression. The next section overviews the equilibrium elastic free energy of a randomly crosslinked network of such filaments, where the main physical effects are the specific strain non-linearity, the negative normal stress usually only encountered in auxetic materials, and the finite range of network stability. The final section examines the constitutive relation in a filament network with dynamic crosslinks, which is a very new topic with much less solid knowledge acquired to date. Still, we were able to demonstrate how the careful analysis of crosslink dynamics allows one to predict the network behavior in stress-relaxation and strain-ramp experiments.