The molecular organization of soft biological networks, and notably the cell cytoskeleton, is at the origin of its rich mechanical behavior, that include nonlinear elasticity, visco-elasticity, activity and remodeling (Burla et al., 2019). Often these networks are made of a variety of filaments with varying persistence lengths, ranging from a few nanometers to several millimeters for microtubules (Gittes et al., 1993). A widely encountered topology consist of hybrid networks made of fairly stiff, or semi-flexible filaments, that are themselves cross-linked by flexible cross-linkers. For instance, in the cytoskeleton of most mammalian cells, stiff actin filaments are connected by flexible cross-linkers such as filamin, spectrin and other proteins that span a wide range in length-scales (Speicher and Marchesi, 1984;Wang et al., 1975). Similarly, in plant and fungal cells, the primary cell wall is a network of stiff cellulose (for plants) or chitin (for fungi) fibrils connected by flexible cross-linking proteins and polysaccharides (Geitmann and Ortega, 2009;Sridhar et al., 2018).
These hybrid networks differ from those made exclusively of flexible chains in that athermal or semi-flexible filaments possess an inherent orientation in their rest state. Thus, even if a network made by these filaments is originally isotropic, its mechanical response will eventually display strain stiffening due to filament alignment over large strains; a mechanism that has been identified as a key contributor to the nonlinear elasticity of semi-flexible networks (Storm et al., 2005). Furthermore, stiff filaments, when preferably aligned in a given direction, are the source of another type of nonlinearity, known as soft elasticity (Warner et al., 1994). The underlying mechanisms here relies on the effortless rotation of the filaments in the direction of principal stretch, with very little increase in elastic energy. At the macroscale, this is manifested by a regime of deformation where the network appears to have a negligible stiffness. While this effect has mostly been discussed for liquid crystal elastomers (LCEs) (Warner and Terentjev, 2007), discrete simulations of hybrid networks of actin filaments predict similar trends (Dalhaimer et al., 2007). A number of theoretical models have been introduced in the literature to describe these nonlinear elastic effects. Strain stiffening effect from filaments alignment was tackled by theoretical models in both the mechanics and physics communities. The mechanics community has mostly been interested in phenomenological models that could describe the nonlinear elasticity of fiber-reinforced biological materials accounting fiber realignment with strain (Gasser et al., 2006), or simply accounting for the presence of aligned fibers via anisotropic elasticity models (Guo et al., 2007). By contrast, the physics community concentrated on the relationship between network and mechanical response (Broedersz and MacKintosh, 2014), motivating models such as the hairy rod model (Kasza et al., 2010;Broedersz et al., 2008). Incorporating the effects of soft elasticity in continuum model has also been challenging as it involves the coupling between rod rotation and the elasticity of cross-linkers. Warner and Terentjev (2007), Biggins et al. (2008) presented such a model for nematic liquid crystal elastomers, via the introduction of an elastic potential expressed in terms of the so-called step length tensor, describing the rod's alignment.
Beyond these complex elastic responses, another important signature of semi-flexible biological networks is their ability to reconfigure over time. These materials are indeed known for their dynamic nature where cross-links are allowed to associate and dissociate over time. This is a prerequisite to their visco-elastic response, self-healing and remodeling over long time scales. A myriad of biopolymer, as well as novel synthetic polymers, such as vitrimers (Imbernon et al., 2016) and covalently adaptable networks (Kloxin and Bowman, 2013) have recently been fabricated with similar capabilities. This was followed by the synthesis of LCEs with dynamic bonds whose peculiar mechanical behavior was discussed in depth in a recent review (Saed et al., 2021). A particular advantage of introducing dynamic bonds in these materials is the possibility to program (and reprocess) rod alignment with mechanical stretch (Wang et al., 2017). From a theoretical standpoint, the time-dependent response of flexible networks was successful described by continuum formulations such as that introduced by Hui et coworkers (Long et al., 2014) and more recently by the transient network theory (Vernerey et al., 2017;Vernerey, 2018). Both frameworks are indeed able to bridge the molecular level (described in terms of the bond's exchange rates and the chain elasticity) and the macroscale, characterized by complex historydependent viscoelastic responses. However, when it comes to hybrid and semi-flexible networks, there is today very little theoretical understanding of the coupling between nonlinear elasticity and bond dynamics during mechanical deformation. A first attempt was described by Sridhar and Vernerey (2020) but concentrated on transient nematic networks, for which rods are strongly aligned and are assumed to rotate with the surrounding continuum. Thus, while this earlier formulation can describe a general class of anisotropic transient networks, it is unable to capture independent rod rotation, the associated soft elasticity, or the strain stiffening from rod realignment.
The objective of this work is therefore to expand on previous theories and introduce a mathematical framework to explain and predict the mechanical response of hybrid network made of elastic rods connected by flexible and dynamic cross-linkers (Fig. 1a). To ensure that the model is deeply rooted in structural mechanisms occurring at the network level, we take a statistical approach that is similar to that used in the transient network theory described in Vernerey et al. (2017). However, in the current work, the network statistics are described by two populations: the rods and the cross-linkers, that mechanically interact to yield the nonlinear effects discussed above. More specifically, the network is decomposed into elementary units (the hairy rods) whose mechanical behavior is captured by considering the non-affine kinematics of a rod deforming within a field of cross-linkers. This allows us to derive a reduced Fokker-Planck equation describing the elastic deformation of the rod and cross-linkers, accounting for the effects of bond dynamics. This model is then utilized to reconstruct a numerical model of the entire network, that is represented by a statistical distribution of hairy rods. We show that the model is not only able to capture the nonlinear elasticity of these hybrid networks (rod rotation, strain stiffening, and soft elasticity), but can also describe its complex interactions with bond dynamics.
This work focuses on percolated molecular network made of two types of filaments (Fig. 1): (a) semi-flexible (or athermal) filaments that are only mildly affected by thermal fluctuations; for simplicity, they are here described here as straight rods with non-zero rest length 𝓁 0 and (b) flexible cross-linkers characterized by their contour length 𝐿 = 𝑁𝑏, where 𝑁 and 𝑏 are the number and length of Kuhn segments, respectively. Due to their relatively small persistence length, these filaments take on average, a coiled conformation with mean end-to-end distance 𝑟 0 = 𝑏 √ 𝑁. It is assumed that many of these cross-linkers can bind to arbitrary locations along the rod as depicted in Fig. 1a. The bonds between the rod and the cross-linkers are assumed to be reversible, with association and dissociation rates 𝑘 𝑎 and 𝑘 𝑑 , respectively (Fig. 1d).
From a macroscopic (continuum) viewpoint, this class of networks may be described by a mean field theory, with descriptors including the concentration of the rods and cross-linkers as well as a statistical measure of their respective configuration, or end-toend vectors 𝒖 and 𝒓, respectively (Fig. 1c). We here assume that the concentration of cross-linkers greatly exceeds that of the rods, such that many cross-linkers can attach to a single rod. In this condition, it is possible to conveniently decompose the network into a collection of small units, the hairy rod, consisting of a single rod and its connected cross-linkers (Fig. 1b). The volume 𝑉 𝐻𝑅 of these units can be estimated by idealizing the hairy rod as a cylinder (Fig. 1b) of length 𝓁 0 +2𝑟 0 and radius 𝑟 0 , yielding 𝑉 𝐻𝑅 = 𝜋𝑟 2 0 (𝓁 0 +2𝑟 0 ). In the remainder of this work, the network is assumed incompressible so that the rod concentration 𝑐 𝑟 = 1∕𝑉 𝐻𝑅 remains constant and equal to:
Based on this hierarchical point of view, the full network can be statistically described by two main distributions.
• The full network is represented by the probability density function 𝑃 𝑢 (𝒖, 𝑡) of the rods' end-to-end vectors 𝒖 (Fig. 1a). The rods are here assumed to be large enough to not be affected by thermal noise. This means that any change in their conformation arises from their mechanical interactions with cross-linkers. In a stress-free configuration, one therefore expect rods to be in a stress-free conformation (i.e. their length is 𝓁 0 ) and a direction determined by the initial condition 𝑃 𝑢 (𝒖, 0) = 𝑃 0 𝑢 (𝒖) and the history of the network deformation (i.e. rod alignment during directional stretch).
• Zooming in to the level of a single rod, the hairy rod is described by the configuration of its connected cross-linkers as shown in Fig. 1b. Statistically, this is captured by (a) the concentration 𝑐 of connected cross-linkers and (b) a probability density function 𝑃 𝑟 (𝒓, 𝑡) of their end-to-end vectors 𝒓. A more thorough discussion of these quantities is provided next.
In what follows, we first concentrate on the mechanical response of a hairy rod with transiently bonded cross-linkers. Since a hairy rod has a well-defined orientation, it may be thought of as a network of perfectly aligned rods oriented along the vector 𝒖. Based on this formulation, a full network model can be constructed as a linear combination of hairy rods whose initial orientations follow the distribution 𝑃 𝑢 . This is presented in Section 5.
Within a hairy rod, the end-to-end vectors 𝒓 of connected cross-linkers' depends on two quantities: (a) the rod orientation 𝒖 = 𝓁𝒏 where 𝓁 is the length of the rod and 𝒏 is the unit vector pointing in the rod's direction and (b) the location 𝑦 of a chain measured from the center of the rod. This local coordinate is normalized with respect to the rod's length, such that 𝑦 ∈ [-1∕2, 1∕2] (Fig. 2). For the moment, we concentrate on a specific hairy rod for which the vector 𝒖 is known; the case of a random rod network characterized by a distribution of the vector 𝒖 is discussed in Section 5. A statistical description of the cross-linkers along the rod may thus be described by the distribution:
where 𝑐 is the concentration of connected cross-linkers per volume and the function 𝑃 𝑟 (𝒓, 𝑦) is the joint probability density function (pdf ) 𝑃 (𝒓|𝑦)𝑃 (𝑦), where the first term is the probability of finding a cross-linker with vector 𝒓 given its location 𝑦 along the rod, while the second term is the probability of finding a chain located at 𝑦. For simplicity, we assume here that the cross-linkers are uniformly distributed along the length of the rod such that 𝑃 (𝑦) = 1. Two remarks can be made regarding the distribution 𝜙(𝒓, 𝑦). First, it is a function of the concentration of connected crosslinkers, which may change over time due to the dynamic nature of the bonds. If one considers that the total concentration of these cross-linkers in solution is 𝑐 𝑡 , then, one can write:
where 𝑐 𝑑 (𝑡) is the concentration of cross-linkers that are detached (and that can potentially associate with the rod in time). Note that we explicitly expressed the dependency of 𝑐 and 𝑐 𝑡 on time to emphasize there transient nature. Furthermore, we assume here that flexible cross-linkers, which are strongly affected by thermal noise, tend to return to a quasi-isotropic distribution with a mean end-to-end distance 𝑟 0 = 𝑏 √ 𝑁 in a natural state (i.e. force free). The natural distribution 𝑃 0 𝑟 (𝒓, 𝑦) of cross-linkers may then be written as a Gaussian with variance 𝑁𝑏 2 ∕3 (Doi, 2013), i.e.:
Note that this expression implies that we are neglecting here the influence of the rod orientation on any bias in the orientation of the cross-linkers. More sophisticated models may be introduced to account for such interactions (see Sridhar and Vernerey, 2020 for instance), but would come at the cost of an increased model complexity.
From a statistical viewpoint, the conformation of the cross-linkers in a hairy rod is entirely described by the pdf 𝑃 𝑟 (𝒓, 𝑦, 𝑡). For the sake of practicality, it is appropriate to reduce our description by only considering moments of this distribution, rather than its full expression. To define these moments, we first introduce an averaging operation over the hair rod. For this, if we consider an arbitrary field ▵, that is an integrable function of the variables 𝒓 and 𝑦, its average ⟨▵⟩ over the hairy rod is expressed by:
Following and generalizing the transient network theory proposed in Vernerey et al. (2017), we introduce three distribution moments, that are expressed by the following second order tensors (denoted by conformation tensors in the remainder of the manuscript):
The first is the cross-linker conformation tensor which characterizes the average mean square stretch of the cross-linkers. The second is the rod conformation tensor, providing a tensorial representation of the rod direction and stretch. The third is the cross-linker/rod interaction tensor, whose physical interpretation becomes clear when the tensor is expressed in a local coordinate system aligned with the rod (Fig. 2). Let us therefore define such a coordinate system spanned by the unit basis vectors (𝒏, 𝒎, 𝒌) where 𝒏 is the unit vector along the rod direction, 𝒌 is normal to the page and 𝒎 is defined by 𝒎 = 𝒌 × 𝒏 = -𝝐 ⋅ 𝒏. Here, 𝝐 is the skew-symmetric
Levi-Civita tensor. Expressing the end-to-end vector 𝒓 in this coordinate system yields 𝒓 = 𝑟 ∥ 𝒏 + 𝑟 ⟂ 𝒎 where 𝑟 ∥ = 𝒓 ⋅ 𝒏 and 𝑟 ⟂ = 𝒓 ⋅ 𝒎.
One can then substitute this expression into (8) and use the fact that 𝒎 ⊗ 𝒏 = -𝝐 ⋅ 𝝂 to obtain:
where the parallel and perpendicular components of the cross-linker/rod interaction tensor are:
A graphical representation of these two components is given in Fig. 2b. It is clear here that 𝜂 ∥ is a measure of extensional forces applied to the rod by the cross-linkers while 𝜂 ⟂ measure the shear force applied in a perpendicular direction with respect to the rod, The latter is therefore prone to induce rotation of the rod in order to balance out the shear deformation of the cross-linkers.
The overall deformation history of an elementary volume (such as that describing the hairy rod) may be captured by the evolution of the deformation gradient 𝑭 (𝑡) in time. Alternatively, this history can be expressed in a rate form by the velocity gradient:
The use of the velocity gradient, rather than the deformation gradient is here justified by the fact that we consider transient networks, that are known to change their reference (or natural) configuration over time (Vernerey et al., 2017). In this section, we discuss the relationship between the overall deformation rate (described by 𝑳) and changes in the conformation tensors. This starts by the introduction of a Fokker-Planck equation that describes the evolution of the cross-linkers' configuration over time.
The evolution of the cross-linkers conformation over time is described by the rate of change of the distribution 𝜙. This rate is provided by the Fokker-Planck equation (Vernerey et al., 2017):
where 𝐷∕𝐷𝑡 is the material time derivative. We see here that the cross-linker conformation is dictated by three distinct ratedependent processes. The first contribution (first term on the right hand side) arises from their rate of elastic distortion ̇𝒓. The second term (second term on the right hand side) originates from the detachment of cross-linkers at rate 𝑘 𝑑 𝑐 in their stretched conformation (i.e. following the pdf 𝑃 𝑟 ). The last contribution (third term on the right hand side) describes the rate of crosslinker association at rate 𝑘 𝑎 𝑐 𝑑 from the pool of detached chains. Cross-linker association is here assumed to occur in the natural configuration (i.e. following the pdf 𝑃 0 𝑟 ). Following the decomposition (2) of the distribution, the above equation can be integrated over the hairy rod to derive an kinetic equation for the concentration 𝑐 (Vernerey et al., 2017):
Since 𝑘 𝑎 and 𝑘 𝑑 are here assumed to be independent of deformation, they remain constant throughout the network's loading history. This implies that the concentration 𝑐 quickly settles in a steady state 𝐷𝑐∕𝐷𝑡 = 0 that does not change over time. Setting (13) to zero, and invoking (3) yields the following relation between the concentrations of attached and detached cross-linkers:
Using this relation in ( 12) together with decomposition (2) yields a Fokker-Planck equation for the pdf 𝑃 𝑟 :
which can subsequently be used to derive evolution equations for the tensors 𝝁, 𝝆 and 𝜼. Before that can be done, however, we need to specify how the rate ̇𝒓 of cross-linker stretch relates to the overall deformation of the hairy rod. This is addressed next.
Let us now attempt to relate the local rate of deformation ̇𝒓 of the cross-linkers appearing in (15) to the velocity gradient 𝑳. For simple networks of flexible chains, this relation is usually stated by the affine assumption ̇𝒓 = 𝑳𝒓 (Vernerey et al., 2017). For the hybrid network considered here, however, such a relationship may be inappropriate since the flexible chains and the rods are usually characterized by very different mechanical properties. This means that the local deformation within the elementary unit will most likely be inhomogeneous. To capture this effect, a correction to the affine approximation can be adopted as follows. Let us first introduce 𝒂 as the vector that links the center of the rod to the end of one of its connected cross-linkers, i.e (Fig. 2a): where the vector 𝒖 points in the direction 𝒏 of the rod and has a length |𝒖| = 𝜆𝓁 0 where 𝓁 0 is the reference length of the rod (i.e. in the absence of applied force). The stretch ratio 𝜆 is defined by 𝜆 = 𝓁∕𝓁 0 , with 𝓁 the deformed length. We now postulate that the deformation is affine at the level of the elementary unit. This is expressed by the relation:
Using the definition ( 16), this yields an equation for the change in the cross-linkers' end-to-end vectors:
We see here that the rate ̇𝒓 is not only a function of the macroscopic velocity gradient and the rod's direction, but also a function of the change u in the rods direction and dimensions.
The change in rod orientation and stretch does not necessary follow the macroscopic deformation. We choose here to describe its kinematics with the micro-velocity gradient 𝑳 * defined by the relation:
This velocity gradient expresses the existence of a variety of deformation modes as depicted in Fig. 3. In particular, the effective rod rotation arises from two effects: the internal shear stress transmitted by the cross-linkers and the overall rotation 𝑾 * of the material point. We here postulate that the rotation of the frame is dictated by the macroscopic deformation, i.e, the angular rotation is directly related to the macroscopic velocity gradient 𝑳. The micro-velocity gradient can therefore be additively decomposed into a symmetric rate of micro-deformation tensor 𝑫 * and an skew-symmetric spin tensor 𝑾 as follows:
As discussed above, this strain rate is subjected to a number of restrictions, owing to the fact that rods can be considered as onedimensional elements, oriented along the direction 𝒏. To apply these conditions, it is first convenient to introduce the projection operators along and perpendicular to the rod, respectively:
With these operators, the axial strains rates along and perpendicular to the rod are:
Owing to its one-dimensional representation, we postulate that the rod can only increase in length, but not change its lateral dimensions; it can only undergo a longitudinal stretch rate ė, but no perpendicular stretch rate, i.e. ė⟂ = 0. Combining the two equations in (21), this condition can be rewritten in the simple form:
Using this relationship, the rate of micro-deformation can be decomposed into a stretch ( ė) and a traceless shear (or deviatoric) component 𝑫 * 𝑑 as:
It is now possible to completely specify the deformation of both rod and cross-linkers through the definition of two tensors: 𝑳 and 𝑫 * . The rod deformation is directly given by ( 19), while the cross-linker deformation is obtained by substituting ( 19) into (18), to yield:
where the rate of macro-deformation is defined by 𝑫 = 𝑳 -𝑾 . In the remainder of this work, we consider the full network to be incompressible, which finally imposes the condition:
Substituting the affine approximation (25) in Eq. ( 15) then leads to:
This constitutes the Fokker-Planck equation for the probability density 𝑃 𝑟 (𝒓). However, as discussed earlier, we here concentrate on the evolution of the moments of this distribution, rather than the distribution itself. We show in Appendices A.1 and A.2 that combining Eq. ( 27) and the definition of the tensors 𝝁, 𝜼 and 𝝆 in ( 6)-( 8), one can derive evolution equations for these tensors.
Accounting for the incompressibility condition 𝑇 𝑟(𝑳) = 0, the cross-linker conformation tensors obeys:
the rod conformation evolves according to
and cross-linker/rod interaction tensor is given by
where the geometric moment of rod relative to the cross-linkers is
These coupled ordinary differential equations are associated with initial conditions for 𝝁, 𝜼 and 𝝆. If we consider that the network is initially at rest, these take the simple form:
where 𝒏 0 is the unit vector that points in the initial rod orientation.
Note 1: In the limit when the rod reduces to a point (𝜒 → 0), the differential equation for 𝜼 (30) is decoupled from the other two equations. In this case, it can be shown that if 𝜼 obeys the initial condition (32), 𝜼 = 𝟎 at all times, regardless of the deformation history. In this case, the evolution equation for the cross-linker's conformation tensor becomes:
This model therefore degenerates to the standard transient network theory for networks of flexible chains (Vernerey et al., 2017).
Note 2: Another interesting limit of this model is given by the case of transiently cross-link nematic networks described by Sridhar and Vernerey in Sridhar and Vernerey (2020). In this model, the rod rotation is assumed rigid and follows the macroscopic motion. Its rate of deformation is therefore given by
In this situation, the rod rotation and stretch are fully determined by the macroscopic deformation and the only component of the cross-linker/rod interaction tensor is parallel to the rod, i.e. 𝜼 = 𝜂 ∥ 𝝂. It is then straightforward to show that the evolution Eqs. ( 28)-( 30) degenerate to the earlier model, albeit a change in notation (𝜂 ∥ is denoted as 𝜌 in Ref. Sridhar and Vernerey, 2020). The two models therefore converge when the principal strain directions coincide with those of the rod and the model can be thought of as a standard fiber-reinforced elastomer, with bond dynamics. When strains are not aligned with the rod, we will however see that the independent rotation of the rods yields nonlinear effects that can dominate the network's response.
To close the model, we now need to connect the concepts of deformation, conformation tensors and network stress. With this in mind, we introduce an energetic approach that is used for two purposes: (a) derive the governing equations for stress measures that are energy conjugates to the deformation tensors introduced previously and (b) express the Clausius-Duhem inequality, that is further used to derive a thermodynamically consistent constitutive relation relating stress to conformation tensors.
The proposed model enhances the kinetics of a classical solid by introducing the micro-deformation (rate) tensor 𝑫 * in addition to the velocity field 𝒗. Generally, one can construct a first-order theory, where the solid's internal energy is expressed in terms of the independent kinematic variables and their first gradient. Using the operators ∶ and ⋮ as the tensor's double contraction and triple contraction, respectively, one can therefore write the power P 𝑖𝑛𝑡 of internal force in a macroscopic domain 𝛺 as:
where the stresses 𝝈, 𝜷 and 𝝉 are energy conjugates of 𝑫, 𝑫 * and ∇𝑫 * , respectively. The tensor 𝝈 is to be interpreted as the Cauchy stress, 𝜷 as the rod stress, and the third-order tensor 𝝉 as a couple-stress tensor that originates from spatial gradients in rod deformation. Note that P 𝑖𝑛𝑡 does not depend on the velocity vector 𝒗 as required by the principle of frame indifference. Similarly, one can derive an expression for the power of external forces applied within and on the boundary 𝛤 of the body 𝛺:
We introduced here the body force density 𝒃 and couple-forces density 𝑩 acting throughout 𝛺. To account for surface forces, we also split the boundary as 𝛤 = 𝛺 𝑣 ∪ 𝛤 𝑡 where subset 𝛺 𝑣 is subjected to a fixed velocity and subset 𝛤 𝑡 is subjected to surface tractions 𝒕. A similar decomposition 𝛤 = 𝛺 𝐷 ∪ 𝛤 𝑇 is introduced for subsets where a micro-deformation rate D * or a surface couple-tractions 𝑻 are applied, respectively. The body force couple 𝑩 (and surface force couple 𝑻 ) are to be interpreted as force-couple originating from an applied external field (magnetic, for instance) that would induce a rod deformation illustrated in Fig. 3(b,c,d).
Using forms ( 35) and ( 36), one can invoke the principle of virtual power (i.e. following the procedure used by Germain in Germain, 1973), to derive the following governing equations for the three stress measures 𝝈, 𝜷 and 𝝉:
We here used the indicial notation for the sake of transparency. These coupled differential equations are subjected to boundary conditions:
In the present work, we focus on the mechanical response of a network that is subjected to a uniform deformation, i.e. the gradient ∇𝑫 * (and the associated couple stress 𝝉) identically vanishes at all times. For clarity, we also do not consider the effect of body forces. In this situation, the governing equations become:
𝜎 𝑗𝑖,𝑗 = 0 (39)
The first equation implies that the Cauchy stress is constant throughout the specimen, while the internal stress 𝜷 vanishes. Note that the latter stress may be thought of as an unbalanced stress on the rod, arising from the combined action of the cross-linkers and the rod deformation. When 𝜷 = 0, these forces are balanced and we reach internal equilibrium. Let us now explore how these stresses relate to the conformation tensors 𝝁, 𝜼 and 𝝆.
The mechanical response of a solid ultimately depends on its ability to store elastic energy; this is captured here by the introduction of an elastic potential 𝜓. In the hybrid networks considered in this work, elastic energy is stored in two ways (a) the entropic stretch of the cross-linkers, described by the conformation tensor 𝝁, and (b) and the deformation of the rod, captured by the tensor 𝝆. This allows us to postulate an elastic energy of the general form:
where the dependency of the energy on 𝝆 is only through its parallel, or stretch component (other components are associated with rod rotation). The definition of this potential must satisfy the second principle of thermodynamics. For simplicity here, we ignore here thermal effects by considering an adiabatic and isothermal process. In these conditions, the Clausius-Duhem inequality (expressing the second principle) may be simply expressed in terms of the rate of dissipation (per unit volume) = P 𝑖𝑛𝑡ψ ≥ 0, where the superimposed dot is used for the material time derivative. Using the internal energy in (35) and neglecting the rod deformation gradient for uniform deformations, we obtain:
The change in energy can be computed from (41) using the chain rule:
where the last terms 𝑝 𝑛 and 𝑝 𝑚 act as a Lagrange multipliers to enforce the inextensibility of the mesogen and the incompressibility of the network. Using ( 28) and ( 29) and substituting into (42) yields:
Since no dissipation arises from elastically deforming the cross-linkers or the rods, the first two terms must identically vanish. This allows us to express the stresses 𝝈 and 𝜷 in terms of the conformation tensors as follows:
The dissipation therefore takes the simple form:
This results implies that the rate of energy dissipation originates solely from bond dynamics due to its explicit dependency on the rate 𝑘 𝑑 . In other words, if cross-links are covalent, the rates of bond exchange vanish (𝑘 𝑑 = 𝑘 𝑎 = 0) and so does the dissipation. By contrast, increasing bond dynamics 𝑘 𝑑 renders the network more fluid-like and increases dissipation. We note here that the dissipation has a similar form as in dynamic network of flexible chains (Vernerey et al., 2017) and has to be interpreted as the elastic energy release rate (per unit volume) from bond dissociation.
In the remainder of this manuscript, we limit ourselves to the situation where the cross-linkers are modeled as Gaussian chains while rods are assumed to be made of a linear elastic material with Young's modulus 𝐸 𝑟𝑜𝑑 . We therefore postulate an elastic energy density of the form:
where 𝝂 = 𝒏⊗𝒏 with 𝒏 = 𝒖∕|𝒖| as the unit vector in the rod direction. The first term on the right-end side is the standard Neo-Hookean model for the cross-linkers, where the shear modulus 𝐺 = 𝑐𝑘𝑇 depends on the concentration 𝑐 of attached chains and the thermal energy 𝑘𝑇 (here 𝑘 is the Boltzmann constant and 𝑇 the absolute temperature). The last term captures the elasticity of the rods where the effective Young modulus 𝐸 depends on rod density 𝑐 𝐻𝑅 = 1∕𝑉 𝐻𝑅 , rod length 𝓁 0 , cross-sectional area 𝐴 𝑟𝑜𝑑 and stiffness 𝐸 𝑟𝑜𝑑 as:
Note that using the definition of 𝝆 in(7) this last term could also be written as 𝑇 𝑟(𝝆 -𝝂) = 𝜆 2 -1, where 𝜆 = 𝓁∕𝓁 0 is the rod stretch. Using energy density (46) in Eqs. ( 43) and (44) yields the following expression for the Cauchy stress 𝝈 and the rod stress 𝜷, respectively:
We have seen in ( 40) that for uniform deformations, 𝜷 = 𝟎. Using this equality in ( 49) and ( 48) quickly leads a simple expression for the Cauchy stress:
In the remainder of this work, we assume that a uniform deformation is applied to the network by enforcing the deformation history 𝑭 (𝑡) (or alternatively the velocity gradient 𝑳(𝑡)). The solution procedure to determine the rod deformation rate 𝑫 * , the conformation tensors and the stress is as follows. The first step is to express Eqs. ( 48) and ( 49) in a rate form so that the rates of 𝝁, 𝜼 and 𝝆 can be substituted from ( 28),( 29) and (30). Enforcing the fact that β = 𝟎 allows us to determine the rate of rod deformation 𝑫 * and estimate the rates μ, η and ρ. Using a classical integration scheme (forward Euler in the present work), one might then determine the tensors 𝝁, 𝜼 and 𝝆 at all times and thus explicitly determine the stress tensor 𝝈 using either (48) or (50). We note that the pressure 𝑝 𝑛 and 𝑝 𝑚 are determined by enforcing that the overall deformation is incompressible, i.e. 𝑡𝑟(𝑫) = 0 and that the rod deformation verifies 𝑫 * ∶ (𝑰 -𝝂) = 0. We next use this numerical approach to explore the model prediction for the mechanical response of a single hairy rod subjected to tensile loading. We first explore its elastic response by setting 𝑘 𝑑 = 0 and extend our analysis to its viscoelastic response (i.e. 𝑘 𝑑 > 0). 𝜆. The model predicts a soft elasticity regime, with a stress plateau followed by a stiffening phase when the rod is not initially aligned with the stretch direction.
The response is however monotonic when the rod is initially aligned with the stretch directions. These results are for a ratio 𝓁 0 ∕𝑟 0 = 20 and 𝐸∕𝐺 = 400.
Response of an elastic hairy rod in tension for different initial rod orientations 𝜃 0 , with respect to the loading direction. The uniaxial stress-stretch relationship and the evolution in rod angle are shown. We identify two regimes, depending on 𝜃 0 : a regime of soft elasticity when the angle 𝜃 0 is above 𝜋∕4 and a regime of pure strain stiffening when the rod is more aligned with the deformation. The transition between these regimes is rather smooth. These results are for a ratio 𝓁 0 ∕𝑟 0 = 20 and 𝐸∕𝐺 = 400.
We start by exploring the mechanical behavior of a single hairy rod with permanent cross-linkers (𝑘 𝑑 = 0). This is equivalent to studying a nematic network, where rods are originally uniformly oriented in a single direction 𝒏 0 , and are subjected to a uniform deformation field. For clarity, we further focus on the response of the network in uniaxial tension, and when the rod stiffness remains significantly larger than the stiffness of the polymer chains. We therefore consider the case where 𝐸∕𝐺 = 400, and verify that increasing this ratio had very little impact on the presented results. In the remainder of the manuscript, all stresses are normalized by the shear modulus 𝐺 of an equivalent network of flexible chains.
The elastic response of a hairy rod is controlled by two key parameters: (a) the relative rod length defined by 𝓁 0 ∕𝑟 0 and (b) the initial orientation 𝜃 0 of the rod relative to the stretch axis. Fig. 4 shows the response of the hairy rod, along with a graphical illustration of the rod and cross-linker deformation as it is subjected to a stretch ratio 𝜆 𝑥 = 2.2 in an orthonormal coordinate system with unit vectors { 𝒆 𝑥 , 𝒆 𝑦 } as shown. We focus here on cases where the rod is (a) strongly misaligned (𝜃 0 = 0.35𝜋) and (b) when it is fully aligned (𝜃 0 = 0) with the stretch direction. In the first situation, the model predicts a soft response distinguished by the existence of a stress plateau of zero stiffness before the network stiffens. This soft plateau is associated with a regime where rod rotation is dominant, while network stiffening starts once the rod has realigned with the stretch direction. We note here that the soft elasticity regime is usually associated with a non-convex strain energy function. When assessing the elastic energy stored in the cross-linker 𝜓 𝑐 = 𝐺∕2[𝑇 𝑟(𝝁) -3] as a function of axial strain 𝜆 -1 (Fig. 6c andd), we indeed observe the loss of convexity and how it is associated with the non-monotony of the stress-strain response. Interestingly, this feature appears despite the apparent convexity of the energy function of Eq. ( 46). This can be explained by that the energy is written as a function of the conformation tensor 𝝁 (rather than the strain) which is itself a non-monotonic function of strain as dictated by the rod/cross-linker interactions captured by evolution Eqs. ( 28)-( 30). In other words, the nonconvexity of the elastic energy naturally emerges from internal interactions, and does not need to be specified as such. This non-convexity is known to be at the origin of strain localization and pattern formation which includes the formation of stripe domains in LCEs subjected to uniaxial stretch (Warner and Terentjev, 2007). When the rod is initially aligned with stretch however, the rod does not need to realign and the mechanical response exhibits a monotonic (mildly stiffening) stress-strain curve that usually characterizes flexible networks. We show later (Fig. 6) that the response is in fact stiffer than that of an isotropic network due to the presence of the stiff rod, frustrating the chain deformation. A notable advantage of the proposed model is its ability to follow the deformation of the rod and the cross-linkers individually. This is illustrated in Fig. 4 by a rod vector whose length and direction represent the rod stretch 𝜆 and orientation, respectively. Similarly, the cross-linker conformation tensor 𝝁 can be represented by an ellipse from which the lengths and directions of the principal axes indicate the eigenvalues and eigenvector of 𝝁, respectively (Vernerey et al., 2017). It is therefore possible to visualize the large stretch experienced by cross-linkers during deformation, despite the very small rod deformation.
The soft elastic response observed above is typically observed in liquid crystal elastomers (LCEs) in a nematic phase, and is usually associated with mechanical instability (Warner and Terentjev, 2007). To better understand how microstructure affects this response, let us now explore the effects of initial rod orientation on the hairy rod behavior. In this context, Fig. 5 depicts the stress-strain response and rod rotation for a range of initial rod orientation ranging from 𝜃 0 = 0.45𝜋 (nearly perpendicular to the stretch direction) to 𝜃 0 = 0 (aligned with stretch). We observe here a smooth transition from a regime of pronounced soft elasticity, to a regime of pure strain stiffening as 𝜃 0 ≈ 0.3𝜋. In the soft elasticity regime, we see a clear correlation between the range of deformation associated with rod rotation (leading to the stress plateau) and the range of deformation when the rod becomes aligned with stretch (leading to a monotonic increase of stress). Thus, when 𝜃 0 decreases, the effect of rod rotation is weaker, leading to the disappearance of the stress-plateau, and the emergence of a purely stress-stiffening response. The relative length 𝓁 0 ∕𝑟 0 of the rod also has a noticeable influence in mediating these responses. To illustrate this, Fig. 6 depicts again the cases of an initially misaligned rod (𝜃 0 = 0.35𝜋) and an exactly aligned rod (𝜃 0 = 0) but now shows the response when the rod length ranges from 0 to 30𝑟 0 . We first observe that when 𝓁 0 = 0, the network displays the response of a Neo-Hookean solid. Indeed, in this case, the evolution of the conformation tensor is given by (33) (and subjected here to 𝑘 𝑑 = 0). In this case, it was shown in Vernerey et al. (2017) that the conformation tensor becomes 𝝁 = 𝑭 𝑭 𝑇 and the free energy (46) degenerates to that associated with a traditional incompressible Neo-Hookean solid. Results with a finite rod length suggest that 𝓁 0 plays a role in mediating two important characteristics of the hairy rod: (a) First, increasing rod length leads to a soft elastic behavior that becomes more and more severe (Fig. 6a). Indeed, not only does the plateau becomes more apparent with increasing 𝓁 0 , but we observe a temporary softening response when 𝓁 0 > 20𝑟 0 , indicating an unstable behavior where the rod would snap between two solutions during stretch. (b) Second, when the stretch is aligned with the rod, we observe that a longer rod yields a stiffer axial response, a response that is characteristic to most fiber-reinforced polymers (Fig. 6b). In fact, it can be shown (see Sridhar and Vernerey, 2020) that for these loading conditions, the model degenerates exactly to the standard formulation for a fiber-reinforced Neo-Hookean rubber, developed by Guo et al. (2007) and Spencer et al. (1984). b) a strain stiffening regime followed by a stress relaxation when 𝜃 0 = 0. These results are for a ratio 𝐸∕𝐺 = 400. Fig. 8. Effect of bond dynamics on the loading and relaxation response of a hairy rod that is not aligned with stretch direction (𝜃 0 = 0.35𝜋). The stretch history is depicted in the inset of subfigure (a). When 𝑘 𝑑 is relatively low (a), we still observe a soft elastic regime, and the reorientation of the rod with deformation is not affected by bond dynamics. When 𝑘 𝑑 becomes large (b), the deformation is dominated by creep, leading to a very weak stress response and a dependence of rod rotation with bond dynamics. For extremely slow bond dynamics (𝑘 𝑑 = 1000 𝑠 -1 ), the rod does realign with the applied stretch. These results are for a ratio 𝐸∕𝐺 = 400.
Let us now explore the effect of bond dynamics (characterized by a finite value of 𝑘 𝑑 ) on the hairy rod behavior. In this case, we use the Weissenberg number 𝑊 = |𝑳|∕𝑘 𝑑 to monitor the competition between the rate of bond exchange 𝑘 𝑑 and the loading rate |𝑳| (where | ⋅ | is taken as the 𝐿 2 norm of a tensor). In the analysis below, we choose to fix the loading rate as ε𝑥 = 1, such that 𝑊 = 1∕𝑘 𝑑 in all presented results.
To illustrate the response predicted by the transient hairy rod model, we consider a loading path that consists of two stages: (a) loading at constant strain rate (with 𝑘 𝑑 = 1) up to a stretch of 1.5, and (b) holding at constant strain for a long period of time to allow stress relaxation. Again, model predictions are shown for the cases of a misaligned rod (𝜃 = 0.35𝜋) and an aligned rod (𝜃 = 0). In the first situation (Fig. 7a), we observe that despite the moderate bond dynamics, the hairy rod still exhibits a soft elasticity regime in the loading stage when the rod is misaligned with the stretch direction. However, by contrast to the elastic network, we predict a relaxation of the stress (until is eventually vanishes at infinite times) during the holding stage. The origin of this behavior can readily be visualized by the elliptical representation of the cross-linker conformation tensor 𝝁, that slowly returns to an circular shape over time. This occurs as the most stretched cross-linkers detach from the rod and reattach in an isotropic fashion and in a state close to their stress-free state. Noticeably, while the rod slowly returns to its stress-free length during this process, it does not return to its initial orientation. In other words, the model predicts a permanent realignment of the rod in the maximum stretch directions. A very similar response is observed when the rod is initially aligned with stretch 7b), but with no realignment. In this loading scenario, the response is that of a viscoelastic fiber-reinforced elastomer (Sridhar and Vernerey, 2020). Now that the effects of bond dynamics are identified, Fig. 8 illustrates its effect on the general stress-strain signature of the hairy rod, and particularly, on its ability to realign rods along the stretch direction. We shows here the evolution of stress and rod orientation over time (starting with a rod orientation 𝜃 0 = 0.35𝜋), for a network with (a) moderate bond dynamics (Fig. 8a) and (b) fast bond dynamics (Fig. 8b). In the case of moderate bond dynamics (0 ≤ 𝑘 𝑑 ≤ 20), the hairy rod retains some of its elasticity (and particularly its soft elastic behavior). We note here that while an increase in 𝑘 𝑑 decreases the stress sustained by the network, it does not significantly affect the rotation of the rod over time. This may be explained by the fact that during loading, most cross-linkers F.J. Vernerey Fig. 9. Deformation and tensile response of an isotropic network. (a) Evolution of the rod stretch and orientation (angular distribution plots) and network order (curve of 𝑄 versus 𝜆 -1). Corresponding uniaxial stress-strain response for three different rod length 𝓁 0 = 0, 10𝑟 0 and 20𝑟 0 (the stress is normalized by the shear modulus 𝐺 = 𝑐𝑘𝑇 ). Compared to a standard network (𝓁 0 = 0), we see that the existence of stiff rod and their alignment with stretch induces a significant stiffening of the network response. We note that within this range of parameters, the order parameter is insensitive to rod length. These results are for a ratio 𝐸∕𝐺 = 400.
remain attached to the rod, and may thus efficiently apply the internal couple moments necessary for its rotation. This is in contrast to the case of fast bond dynamics, where 𝑘 𝑑 is several orders of magnitude larger than the strain rate. In this case, the network effectively behaves like a viscous fluid which produces very small stresses. While this was expected from previous work on transient networks, we further observe that fast bond dynamics may decrease and even completely suppress rod rotation during loading (Fig. 8b). In other words, the model predicts that such a dynamic network can be subjected to very large deformation without rod reorientation, provided that the loading is slow enough.
Now that the response of a single hairy rod unit is characterized, we turn to the length-scale of a full network depicted in Fig. 1a. At this level, the network is described by the rod pdf 𝑃 𝑢 (𝒖, 𝑡), that evolves over time from an initial state characterized by
Since the evolution of a hairy rod conformation 𝒖 is entirely described by the trio {𝝁, 𝝆, 𝜼}, the full network is here modeled by the discrete 𝑝𝑑𝑓 𝑃 𝑢 (𝒖 𝑖 ) where the vector 𝒖 𝑖 are 𝑛 vectors (1 ≤ 𝑖 ≤ 𝑛) that span the space of rods' end-to-end vectors. Consider an initial network characterized by the 𝑝𝑑𝑓 𝑃 𝑢 (𝒖 0 𝑖 ) for which the collection of 𝑁 vectors 𝒖 0 𝐼 (𝐼 = 1...𝑁) span the space of initial rod's end-to-end vectors, each of them being assigned the probability 𝑃 𝑢 (𝒖 0 𝑖 ). This discrete distribution therefore verifies:
where we used the notation 𝑃 𝐼 = 𝑃 𝑢 (𝒖 0 𝐼 ). Upon deformation, each original vector 𝒖 0 𝐼 is transformed into a deformed vector 𝒖 𝐼 as dictated by the hairy rod model (and particularly Eqs. ( 28)-( 30)). For each hairy rod, one can also compute the stress vector 𝝈 𝐼 , such that the product 𝑃 𝐼 𝝈 𝐼 may be thought of as the partial stress associated with a rod that was originally defined by the vector 𝒖 0 𝐼 . Let us now consider a arbitrary measures 𝛥 𝐼 that can be defined for the 𝐼 𝑡 ℎ hairy rod. The statistical average ⟨𝛥⟩ of this field over the entire network can then be defined as the weighted sum:
In what follows, we are particularly interested in the network stress and the (two-dimension) rods' order parameter defined by:
The traceless order parameter tensor maybe better represented in its principal directions as 𝑸 = 𝑄 𝑑𝑖𝑎𝑔(1, -1) where 𝑄 is a scalar order parameter. With this definition, 𝑄 = 1 corresponds to a perfect nematic order in the 𝑥-direction, 𝑄 = -1 to a perfect nematic order in the 𝑦-direction, while 𝑄 = 0 represents a fully isotropic network.
Before exploring the role of dynamic bonds on anisotropic networks, we first give in Fig. 9 a glimpse of the model prediction for an elastic network. For convenience, let us consider an isotropic network of equal length rod. In this case, the rod orientation can be represented by a uniform angular distribution shown in Fig. 9a yielding an initial order 𝑄 = 0 in its undeformed state. Again, all results in this section are shown for a ratio 𝐸∕𝐺 = 400. We now examine the stress-strain curve and change in network orientation following the application of a uniaxial stretch in the horizontal direction, up to a stretch ratio 𝜆 ≈ 4. For this, each direction of the network corresponds to a representative hairy rod whose evolution with deformation is described in Section 4.2. As the deformation is applied, one can then follow the evolution in network anisotropy as represented by the change in the rod's direction and the order parameter 𝑄 (Fig. 9a). We observe a strong rod alignment with stretch, until the order parameter asymptotes to 1(that would correspond to a fully aligned network).
To visualize the effect of rod alignment on the network response, Fig. 9b presents its uniaxial stress-strain response for three different initial rod lengths ranging from 𝓁 0 = 1 to 𝓁 0 = 20𝑟 0 . We find that while alignment is not sensitive to rod length, the stress increases quadratically with 𝓁 0 ; this is indeed a direct consequence of the dependency of the effective rod modulus on length as expressed in (47). We finally note that for this elastic network, the deformation is entirely reversible and the loading-unloading curves are identical.
The mechanical response of a dynamic network 𝑘 𝑑 > 0 is now examined by subjecting it to various values of constant true strain rate (represented by their normalized value 𝑊 = ε∕𝑘 𝑑 ) where 𝑘 𝑑 remains constant for all cases. Again, the deformation history is split into two regimes:
(a) A creep-like regime in which the stress is highly strain-rate dependent. The Weissenberg number measures the ratio of elastic to viscous deformation (Vernerey et al., 2017): a large Weissenberg number (𝑊 ≫ 0.5) yields an elastic dominant response, characterized by a measurable stiffness and strain stiffening at large strains (see the curve corresponding to 𝑊 = 0.5 in Fig. 10a). By contrast, at small Weissenberg numbers, the network dynamics occurs faster than the loading rate and the network behaves more like a Maxwell fluid. In this situation, the mechanical response is soft and reminiscent of a creeping flow (see the curves associated with 𝑊 = 0.2 and 𝑊 = 0.1 in Fig. 10a). Despite these large differences in mechanical response, the model predicts that rod alignment is nearly insensitive to loading rates in these regimes (Fig. 10c).
(b) A relaxation regime where the stress exponentially decays to zero at a characteristic time 1∕𝑘 𝑑 . While this exponential stress decay is expected for dynamic networks with constant 𝑘 𝑑 (Vernerey et al., 2017), a notable observation here is that rods do not revert back to their original orientation during relaxation, but instead remain aligned with the previous deformation history. The rod stretch 𝜆 however relaxes back to its original value of 1 over time as can be observed in the elastic deformation of the rods in Fig. 10c. The irreversibility in rod orientation is a particularly interesting feature of these dynamic networks, that could be used to program the anisotropy of a network via prestretch and relaxation. To further explore the required condition for network alignment through stretching, we next decrease the relative loading rate 𝑊 to extremely low values. For the sake of illustration, Fig. 11 Fig. 12. Illustration the response of an initially isotropic network to a loading path history that is composed of four phases (a) a stretch at constant true strain rate (𝑊 = 0.2) along the 𝑥-direction for a time 𝑡 = 1.5𝑘 𝑑 , (b) and relaxation state until the stress fully vanishes, (c) a stretch at constant true strain rate (𝑊 = 0.2) along the 𝑥 ′ -direction for a time 𝑡 = 1.5𝑘 𝑑 , followed by another stress relaxation at constant strain. We show here the normalized tensile stress (measured in units of 𝐺) along the principal stretch directions. The first phase shows a typical elastic-creep response with rod realignment, while the second phase shows a general realignment of the anisotropic network, with a soft elastic regime occurring early during the deformation history.
presents relaxation results after loading an isotropic network at Weissenberg number 𝑊 = 0.01, 0.001 and 0.0001, respectively. The change in order parameter with time and the final rod distribution long after relaxation are presented in Fig. 11a. These results suggest that network reprogramming can only be achieved if the loading rate is fast enough. Indeed, at low Weissenberg numbers (i.e. 𝑊 = 0.0001), the loading rate is so low that cross-linkers dissociate before they can elastically interact with the rods. In this scenario, rods behave as if they were surrounded by a viscous fluid that does not generate enough of a couple moment to trigger their rotation. In this regime, the stress response can be characterized by that of creeping flow (Fig. 11b), where the stress is only a function of strain rate, and thus remains constant throughout the loading history.
To finally illustrate the importance of loading history in these hybrid dynamic networks, we here examine a loading path that consists of two main stages. Starting from a fully isotropic network, we first consider a programming phase where the network is stretched relatively quickly (𝑊 = 0.2) in a direction shown by the 𝑥-axis in Fig. 12. In these conditions, the rods strongly align with the stretch, such that after imposing a relaxation stage (i.e. keeping a fixed deformation for a long enough time), the rods remain ordered in a stress-free configuration. Following this stage, the network becomes nearly nematic, and thus exhibits anisotropic properties as discussed in Sridhar and Vernerey (2020). If the network represents a LCE, this step can also be used to uniformly align mesogens across a macroscopic sample (i.e. obtain a monodomain) and overcome the traditional limitations that arise with multi-domains LCEs (Saed et al., 2021). In this example, we particularly emphasize on the soft elastic regime that can originate from these aligned networks. For this, we consider a second deformation phase where a purely tensile strain is applied along the 𝑥 ′ -axis, oriented at angle of 𝜋∕3 with respect to the original direction -𝑥. In these conditions, the model predicts a general realignment of the rod along the new principal directions. However, because the loading rate is large enough (𝑊 = 0.2), the stress-strain response displays a soft elastic regime that is similar to that previously examined in Fig. 4a. Generally, this example shows that a variety of nonlinear effects can result from the combined structural rearrangement of the network and its inherent and complex elastic response.
To summarize, this paper presents a physically-derived model of the time-dependent response of hybrid networks made of athermal (or semi-flexible) rods dynamically linked by flexible cross-linkers. Stemming from statistical mechanics, the proposed model directly relates the geometrical and physical description of the network -expressed by the concentration, length, and elasticity of the cross-linkers and rods as well as the orientational rod distribution and the bond dynamics -to its complex emerging response. This response includes a combination of nonlinear elasticity and viscoelasticity, along with the network's structural and orientational evolution. In a nutshell, the theory rests on three main pillars:
• A set of governing equations, expressed by ( 37) and ( 38), that enforce mechanical equilibrium. While the current work concentrates on the case of uniform deformations, without strain gradient, the theory allows for the incorporation of gradients in rod deformation, and thus incorporate size-effects as commonly observed in liquid crystals (i.e. Frank elasticity arising from distortions in mesogen alignment Frank, 1958). In this case, the theory takes a form that may be compared with the micromorphic theory (Germain, 1973;Vernerey et al., 2007).
• The definition of a Helmholtz free energy functional 𝜓, expressed in terms of the three conformation tensors 𝝁, 𝜼 and 𝝆 (Eq. ( 41) in this work). This may generally be achieved by recognizing that these tensors represent the average elastic deformation of the cross-linkers (𝝁), the elastic interaction forces between rods and cross-linkers (𝜼), and the average elastic deformation of the rod (𝝆). The stresses may then be derived as gradients of the potential 𝜓 with respect to the conformation tensors as expressed in ( 43) and ( 44).
• The introduction of kinetic Eqs. ( 28), ( 29), (30) that describe the evolution of the three conformation tensors in terms of the velocity gradient 𝑳, the rate of rod deformation 𝑫 * , and bond kinetics (𝑘 𝑑 in this work). The form of these equations originates from a Fokker-Planck equation for the rod and cross-linker distributions and the non-affine kinematics of the hairy rod discussed in Section 3.
Thus, in general, the network elasticity is captured by the definition of the free energy 𝜓, while its viscous response is captured by bond dynamics appearing in the evolution equations for 𝝁, 𝜼 and 𝝆. We have shown here that this model naturally captures important signatures of the nonlinear response of reversibly cross-linked filamentous networks including:
• The nonlinear, strain stiffening behavior of an isotropic network due to filament realignment in the direction of stretch. This effect is commonly observed in most biopolymers and semi-flexible networks such as the actin cytoskeleton (Erk et al., 2010) or collagen gels (Vader et al., 2009;Motte and Kaufman, 2013).
• The soft elasticity plateau of a network of aligned filaments (nematic) when stretched perpendicular to the average filament direction. This phenomenon occurs when filament size is on the same order of magnitude as cross-linkers. Most experimental observations of soft elasticity have therefore pertain to LCE networks (Herbert et al., 2021), but have also been predicted by numerical simulations of networks of short actin fibers (Dalhaimer et al., 2007).
• The creep and stress-relaxation response of the network due to bond dynamics. Depending on the relative length of the filament, we recover the nonlinear visco-elasticity of polymers with transient bonds when 𝓁 ≪ 𝑟 (Hui et al., 2021;Winne et al., 2019) and the behavior of dynamic semi-flexible polymers, including the creep and stress relaxation when 𝓁 ≫ 𝑟, Fernandez-Castano Romera et al. (2018). When fiber are aligned, fiber rotation follows that of the continuum and the model describes the response of anisotropic polymers (Sridhar and Vernerey, 2020) ranging from anisotropic fluids (Ericksen, 1960) to fiber-reinforced rubbers (Spencer et al., 1984).
• The combined nonlinear elastic and time-dependent response of semi-flexible networks and LCEs when 𝓁 ≈ 𝑟. In this case, we have shown that the model can capture the interplay between soft-elasticity, stress relaxation and energy dissipation observed in these materials (Saed et al., 2020). We also illustrated the possibility of programming the network's alignment into a monodomain by purposefully subjecting it to a series of stress states over time as experimentally observed in Pei et al. (2014).
This theory represents an extension of the standard transient network theory to describe a variety of networks, that includes fibrous connective tissues (Pritchard et al., 2014), the cell cytoskeleton (Fletcher and Mullins, 2010), the cell cortex (Dalhaimer et al., 2007), the wall of plant/fungi cells (Sridhar et al., 2018), or LCEs (Warner and Terentjev, 2007). The goal of this work was however to lay out the fundamental structure of the theory without incorporating some complex and specific phenomena occurring in the above network. The presented model therefore suffers from limitations, that could be addressed in future work. Among them is (a) the consideration of non-uniform deformations and size effects. This can be addressed by adding a dependency of the free energy 𝜓 on the gradient of rod deformation, similarly to the Frank elasticity terms in the theory of LCEs (Frank, 1958). (b) The nonlinear elasticity of cross-linkers and rods, including rupture. Network of semi-flexible rods may incorporate the nonlinear elastic response of the rod, as well as the stiffening of the cross-linkers as they reach their contour length. These effects can potentially be incorporated into the elastic potential 𝜓 following a similar approach as in Vernerey (2018). (c) Thermal effects on rod entropy. In this work, rods were assumed to be athermal, allowing us to neglect their entropy. If these effects are included, an increase in temperature could induce a gain in the rod's entropy, forcing them to lose their order and regain a random configuration. The incorporation of this effect would thus be essential in understanding the effect of bond dynamics on the isotropic-nematic phase transition in LCEs. We note here that the extension of this theory to LCEs will also necessitate to revisit the kinematic assumptions for different network topologies (i.e. main chain, side chain, . . . ). (d) The nonlinearity of bond dynamics. Most reversible bonds are known to decrease their lifetime with applied force, as predicted by Bell's law (Bell, 1978). This can induce a shear-thinning response of the network, or vice-versa, a shear thickening in the case of catch bonds. Such effects can be induced by explicitly writing the rate 𝑘 𝑑 as a function of cross-linker stress (Vernerey et al., 2018). (e) Finally, the dynamic networks studied here typically enable a Rouse-type diffusion of rods and filaments over time, a process that could play an important role in the self-healing capacity of the network (Stukalin et al., 2013). Recent work on the diffusion of rods enabled by dynamic cross-links (Lalitha Sridhar et al., 2021;Koo et al., 2021) show that diffusion is favored by a smaller number of cross-linkers per rod. Such feature may be incorporated into the present theory to describe the mechanics of self-healing.
Franck J. Vernerey: Conceptualization, Methodology, Visualization, Investigation, Writing -review & editing.
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.