The role of topology and thermal backbone fluctuations on sacrificial bond efficacy in mechanical metalloproteins

Sacrificial bonding is a ubiquitous cross-linking strategy for increasing toughness that is found throughout nature in various biological materials such as bone, wood, silk and mussel byssal threads. However, the molecular mechanism of sacrificial bonding remains only poorly understood. Molecular modeling possesses a strong potential to provide insights into the behavior of these cross-links. Here we use Monte Carlo simulations to investigate the effect of reversible sacrificial binding sites on the mechanical properties of single linear polymer chains based on load-bearing metalloproteins found in the mussel byssus. It is shown that the topology of the bonds determines the position and spacing of sacrificial force peaks, while the height of these peaks is intimately tied to the magnitude of thermal fluctuations in the chain that are dependent on effective chain length. These results bear important implications for understanding natural systems and for the generation of strong and ductile biomimetic polymers.


Introduction
Evolution via natural selection has provided an effective means for achieving practical and economic solutions to various physical challenges encountered by load-bearing structures in nature. Natural materials, therefore, are a potentially valuable source of inspiration for the design of novel man-made materials. However, this first requires a thorough understanding of the underlying structure-function relationships that determine the material behavior. For example, in many natural materials increased toughness is achieved via careful tuning of molecular interactions at numerous levels of hierarchy [1]. One particularly effective strategy for increasing toughness found in bone [2,3], wood [4] and some softer biological fibers, such as silk [5][6][7], mussel byssus [8][9][10] and whelk egg capsule [11], is the use of sacrificial bonds (SBs). SBs are load-bearing cross-links which are weaker than covalent bonds making up the backbone structure of the materials that rupture 'sacrificially' in order that the overall material integrity survives [2]. In doing so, these bonds dissipate mechanical energy by unraveling of 'hidden' folded protein length. These bonds can often be reformed after rupture leading to a kind of molecular repair and possibly self-healing behavior in the material. This combination of hidden length unraveling and self-repair capability makes SBs a very effective energy dissipation mechanism and increases the toughness of a material dramatically [2,12].
Hydrogen bonding, electrostatic interactions and most recently, protein-metal coordination bonds have been identified as SBs in natural materials. In metal coordination complexes, several ligands (amino acid side chains such as 3,4-dihydroxyphenylalanine (DOPA) and histidine in protein based systems) contribute lone pairs of electrons to the outer orbitals of transition metal ions. These bonds are ideal as reversible SBs because they require large energy inputs to rupture, but, unlike typical covalent cross-links, will reform rapidly afterwards [13,14]. For example, in mussel byssal threads that are produced by marine mussels to adhere to rocky substrates, loadbearing protein-metal complexes based on histidine-Zn 2+ and DOPA-Fe 3+ interactions have been implicated in increased toughness, hardness and self-healing behavior in various structures of the threads [8][9][10].
While SBs have been correlated with increased toughness in biological materials, as well as in some recently developed supramolecular polymers [15] and bio-inspired metallopolymers [16,17], there is not a solid understanding of how the molecular environment of the SBs influences mechanics. Cao et al [18] demonstrated at the molecular level that recombinantly engineered metal bonds can increase the force to unfold a protein and the total energy dissipated via unfolding. However, the magnitude of this increase was found to be highly dependent upon the location of the cross-link in the folded structure. In other words, the mechanical behavior of the system is influenced by how the SBs are situated within the rest of the chain. Along these lines, it was suggested from molecular simulation that there is an ideal size of β-sheet nanocrystals for maximizing toughness in spider silk, which reflects the most effective use of sacrificial hydrogen bonds [7]. Clearly, there is more to be understood regarding the molecular mechanism of toughening inherent to SBs. Here, we take a step toward unraveling this phenomenon by using Monte Carlo (MC) simulations to examine the influence of topology, chain length and thermal fluctuations on the mechanical behavior of SB reinforced fibrous systems. As a starting point for our model, we consider the histidine-and DOPAcontaining proteins from the mussel byssus; however, the extracted concepts can be abstracted and applied more broadly, such as in the design of mechanical metallopolymers. The results of the simulations show that the efficacy of these SBs is closely related to the magnitude of thermal fluctuations in the backbone of the chain, which depends upon chain length and temperature.

The model
The biological system inspiring this work are the SB-forming histidine-rich regions of the loadbearing proteins in the tough fibrous core of mussel byssal threads [8] and DOPA-rich proteins found in the plaque [10] and the protective cuticle [9]. The his-rich sequences from the core proteins are relatively short segments that consist typically of 30-80 amino acids [19] and are clamped between collagen and flanking domains [8]. This clamping defines an effective length of the polymer over which it is free to fluctuate. The DOPA-rich proteins from the thread plaque and cuticle consist of short tandemly repeated amino acid sequences containing between 5 and 25 mol% DOPA [20]. In the present paper this situation was modeled by a linear chain of N hard spheres with radius R (which we set as the unit of length). The covalent interaction between neighboring beads was described by a Morse potential with the depth E 0 = 5 eV, the width β −1 = 0.5 R, the equilibrium bond length r 0 = 3 R and the actual distance between two neighboring beads i and j, r i j . The contour length of the chain is given by L 0 = (N − 1)r 0 . To model the presence of SBs, N s of the beads were defined as sticky, two of which could form a SB (ρ s = N s /N is the sticky site density). The sticky sites were introduced regularly: the same number of non-sticky sites separating them. The protein backbone is a linear polymer consisting of C-C and C-N bonds which possess average bond energies of several eV. The energetics of the SBs in proteins have not been well characterized; however, several studies on DOPA and histidine-based metal complexes suggest they possess binding energies that range between 20-30% of a C-N bond [21,22]. Thus a value of E SB 0 = 1.25 eV was chosen for the SB in this study, corresponding to 25% of the bond energy of the covalent backbone bond in our simulated polymer. Additionally, the SBs were allowed to open and close reversibly. Simulations mimicking loading experiments were performed by starting from a small end-to-end distance L that was defined by pinning the two outer beads. The positions of the inner beads were updated according to a standard MC procedure using the Metropolis algorithm [23]. For SB updating, one sticky site was chosen randomly and the probability for bond breaking or bond formation was calculated using the Metropolis algorithm, depending on whether the bond was intact or already broken, respectively. During the simulation the forces on the outer beads were recorded and averaged. Subsequently, L was increased and the simulation re-run until the contour length L 0 was reached. Up to 90 million MC step (i.e. jump trials per bead) were performed for the single step. For each simulation an independent starting configuration was produced by relaxing a fully stretched chain L/L 0 = 1 without sticky sites. When the initial end-to-end distance was reached the sticky sites were added and allowed to form SBs. Subsequently this structure was stretched. The model presented bears some resemblance to the dynamic loop model used to describe mitotic chromosomes [24]. The main differences are, firstly, in our model SBs can't form between all monomers but only between sticky sites. Secondly, SB formation and rupture are determined by a distance dependent potential, rather than by rate constants and lifetimes drawn from a Poisson distribution. Thirdly, the monomers can move freely in space, because no lattice model is used.  Figure 1 shows an averaged load-displacement curve for a chain with N s = 4 sticky sites separated by ten monomers (sticky site density ρ s = 0.08) at ambient temperature of k B T = 25 meV. The curve shows four distinct peaks (the second being a double peak IIA and IIB) approximately equally spaced but with significantly different heights. This behavior closely resembles the experimentally found sawtooth patterns, which were reported in a variety of biological systems, e.g. in nacre [25], in single molecule measurements of titin [26,27], tenascin [28], DOPA [22] and modular proteins [29] as well as in the adhesive mucilage pads of diatoms [30,31] and polymer brushes from rat tail tendon [32]. In the following we will focus at explaining the number, position and especially the height of the observed peaks. Peak IV is the trivial contribution corresponding to the onset of backbone stretching and will not be considered here.

Results and discussion
First, we focus on the number and position of the peaks. In the crumpled starting configuration the four sticky sites form two SBs out of three possible topologies: the independent, pseudoknotted and nested configuration [33]. Figure 2 shows a sketch of the different topologies together with the expected load-displacement curve and a simulated one. For the case of the independent configuration the SBs form between sticky sites (1,2) and (3,4), respectively (see figure 2(a)). Due to symmetry the inner part of the chain is now divided into three parts of equal length d = 11r 0 , the outer part of the chain having a total length of d 0 = 16r 0 . The first bond starts stretching when L/L 0 = (d 0 + d + 2r 0 )/L 0 ≈ 0.59 (the additional term 2r 0 taking into account the two SBs). The opening of this loop reveals an additional length of d giving rise to another force peak at L/L 0 ≈ 0.8. When this second SB has broken backbone stretching sets in at L/L 0 ≈ 1. These two bond stretching events correspond to peak (IIB) and (III) in the averaged load-displacement curves shown in figure 1. In the pseudoknotted configuration (figure 2(b)) sticky sites (1,3) and (2,4) form SBs. Ideally the force peak corresponding to the stretching of this structure should be found at the same value of the  chain extension as the first peak in the independent configuration L/L 0 = (d 0 + d + 2r 0 )/L 0 . Nevertheless, it is found at a slightly lower chain extension giving rise to the peak (IIA) in figure 1. This decrease can be attributed to the hard core repulsion of beads that forbid that the chain can exactly align in the direction of the load, which is possible for the independent topology. Whenever the pseudoknotted configuration ruptures the reversibility of SBs assures that another SB is formed between sticky sites (2,3). Like in the independent configuration this SB starts stretching at L/L 0 ≈ 0.8 corresponding to peak (III). Finally, in the nested configuration (see figure 2(c)) the SBs form between sticky sites (1,4) and (2,3), respectively. Thus, the first SB stretching occurs already at extensions L/L 0 = (16r 0 + r 0 )/L 0 ≈ 0.35 corresponding to peak (I) in figure 1. As in the preceding case the second stretching takes place whenever the chain is elongated to L/L 0 ≈ 0.8 once again contributing to peak (III).
While the position of the peaks occurring in figures 1 and 2 can well be explained with the different topologies of the involved SBs, their height can not be understood that easily. Partly the height of the peaks in the averaged load-displacement curve shown in figure 1 can be attributed to the different occurrences of the different topologies. Out of 20 simulation runs there were 12 independent, 2 pseudoknotted and 6 nested configurations. It is due to this different frequency that peak (I) and (IIA) belonging to the less probable pseudoknotted and nested configuration are reduced in height compared to peak (IIB) and (III). Nevertheless, the different probability of topologies can explain only partly the different heights. Also in the single runs shown in figure 2 not a single SB in any topology attains its expected strength (which we define as the maximum force that can be generated by the used potential) shown by the solid lines. For all bond rupture events except for peak (IIA) this expected peak height coincides with the theoretical strength of one SB F max = β E SB 0 /2 = 1.25 eV R −1 . In the case of the pseudoknotted configuration the rupture force should be twice the strength of one SB, because the force can be simultaneously distributed between the two SBs. Nevertheless, the simulation results show that all forces are considerably reduced compared to the expected value of F max . This is surprising because the energy of SBs at the distance of maximum force r max = ln 2/β + r 0 compared to temperature is E(r max )/k B T = 37.5. Thus, the probability of bond failure at this point is extremely small. Therefore it is not expected that the SBs fail, before they are stretched beyond r max corresponding to a peak with height F max in the load-displacement curve.
To understand the reduced height of the forces, a very simple distribution of sticky sites in the chain was investigated. The number of sticky sites was set to N s = 2 (i.e. the formation of only one SB was possible). The position of the sticky sites was in the middle of the chain, one non-sticky site separating them. The starting configuration was chosen such that the system was in its energetic ground state with L/L 0 = (N − 2)/(N − 1) ≡ L * (see also figure 3(a)). Figure 3(b) shows the resulting load-displacement curve in stretching for a chain with N = 5 at k B T = 25 meV. When the chain is stretched from its equilibrium configuration, the SB starts stretching, resulting in an elevated force that reaches its maximum around F = 1.25 eV R −1 , in good agreement with the theoretical strength of a SB. Eventually, the SB fails and the chain is further stretched until backbone elasticity sets in around L/L 0 ≈ 1. The situation changes, however, when the same simulation is repeated on the longer chain with N = 50 (see figure 3(b)). In the corresponding length region no additional load on SB can be observed. Only elevated backbone stretching takes place, being due to a higher entropic contribution from the higher number of monomers. The reason for this unexpected absence of the sacrificial force peak is due to SB opening during the very first simulation steps of the unstretched chain. As discussed before, due to the high ratio of binding energy to thermal energy from a purely energetic point of view the SBs should be thermally stable producing an additional peak in the load-displacement curve as observed for the short chain. The only reasonable physical explanation for SB bond breaking at these low temperatures is that for longer chains, fluctuations of the covalent backbone bonds lead to such large loads on the SB that it fails already before the external load sets in. In other words, it is backbone entropy that leads to SB failure. To validate this explanation the temperature in the simulation was changed. Figure 4 shows load-displacement curves for chains with N = 50 for different temperatures and for a larger range of L/L 0 . A sacrificial force peak is observed for all temperatures and comes close to the theoretical strength of SBs for low temperatures. When the temperature is increased this peak shifts to smaller values of L/L 0 and decreases in height. Simultaneously the entropic background (the baseline of the curves) is rising. Finally, for k B T = 25 meV the sacrificial force peak is shifted below L * . This means that for this temperature, the elevated loading is a purely entropic effect stemming from internal thermal fluctuations and not from external loading. For L/L 0 < L * no load larger than zero can prevail when the temperature is reduced to zero.
To investigate the effect of chain fluctuations on SBs further, the force experienced by SBs during the simulation was recorded. Figure 5 shows the distribution of bond lengths (a) and the loading of SBs (b) for two different chain lengths: N = 5 and N = 33, the latter one corresponding to the longest chain where no SB rupture was observed below L * . For all chains longer than 33 monomers backbone fluctuations are sufficiently large to rupture SBs without external loading. An increase in monomer number leads to an increase of the mean bond length and a corresponding shift toward higher forces. These effects correspond to an effective weakening of the SB and are the reason for premature SB rupture. For increasing chain length, the most probable bond length is shifted to higher values and the distribution is broadened significantly. These two effects lead to a non-zero probability that the SB is strained beyond r max . Subsequently, the bond is brought close to its maximum load, leading to bond rupture. Both distributions show a pronounced asymmetric shape, but-in contrast to the bond length distribution-the force distribution shows a pronounced drop at its right flank. This drop takes place at the maximum load, since above a value of F = 1.25 eV R −1 the distribution is exactly zero by definition. This is also the reason for the slight narrowing of the force histogram for larger N . Therefore, whenever the number of monomers is high enough that the SB length fluctuations considerably exceed r max or-equivalently-the position of the force histograms coincides with the theoretical maximum force, the SB spontaneously fails. The inset in figure 5(b) shows the width and position of the distributions as a function of N . The position was evaluated at the maximum of the curve, while the width was defined to be the full width at half maximum. The temperature, although small compared to the involved binding energies, is a crucial parameter in determining SB rupture. Figure 5(c) shows sacrificial force probability distributions for N = 33 at two different temperatures. An increased temperature leads to a broadening of the distribution and to a shift of the maximum to higher loads. While the shift of the maximum with increasing temperature is similar to the behavior of the probability distributions for increasing N ( figure 5(b)), the width of the distributions behaves differently. On first sight, the observed narrowing of the curve with increasing N is counter-intuitive since one would expect stronger fluctuations with larger N . However, it is the difference in relative length that explains this effect. For large N , the relative length of the chain in the starting configuration reaches L * ≈ 1, while it is L * = 0.75 for N = 5. A short relative length gives the SB enough freedom to follow the covalent bond fluctuations and transmit additional load through the system (see figure 3(b)). On the other hand, a relative length close to 1 reduces the number of possible configurations of the system (in the limit of T = 0 there is only one single state), thus effectively reducing its fluctuations. Figure 5(d) summarizes this effect. Force histograms for N = 33 at two different starting lengths are shown. Clearly, the distribution shifts to larger loads and narrows upon increasing L. For chains longer than N = 33 the SB already fails, when the relative length of the chain is below L * . Thus, only a reduced sacrificial load of ≈ 0.8 eV R −1 that is purely entropic in nature can be transmitted (see figure 4). Such entropic loads have been found experimentally [32,34] and are often described using the freely jointed or -its variant for small bond angles-the worm like chain (WLC) model [35]. Also the shape of the single rupture events in figure 2 can be well fitted with the analytical load-displacement expression obtained for the WLC model [36] with the persistence length L p and an effective contour length L eff 0 as the two fitting parameters. The fitting results yielded a similar persistence length of approximately L p = 5.63 for all peaks except the first rupture of the pseudoknotted configuration that gives L p = 3.1. The effective contour lengths obtained are smaller than the actual contour lengths L 0 . This effective contour length is given by L 0 minus the hidden length shielded by the SBs. The rather small persistence length of the chain reflects its rotational freedom resulting in the high entropic forces. Thus, close to bond rupture the chain behaves as a WLC that is effectively stretched close to its contour length and whose elasticity is described as an entropic spring. Although the WLC equation can well reproduce the found curves, it should be noted that care has to be taken when analyzing polymer stretching using this simple model. Recently it has been shown that the WLC model breaks down for chains too short due to excluded volume effects [37,38]. Nevertheless, because in the present case the chain is rather extended close to bond rupture it can be expected that these effects can be neglected.
These results explain the reduced peak forces observed in figures 1 and 2. Furthermore, the results also explain why SBs that rupture at lower chain extensions are effectively stronger than those that rupture at larger chain lengths (compare e.g. the first and second peak in figure 2(c)). This is because at lower chain extensions the effective length of the chain is reduced which effectively stabilizes SBs. This shows that a reduction in effective bond strength upon loading can also be found for identical SBs, which is an alternative interpretation to the view that such a reduction can only be explained by a successive decrease in SB strength or by multiple molecules loaded in parallel [39]. These results bear the important conclusion that load-displacement curves measured in single-molecule experiments need not necessarily directly reflect the underlying microscopic potential. Rather, the effective potential may be smeared out and significantly reduced in its strength due to entropic contributions. Single-molecule measurements with different chain length might clarify this question. Furthermore, our results indicate that the effective breaking forces of ligand-metal bonds could be considerable lower than the previously measured breaking forces of single bonds (e.g. His-metal complexes [40,41]), if they were embedded in a chain that is free to vibrate. This has important consequences for our understanding of the biological materials and for the design of biomimetic metallopolymers.
The idealized model presented in this paper is a simplified representation of the real situation. Firstly, one isolated polymeric chain is investigated here while in biological systems many protein chains are assembled into larger structures into which the SBs are embedded. This confines the chain into a reptation tube hindering its transverse fluctuations [35,42]. We studied the magnitude of the root mean-square transversal fluctuations in our model and found them to be less than a monomer-monomer distance. This is a very short length compared to the reptation tube diameters of several nanometers found for polymer melts, for instance [43]. This short length indicates that the hinderance of transversal fluctuations due to local reptation is not sufficient to prevent SB rupture. Secondly, the polymer model chosen in this letter corresponds to a freely jointed chain without considering any bending or torsional contributions to the energy. It is likely that additional interactions resulting in special folding patterns may effectively reduce monomer fluctuations, thus stabilizing SBs. Third, in the presented model a SB is always formed between two monomers, while real metal coordination bonds involve three or more partners that are cross linked. Thus, even when one partner detaches, there remain ligands that hold the structure together. However, since the effect of bond weakening is significant for SBs with a rather high but realistic binding energy compared to ambient temperature, it can be concluded that despite the obvious simplifications of the model, protein backbone fluctuations will definitely reduce SB strength in real biological materials. These fluctuations may even prove beneficial for the self-healing process. Thus, the interplay of thermal fluctuations, backbone rigidity, covalent bond strength and SB strength dictate the overall mechanical behavior of these structures. An important question to answer in this context is, whether there is an optimum value of the strength of SBs compared to covalent bonds and if there is an optimum length of the protein chains? Our results indicate that the relatively short segment length of the SB rich domains of only 30-80 amino acids is not incidental, but rather due to SB weakening for longer segments. For hydrogen bonded beta sheet configurations it has already been shown that thermodynamic stability determines the optimum length of the involved beta-strands [44]. Our results in addition indicate that there is a strong temperature dependence of the mechanical properties of materials relying on SBs. Thus, mechanical tests on e.g. byssus fibers at different temperatures seem promising to further decipher the secrets of their mechanical performance. First experiments confirming the importance of temperature have shown that self-healing proceeds faster at higher temperatures [8].

Conclusion
Employing a simple model with reversible crosslinks mimicking SBs between proteins and metal ions in the mussel byssus makes it possible to reproduce characteristic features found in experimental load-displacement curves in natural materials. Characteristic sawtooth patterns corresponding to the rupture of single bonds were observed. The distance between two peaks (the hidden length revealed) is directly linked to the topology of the bonds and corresponds to the length of the loops defined by the SBs. The height of the peak force is considerably lower than the theoretical strength of a SB. It was shown that this reduction is of entropic origin. The capability of SBs to transmit load are drastically reduced at ambient temperature due to thermal fluctuations in the backbone of the chain.