Molecular origin of viscoelasticity in mineralized collagen fibrils

Abstract

Bone is mineralized tissue constituting the skeletal system, supporting and protecting the body's organs and tissues. In addition to such fundamental mechanical functions, bone also plays a remarkable role in sound conduction. From a mechanical standpoint, bone is a composite material consisting of minerals and collagen arranged in multiple hierarchical structures, with a complex anisotropic viscoelastic response, capable of transmitting and dissipating energy. At the molecular level, mineralized collagen fibrils are the basic building blocks of bone tissue, and hence, understanding bone properties down to fundamental tissue structures enables better identification of the mechanisms of structural failures and damage. While efforts have focused on the study of micro- and macro-scale viscoelasticity related to bone damage and healing based on creep, mineralized collagen has not been explored at the molecular level. We report a study that aims at systematically exploring the viscoelasticity of collagenous fibrils with different mineralization levels. We investigate the dynamic mechanical response upon cyclic and impulsive loads to observe the viscoelastic phenomena from either shear or extensional strains via molecular dynamics. We perform a sensitivity analysis with several key benchmarks: intrafibrillar mineralization percentage, hydration state, and external load amplitude. Our results show an increase of the dynamic moduli with an increase of the mineral percentage, pronounced at low strains. When intrafibrillar water is present, the material softens the elastic component, but considerably increases its viscosity, especially at high frequencies. This behavior is confirmed from the material response upon impulsive loads, in which water drastically reduces the relaxation times throughout the input velocity range by one order of magnitude, with respect to the dehydrated counterparts. We find that, upon transient loads, water has a major impact on the mechanics of mineralized fibrillar collagen, being able to improve the capability of the tissue to passively and effectively dissipate energy, especially after fast and high-amplitude external loads. Our study provides knowledge of bone mechanics in relation to pathologies deriving from dehydration or traumas. Moreover, these findings show the potential for being used in designing new bioinspired materials not limited to tissue engineering applications, in which passive mechanisms for dissipating energy can prevent structural failures.

---

Introduction

Bone is the main constituent of the skeletal system, providing shelter to the organs and soft tissues, and allowing sound transmission and conduction for hearing.^1–3 It possesses a special hierarchical structure that has been deeply studied both with experimental and modeling tools across different length scales.^4–12 Among them, the molecular scale is the observation level that is still open to new research avenues, since the material composition and its microstructural factors determine bone properties and dynamic behavior.^13–15 Bone presents a staggered periodicity length, derived from pure collagenous tissues, called the D-period, which is of ≈67 nm, which is composed of two parts named overlap and gap regions.^10,16 The latter is the space where the mineral component, hydroxyapatite (HA), fills most of the interfibrillar voids.¹⁷ As a biopolymer material working at 37 °C, collagen is subject to viscoelastic phenomena.

Viscoelastic properties of bone are crucial to understand how the tissue behaves in terms of fracture toughness,¹⁸ failure,¹⁸ and fatigue, among others,¹⁹ since our daily activities (e.g., walking and running) dynamically load our bone tissues with different intensities and frequencies. Since bone has hierarchical structures, a deep knowledge of the origin of molecular level viscoelasticity may further elucidate macroscale mechanisms observed in previous studies. Specifically, bone exhibits creep and stress relaxation and is able to efficiently dissipate high-energy inputs.²⁰ Crack propagation and catastrophic failure of bones and tendons are the most predominant topics that have been widely investigated, but the origins of micro-damage are still under debate.^21,22 The most accredited hypotheses were related to the heterogeneity of the microstructural composition that is affected by several factors (e.g., gender),²³ and atomic-scale structural changes due to shearing by accumulating the damage from cyclic loadings on tendon.²⁴

A number of studies have provided a description of the mechanics of dry mineralized fibrils, highlighting the important role of HA in strengthening and toughening the collagen macromolecular chains.^10,25,26 At the same time, a previous study²⁶ found that collagen is the component primarily responsible for dissipating energy upon transient loads, with relaxation times that tend to decrease as fast as the increase of the input load velocity, independent of the mineralization percentage and load directionality. This feature is in good agreement with the mechanics of the tissue on the macroscale, thus giving new insights into how the physiology of the macro tissue works. Another component that influences the mechanics of bone is water, which represents about 10% of the total weight of bone, in addition to 30% of collagen and 60% of HA (Fig. 1A).^25,27,28 Researchers have demonstrated the key role of water in providing bone quality in terms of toughness²⁹ and in mitigating the deterioration of aging bone.^30–32

Fig. 1 (A) Hierarchical structure of bone: from the macro tissue to the atomistic scale at which level we find tropocollagen (≈30% w/w), hydroxyapatite (≈60% w/w), and water (≈10% w/w).²⁶ (B) Dry and hydrated topologies employed in our studies.

Water is present in different forms: (i) bulk water fills the spaces of the vascular network to deliver nutrients and evacuate waste. It plays a negligible role in stiffening the material,^33–35 but has effects on the viscoelasticity at the macroscale;³⁶ (ii) structural water is present between the staggered mineralized fibrils to ensure the stabilization and assembly of bone ultrastructure;^37,38 and (iii) bound water fills the voids of the mineralized fibrils,^39,40 and is observable with nuclear magnetic resonance procedures.³⁴ Recently, a 2D finite element model has been proposed to investigate the effect of bound and structural water on bone mechanics, highlighting the ability of the material to tolerate external loads and microcracks, but without any specific speculation on intrinsic viscoelastic phenomena.⁴¹ Moreover, the stiffness of the mineralized fibrils, which is monotonic with the %HA, has a linearly decreasing trend with the increase of the intrafibrillar water content, independent of the presence of minerals.²⁵This interesting behavior gives the opportunity for a research question related to the effect of bound water on the viscoelastic properties of bone that, currently, have been studied mainly at the macroscale.^{20,36,42–44}

In this paper, we explore the viscoelastic properties of mineralized collagen fibrils based on the different mineral and water contents at the atomic scale. Full atomistic models are utilized to understand the atomic scale behavior of dry and hydrated (also later referred to as ‘wet’) mineralized fibrils under cyclic shearing and impulsive loads by using molecular dynamics simulations. Understanding the mechanics of the materials upon cyclic angular loads will complement the studies on bone failure and damage in which shear deformation is more prominent than homogeneous uniaxial loads, especially at the molecular scale.

Previous works have speculated on the viscoelasticity of bone at the macroscale with experiments,^20,42,43 or with simulation tools focusing on a single collagen peptide^45,46 or the mechanics of dry mineralized fibrils upon quasi-static and impulsive loads.²⁶ Recently, Fielder et al. have presented a first study investigating the effect of water on the mechanical behavior of the D-period varying the %HA. They observed a decrease of Young's moduli with respect to similar dry topologies.²⁵ Previous studies have highlighted the role of collagen as the main responsible component for the creep phenomenon in the bone matrix at the macroscale, due to its polymeric structure that is loading-rate dependent; however, its description does not take into account the synergies with the mineral component and water.⁴⁷ Also, other studies have speculated on a possible role of water in improving the viscosity of bone, comparing the viscoelastic response of dry/wet bone upon experimental extensional cyclic loads,⁴³ and confirming an increase of tan δ.^20,42 To the best of our knowledge, a systematic investigation on the viscoelastic properties of mineralized fibrils, considering the role of water and %HA, has been firstly explored in this study.

Our work can also be exploited to create efficient bioinspired synthetic materials that may employ energy transmission and dissipation mechanisms, overcoming the limitations of the biological material. Moreover, from the clinical standpoint, an improved knowledge on the capability of collagen-based tissues to dissipate energy will pave the way to new treatments of pathologies or traumas, which mostly affect the elderly population,^48–50 by actively exploiting the intrinsic features of bone damage.⁵¹

Results and discussion

We first conduct a study considering the D-period of the mineralized collagenous fibrils (Fig. 1B), and perform a sensitivity analysis on the elastic and viscous components of such structures under axial and tangential transient loads. Based on earlier studies on bone-like materials at the macroscale, we take advantage of the dynamic mechanical analysis (DMA, see Materials and methods section) to observe the shear viscoelastic mechanisms⁴² by performing a sensitivity analysis on three different parameters: %HA, angular strain amplitude, and hydration state.

First, we investigate the properties of the dried mineralized collagen fibrils by varying the %HA. Fig. 2A–C present the viscoelastic behavior of the dry topologies with a mineral percentage ranging from 0% to 40% upon angular strain amplitudes (γ₀) of 0.017 (γ_{0 min}) and 0.17 (γ_{0 max}). We obtain G′ (storage modulus) and G′′ as a function of γ₀ (see details in the Materials and Methods section) and describe through mathematical interpolation of the dependency of the moduli on the %HA. G′ appears almost constant over the frequency, independent of γ₀ and the mineral percentage. The frequency averaged value of G′ shows a clear linear trend with %HA for both γ₀s, with a higher angular coefficient for γ_{0 min} (Fig. 2D):

Fig. 2 Results from the shear DMA varying the %HA and shear strain amplitude γ₀ on dry topologies. (A–C) Elastic component (G′), loss modulus, or viscous component (G′′), and the lag between shear and strain (tan δ) over frequency. We observe a G′ that is almost constant independent of the frequency for each strain amplitude. The increase of %HA leads to a stiffening of the compound that is quantitatively higher at γ_{0 min}. In contrast, the increase of γ₀ implies a reduction of the elastic component and a simultaneous increase of the viscosity that is more pronounced above 5 GHz. (D–F) G′ and G′′ over %HA. (D) G′ presents linear trends over %HA with the angular coefficient that is higher at γ_{0 min}. (E) and (F) show the angular coefficient (G′′_m) and the intercept (G′′₀) for the viscous component. Also in this case, we denote linear trends in both cases. The angular coefficient at γ_{0 min} is very low, almost constant, while the trend for γ_{0 max} is markedly increasing, confirming the triggering of viscosity at high strains and frequencies. The intercepts have monotonic increasing values due to the increase of %HA.

At γ_{0 max}, we estimate an increase of +66% from 0%HA to 20%HA, and +25% from 20%HA to 40%HA. In contrast, these increments are more pronounced at γ_{0 min} with +72% from 0%HA to 20%HA, and +53% from 20%HA to 40%HA (Fig. 2A and D).

The values of G′′ show a constant behavior over the frequency range for γ_{0 min} (between 100 MPa and 200 MPa), but show a linear increase for γ_{0 max}, especially for 40%HA (Fig. 2B). Therefore, we linearly fit every data series and, subsequently, linearly fit the extracted angular coefficients (G′′_m) and intercepts (G′′₀) of the loss modulus over the %HA (Fig. 2E and F) to evaluate the dependency of such modulus against the mineral content in the topology:

From the intercepts, we deduce that the frequency-independent component of G′′ is linearly dependent on the %HA, but the increase is stronger for γ_{0 max}, which also registers higher G′′ values throughout the whole frequency range. The frequency dependence of G′′, deduced from the analysis of the angular coefficients, shows a clear linear trend for γ_{0 max}, which is almost absent for γ_{0 min}, as formalized by the low R² value. The above results are consolidated through the analysis of the loss factor tan δ (Fig. 2C), which is lower and frequency-independent at γ_{0 min} (∼0.1), while higher and frequency-dependent for γ_{0 max}. This shows that the tissue hardens when increasing the %HA but, at the same time, it does not show significant variations across input frequencies. In contrast, the viscous behavior is activated at high strains and is more pronounced at high frequencies, confirming the natural propensity of the tissue to efficiently dissipate high-energy inputs.

Then, we explore the difference between wet and dried systems through the same methodology. Fig. 3 shows a synoptic table comparing the DMA results for the dry/wet systems at 0%HA and 20%HA. By looking at the first column, we observe the growth for G′ over the %HA that is enhanced in the dry state rather than in the hydrated one (≈+65% vs. 48%). G′ shows a frequency-independent behavior in the dry state, while we notice a marked frequency-dependent behavior at low frequency in the wet state. For both dry and wet systems, the numerical values for G′ at γ_{0 min} are about double the ones achieved at γ_{0 max}. In contrast, the viscous component of the material is more pronounced for the hydrated states, especially for 0%HA at high frequencies above 5 GHz with γ_{0 max} (Fig. 3C). G′′ shows a frequency-dependent behavior at γ_{0 max} with a numerical value for tan δ up to ≈1.4 (Fig. 3B and C). These features are almost absent for γ_{0 min}, in which G′′ and tan δ are equal to 600 and 0.5, respectively (Fig. 3E and F). Viscous behaviors are, thus, triggered by a stronger shear deformation and hydration, a trend which is also observable in Fig. 2B.

Fig. 3 Results from the shear DMA varying the %HA and shear strain amplitude γ₀ on dry/wet topologies. Panels A–C show the elastic component (G′), viscous component (G′′), and the lag between shear and strain (tan δ) over frequency, at γ₀ = 0.17, while panels D–F present the results at γ₀ = 0.017. We notice a different ratio of the G′ between 0%HA and 20%HA for the dry and wet topologies. Specifically, hydration softens the materials and induces a monotonic increase of G′ up to a saturation above 3 GHz. The role of the hydration is more pronounced at 0%HA since the intrafibrillar voids are completely filled with water. It enhances the viscosity of the material, especially at high frequencies, with a growing trend above 5 GHz at high angular strains. In contrast, at lower γ₀, we observe a high viscosity that tends to a saturation over frequency.

To understand the behavior under the impact of creep loading, we utilize the characterization based on the response upon axial impulses (later referred to as wave propagation – “WP”). This approach is a validated technique, exploited across different scales in earlier works, able to deliver a complete description of a material or structure with a simple broadband impulsive load.^26,45,52

Next, we analyse the propagation and dissipation of energy along the hydrated 0%HA and 20%HA topologies from WP, and later compare the results with the outcomes from the dry materials studied in earlier works.²⁶ We use impulsive loads delivered with compressive displacements, investigating the material behavior across input velocities from 100 m s⁻¹ to 1000 m s⁻¹. We use wave speed and relaxation time as benchmarks to estimate the energy propagation and dissipation. Fig. 4 shows the normalized and average displacements of the C_α atoms over time and position from the loading interface to the fixed constraint, while Table 1 summarizes the main results after post-processing the data, comparing the outcomes with the dry topologies from earlier works.²⁶

Fig. 4 Wave propagation in hydrated mineralized fibrils. Each plot shows the average displacement of the C_α atoms along the topologies from the loading interface to the fixed constraint over time. Panels A and B, panels C and D, and panels E and F relate to 100 m s⁻¹, 500 m s⁻¹ or 1000 m s⁻¹ input velocities (v_i), respectively, for the 0%HA and 20%HA wet topologies. We observe that the slope of displacement over time gets steeper (viz., wave speed increases) with the input velocity and %HA. Moreover, we observe a cleaner energy transmission and faster energy at a high v_i, which are less disturbed by the water molecules.

Table 1 Benchmarks for energy propagation and dissipation. A. Wave speeds estimated for the 0%HA and 20%HA dry and hydrated topologies at different input velocities, along with the average increase introduced by water in the compound. B. Average relaxation time estimated for the 0%HA and 20%HA dry and hydrated topologies at different input velocities, along with the average increase introduced by water in the compound

A	Wave speed [m s^−1]		Input velocity [m s^−1]						Average increase [%]
			100		500		1000
			Dry^26	Wet	Dry^26	Wet	Dry^26	Wet
	% HA	0	2255	3640	1817	3244	1462	3039	≈+83%
		20	3224	3781	2654	3301	2528	2974	≈+20%

---

B	Average relaxation time [ps]	Input velocity [m s^−1]	Dry^26	Wet	Ratio [%]
		100	34.30	3.78	11
		500	14.04	2.60	19
		1000	12.91	2.55	20

---

We observe a decreasing trend of the wave speed numerical values with the increase of the input velocity (Table 1A) in the wet fibril. Specifically, 0%HA at 100 m s⁻¹ presents a wave speed of 3640 m s⁻¹ that decreases to 3039 m s⁻¹ at 1000 m s⁻¹. In contrast, 20%HA shows, on average, higher velocities from 3781 at 100 m s⁻¹ to 2974 m s⁻¹ at 1000 m s⁻¹. The same topologies without water molecules in the intrafibrillar voids show a general increase of the wave speeds that is more pronounced for 0%HA (≈+83%) than 20%HA (≈+20%).

We find the relaxation time decreases over the input velocity, although more slowly than in the dry material results, as shown in Table 1B. We observed the same trend in the hydrated one. Quantitatively, comparing the dry/hydrated states, we observe a noticeable reduction of the numerical values of about one order of magnitude, from 3.78 ps at 100 m s⁻¹ to 2.55 ps at 1000 m s⁻¹ (between 11% and 20% of the results estimated in ref. 26). These results highlight the prominent role of water in dissipating energy with relevant implications on the toughness of the tissue in the case of specific pathologies that are characterized by a reduced presence of water in bones. Our results on the relaxation times from the WP are comparable with those of the DMA upon impulsive loads. It is important to note that the definition used for this parameter is expressed by eqn (6), also used in ref. 26 and 45, and is different from the time required for the convergence of the strain upon external loads to investigate the creep in the material.^46,53 A previous study has observed a decreasing trend for the relaxation time that is stabilized at high-velocity loads in the order of 10 ps.²⁶ The addition of water drastically flattens the curve of one order of magnitude. This also appears in the outcomes from the DMA, in which the addition of water marks a high G′′, especially at a frequency above 5 GHz. 20%HA is the topology that is less affected than 0%HA, since part of the intrafibrillar voids are filled by HA.¹⁰ This outcome is mainly due to the molecular movement within the material that is enhanced by water that acts as a lubricant (i.e., plasticizer) between the collagen and HA.⁵⁴ Looking at the bonds (data not shown), we observe an increment of the H-bonds for the wet topologies (∼+456% at 0%HA and ∼+290% at 20%HA) and their cyclic breaking and recovery induces a stick-slip motion that significantly affects energy dissipation.^46,55–57

Lastly, we estimate the mechanical properties along the fibrils by using the front-wave speed from the WP. In this case, compared to ref. 26, water fills the intrafibrillar voids (bound water), ensuring an almost full continuity to the material. From the WP, we also observe an increase of the wave speeds ≈83% and 20% for 0%HA and 20%HA, respectively. A direct measurement of Young's modulus is not possible using the equation:

E = v²ρ, (3)

which relates it with the wave speed (v) and density (ρ) since the material is not homogeneous and presents different phases.⁵⁸ However, the effect of water on the elasticity behavior is manifested through the DMA, in which we observe the numerical values of G′ that are always lower when the topologies are hydrated. This is in good agreement with the outcomes presented in ref. 25, which showed a softening due to the water component upon quasi-static loads. Moreover, the effect of water is dependent on the load frequency since it gives a monotonic growth to G′ up to a substantial saturation above 3 GHz. In contrast, the elasticity of the dry topologies presents a constant behavior over frequencies.

A major impact is the applicability of our study to auditory apparatus, whose components are collagen-based structures in which the frequency-dependent mechanical response is critical for the physiological functions. The tympanic membrane is a thin hydrated membrane mainly composed of collagen with fiber bundles that have different arrangements across the thickness.⁵⁹ Due to its conformation and structural organization, it filters the acoustic energy collected by the pinna and delivers a selection of the energy content to the three ossicles of the middle ear. Although the macro functions of the eardrum have been deeply studied (see for instance: ref. 52 and 60), its nanostructural behavior is still being investigated.⁶¹ A preliminary study explained the role of water and load directionality in energy transmission and dissipation on a single collagen peptide,⁶² but our outcomes on the dry/wet 0%HA topologies have the chance to elucidate the viscoelastic mechanisms at a higher hierarchy. We observe a reduced elasticity and viscosity in pure dry collagenous structures that are independent of the amplitude of the angular displacement, with numerical values almost constant across frequencies. In contrast, the addition of water promotes the viscous component of the material with a monotonic increasing trend across frequencies with numerical values that are much higher than those of the mineralized topologies. Water, indeed, enables the sliding mechanisms between collagenous fibers that are not constrained by the presence of rigid HA molecules. A similar scenario is also highlighted by the wave transmission study, in which we observe a significant reduction of the relaxation time. Hydrated collagen is the material that mostly constitutes tissues like skin, muscles and tendons that primarily receive external loads and are responsible for the first reaction or response upon external energetic inputs that may damage living bodies. Understanding the capabilities and limits of such structures may really help designing efficient tissue replacements to treat pathologies or traumas.

Furthermore, the results on mineralized fibrils help in understanding the viscoelastic mechanisms of the ossicular chain and the temporal bone, which with different mechanisms provide conductive hearing.^2,3 A deeper knowledge of healthy collagen-based structures can also pave the way to new bioinspired materials to replace and treat native tissues after traumas or pathologies.^51,63–65 With the advent of additive manufacturing and machine learning, a new generation of optimized prostheses and structures can take advantage of this study to mimic, and even improve, the behavior of native tissues for the eventual benefit of patients suffering different bone-related diseases (e.g., osteogenesis imperfecta and osteoporosis) or injuries.^66–69 In view of this, new research avenues could be pursued to develop novel inks for replacing load-bearing bone parts. Hydrated constructs with peculiar viscoelastic properties could be fabricated by tuning the amount of water present in the ink.

The use of biocompatible materials has some major advantages compared to that of synthetic materials, as it promotes the inclusion of the implant in the host site by lowering the probability of rejection. Although this approach is promising, several challenges have to be considered. Firstly, the rheology of the ink highly influences the fabrication process. In fact, high viscosity promoted by water facilitates extrusion through the nozzle, but can also jeopardize the stability of the final structure. On the opposite site, the dehydration of the constructs that occurs over time may change the mechanical properties of the implant and should be prevented by efficiently incorporating water into the constructs.

Conclusions

Bone is indisputably one of the most interesting and peculiar tissues of the human body due to its excellent mechanical properties originating from its hierarchical structures. It accomplishes different mechanical functions related to the body's support and movement, containment, and protection of vital organs, as well as sound conduction and transmission.⁷⁰ Despite the different tasks, all these applications take advantage of bone's capability in managing energy from external sources and conveying/dissipating it to prevent damage to the living tissues or to accomplish specific functions. Based on this premise, we design this work to understand bone's viscoelastic behavior from the molecular perspective, delivering a description of this specific aspect of mechanics of the material, starting from its building blocks: collagen, HA, and water.

In this study, we investigated the viscoelastic properties of dry versus hydrated mineralized fibrils under shear and longitudinal transient loads. From our approach, we achieve two complementary sets of results that provide the understanding of the origin of viscoelasticity in bones, and broadly, of collagen-based materials at the microscale. Even a simple axial loading is transferred as a combination of extensional and shear strain.^71–74 We unveil the origin of the tissue's capability for efficiently transferring energy with its elastic features and, simultaneously, dissipating high energy quickly from loading at its fibrillar levels.

Our study provides new knowledge on bone mechanics, highlighting new aspects at the molecular scale. The viscoelastic behaviors of mineralized fibrils help in understanding the responses of the nanostructures upon transient loads and give a new perspective on how microcracks are tolerated without causing a macroscopic failure of the tissue.^75,76 Moreover, hydration is confirmed as a key player in tissue biomechanics since it severely affects the mechanical properties, creating a synergistic effect with the solid phases of the materials and enhancing the passive mechanisms to prevent structural failure.^25,43,45,46

Materials and methods

Topologies

We use the dry topologies employed in previous works to investigate the effect of mineralization on collagen fibrils with three different percentages of %HA (w%/w%): COL/HA 100/0–0%HA, representing pure collagenous structures (e.g., tendons and the tympanic membrane), COL/HA 80/20–20%HA and COL/HA 60/40–40%HA. Our model is able to properly describe the intrafibrillar mineralization with the mineral component that fills the voids between collagen fibrils, mostly in the gap region.^10,62 We hydrate and randomly ionize 0%HA and 20%HA using VMD software, filling the vacuum in the simulation boxes. We exclude the 40%HA topology since the intrafibrillar voids are almost fully filled with HA.

Force field for molecular dynamics

We perform molecular dynamics simulation using LAMMPS software (https://lammps.sandia.gov/index.html).⁷⁷ We employ an improved version of the CHARMM force field in all MD simulations, consistent with previous works.^78–80 Our parameters include hydroxyproline (HYP), which is one of the building blocks of the collagen protein structure. We use previous studies on atomistic simulations matching quantum mechanical models to provide a suitable parameter set for HYP.⁸¹

Equilibration

Full periodic topologies. The following procedure is applied to the topologies with full periodic boundary conditions. The equilibration of the full periodic topologies includes in first instance 11 cycles, 0.5 ns each, of NVT integration at 310 K with a 1 fs time step, fixing the net momentum of all the atoms as zero. Thus, achieving the energy convergence after 5.5 ns, we use an anisotropic NPT ensemble for 6 cycles, 0.5 ns each, at 310 K and 1 atm. The final energy convergence was reached at 8.5 ns.

Chopped topologies. In order to study the wave propagation and energy dissipation along the structures, we chopped the bonds at the edges of the simulation box from the reference geometry, leaving the periodic conditions only along the y- and z-directions.

We preliminarily define two groups of atoms. The first group is named “bound” – BD, including the C_α atoms at the end of the triple helices, and the second group is labelled “mobile” – MB, including all the remaining atoms. We relax the structures with several steps to prevent abrupt structural changes from the reference. The duration of the first step is 0.1 ns with an NVE ensemble fixing the momentum of the atoms belonging to the MB group. Then, we perform two cycles of 50 ps NVT ensemble each with a Langevin thermostat: in the first one, we change the temperature from 10 K to 310 K, and in the second one, the temperature is kept fixed at 310 K. Afterwards, we fix the temperature of the MB group at 310 K with a Berendsen thermostat for 0.2 ns. During the equilibration, we fix the displacements and velocities of the MB group along the y-/z-axis, without any constraint along the x-axis.

We stabilize the system without this last constraint for 0.2 ns and, then, we use an NVT ensemble for all the atoms for 0.2 ns, fixing the temperature at 310 K. Finally, we employ an NPT ensemble for 2 ns at 310 K with a pressure constraint at 1 atm only along the y-/z-direction. A final energy convergence is reached at 2.8 ns based on RMSD, energy and pressure.

Dynamic mechanical analysis

Dynamic mechanical analysis (DMA) is an experimental technique that aims at characterizing the viscoelastic behavior of a material upon cyclic strains.⁸² In this work, we implement a computational version of the DMA by imposing a periodic angular deformation (γ) to the simulation box:

γ = γ₀ sin(2πft), (4)

in which γ₀ is the deformation amplitude and f is the frequency of oscillation (Fig. 5A and B).

Fig. 5 Methodologies to study the viscoelastic behavior of the material. (A and B) Shear DMA of the material. We apply a cyclic angular deformation with an amplitude equal to γ₀ to the full periodic topologies, varying the %HA, amplitude, frequency, and hydration state of the material. (C and D) Wave propagation in the chopped topologies. (C) Topology fixed at one end and loaded by an axial compressive displacement-based load at the opposite end. (D) Impulsive load evolution over time. The load is applied for T = 2 ps with an input velocity (v_i) equal to 100 m s⁻¹, 500 m s⁻¹ or 1000 m s⁻¹. The observation time t_f is equal to 140 ps.

We study the behaviour of the full periodic dry/hydrated topologies as a function of different features: the deformation frequency (viz., 1 GHz–10 GHz), amplitude (0.017 rad and 0.17 rad, corresponding to the ε_xy strain of the top face of the box of 1% and 10% of the total length of the box), mineralization percentage (0%HA–40%HA) and hydration. Following earlier studies, we also investigate the behavior of wet structures at lower frequencies (0.01 GHz and 0.1 GHz). We do not extend the simulation for the dry state at such low frequencies, since the obtained curve shows no trend. This choice is motivated by the computational effort that these simulations require. Simulations are carried out in the NVT ensemble at 300 K. We record stress data every 1 ps. The output data are averaged over the values calculated at every time step. We fit the shear stress (τ_xy) with a sinusoidal function:

τ_xy = τ₀ sin(ωt + δ), (5)

in which δ is the delay between the stress and the strain.

According to ref. 83, we define the storage modulus (G′) and loss modulus (G′′) as parameters that give the capability of the material to store (elastic component) and dissipate (viscous component) energy, respectively:

From eqn (3), we estimate tan δ as:

We organize our outcomes in order to provide a comparison among related systems to highlight how the viscoelastic behavior is affected by: (1) the %HA in the compound; (2) the presence of water; and (3) the amplitude of the angular strain. For each frequency of deformation, we simulate three shear loading cycles and report the average value along with the error bars to highlight the main tendencies.

Wave propagation analysis (WP)

We fix the atoms belonging to the BD group at one end and use the other to apply a displacement-based load. Based on a previous study that demonstrates the similar outcomes from tensile/compressive loads,²⁶ we carry out our study using impulsive compressive loads with three reference input velocities (v_i): 100 m s⁻¹, 500 m s⁻¹ and 1000 m s⁻¹. The load is delivered in 2 ps, while the total observation time (t_f) is 140 ps (Fig. 5C and D).

We post-process the displacements of the C_α atoms along the x-axis and estimate the wave speed by taking into account the evolution of the average displacement peak of the traveling wave.

We study the dissipative behavior of the material through the evolution of the square value of the velocity of the wave front (v_max²), which resembles the normalized kinetic energy of the system. By assuming the following approximation:

we define τ as the relaxation time of the system, which is the amount of time needed to reduce the amplitude of the function of a factor equal to e⁻¹.^26,55,62

Author contributions

Conceptualization, M.M., A.D., and M.J.B.; data curation, M.M. and A.D.; formal analysis, A.D.; funding acquisition, M.M. and M.J.B.; investigation, M.M., A.D., S.D., and M.J.B.; methodology, M.M., A.D., G.S.J. and M.J.B.; project administration, M.M. and M.J.B.; resources, M.J.B.; software, M.J.B.; supervision, S.D. and M.J.B.; validation, M.M., A.D., and G.S.J.; visualization, M.M. and A.D.; roles/writing – original draft, M.M.; and writing – review and editing, M.M., A.D., G.S.J., S.D., and M.J.B. All authors intellectually contributed and provided approval for publication.

Conflicts of interest

All the authors declare no conflict of interest.

Acknowledgements

This work was supported by the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement COLLHEAR no. 794614. G. S. J. and M. J. B. acknowledge additional support from ONR (N000141612333) and AFOSR (FATE MURI FA9550-15-1-0514), as well as NIH U01HH4977, U01EB014976, and U01EB016422. A. D. was financially supported by the Progetto Roberta Rocca (Doctoral at MIT) Fellowship. G. S. J. acknowledges support from the Laboratory Directed Research and Development (LDRD) Program of Oak Ridge National Laboratory. ORNL is managed by UT-Battelle, LLC, for DOE under contract DE-AC05-00OR22725.

References

1. G. S. Jung and M. J. Buehler, Annu. Rev. Biomed. Eng., 2017, 19, 435–457.
2. M. Milazzo, S. Danti, F. Inglese, G. Jansen van Vuuren, V. Gramigna, G. Bonsignori, A. De Vito, L. Bruschini, C. Stefanini and S. Berrettini, J. Biomed. Mater. Res., Part B, 2016, 105, 2495–2506.
3. S. Stenfelt and R. L. Goode, Otol. Neurotol., 2005, 26, 1245–1261.
4. F. Hang and A. H. Barber, J. R. Soc., Interface, 2011, 8, 500–505.
5. J. Samuel, J.-S. Park, J. Almer and X. Wang, J. Mech. Behav. Biomed. Mater., 2016, 57, 128–138.
6. H. C. Wells, K. H. Sizeland, H. R. Kayed, N. Kirby, A. Hawley, S. T. Mudie and R. G. Haverkamp, J. Appl. Phys., 2015, 117, 44701.
7. B. Depalle, Z. Qin, S. J. Shefelbine and M. J. Buehler, J. Bone Miner. Res., 2016, 31, 380–390.
8. D. K. Dubey and V. Tomar, J. Eng. Mater. Technol., 2013, 135, 021015.
9. M. C. Tourell and K. I. Momot, J. Phys. Chem. B, 2016, 120, 12432–12443.
10. A. K. Nair, A. Gautieri, S. W. Chang and M. J. Buehler, Nat. Commun., 2013, 4, 1–9.
11. K. Tai, F.-J. Ulm and C. Ortiz, Nano Lett., 2006, 6, 2520–2525.
12. H. S. Gupta, W. Wagermaier, G. A. Zickler, D. Raz-Ben Aroush, S. S. Funari, P. Roschger, H. D. Wagner and P. Fratzl, Nano Lett., 2005, 5, 2108–2111.
13. L. M. McNamara and P. J. Prendergast, J. Biomech., 2007, 40, 1381–1391.
14. S. Askarinejad, J. E. Johnson, N. Rahbar and K. L. Troy, J. Mech. Behav. Biomed. Mater., 2019, 96, 269–278.
15. C. J. Caitlyn, O. G. Andriotis, V. Nedelkovski, M. Frank, O. L. Katsamenis and P. J. Thurner, Bone micro- and nanomechanics, in Encyclopedia of Biomedical Engineering, ed. R. Narayan, Elsevier, 2019, vol. 2, pp. 22–44.
16. J. P. R. O. Orgel, T. C. Irving, A. Miller and T. J. Wess, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 9001–9005.
17. A. Gautieri, S. Vesentini, A. Redaelli and M. J. Buehler, Nano Lett., 2011, 11, 757–766.
18. R. M. Kulin, F. Jiang and K. S. Vecchio, J. Mech. Behav. Biomed. Mater., 2011, 4, 57–75.
19. T. M. Wright and W. C. Hayes, J. Biomed. Mater. Res., 1976, 10, 637–648.
20. R. S. Lakes, J. L. Katz and S. S. Sternstein, J. Biomech., 1979, 12, 657–678.
21. T. L. Norman and Z. Wang, Bone, 1997, 20, 375–379.
22. R. A. Schapery, Int. J. Fract., 1975, 11, 141–159.
23. Z. Wu, T. C. Ovaert and G. L. Niebur, J. Orthop. Res., 2012, 30, 693–699.
24. J. L. Zitnay, G. S. Jung, A. H. Lin, Z. Qin, Y. Li, S. M. Yu, M. J. Buehler and J. A. Weiss, Sci. Adv., 2020, 6, eaba2795.
25. M. Fielder and A. K. Nair, Biomech. Model. Mechanobiol., 2019, 18, 57–68.
26. M. Milazzo, G. S. Jung, S. Danti and M. J. Buehler, ACS Nano, 2020, 14, 8307–8316.
27. J. D. Currey, J. Mater. Sci., 2012, 47, 41–54.
28. S. V. Dorozhkin and M. Epple, Angew. Chem., Int. Ed., 2002, 41, 3130–3146.
29. J. Samuel, D. Sinha, J. C.-G. Zhao and X. Wang, Bone, 2014, 59, 199–206.
30. X. Wang, H. Xu, Y. Huang, S. Gu and J. X. Jiang, J. Bone Miner. Res., 2016, 31, 1026–1029.
31. X. Wang, R. Hua, A. Ahsan, Q. Ni, Y. Huang, S. Gu and J. X. Jiang, JBMR Plus, 2018, 2, 164–173.
32. Y. N. Yeni, C. U. Brown and T. L. Norman, Bone, 1998, 22, 79–84.
33. E. E. Wilson, A. Awonusi, M. D. Morris, D. H. Kohn, M. M. J. Tecklenburg and L. W. Beck, J. Bone Miner. Res., 2005, 20, 625–634.
34. E. E. Wilson, A. Awonusi, M. D. Morris, D. H. Kohn, M. M. J. Tecklenburg and L. W. Beck, Biophys. J., 2006, 90, 3722–3731.
35. F. W. Wehrli and M. A. Fernández-Seara, Ann. Biomed. Eng., 2005, 33, 79–86.
36. N. Sasaki and A. Enyo, J. Biomech., 1995, 28, 809–815.
37. J. Bella, B. Brodsky and H. M. Berman, Structure, 1995, 3, 893–906.
38. P. Fratzl, Collagen: Structure and mechanics, 2008.
39. C. Jaeger, N. S. Groom, E. A. Bowe, A. Horner, M. E. Davies, R. C. Murray and M. J. Duer, Chem. Mater., 2005, 17, 3059–3061.
40. T. T. Pham, T. Lemaire, E. Capiez-Lernout, M. Lewerenz, Q.-D. To, J. K. Christie, D. Di Tommaso, N. H. de Leeuw and S. Naili, Theor. Chem. Acc., 2015, 134, 59.
41. M. Maghsoudi-Ganjeh, X. Wang and X. Zeng, J. Mech. Behav. Biomed. Mater., 2020, 101, 103454.
42. J. Yamashita, B. R. Furman, H. R. Rawls, X. Wang and C. M. Agrawal, J. Biomed. Mater. Res., 2001, 58, 47–53.
43. J. Yamashita, X. Li, B. R. Furman, H. R. Rawls, X. Wang and C. M. Agrawal, J. Biomed. Mater. Res., 2002, 63, 31–36.
44. J. Black and E. Korostoff, J. Biomech., 1973, 6, 435–438.
45. M. Milazzo, G. S. Jung, S. Danti and M. J. Buehler, ACS Biomater. Sci. Eng., 2020, 6, 1367–1374.
46. A. Gautieri, S. Vesentini, A. Redaelli and M. J. Buehler, Matrix Biol., 2012, 31, 141–149.
47. S. M. Bowman, L. J. Gibson, W. C. Hayes and T. A. McMahon, J. Biomech. Eng., 1999, 121(2), 253–258.
48. J. Parkkari, P. Kannus, M. Palvanen, A. Natri, J. Vainio, H. Aho, I. Vuori and M. Järvinen, Calcif. Tissue Int., 1999, 65, 183–187.
49. S. N. Robinovitch, W. C. Hayes and T. A. McMahon, J. Biomech. Eng., 1991, 113(4), 366–374.
50. M. Beunza-Sola, Á. M. Hidalgo-Ovejero, J. Martí-Ayerdi, J. G. Sánchez-Hernández, M. Menéndez-García and S. García-Mata, Postgrad. Med. J., 2018, 94, 76–80.
51. J. G. Snedeker and J. Foolen, Acta Biomater., 2017, 63, 18–36.
52. M. Milazzo, E. Fallah, M. Carapezza, N. S. Kumar, J. H. Lei and E. S. Olson, Hear. Res., 2016, 346, 1–13.
53. N. Karathanasopoulos, G. Arampatzis and J.-F. Ganghoffer, J. Mech. Behav. Biomed. Mater., 2019, 90, 256–263.
54. W. D. Callister, et al., Fundamentals of materials science and engineering, Wiley, London, 2000, vol. 471660817.
55. Z. Xu and M. J. Buehler, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 061910.
56. S. Keten and M. J. Buehler, Nano Lett., 2008, 8, 743–748.
57. G. E. Fantner, J. Adams, P. Turner, P. J. Thurner, L. W. Fisher and P. K. Hansma, Nano Lett., 2007, 7, 2491–2498.
58. W. Kausel, Vibrations and Waves, in Springer Handbook of Systematic Musicology, ed. R. Bader, Springer Handbooks, Springer Berlin Heidelberg, Berlin, Heidelberg, 2018, pp. 29–47.
59. D. J. Lim, Acta Oto-Rhino-Laryngol. Belg., 1995, 49, 101–115.
60. J. P. Fay, S. Puria and C. R. Steele, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 19743–19748.
61. S. Moscato, A. Rocca, D. D'Alessandro, D. Puppi, V. Gramigna, M. Milazzo, C. Stefanini, F. Chiellini, M. Petrini and S. Berrettini, et al. , Appl. Sci., 2020, 10, 3043.
62. M. Milazzo, G. S. Jung, S. Danti and M. J. Buehler, ACS Biomater. Sci. Eng., 2020, 6, 1367–1374.
63. M. Milazzo, S. Danti, F. Inglese, S. Berrettini and C. Stefanini, in European Society for Precision Engineering and Nanotechnology, Conference Proceedings - 18th International Conference and Exhibition, EUSPEN 2018, ed. O. Riemer, E. Savio, D. Billington, R. K. Leach and D. Phillips, euspen, 2018, pp. 375–376.
64. C. Mota, S. Danti, D. D'Alessandro, L. Trombi, C. Ricci, D. Puppi, D. Dinucci, M. Milazzo, C. Stefanini, F. Chiellini, L. Moroni and S. Berrettini, Biofabrication, 2015, 7, 025005.
65. C. Mota, M. Milazzo, D. Panetta, L. Trombi, V. Gramigna, P. A. Salvadori, S. Giannotti, L. Bruschini, C. Stefanini, L. Moroni, S. Berrettini and S. Danti, J. Mater. Sci.: Mater. Med., 2018, 29(5), 1–14.
66. M. Milazzo, N. Contessi Negrini, S. Scialla, B. Marelli, S. Farè, S. Danti and M. J. Buehler, Adv. Funct. Mater., 2019, 29, 1903055.
67. M. Milazzo, P. G. G. Muyshondt, J. Carstensen, J. J. J. Dirckx, S. Danti and M. J. Buehler, J. Mech. Behav. Biomed. Mater., 2020, 103, 103541.
68. G. X. Gu, C. T. Chen and M. J. Buehler, Extreme Mech. Lett., 2018, 18, 19–28.
69. Y.-T. Wang, S.-Y. Chang, Y.-C. Huang, T.-C. Tsai, C.-M. Chen and C. T. Lim, Nano Lett., 2013, 13, 5247–5254.
70. E. C. Ekwueme, J. M. Patel, J. W. Freeman and S. Danti, in Bioresorbable Polymers for Biomedical Applications, Elsevier, 2017, pp. 391–422.
71. T. Stauber, U. Blache and J. G. Snedeker, Matrix Biol., 2020, 85, 68–79.
72. S. E. Szczesny, J. L. Caplan, P. Pedersen and D. M. Elliott, Sci. Rep., 2015, 5, 14649.
73. H. R. C. Screen, J. Biomech. Sci. Eng., 2009, 4, 481–499.
74. S. E. Szczesny, K. L. Fetchko, G. R. Dodge and D. M. Elliott, J. Orthop. Res., 2017, 35, 2127–2134.
75. R. O. Ritchie, M. J. Buehler and P. Hansma, Phys. Today, 2009, 62, 41–47.
76. M. J. Buehler, Nanotechnology, 2007, 18, 295102.
77. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
78. B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. a. Swaminathan and M. Karplus, J. Comput. Chem., 1983, 4, 187–217.
79. Z. Qin, A. Gautieri, A. K. Nair, H. Inbar and M. J. Buehler, Langmuir, 2012, 28, 1982–1992.
80. A. K. Nair, A. Gautieri and M. J. Buehler, Biomacromolecules, 2014, 15, 2494–2500.
81. S. Park, R. J. Radmer, T. E. Klein and V. S. Pande, J. Comput. Chem., 2005, 26, 1612–1616.
82. R. P. Chartoff, J. D. Menczel and S. H. Dillman, Thermal Anal. Polym., 2009, 387–495.
83. S. S. Bair, in Handbook of Lubrication and Tribology, Volume II: Theory and Design, 2nd edn, 2012.

---