Physics-Based Computational Approaches to Compute the Viscoelasticity of Semiflexible Filamentous Biomaterials

Rizzi, L. G.

doi:10.3389/fphy.2022.893613

MINI REVIEW article

Front. Phys., 20 June 2022

Sec. Soft Matter Physics

Volume 10 - 2022 | https://doi.org/10.3389/fphy.2022.893613

This article is part of the Research Topic Exploring the Mechanical Properties of Soft Materials at Multiple Length Scales View all 5 articles

Physics-Based Computational Approaches to Compute the Viscoelasticity of Semiflexible Filamentous Biomaterials

L. G. Rizzi*

Departamento de Física, Universidade Federal de Viçosa (UFV), Viçosa, Brazil

This mini-review highlights recent advances on computational approaches that have been used in the characterisation of the viscoelastic response of semiflexible filamentous biomaterials. Special attention is given to the multiscale and coarse-grained approaches that might be used to model the mechanical properties of systems which involve biopolymer assemblies, for instance, actin, collagen, vimentin, microtubules, DNA, viruses, silk, amyloid fibrils, and other protein-based filaments. Besides the basic features of the most commonly used models for semiflexible filaments, I present a brief overview of the numerical approaches that can be used to extract the viscoelasticity of dilute and concentrated solutions, as well as systems with cross-linked networks. Selected examples of simulations that attempt to retrieve the complex shear moduli at experimentally relevant time and length scales, i.e., including not only the fully formed filaments and networks but also their self-assembly kinetics, are also considered.

1 Introduction

Solutions of semiflexible filaments formed from the self-assembly of biomolecules are ubiquitous in living organisms [1, 2]. Understanding how their viscoelastic properties emerge is crucial not only for a better comprehension about the transport and structural properties of fluids and hydrogels at the cellular level [3–6], but also because they seem to play a significant role in many disruptive processes, like cell invasion in several types of cancer [7] and protein aggregation in tens of proteinopathies [8].

Besides being a highly interdisciplinary problem [9], the characterisation of the viscoelastic response of self-assembled molecular systems involves time and length scales that span several orders of magnitude [10], with the typical building blocks at length scales of a few nanometers (10⁻⁹ m) forming structures of micrometers (10⁻⁶ m) to centimeters (10⁻² m) along time scales that range from nanoseconds (10⁻⁹s) to hours (10⁴ s). Experimentally, the mechanical properties of biomolecular systems and their building constituents have been probed at different scales mostly with aid of single-molecule [11] and microrheological techniques [12], and now, more than ever, multiscale and coarse-grained computational simulations [13] are becoming also a valuable tool in testing and validating modelling concepts in order to both understand and predict the viscoelastic behaviour of solutions of semiflexible filaments.

From the practical point-of-view, one aims to understand how the molecular information can be used to design the kind of response the biomaterial will display, e.g., liquid-like or solid-like [14], as well as to estimate the characteristic time scales that determine their viscoelastic behaviour. Figure 1 illustrates the typical viscoelastic responses that are obtained from rheology and microrheology experiments for three different types of solutions of semiflexible filaments. Liquid-like solutions, for instance, are primarily characterised by the value of their viscosity at low frequencies, i.e., η₀ = lim_ω→0η′(ω), with the frequency-dependent viscosity η′(ω) being associated to the loss modulus G′′(ω) as η′(ω) = G′′(ω)/ω. As shown in Figure 1A, the loss modulus displays a characteristic linear dependence on the frequency, i.e., G′′(ω) ≈ η₀ω, while the storage modulus displays a quadratic behaviour, i.e., G′(ω) ∝ ω², which are the expected low-frequency behaviours that one would obtain theoretically from the constitutive Maxwell and Rouse models [15]. Figures 1B,C illustrate the typical viscoelastic responses observed for solutions containing entangled and cross-linked semiflexible filaments, respectively. In both cases the solutions will display a semisolid/gel-like behaviour, and the interesting quantities are the entanglement modulus G_e, i.e., where G′(ω) displays a plateau-like regime (Figure 1B), and the low-frequency storage modulus G₀ = lim_ω→0G′(ω) (Figure 1C). In all cases one might want to predict both the exponents α and the corresponding frequency ranges of the power law regimes, i.e., where G′(ω) ∝ ω^α and/or G′′(ω) ∝ ω^α.

FIGURE 1

FIGURE 1. Typical viscoelastic responses characterised by the shear moduli G′(ω) and G′′(ω) for (A) complex fluid/liquid-like, (B) entangled, and (C) semisolid/gel-like solutions of semiflexible filaments. Inset drawings illustrate the configurations of the semiflexible filaments and their interactions in solution. Besides the identification of the distinct viscoelastic behaviours, the frequency-dependent shear moduli are often determined by experimentalists to provide a lower bound for the stiffness of the polymer biomaterials at different time scales [12].

In this mini-review I will focus mainly on simulations that have been used to study the aforementioned behaviours illustrated in Figure 1, including the effective modelling approaches of single semiflexible filaments, and the numerical methods used to describe entangled and cross-linked filament networks. Also, whenever it is pertinent, I will include information about the related self-assembly processes.

2 Modelling Approaches

2.1 Viscoelasticity and Relaxation Spectrum

As discussed in Ref. [16], polymer solutions are generally composed of structures that span several length scales so that they should contain many relaxation modes that can be characterised by a distribution of characteristic times τ, which is also known as the relaxation spectrum [17] H(τ). As extensively discussed in Ref. [17], one can relate the different mechanical responses illustrated in Figure 1 to different spectra H(τ) by assuming that the system is in the linear viscoelastic (LVE) response regime. Usually, the LVE regime is attained for small strains γ where the integrity of the biomolecular filaments (and, possibly, of the network structures) in the system is maintained during the whole measurement [12]. Although this might exclude drastic phenomena as those involved in shear banding and fracture experiments, situations that include the self-assembly of filaments (and networks) can be still studied if the time scales involved in the macroscopic reorganisation of the structures are greater than the time scales probed by the oscillatory or microrheology experiments. In practice, this means that the relaxation spectrum H(τ) or, equivalently, the stress relaxation modulus [17] G(t), does not change during the observation time [18], and one can evaluate the shear moduli in the LVE through a one-sided Fourier-Laplace transform [15, 17] as $G^{*} (ω) = i ω \int_{0}^{\infty} G (t^{'}) e^{- i ω t^{'}} d t^{'} = G^{'} (ω) + i G^{''} (ω)$ .

For solutions, one might consider to perform nonequilibrium simulations and implement shear flow conditions through driven, e.g., Lees-Edwards [19], boundary conditions and similar approaches (see Figure 2A). Alternatively, in order to obtain the response of the system in the LVE, one may avoid working with transient behaviours by considering simulations at equilibrium [20–22], and evaluate the relaxation modulus G(t) from the stress-stress autocorrelation (see Figure 2B). This kind of approach has been used not only to demonstrate the characteristic high-frequency regime where G′ ∝ G′′ ∝ ω^3/4 for dilute solutions of semiflexible filaments, but also in the evaluation of the plateau (i.e., entanglement) modulus G_e in entangled solutions [23–26]. Although in the highly entangled regime it is not always easy to determine the relevant characteristic length scales (i.e, the entanglement length L_e) and the interactions that significantly contribute to the stress tensor [27], an alternative modelling approach based on primitive path analysis [28] (see Figure 2B) has been successfully used to obtain values for G_e [29, 30].

FIGURE 2

FIGURE 2. Examples of computational particle-based methods used to estimate the mechanical properties of biopolymer materials. (A) Shear flow conditions can be implemented through Lees-Edwards [19] driven boundary conditions [40], and other related methods, e.g., the SLLOD algorithm [40, 41], and the reverse perturbation method [42, 43], in order to evaluate the viscosity η₀ at low shear rates $\dot{γ}$ , i.e., $η_{0} \approx \lim_{\dot{γ} \to 0} η (\dot{γ})$ . (B) Bead-spring and bead-rod coarse-grained filament models [44] used to simulate systems with a single or multiple chains. Full simulations at isothermal conditions can be carried out either through molecular dynamics (MD) simulations [45], where particles have thermostated degrees of freedom (e.g., peculiar momenta ${\vec{p}}_{i}$ ), or through stochastic simulations, which include mesoscale approaches like [40], e.g., multi-particle collision dynamics [46] (MPCD), dissipative particle dynamics (DPD), and Langevin/brownian dynamics (BD), where solvent effects are implemented through stochastic forces ${\vec{F}}_{s}$ . In any case the relaxation modulus for both dilute and entangled solutions can be evaluated as G(t) = ⟨σ_xy(t)σ_xy (0)⟩/k_BT, where σ_xy(t) denote entries of the stress tensor [15] which are given in terms of relative positions ${\vec{r}}_{i j} (t)$ and forces ${\vec{F}}_{i j} (t)$ between all particles in the system [47]. Alternatively, the entanglement modulus G_e can be effectively determined from the properties of primitive paths [28] defined by tubelike regions of length L_e that are obtained by considering that all the other chains are frozen [29]. (C) As discussed in Ref. [31], one can implement an oscillatory shear setup by imposing a time-dependent strain γ(t) = γ₀ sin (ωt), and, after a transient period, fit the numerically obtained stress σ(t) ≈ (1/V) (∂E/∂γ) to a sinusoidal function σ(t) = σ₀ sin (ωt + δ), which yields G′(ω) = (σ₀/γ₀) cos(δ) and G′′(ω) = (σ₀/γ₀) sin(δ). (D) Effective elastic network formed from self-assembled filaments with different thicknesses [38]. Displacements $\vec{u} (t)$ of the crosslinks in the steady-state regime are obtained from the overdamped equations so that the eigenvalues of the Hessian matrix H can be related to the relaxation times τ_i = 6πaη_s/λ_i. The different crosslinking connectivities z and the “three-particle” (ikj) bending terms in E_stretching lead to wide distributions of relaxation times H(τ) [34]. (E) By considering energy minimization schemes (see, e.g., Ref. [48] and references therein) one can obtain the low-frequency value of the storage modulus G₀ for either self-assembled [38] or arbitrarily defined networks [35–37]. Network configurations of athermal filaments [33] have been obtained in several ways [49], e.g., from the random deposition of line segments (i.e., Mikado networks), or by erasing a fraction p of bonds in a pre-established regular network. (F) In microrheology-based approaches the compliance J(t) of a material is related to the MSD ⟨Δr²(t)⟩ of probe particles [18] obtained either from passive [50] or active [39], i.e., externally driven, simulations. In both cases the complex modulus G*(ω) can be obtained from the Fourier-Laplace transform of a convolution integral involving G(t) and J(t) [12].

When the system present a percolating network one can, in principle, implement simulations based on oscillatory setups which are similar to experiments in rheology [14], and evaluate G*(ω) by direct means. For instance, in Ref. [31] the authors studied a network of cross-linked semiflexible filaments placed in a finite volume V by imposing a strain γ(t) = γ₀ sin (ωt) and measuring the shear stress as¹ $σ \approx (1 / V) (\partial E / \partial γ)$ , with E being the total elastic energy of the system (see Figure 2C). More recently, a numerical method to compute G*(ω) which avoids transient periods have been developed in Ref. [34] for systems with cross-linked filament networks. In this method one consider the overdamped equations for the displacements $\vec{u}$ of the crosslinks so that the nonzero eigenvalues λ_i of the Hessian matrix H, which elements are computed from the linearization of the elastic energy $E (\vec{u}) = E_{stretching} (\vec{u}) + E_{bending} (\vec{u})$ , can be associated to relaxation times τ_i (see Figure 2D), and the shear moduli G′(ω) and G′′(ω) can be determined from the measured strains at the boundaries of the system, which are given in terms of the eigenvectors of H. Alternatively, if one is interested in evaluating only the value G₀ of the storage modulus G′(ω) of cross-linked networks at the low-frequency regime, it might be useful to consider optimization schemes (see Figure 2E), e.g., the conjugate gradient algorithm [35, 36], to find the energy difference ΔE = E(γ) − E(0) so that [37, 38] G₀ = 2ΔE/Vγ².

It is worth mentioning that one can also consider the compliance function J(t) that is usually obtained from creep experiments to evaluate G*(ω), since G(t) is also related to J(t) through a convolution integral [12, 17]. In fact, since the compliance can be related to the mean-squared displacement (MSD) ⟨Δr²(t)⟩ of probe particles with radius a in d dimensions as [18] J(t) = (3πa/dk_BT)⟨Δr²(t)⟩, one can also use approaches based on microrheology (see Figure 2F). In some cases it might be even possible to speed up the simulations by considering active, i.e., externally driven, approaches that are based on fluctuation-dissipation relations [39], where the equilibrium fluctuations in the position of the probe particles can be estimated from their displacement Δz in the direction of the external force ${\vec{F}}_{ext} = F_{0} \hat{z}$ .

2.2 Coarse-Grained Models

Ideally, a full bottom-up modelling approach would have to incorporate all information about the molecular structures of the system, including not only the chemically specific features of the building blocks of the filaments but also additional solvent-specific details (Figure 2B). However, due to the intrinsic multiscale character of the viscoelastic behaviour, such atomistic-based approaches are only considered in a complementary manner, and mesoscopic (i.e., coarse-grained) modelling approaches are usually inevitable [26, 51, 52].

2.2.1 Self-Assembly of Filaments

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).

2.2.2 Models for Single Semiflexible Filaments

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) \sum_{j = 1}^{N - 1} {({\vec{r}}_{j + 1} - {\vec{r}}_{j})}^{2}$ , with κ being the elastic constant and ${\vec{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 \approx \sqrt{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} \sum_{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) \sum_{j = 2}^{N - 1} {({\vec{r}}_{j + 1} - 2 {\vec{r}}_{j} + {\vec{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_BT, 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].

2.3 Numerical Simulations

In the following I will describe additional approaches that are generally used in computational simulations, including a few selected examples that illustrate how the methods and models mentioned in the previous sections can be used to extract the viscoelastic responses of solutions like those displayed in Figure 1.

2.3.1 Dilute Solutions of Semiflexible Filaments

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 $η (\dot{γ})$ 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.

2.3.2 Solutions of Entangled Semiflexible Filaments

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).

2.3.3 Cross-Linked Networks

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.

3 Outlook and Challenges

There are still many challenges to the physics-based computational approaches involving multiscale simulations that attempt to evaluate the viscoelastic response of solutions of semiflexible filamentous biomaterials. Although generic coarse-grained polymer models have been developed to describe the self-assembly processes of filaments [54, 55], there are only a few computational studies on the association of fully formed semiflexible filaments [61, 76], indicating the feasibility of large scale simulations using mesoscopic models to compute the LVE response. Coarse-grained models seem to be unavoidable when performing simulations in those cases, and besides systematic coarsening modelling approaches [26, 51, 56], one could also explore simple heuristic models which take into account specific details of real biomolecules [77]. Additionally, one can consider the dynamics of probe particles obtained from simulations like the one presented in Ref. [78] to estimate the shear moduli from microrheological approaches [18]. While the entanglement modulus have been successfully determined from simulations of models based on primitive path analysis [29], it might be interesting to verify whether this and other approaches can be used to investigate issues related to dependence of G_e on the persistence length L_p for solutions in the tightly entangled regime [79]. As discussed in Ref. [27], it might be important to assess how significant are the correlations between different chains in entangled solutions, but only recently such large scale simulations have been reported for semiflexible filaments [80], even though their viscoelastic properties were not computed. Finally, it is worth mentioning that this mini-review focused on isotropic disordered biomaterials, but one can further explore solutions of semiflexible filaments which display nematic phases [81], and where anisotropic viscoelastic responses are expected. Also, it should be noted that, although configurations of filaments and their cross-linked networks have been mostly defined in an arbitrary manner (see, e.g., Figure 2E), experimentally-revelant mesoscopic information about semiflexible filamentous biomaterials are now becoming more available [82–84], and those may provide a strong driven-force in the implementation of novel physics-based computational simulations.

Author Contributions

The author confirms being the sole contributor of this work and has approved it for publication.

Funding

The author acknowledges the financial support of the Brazilian agencies CNPq (Grant No. 426570/2018-9 and 312999/2021-6) and FAPEMIG (Process APQ-02783-18).

Conflict of Interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

¹It is worth noting that, in general, there should be also an entropic contribution to the shear stress [32], i.e., −(T/V)(∂S/∂γ). Even so, semiflexible biopolymers are usually described as athermal filaments [33], which mean that their shapes are practically not affected by thermal fluctuations.

²Generally, Monte Carlo (MC) methods provide the most efficient ways to equilibrate complex polymer systems [67].

References

1. Burla F, Mulla Y, Vos BE, Aufderhorst-Roberts A, Koenderink GH. From Mechanical Resilience to Active Material Properties in Biopolymer Networks. Nat Rev Phys (2019) 1:249–63. doi:10.1038/s42254-019-0036-4