Thermodynamic modeling of a phase transformation in protein filaments with mechanical function

Protective egg capsules from whelks (intertidal marine gastropods) were recently demonstrated to derive their impressive mechanical behavior—reminiscent of the pseudoelastic behavior in some alloy systems from a reversible phase transition of component protein building blocks from a compact α-helical conformation to a more extended softer conformation, called β*. This behavior was previously modeled under equilibrium conditions, demonstrating that the transition from the α- to β*-phase could account for the pronounced yield plateau and reversibility; however, a theoretical understanding of the non-equilibrium behaviors of the egg capsule (e.g. strain rate dependence and hysteresis) requires a new approach. Here, we modify the previously proposed model in order to address the time-dependent behaviors of the whelk egg capsule biopolymer. Our results indicate that hysteresis during cyclic loading originates from a mismatch between the speed of the mechanical driving force and the rate at which the phase transition occurs. Furthermore, the characteristic curved shape of the stress–strain plot arises from a nonlinear relationship between the transformation rate and the amount of applied load. These results have important implications for our understanding of the mechanics of biological polymers and may have implications for the design of biomimetic pseudoelastic polymers.


Introduction
The forces experienced by organisms living on the rocky seashore among crashing waves have been compared in the literature to being strapped to the nose of a supersonic jet [1]. Therefore, it is not altogether surprising that species inhabiting the intertidal zone dedicate significant resources to strategies and materials for protecting soft and vulnerable tissues. Along these lines, many marine organisms construct hard and tough mineralized parts such as the nacreous shells of abalone or the carapace of crabs; however, this strategy is not ideal for all applications, for example, adhesive attachment or encapsulation of embryos. An alternate strategy is to create softer organic structures capable of dissipating large amounts of mechanical energy via reversible deformation (e.g. the byssal threads of marine mussels [2] or various marine egg capsules [3,4]). The egg capsules of marine whelks are particularly interesting due to their characteristic pseudoelastic mechanical behavior.
The marine whelk Busycotypus canaliculatus lays its eggs in disc-like egg capsules arranged in a string, which can be found lining beaches along the east coast of North America (figure 1) [5]. Mechanical studies on the whelk egg capsule (WEC) material reveal unusual properties for a biopolymer, including pseudoelastic behavior and high hysteresis, which is completely reversible [3,[6][7][8]. Figure 1 shows a stress-strain curve of a cyclically loaded strip of egg case material. The material is initially stiff; however, at ∼5% engineering strain, it yields and then proceeds to elongate without a marked increase in stress, after which it eventually stiffens until it ruptures and fails [3,8]. If the egg case polymer is, instead, unloaded prior to failure, a large hysteresis loop is observed as the material returns to its initial length. Immediately upon unloading, it can be reloaded and will exhibit a practically identical stress-strain curve as the first cycle, and this can be repeated over numerous loading cycles.
As mentioned, this behavior is reminiscent of the pseudoelastic mechanical behavior of certain metal alloys, albeit with higher strain values and lower stiffness. Aside from certain metal alloys, similar mechanical behavior is also observed in some rubber-like materials, sometimes referred to as the 'Mullins effect' [9][10][11][12]. This is described within the framework of network models, assuming two interpenetrating networks of different strength. Additionally, hysteretic mechanical behaviors have been observed in other biological systems and described as diffusionless or displacive transformations (often termed 'martensitic' if strain energy is involved in the transformation). A well-characterized example is the length change within the tail sheath of a T4-bacteriophage, in which the cylindrical protein crystal comprising the capsid Protein structural transformations dictate WEC mechanics. Marine whelks lay long chains of disc-like egg capsules. When strips of the WEC tissue are cyclically loaded, there is a large hysteresis, as observed in the above stress-strain curve. The mechanical performance of the material has been correlated with the reversible transformation of α-helical protein subunits (lowstrain phase, LS) into an unstructured, extended conformation (high-strain phase, HS). of the virus is reduced to ∼35% of its initial length during injection of viral DNA into a bacterial host [13,14]. The associated change in the helical arrangement of proteins on the cylinder surface was described as a martensitic phase transformation [13,14]. The transformation from the extended, initial configuration to the contracted, final configuration is irreversible and driven only by the surplus of chemical free energy in the parent phase, which is partially consumed by the deformation of the cylindrical crystal. The transformation in the WEC material during cyclic loading, on the other hand, is completely reversible, leading researchers to investigate the molecular level structure-property relationships in more detail.
It was recently determined that the notable mechanical behavior of the WEC arises from load-bearing α-helical proteins, which undergo a reversible structural transformation during cyclic loading [6,7,15]. Transitions of α-helical proteins to a β-sheet conformation via mechanical stretching have also been observed in materials such as hair, intermediate filaments and in hagfish slime fibers [15][16][17][18][19]; however, in these cases the complete reversibility that characterizes the pseudoelastic behavior of the WEC is not observed. An equilibrium-state model of the WEC mechanical behavior was developed based on pseudoelastic theory, in which the structural changes of the protein building blocks were modeled as a diffusionless or displacive transformation between a low-strain and a high-strain phase, with the coexistence of both phases at intermediate strains [6]. The low-strain phase (α-phase) was determined to behave elastically, whereas the high-strain phase (extended β * -phase) was modeled using a worm-like chain equation [20][21][22] (figure 1). Consistent with experimental observation, the transition from α to β * structure was predicted to occur at a constant stress, following a classical phase transition. While this model fits well with much of the observed structural and mechanical data, it is an equilibrium model and thus it is unable to account for the viscoelastic properties of the material, including the notable hysteresis. In order to address the hysteresis of the WEC material, we present here a modified version of the earlier model [6], which considers specifically non-equilibrium states. This strategy has proven valuable in modeling the mechanical behaviors of other protein-based fibrous materials such as silk [23] and in shape memory polymers [24] or even shape memory alloys [25], such as Ni-Ti [26]. Moreover, models have been developed [27] to describe protein unfolding based on single-molecule experiments including titin [28], tenascin [29] or ubiquitin [30]. Previously, the nonlinear mechanical behavior of the WEC was modeled using a continuum model in which a homogeneous nucleation of the new phase is assumed to occur along the whole chain [31]; however, this macroscopic approach does not incorporate recent biochemical and biophysical evidence for the discontinuous behavior of the protein building blocks during cyclic loading. Based on our new model, our conclusions are two-fold. Firstly, the origin of hysteresis in the WEC appears to stem from a mismatch between the speed at which the molecular phase transformation occurs and the rate at which the material is loaded. Secondly, the curved shape of the WEC stress-strain loading cycle likely originates from a nonlinear relationship between the driving force and the rate of transformation between the two phases (i.e. transformation is slower at higher stress levels). The results provide several experimentally testable hypotheses that should be further pursued.

Deformation and stress state
The transformation takes place within protein fibers, which can be viewed as parallel stringlike arrangements of molecular segments, which-in the low-strain phase-correspond to α-helices, and in the high-strain phase as an extended polypeptide chain. A bundle of these parallel molecular filaments (figure 2) constitute a parallel fiber bundle whose cross-section does not depend on the amount of longitudinal strain. We, therefore, represent the material by a one-dimensional (1D) continuum. We call p the force applied to a single molecular filament, which relates to the stress σ applied to the fiber through the relation where ρ is the number of molecular filaments per cross-sectional area of the fiber (see figure 2). The length of a molecular segment in the reference configuration without load is denoted as s 0 , its length in the actual (deformed) configuration as s and its maximum length after complete extension of the protein (contour length) as L. As a simple measure of strain, we introduce the local strain ε (often denoted as engineering strain or shortly 'strain') as Sketch of a WEC protein fiber during mechanical deformation. It consists of a number ρ of parallel molecular filaments per unit area of fiber crosssection (green dots). Each molecular filament consists of a series of segments that might be either in the low-strain (L) or in the high-strain (H) conformation dependent on tissue deformation.
Given that the largest possible length of the protein is L, the maximum possible strain in this model is given by Of course, one could work also with the natural strain ε n = ln (1 + ε). However, since the following stress and energy terms are measured and expressed in ε and not in ε n , we opt to work with ε. Two phases are distinguished, the 'low-strain phase' (L) and the 'high-strain phase' (H), when the fiber is subjected to a stress σ withσ being the corresponding dimension-free stress, σ = σ/C, and C being a (constant) reference stress used to make the variable dimensionless. The stress σ is not influenced by the cross-sectional area, since it depends only on the number Parameters defining the two models [33,34] used to describe the highstrain phase. The tensile force per molecule is denoted as f. In both models the force is proportional to k B T , the product of the Boltzmann constant and the absolute temperature. The Z-model has bending fluctuations according to a persistence length l p and the M-model has kinking fluctuations with a total number ν of kinks for an extended (contour) length of L.
of molecules in the fiber cross-section (according to (1a)). Consequently, σ is also denoted as engineering stress, which we hereafter refer to just as 'stress'. The string can only transfer a tensile stress, p > 0 (or σ > 0). The following stress-(local) strain relations can be observed for constant temperature, expressed by ε as The quantityk = k/C is a dimension-free elasticity coefficient; F (ε) refers to the mechanical behavior of the high-strain phase, which corresponds to an extended conformation of the molecule [6]. Different models are conceivable to describe the force-extension curve of a fluctuating filament [32].
Here we consider two very simple analytical models yielding nonlinear elastic behaviors for F (ε), namely that according to Misof et al [34], denominated as F M and listed below as F M (ε), and to Zhang et al [33], denominated as F Z and listed below as F Z (ε), as (see figure 3) where the indices Z and M refer to the worm-like chain model [20] in the formulation by Zhang et al [33] and the kinked molecule model by Misof et al [34]. Figure 3 shows sketches of the underlying assumptions and gives a definition of the physical parameters involved. Note that both models have the property that the stress diverges when the strain ε is approaching its limiting value ε * . The reason for this divergence is that we consider the covalent bonds in the polypeptide chain to be inextensible. This implies that the models cannot be used for strains coming too close to ε * . The reference stress C is in both models just the product of the thermal energy k B T (k B being the Boltzmann constant and T the temperature), the molecular density ρ in the fiber (as defined in figure 2) and an inverse length scale characterizing the rigidity of the polypeptide chain. This length scale is the persistence length l p in the worm-like chain model or the contour length divided by the number of kinks L/ν in the other model. Then we have for C z and C M where the labels 'Z' and 'M' point to the material types, as defined in figure 3. The numerical value of ρ has not been confidently determined experimentally for the WEC and was set in the calculations so that the calculated stress values matched the experimentally determined values.

Energetics
We define a Helmholtz free energyW (ε), a dimension-free function of ε, describing the sum of chemical and elastic energy of the molecular segments (see figure 2) in both phases: According to equation (2), we have for the low-strain L-phase Note that with respect to chemistry we address a zero reference chemical energy to the (lowstrain) L-phase. Considering equation (5), one finds forW (ε,s 0 ) in the (high-strain) H-phase HereW 0 stands for the sum of the chemical contribution and the contribution due to the prestress. The integration on the right side of (7) can be performed analytically and we obtain, using (4a) and (4b), respectively, The quantitiesW 0M ,W 0Z are constants and their values are reported below.
Here it should be mentioned that the formulations ofW M ,W Z are in accordance with continuum mechanics. The equilibrium process was discussed previously by Harrington et al [6] and is only summarized here. The equilibrium energy for the complete system is the convex envelope of the energiesW H andW L as a function of the strain ε. At small strains, the total energy isW L (ε). Starting from a lower critical strain ε * L and straining up to a higher critical strain ε * H , the energỹ W T (ε) is given as the common tangent to the two energy curvesW L (ε) andW H (ε) (see figure 4) asW The critical strains ε * L and ε * H can be found by solving the set of nonlinear equations 1 2k Finally, above ε * H , the energy isW H (ε).

Transformation condition and kinetics
We now study the kinetics of the transformation from the L to the H phase (on increasing strain) or from the H to the L phase (on decreasing strain), which was first worked out by Eshelby [35] and later by several researchers (for references see [36]). Our model considers a sharp interface between the regions in the H and L phase, with a thermodynamic force acting on it. The experimental data do not allow us to decide whether there is actually just one or several interfaces as one would expect if there were multiple nucleation sites. Hence, we stick to the simplest hypothesis here. It is also worth mentioning that other theoretical approaches study the movement of an interfacial zone using either a variational principle [37] or a levelset concept [38]. In the current treatment, we prefer the simpler sharp-interface model. If we consider a likely mechanical challenge of the egg capsule such as crashing waves or perhaps biting of a predator, then a stress-controlled process is more likely than a strain-controlled process. We thus consider a stress that increases linearly with time t from t = 0 to t m and then decreases again linearly to zero until time t f : whereσ m is the maximum dimension-free applied stress at t = t m .
From the point of view of thermodynamics, the dimension-free thermodynamic forcef (i.e. the force per unit area of interface, normalized by the constant C Z or C M , with labels LH and HL depending on the direction of the transformation) is acting at the transformation front or, in other words, at the interface between the L and H phases. The thermodynamic forcẽ f LH refers to a transformation from L to H,f HL from H to L. Eshelby [35] was the first who introduced the thermodynamic force as the projection of the so-called energy momentum tensor in the normal direction to the interface. The difference of the thermodynamic forces yields the resultant (driving) thermodynamic force. A very detailed derivation of the thermodynamic force on an interface can be found in Simha et al [39]. We follow the derivation by Abeyaratne and Knowles [36], who previously worked out an expression forf in a small-strain 1D setting formulated forf LH as in which ε L and ε H are the strain values in the L and H phases, respectively, corresponding to a given value of the applied stressσ , according to the relations given in equations (2) and (3), respectively. This makes an inversion of F (ε) necessary, which can be performed analytically for both F M and F Z , but for the sake of convenience, can also be done numerically. Note that no physical interface is formed, and therefore it is unnecessary to consider the interface energy. The advancement of the transformation is characterized by the volume fraction of transformed phase ξ , 0 ξ 1, and its time rateξ . Note that ξ takes the role of an internal variable. Then the rateξ is the conjugate variable to the thermodynamic forcef yielding, according to the thermodynamics of continua, the local dissipation fξ . However, since thermodynamics does not provide a kinetic (or evolution) equation for the internal variable (in our caseξ ), we assume, as a simple and standard assumption in non-equilibrium thermodynamics, a linear kinetic equation such aṡ The quantity κ LH is proportional to the speed of the transformation from L to H with the dimension [κ LH ] = 1 s −1 . Since we assume that no dissipation occurs both in the L phase and in the H phase, the only dissipative term arises from the interface motion, which corresponds, as already mentioned, to the product of interface velocity (represented byξ LH ) and thermodynamic force (f LH ). The dissipation term per unit area of interfaced LH , [d LH ] = 1 s −1 , which is again normalized by the constants C Z or C M , follows as To ensure positive dissipation, the mobility κ LH must be positive. Of course, it does not need to be constant (as considered here for simplicity) but may depend onf LH or its rateḟ LH , see e.g. the discussion in [36]. The dissipation can be considered as the power necessary to break bonds during loading or, alternatively, the energy per time unit extracted from the system to form bonds during unloading.
The total strain ε tot of the molecular filament depends on ξ LH in the form which makes it necessary to integrate equation (10) with respect to time. This is easiest performed by replacing time t with the instantaneous stressσ (dt = dσ /σ LH or dt = dσ /σ HL withσ LH =σ m /t m ,σ HL = −σ m /(t f − t m ), respectively, with the definitions at the beginning of section 2.3), yielding The process starts withf LH = 0 atσ =σ start and continues withf LH > 0. The quantityσ start follows from setting equation (10) to zero and solving the equation with respect toσ start with ε * L = ε L (σ start ) and ε * H = ε H (σ start ). This amounts to solving equation (9b), whereσ start = kε * L = ∂W H ∂ε (ε * H ). The strain at the interface jumps at the stress valueσ start from the Lequilibrium configuration ε * L to the H-equilibrium configuration ε * H . The difference is denoted as transformation strain ε T = ε * H − ε * L . The correspondingσ − ε tot diagram shows a horizontal line, the so-called Maxwell line according toσ start , which is called Maxwell stress. For a further interpretation, see [6].
The above set of equations (9)-(13) can be also used to describe the transformation path from the H to the L phase (in our case the 'unloading' process). Only the labels L and H have to be exchanged to H and L. The mobility κ HL can be different from κ LH . The back-transformation starts atσ =σ start and stops when ξ HL has reached the value 1.
We would like to mention that the so-called Jarzynski equation [40] is often used to describe the pulling of proteins. If we reformulate the thermodynamic force asf LH =g L −g H withg =W −σ ε L , whereg is the (specific) Gibbs (free) energy, then this transformation proceeds forf LH = g =g L −g H > 0. This formulation is in fact equivalent to the Jarzynski equation.

Modeling and experiments
The model introduced in section 2 is formulated in a dimension-free way. All kinetic quantities, such as energies, dissipation and thermodynamic forces, as well as stress, are related to a reference stress C, namely C M or C Z , which means that all these dimension-free quantities must be multiplied by C M or C Z to get to physical quantities.
We now model the transformation process in WEC protein reported in an earlier experimental paper [6]. The relevant parameters are taken from this paper: we consider a molecular segment with the initial length As the free parameter, we use the kinetic parametersσ LH /κ LH andσ HL /κ HL , respectively. Figure 5 demonstrates the stress-strain hysteresis for three different ratiosσ LH /κ LH for material types Z according to Zhang and Evans [33] and M according to Misof et al [34] in comparison with a stress-strain hysteresis measured by Harrington et al [6], figure 1(b).
One can easily recognize that at least semiquantitative agreement of the hysteresis from the model with the experimental hysteresis can be obtained. It should be mentioned that the rather strange behavior of material type Z leading to negative values of ε on the unloading path has Hysteresis for material type Z using the enhanced kinetics equation (14). Numbers 1, 2, 3 designate curves corresponding to the kinetic parameteṙ σ LH /((1 + α) κ LH ) = 5 together with α = 9, 3, 0, respectively. The Maxwell stress is marked by a dashed line.
its origin in the fact that the H-phase covers a certain pre-stress in the amount of 0.78 MPa at strain ε = 0. This part of the curve is not depicted in figure 5(b), since this part is of no physical relevance.
As a further possibility of modeling, we study a so-called enhanced kinetics by applying instead off LH on the rhs of equation (11a) the quantityf LH · ξ −α , α 0, ensuring also the second law (see the text passage in context with equations (11a) and (11b)) and yielding Obviously, the significantly increased transformation rate in the start phase influences the hysteresis remarkably. The results are depicted in figure 6. Curves 1 and 2 are calculated for the enhanced kinetics (see the figure caption). Note that the unloading path above the Maxwell line is nearly the same for all three curves. It is very interesting to note that the enhanced kinetics leads to a hysteresis, whose shape is much closer to the observed one, see figure 5(c).

Discussion
From a biological perspective, the properties of materials such as the WEC are shaped via natural selection to handle cyclic mechanical perturbances by undergoing large extensions and dissipating mechanical energy, observable at the macroscale in the large hysteresis loop. Experimental evidence informs us that mechanical hysteresis stems from reversible transitions in the conformational structure of the protein building blocks, which manifests itself during material yield as a classical phase transition between an ordered α-helical protein structure (α) and an unstructured protein structure (β * ) [6,7]. The current work provides a theoretical framework for understanding the kinetics of this process and how this is manifested in timedependent viscoelastic mechanical behavior at the macroscopic level.
Based on this work, three important conclusions can be drawn that add to the understanding of WEC mechanics and of martensitic behaviors in biopolymers, in general. The first is that modeling the biomolecular subunits as a worm-like chain versus a kinked molecule has only a small influence on the outcome of the model. Only the relative magnitude of the changes differed, but the shapes of the resultant stress-strain curves were preserved. The second observation is that development of hysteresis in the cyclic loading curve, reflecting a dissipative process, arises when the loading rate in the biomolecules is larger than the rate at which the LS phase can transform into the HS phase. Only at very low loading rates (i.e. equilibrium) would one expect no dissipation and, consequently, that the course of the stress-strain curve would follow along the Maxwell line for both the loading and unloading paths. The third observation from our study is that while the simple assumption of a linear relation between thermodynamic driving force and rate of phase transformation (thermodynamic flux) generates hysteresis in the model, it cannot account for the shape of the WEC deformation curve. By considering an 'enhanced' kinetics, where the proportionality constant in the linear relation is replaced by a positive function of the internal variable (i.e. the volume fraction of transformed matter), we can modify the shape of the hysteresis and arrive at some of the experimentally observed features. In the current case, one learns from the nonlinear kinetics that the transformation process occurs much more intensively in the starting phase than in the finishing phase, an observation that is also well known for some shape memory materials [25,26]. These collective observations provide testable hypotheses, namely that there should be a temperature and strain rate dependence in the amount of hysteresis observed in the WEC. In fact, Miserez et al [7] observed that increasing temperature does result in decreased stiffness, yield stress and hysteresis; however, other experiments show that the hysteresis does not depend on stress rate [8]. This second observation is in contrast to our modeling results where the hysteresis is predicted to increase with stress rate. This disparity may have several causes: one is the comparatively small-strain amplitude in the experiment [8] or the frequency range of the loading in comparison to the mobility. Moreover, our model assumes a constant value of the mobility κ LH , which shows the essence of the physical behavior of the fiber. In a real material, κ LH may depend onf LH and, thus, indirectly on the loading speedσ LH . This may influence the way in which hysteresis depends on frequency in cyclic loading.
Experimental measurements on α-helical proteins with single-molecule force spectroscopy provide a biophysical basis for understanding some of the observations of the models in this study. Single-molecule load-extension curves of fibrinogen, a coiled-coil α-helical protein, reveal a characteristic sawtooth pattern indicative of molecular energy dissipation [41,42], whereas myosin exhibits elastic loading behavior lacking notable hysteresis [43]. It was suggested that the relative differences in mechanical behavior between these two α-helical proteins stem from variation in the refolding rate based on the relative structural complexity of the two proteins, i.e. fibrinogen has many disruptions in the α-helical motif sequence, whereas the myosin helix is relatively uninterrupted. Implicit in this result is the fact that the rate of transformation between the L phase and the H phase is dependent on the protein sequence and structure, and therefore the mechanical behavior (and therefore hysteresis) of materials based on α-helical proteins is tunable by molecular augmentation processes such as natural selection or genetic engineering [44]. Finally, as discussed above, stress-controlled processes are typical for natural events. However, strain-controlled processes are usually performed in the laboratory. For the sake of completeness, a corresponding algorithm is presented in the appendix.

Appendix. Strain-controlled process
We start with the strain-controlled process from the L phase to the H phase with given quantities ε tot andε tot . Considering equation (12) one finds and for the rateε tot using (A.1) Sinceσ acts in both phases, one can express ε L by ε H with equation (3) as With F (ε H ) the derivative dF (ε H )/dε H is denoted. We can now expressf LH , equation (10), as a function of ε H by using the first equation of (A.3) and equation (10) as Note that ε tot andε tot are known quantities, so equation (A.5a) represents a first-order ordinary differential equation for ε H (t). As the initial condition we have for t = t start the stressσ start , see the discussion below equation (13), yielding with equation (2) ε H,start = F −1 (σ start ) for t = t start , (A.5b) which enforces inverting F (ε) to F −1 (ε). Then we have to solve the initial value problem for ε H = ε H (t) until ξ LH becomes 1, yielding ε H = ε H,final . Then further straining occurs only in the H phase. The backtransformation starts when ε H has reached ε H (σ start ). Then the strain-controlled backtransformation process from the H phase to the L phase can also be described with the differential equation (A.5a), however, by replacing κ LH by −κ HL with κ HL > 0. Note thatε tot is now a negative quantity. The same start condition as (A.5b) is also valid. The process stops when ξ HL has reached the value 1. Then only the L phase exists again.