In fact, even when simulating just the formation of filaments one may need to resort to coarse-grained models, which generally attempt to describe the folding and self-assembly processes of the biomolecules in an implicit solvent using effective interactions [53]. Unfortunately, there are not many studies in the literature that explore coarse-grained approaches to describe the self-assembly processes of semiflexible filaments [54, 55], and their computational implementation correspond to a challenge itself, as it might involve nucleation pathways that usually requires special rare-event sampling techniques [56]. Alternatively, one may resort to multiscale modelling approaches like the one introduced in Ref. [38], where a lattice model with anisotropic interactions was used to simulate the formation of the fibers, and the resulting network structure was considered as the input configuration for an effective elastic model (see Figure 2D).

Accordingly, in order to obtain the viscoelasticity of solutions at experimentally relevant time and length scales, one has to rely on coarse-grained models even at the single filament level. In that case, the individual filaments are usually described by discrete chains where N beads are connected through springs or rods (see Figure 2B). The simplest potential for the springs is the hookean, or harmonic, potential, U h = ( κ / 2 ) ∑ j = 1 N − 1 ( r ⃗ j + 1 − r ⃗ j ) 2 , with κ being the elastic constant and r ⃗ j the position vector of the jth bead. Such potential is popular because it provides results for pure flexible filaments that can be conveniently compared to the theoretical predictions of the Rouse model [15]. However, if only the hookean potential is included in the model, then the filament will not display a definite contour length L ≈ Nℓ _b like the bead-rod and freely-jointed chain models [13, 47], since it allows the beads to overlap each other and the expectation value for the bond length is only an average value given by [47] b ≈ 3 k B T / κ . Alternatively, the springs between consecutive beads are often modelled by the finitely extensible non-harmonic elastic (FENE) potential [57] U _FENE, which is locally equivalent to U _h for small deformations, but yields more precise values for ℓ _b along the filament. In addition, it is common to consider excluded volume potentials between nonbonded beads such as the shifted and truncated repulsive Lennard-Jones potential also known as the Weeks-Chandler-Andersen (WCA) potential [58] U _WCA. The semiflexibility of the filament can be explicitly considered through bending potentials, which can be written as [29, 30, 43, 59, 60] U b , θ = κ p ∑ k 1 − cos ( θ k ) , with θ _k being the angle between the bonds that connect three successive beads, and κ _p is the bending modulus, or bending constant [33]. Alternatively, the bending potential may be approximated by [23, 39, 47, 61] U b , r = ( κ b / 2 ) ∑ j = 2 N − 1 r ⃗ j + 1 − 2 r ⃗ j + r ⃗ j − 1 2 , where κ _b controls the strength the bending energy which is also somewhat related to the changes in θ _k . Although κ _b (in pN/nm) is expected to be proportional to κ _p (in pN.nm²), the potential U _b,θ is sometimes preferred because its constant can be directly related to the persistence length of the filament, L _p ≈ κ _p /k _B T, as it is formally defined as the correlation length between consecutive segments in the filament [33, 62, 63]. The persistence length L _p is the most basic property of a filament and it can be used to categorize it as flexible (L _p ≪ L), semiflexible (L _p ≤ L), and rod-like (L _p ≫ L). Although the above potentials are assumed in most of the simulations, limitations may occur especially when the filaments approach the rod-like regime, so alternative coasening modelling approaches have been also considered [64].

Finally, it is worth noting that, besides the already mentioned excluded volume and bending effective interactions, implicit effects on the bending rigidity of the filaments may also occur due to other sources. For instance, interactions between charged beads in the filament (and possibly) with ions in solution can be incorporated through bare (or screened) Coulomb potentials [65, 66]. In addition, at the coarse-grained level, hydrodynamics effects might be also modelled as “hydrodynamic interactions” between beads [51, 59].

Since the intrinsic relaxation modulus [20] [G(t)] can be retrieved from the stress-stress autocorrelation function (Figure 2B), one can use the dynamics of a single filament to obtain the intrinsic shear moduli [G*(ω)] for an infinitely dilute solution. In that case one might estimate the actual modulus of dilute solutions by multiplying the intrinsic modulus by the number density of the filaments n _f [15]. Usually, the dynamics of single filaments is obtained either from molecular dynamics [45] (MD) or from stochastic mesoscale approaches [40] (see Figure 2B). In particular, Ref. [43] includes results obtained for the shear-dependent viscosity η ( γ ̇ ) evaluated through the reverse perturbation method [42] (see Figure 2A), indicating that higher the bending constant κ _p higher the value of η ₀. As mentioned in Section 2.1, it is also possible to obtain the shear moduli from approaches based on microrheology (see Figure 2F), and Ref. [39] used a relaxation method based on fluctuation-dissipation theorem to obtain the response of dilute solutions and showed that the range where G′ ∝ G′′ ∝ ω ^3/4 increased for higher values of the bending constant κ _b .

In principle, models for entangled solutions can be obtained simply by including a large number M of filaments in a simulation box with volume V, so that n _f = M/V. In that case, the dynamics of a system with several entangled chains can be also obtained from full simulations [26] (see Figure 2B). However, as detailed in Ref. [47], this kind of approach might face limitations as the molecular weight of the filaments exceeds a “critical” value, and alternative methods may be required ² . As mentioned in Section 2.1, simulations based on primitive path analysis (see Figure 2B) have been successfully used to obtain the plateau moduli G _e of entangled solutions which are consistent with the values that were experimentally determined for many polymer melts [29, 30]. In that case, the semiflexible filaments are described by the so-called Kremer-Grest model, which includes U _FENE, U _WCA, and bending U _b,θ potentials (see Section 2.2.2).

Unfortunately, without many bottom-up approaches that incorporate the self-assembly of filaments (see, e.g., Figure 2D), it is sometimes difficult to generate and equilibrate systems with disordered cross-linked networks. Even so, a few procedures have been developed so that generic features of fully formed networks can be systematically studied. In this context, protocols for constructing ad hoc configurations (Figure 2E) include, e.g., (i) erasing a fraction p of the bonds of pre-established regular networks [37, 68, 69], and (ii) placing line segments in the system at random until the required crosslinking density is reached [35, 36, 70]. It is worth noting that Monte Carlo-based schemes [71–73] have been conveniently used to equilibrate networks generated with the protocol (ii). Besides the number density of filaments n _f , the contour length L, and the persistence length L _p , useful quantities that can be used to characterise cross-linked networks of semiflexible filaments are the mean distance between crosslinks ℓ _c (Figure 2E), which also defines a crosslinking density ρ _c = 1/ℓ _c , and the mean network connectivity ⟨z⟩ (Figure 2D). In particular, the systematic studies presented in Refs. [36, 37] computed the plateau modulus G ₀ to obtain L-vs-ρ _c and z-vs-L _p phase diagrams, respectively. Interestingly, those studies indicated the presence of nonaffine and bending-dominated viscoelastic responses at small values of ⟨z⟩ and ρ _c , which have been also observed for heterogeneous networks [38]. In Ref. [38], the cross-linked networks were generated through a self-assembly process using a lattice model (Figure 2D), and were also explored in Ref. [34]. In the latter reference one can find the method based on the Hessian of the elastic energy (see Section 2.1), which allows one to assess the contributions of both affine and nonaffine deformations (Figure 2E) to the shear moduli G′(ω) and G′′(ω). In addition, Refs. [34, 38] considered a generalisation of the freely-hinged model used in Ref. [37] that incorporates the influence of heterogeneous structures (see Figure 2D), i.e., filaments with thickness-dependent stretching and bending constants, into the effective elastic energy E. As discussed in Ref. [33], semiflexible filaments are less prone to entangle than the flexible ones, so the viscoelastic response of their networks might rely much more on the cross-linker properties. For instance, the role of the flexibility of cross-linkers has been investigated in Ref. [74] in arbitrarily generated networks, with the value G ₀ computed from the derivatives of the total elastic energy of the system. In addition, the study in Ref. [75] investigated the effects of transient cross-linkers on the viscoelasticity of networks of stiff biopolymers, showing that they can lead to a wide distribution H(τ) yielding the power law behaviour observed for the shear moduli at low frequencies.