A computational study of mechanical properties of collagen-based bio-composites

ABSTRACT Studying changes in collagen deformation behavior at the nanoscale due to variations in mineralization and hydration is important for characterizing and developing collagen-based bio-composites. Recent studies also find that carbon nanotubes (CNTs) show promise as a reinforcing material for collagenous bio-composites. Currently, the effects of variation in mineral, water, and CNT content on collagen gap and overlap region mechanics during compression is unexplored. We use molecular dynamics simulations to investigate how variations in mineral, water, and CNT contents of collagen bio-composites in compression change their deformation behavior and thermal properties. Results indicate that variations in mineral and water content affect the collagen structure due to expansion or contraction of the gap and overlap regions. The deformation mechanisms of the gap and overlap regions also change. The presence of CNTs in non-mineralized collagen reduces the deformation of the gap region and increases the bio-composite elastic modulus to ranges comparable to mineralized collagen. The collagen/CNT bio-composites are also determined to have a higher specific heat than the studied mineralized collagen bio-composites, making them more likely to be resistant to thermal damage that could occur during implantation or functional use of a collagen collagen/CNT bio-composite biomaterial.


Introduction
Collagen is a fundamental part of biological tissues such as bone, tendon, muscle, and skin (Fratzl 2008). A hierarchical structure is a primary characteristic of collagenbased bio-composites (Currey 2002;Fratzl 2008). At the nanoscale, collagen in bone is mainly bundles of fibrils embedded in an extrafibrillar mineral matrix. The fibrils are a composite of primarily type I collagen protein, intrafibrillar mineral, and water. The collagen molecules are arranged in a periodic triclinic structure (Wess et al. 1995;Orgel et al. 2006), forming regions where there are gaps and overlaps between collagen molecules. The combined length of a single gap and overlap region forms what is known as D-banding, where D is approximately 67 nm (Nikolov and Raabe 2008;Streeter and de Leeuw 2010). The brittle mineral phase is primarily hydroxyapatite (HAP), and experiments have shown that fibril mineralization initially begins in the gap region (Nudelman et al. 2010;Wang et al. 2012). Water fills the remaining fibril voids and can be either mobile water typically found in pore channels, or structural water found between collagen and mineral molecules. (Timmins and Wall 1977;Nyman et al. 2006;Zhang et al. 2007;Gul-E-Noor et al. 2015).
Inspired by the nanoscale structure and composition of bone, this study investigates collagen-based bio-composites to aid the development of biocompatible materials. Studying the nanoscale structure and mechanics of collagenous tissues and composites is important because their macroscale properties can change depending on many factors. For example, mineral and water content in bone have been found to decrease as a person gets older, making them more brittle (Mueller et al. 1966;Rauch and Glorieux 2004;Saxon et al. 2014). These changes in properties begin with nanoscale changes in bone structure and composition (Hamed et al. 2010;Pradhan et al. 2011;Hamed and Jasiuk 2013;Gul-E-Noor et al. 2015;Depalle et al. 2016).
The effects of the nanoscale mechanics and structure of collagenous tissues on their ultrastructural mechanical properties and the underlying mechanisms responsible for these effects have been studied using experimental techniques such as x-ray diffraction (Sasaki and Odajima 1996), atomic force microscopy (Eppell et al. 2006;van der Rijt et al. 2006;Minary-Jolandan and Yu 2009), and nuclear magnetic resonance spectroscopy (Gul-E-Noor et al. 2015). Computational methods such as molecular dynamics simulations CONTACT Arun K. Nair nair@uark.edu Institute for Nanoscience and Engineering, University of Arkansas, Fayetteville, AR, USA (Zhang et al. 2007;Streeter and de Leeuw 2010;Tang et al. 2010;Gautieri et al. 2011;Nair et al. 2013;Depalle et al. 2016), finite element analysis (Thomopoulos et al. 2006;Hamed and Jasiuk 2013;Ahsan 2017;Lin et al. 2017), and ab initio methods (Dubey and Tomar 2013) have also been used to study the mechanical properties of collagenous tissues.
Studying the thermal properties of collagen-based bio-composite scaffolds is important because a high specific heat capacity and thermal conductivity may prevent damage if the functional use of the bio-composite (or the placement of it in vivo) produces a lot of heat. Experiments by Holmes et al. find that telopeptide crosslink cleavage alters the flow of heat and structure of collagen (Holmes et al. 2017). Qu and Tomar used molecular dynamics simulations to study the thermal properties of collagen-hydroxyapatite nanocomposites and found that the weight fractions, material geometry, and strain all affect the bio-composite's thermal properties (Qu and Tomar 2015).
Carbon nanotubes (CNTs) are being studied as a material to strengthen collagen-based bio-composites (MacDonald et al. 2005;Zanello et al. 2006;Meng et al. 2009;Hirata et al. 2011;Amirian et al. 2012;Jing et al. 2017;Silva et al. 2017;Tanaka et al. 2017). Studies investigating the structure of collagen/CNT bio-composites observed the nanotubes primarily in the collagen gap region. A schematic of the structure of a collagen fibril with CNTs can be seen in Figure 1, with the CNTs between collagen molecules in the gap region. Carbon nanotubes (CNTs) have excellent thermal and electrical conductivity (Balandin, 2011;Charlier et al. 1996), effective bio-molecular adsorption properties (Kang et al. 2008), and a relatively high elastic modulus of about 1 TPa (Zhang et al. 2002). This makes CNTs practical as a stand-alone material or for use in composites in applications that require high stiffness.
Previous studies have investigated the nanomechanics of collagenous bio-composites as a singular function of either water content or mineral content. However, both mineral and water content in collagen can vary simultaneously and it has yet to be investigated how this affects the mechanical and thermal properties of collagenous bio-composites. Previous studies have also focused on the structure and biological viability of collagen/CNT bio-composites, but the nanoscale mechanics and thermal properties of collagen/CNT biocomposites and its comparison to mineralized collagen is still unexplored. We investigate how the mechanical and thermal properties of collagen-based bio-composites, simulated with molecular dynamics, change as a function of both mineral and water content. We also compare the elastic modulus values of the bio-composite to established theoretical models. Finally, we aim to determine if collagen/CNT bio-composites can mimic the structure, mechanics, and thermal properties of mineralized collagen bio-composites.

Materials and methods
Studying the nanoscale mechanics of bio-composites can be challenging using experimental techniques. One method for studying their nanoscale mechanics is through computational methods such as molecular dynamics (MD). In this study, MD simulations are performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) program (Plimpton 1995). Details on the development of the bio-composites are available in the S1 section of the supplementary material. The bio-composites are mineralized at 0, 20, and 40 wt% to model varying degrees of mineralization. The bio-composites are hydrated at 0, 2, 4, 10, or 20 wt% water. Figure 2(a) shows the equilibrated bio-composite gap and overlap regions, and the model unit cell whose length along the x-axis corresponds to  Details of the development of the CNT models for the collagen/CNT bio-composites are available in the S2 and S3 sections of the supplementary material. We vary the water content from 10 wt% to 20 wt% for a bio-composite with 10 wt% single-walled carbon nanotubes (SWCNTs). Higher water contents of 10-20 wt% are chosen for the collagen/CNT composites in order to compare their mechanical properties to the non-mineralized collagen composites with the water content of healthy fibrils (10-20 wt%), while the water contents of the mineralized composites were varied until they reached modulus values comparable to the collagen/CNT composites with 10-20 wt% water content. We also compare collagenous bio-composite with 7 wt% SWCNTs with 11 wt% multi-walled carbon nanotubes (MWCNTs). Figure 2(b) shows an equilibrated bio-composite with 10 wt% SWCNT and 20 wt% water. A zoomed-in view of the unit cell shows a single SWCNT per unit cell. In agreement with experiment, the SWCNT bends slightly during equilibration to fit between collagen molecules in the gap region (MacDonald et al. 2005;Tan et al. 2010). Table 1 shows the material contents for all the bio-composites in our study, as well as the total number of atoms in each model.
The method for applying quasi-static stress is the same as used in a previous study (Fielder and Nair 2018), and further details are available in the S1 section of the supplementary material. In short, the unit cell dimensions are allowed to change until the pressure (which corresponds to the stress) of the unit cell in the x-direction reaches the target stress value. The compressive stress is plotted vs. the linear strain, the Young's modulus is plotted vs. the water content, and we plot the gap/overlap ratio vs. applied stress. The gap and overlap region lengths are determined by the measuring distance between the terminal amino acid residues of the collagen molecules, which corresponds to the   beginning and end of the gap and overlap regions. This can be seen in the schematic of Figure 1 with the red and green points showing the distances between the gap and overlap regions. The thermal properties of the biocomposites are determined by the specific heat capacity, whose method of calculation can be found in section S1 of the supplementary material. Briefly stated, based on the total energy, pressure, and volume of the simulations, the enthalpy of the system can be determined. Subsequently, the specific heat capacity can be calculated utilizing the ideal gas constant and the enthalpy, mass, and temperature of the system.

Stress vs. Strain
The compressive stress vs. strain is plotted in Figure 3 for all samples. We find that as the water content increases the strain increases, and as the mineral content increases from 0 to 20 to 40 wt% ( Figure 3(a-c) respectively), the strain decreases, and the stress-strain relation becomes more linear. Next, we compare the effect of an increase in water content from 2 to 4 wt% at an applied stress of 60 MPa and find that at 0, 20, and 40 wt% mineral there is 30%, 10% and 114% increase in strain, respectively. Comparing the bio-composites with 2 wt% water at 60 MPa applied stress, an increase in mineral content from 0 wt% to 20 wt% results in 44% decrease in strain. A further increase in mineral content from 20 wt% to 40 wt % results in 74% decrease in strain. The stress vs. strain for collagen/CNT bio-composites is shown in Figure 3(d). At 60 MPa stress, an increase in SWCNT content from 4 wt% to 10 wt% results in 45% decrease in strain. The bio-composites with CNTs have strain values comparable to those of mineralized biocomposites without CNTs. We see at 60 MPa applied stress that the maximum strain for bio-composites with CNTs is 3.8% while for bio-composites with 20 wt% mineral the maximum strain is 4.4%. At 60 MPa stress, the bio-composites with 40 wt% mineral still had the overall minimum strain at 0.7%, however, the weight fraction of CNTs reinforcing the bio-composite is lower than the weight fraction of mineral needed to reach a similar strain value.

Elastic modulus
The Young's modulus values for the bio-composites are plotted in Figure 4. The modulus values determined in a previous study by Nair et al. (Nair et al. 2014) using a similar model are also shown in Figure 4 for comparison and are represented by pentagram-shaped data points. We compare the effect of an increase in water content from 2 to 4 wt% on the bio-composites' modulus values. Overall, as the water content increases, the modulus decreases. The trend is nonlinear to account for a saturation point where an increase in water content will have a negligible effect on the composite modulus. For compositions with 0, 20, and 40 wt% mineral, an increase in water content results in decreases in modulus by 26%, 13%, and 45%, respectively. We see that as the mineral content increases from 0 to 20 wt%, and subsequently from 20 to 40 wt %, the magnitude of the decrease in modulus increases. This is likely as the water is having a more significant effect on the dominant deformation mechanism of mineralized bio-composites (shearing of collagen between mineral phases) than that of the non-mineralized bio-composites (collagen bending). Evidence of the deformation mechanisms is discussed further in section 3.3.
An increase in SWCNT concentration from 0 to 10 wt % results in an increase in modulus. Due to small variation, the average modulus of bio-composites with 20 wt % water in Figure 4(a) is plotted in Figure 4(b) as a single value. In Figure 4(b), the modulus values for the collagen/CNT bio-composites are comparable to those of the mineralized bio-composites, ranging from approximately 1.4 to 4 GPa. As the bio-composite mineral content increases, the modulus increases. However, the presence of water reduces the degree to which mineralization changes the modulus. For example, in Figure 4 (b) at 2 wt% water, for an increase in mineral content from 0 to 40 wt% the modulus increases from 1.16 to 7.63 GPa. At 4 wt% water, an increase in mineral content from 0 to 40 wt% results in an increase in modulus from 0.86 to 4.18 GPa.
We compare the bio-composite compressive modulus values from Figure 4 to fibril compressive modulus values determined in other studies, which is shown in Table 2. Studies of highly mineralized fibrils that are dry or fully hydrated are compared to the 40 wt% mineral models in Figure 4, with 2 and 6 wt% water, respectively. Studies of non-mineralized fibrils that are dry or fully hydrated are compared to the 0 wt% mineral models in Figure 4, with 2 and 20 wt% water, respectively. The bio-composite modulus values for this study are in good agreement with the modulus values of the other studies in Table 2. The results of this study also give insight into the modulus values of partially mineralized or hydrated collagenous bio-composites. This is crucial since mineral and water content may vary due to external factors, or the mineral and water distribution within the bio-composite may not be homogeneous. This can be seen in Table 2 where, for the experimental values from other studies, there is often a significant variation in the compressive modulus values at the same mineral and water contents. The variation in compressive modulus is likely due to factors not considered in these experiments such as the amount of mineral substitution, the amount of mobile vs. bound water, and the presence of defects or non-collagenous proteins, all of which can change the bio-composite modulus. It should also be noted when comparing to the MD study by Dubey and Tomar (2009), that their mineralized fibril models were created with idealized geometric structure, which accounts for their variations in modulus values.
In Figure 5, the modulus values from this study are compared to those found using the Gao model (Ji et al. 2004) and the modified Padawer and Beecher model (Padawer and Beecher 1970). The formulation of these models are available in the S4 section of the supplementary material. The experimental values in Table 2 are also plotted in Figure 5 for comparison. As Figure 5 shows, the Young's modulus results from this MD study agree well with the experimental modulus values, as they are within the error bars, and follow a similar trend as the Gao model. However, it should be noted that at low mineral volume fractions the Gao model effective modulus reduces to zero and does not account for the effect of the remaining collagen and water on the bio-composite modulus.
The formulated model for the bio-composites (based on the Padawer and Beecher model) agrees well with the experimental results from other studies in Figure 5. Our  bio-composite model also agrees well with the simulation results, although at high mineral volume fractions the formulated model for bio-composites with 2 wt% water predicts modulus values lower than those from our MD simulations. This is likely because at higher volume fractions of mineral content, the models do not account for extrafibrillar mineralization, where mineral is also located between individual collagen fibrils, which also affects the deformation behavior. While beyond the scope of this paper, further investigation could incorporate the effect of extrafibrillar mineralization into biocomposite model. Linear regression analysis for the biocomposites determines that the linear relation between the formulated bio-composite model and the simulation results is significant and the model can accurately predict the modulus of the bio-composite (see section S4 of the supplementary material).

Gap/Overlap ratio and deformation behavior
Next, we analyze the gap/overlap ratio (shown in Figure  6) and atomic trajectories to determine the deformation mechanisms (shown in the schematics of the atomic trajectories in Figure 7). At 0 wt% mineral and under no stress, an increase in water content from 2 wt% to 4 wt% results in a 13% increase in the gap/overlap ratio, as seen in Figure 6(a). This is because at low water contents the water primarily occupies gap region voids. As the water content increases, the gap region length expands while the overlap region length remains approximately the same, as seen in Figure 7(a). As the water content increases from 4 wt% to 20 wt% the water hydrates both the gap and overlap regions, which leads to expansion in both regions and a decrease in the gap/ overlap ratio, as seen in Figure 6(a). Additionally, as the stress increases the gap/overlap ratios in Figure 6(a) decrease approximately linearly, demonstrating that the gap region deforms more than the overlap region under compressive stress. In Figure 7(a), at 0 wt% mineral and 2 wt% water, the increase in stress results in compression/bending of the collagen molecules in the gap region. As the water content increases to 4 wt% and under no stress, the water causes slight folding of collagen in the gap region. As the applied stress increases on the bio-composite with 4 wt% water, the overall gap length decreases, but the collagen molecules also straighten out, which is characteristic of shearing.
The gap and overlap regions' deformation behavior for bio-composites with 20 wt% mineral is shown in Figure 7(b). At 2 wt% water, the water is primarily in the overlap region voids due to the presence of mineral in the gap region. The gap/overlap ratio increases with increasing stress because the water results in larger compression of the overlap region, while mineral in the gap region resists compression, although this increase in gap/overlap ratio is relatively small. As the water content increases from 2 to 4 wt%, the gap/overlap ratio increases with no applied stress. This is because at 4 wt % water content the water distributes evenly throughout the bio-composite and causes gap region expansion. As stress is applied to the bio-composite with 4 wt% water, the gap/overlap ratio remains approximately constant at 1.3. The results in Figure 6(c) are also within the range of values determined in a previous study by Nair et al. (Nair et al. 2014). An increase in water content from 4 wt% to 8 wt% causes the gap/overlap ratio to decrease, since the overlap region expands due to hydration while the mineral in the gap region resists expansion. This agrees with experiments that found compressive pre-strain in the mineralized portions of hydrated collagenous bio-composites (Samuel et al. 2016). An increase in water content also led to an increase in the D-period of the simulated bio-composites. This agrees with the results of the previously mentioned experiment (Samuel et al. 2016), which also observed an increase in the D-period of mineralized fibrils as the water content increased.
For the 40 wt% mineralized bio-composites, at 2 wt% hydration the water is already uniformly distributed throughout the bio-composite. In addition, as the stress increases on the 40 wt% mineralized bio-composite samples, the gap/overlap ratio remains approximately constant, and while collagen bending due to hydration is observed, the amount of collagen bending is not significant relative to the total bio-composite size as shown in Figure 7(c). For the mineralized bio-composites under no stress an increase in water content from 2 to 4 wt% causes bending/folding of the collagen molecules in the overlap region, while the collagen in the gap region does not bend due to the minerals' resistance to deformation, as shown in Figure 7(b). For both water contents, the increase in stress results in apparent shearing of the collagen molecules. We observe this by the straightening of the collagen molecules, particularly in the overlap region of the bio-composite with 4 wt% water.
The gap/overlap ratios for bio-composites with 40 wt % mineral is presented in Figure 7(c). The results at 2 wt % water and 4 wt% water are in good agreement with the results from a previous study of a fibril with 6 wt% water in (Nair et al. 2014), with a gap/overlap ratio between 1.4 and 1.6. Even at 2 wt% water, the water distributes uniformly throughout the bio-composite. As the water content increases the gap/overlap ratio decreases, since water causes an expansion of the overlap region, while the mineral in the gap region resists expansion of the gap region. Additionally, we again observe an increase in the D-period length with an increase in water content. As the stress increases, the gap/overlap ratio stayed approximately constant.
Next, we analyze the deformation mechanisms of the collagen/CNT bio-composites. Comparing Figure 6(a,b), for 0 wt% mineral and at no applied stress an increase in CNT content from 0 wt% to 4 wt% increases the gap/ overlap ratio by 5%. This is because the CNTs cause the gap region to expand while the overlap region length remains approximately constant. As the CNT concentration increases to 7 wt%, the gap/overlap ratio decreases (Figure 6(b)). This is because, while the CNTs cause expansion of the gap region, they also displace water from the gap region to the overlap regions, which causes the overlap region to expand. At a CNT content of 4 wt%, the CNTs also displace water, but the CNTs only take up 28% of the gap region volume and the displaced water remains primarily in the gap region. At a CNT content of 7 wt%, CNTs take up 56% of the gap region volume so a significant amount of the water is displaced to the overlap region.
For an increase in CNT concentration (see Figure 6(b)) to 10 wt%, there is 11% increase in the gap/overlap ratio. At 10 wt% CNT, the CNTs occupy 83% of the gap region volume, expanding the gap region more than the overlap region expansion (which is displaced by water). Overall, as the stress increases, the gap/overlap ratio decreases. Although this shows the gap region deforms more than the overlap region, as stress increases, the decrease in gap/overlap ratio for bio-composites with CNTs is smaller than for bio-composites without CNTs. We next compare the bio-composites with 7 wt% CNTs (single-walled) and 11 wt% MWCNTs, as both cases contain 20 wt% water and 20 nm long CNTs. At no applied stress, the gap/overlap ratio of the bio-composites with MWCNTs is lower than those with SWCNTs because the larger MWCNT volume displaces additional water from the gap region. As stress increases, the gap/overlap ratio of the bio-composite with 7 wt% SWCNT decreases, while for the bio-composite with 11 wt% MWCNTs stays approximately constant. This demonstrates that the MWCNTs reduce the deformation of the gap region compared to the bio-composite with SWCNTs. At no applied stress (see Figure 7(d)) water causes folding of the collagen in the gap region. As the CNT concentration increases, this folding is reduced. As stress increases, although the overall bio-composite length decreases, the collagen molecules undergo straightening, a characteristic of shearing.
We compare the specific heat capacity for each of our models to the specific heat capacity values determined by Qu and Tomar (Qu and Tomar 2015) in Table 3. We find that the heat capacities for our models were an order of magnitude larger than those determined by Qu and Tomar. This is expected since heat capacity is dependent on the material mass, and the masses of the bio-composites in this study are an order of magnitude larger than the bio-composites studied by Qu and Tomar. For easier comparison we plot the logarithm of the bio-composite heat capacities in Figure 8. We see that an increase in mineral content increases the heat capacity due to an increase in the model mass, as observed by Qu and Tomar. A non-mineralized collagenous bio-composite with 4-10 wt% SWCNT has a 19% difference in heat capacity from a 40 wt% mineralized Table 3. Average values of the bio-composites' heat capacities.

Conclusions
We developed models of collagenous bio-composites with varying mineral, water, and CNT content and tested them under compressive loading using molecular dynamics to determine the viability of CNTs as a substitute for mineral in collagenous bio-composites. In addition, we studied the deformation behavior of the gap and overlap regions of the bio-composites and calculated their heat capacities. An increase in the mineral content increased the stiffness of the bio-composites and reduced the change in elastic modulus caused by the inclusion of water. The inclusion of CNTs in non-mineralized bio-composites also increases the elastic modulus of the bio-composites to ranges comparable to the modulus values of mineralized bio-composites. There is evidence of the shearing of collagen molecules between mineral molecules in the gap region, but at lower mineral contents the bending of collagen molecules is also observed. These observations were used to formulate a model to determine the modulus of mineralized bio-composites. The formulated model agrees with the simulations of this study, as well as modulus values determined in experimental studies by others.
For non-mineralized bio-composites at low water contents, the water occupies gap region voids, causing gap region expansion and overlap region contraction. At higher water contents some water also occupies overlap region voids, which reduces the overlap region contraction, and when stress increases, the gap region deforms more than the overlap region. In mineralized bio-composites, the water primarily occupies overlap region voids because of mineral in the gap region, which causes the overlap region to expand while the gap region contracts. As the stress increases, the gap/overlap ratio remains approximately constant. CNTs in non-mineralized bio-composites occupy gap region voids, with the amount and distribution of CNTs and water determining the relative deformation between gap and overlap regions.
The collagen/CNT bio-composites have a much higher heat capacity than the mineralized bio-composites and a lower total bio-composite mass. Overall, collagen/CNT bio-composites have stress vs. strain behaviors, elastic moduli, and gap/overlap ratios comparable to those of mineralized bio-composites, but with a higher heat capacity that would allow collagen/CNT bio-composites to be less susceptible to thermal damage during implantation or functional use in applications like bone machining, as pointed out by Qu and Tomar (Qu and Tomar 2015). This includes uses such as high-speed drilling (Wiggins and Malkin 1976), laser ablation (Nelson et al. 1988), or the curing of cements in implants (Huiskes 1980). The mechanical and thermal characterization of these collagen-based bio-composites is important for designing biocompatible bio-composite materials to treat conditions that affect collagenous tissues via applications such as wound healing, medical implants, drug delivery systems, or prosthetics.