Model overview

The goal is to develop a full-atomistic model that enables us to perform a systematic study of bone nanomechanics from a fundamental, molecular point of view. This goes beyond earlier theoretical, computational and analytical models (references and discussion see main text) that do not include the full biochemical details, and allows us to ask questions such as ‘what is the deformation mechanism in bone fibrils?’ or ‘what is the load transfer mechanism between collagen and mineral platelets?’. The geometry and composition of the model are summarized in Fig. 1b, indicating varied levels of mineral content. The 0% case corresponds to a non-mineralized collagen microfibril and also clearly shows the gap and overlap region (Fig. 1b) that emerges from the geometrical arrangement of collagen protein molecules. The D-period corresponds to the periodicity typically observed in collagen fibrils. The 20 and 40% mineral-density cases correspond to the two different mineral concentrations in the collagen microfibril. A recent study30 conducted of a pure collagen fibril has shown good agreement with experimental results. This prior work forms the basis for the model of bone, accomplished here by incorporating minerals in collagen in increasing densities (as it occurs naturally during bone formation), whereas accounting for the complexity of the chemical interactions26,45,46,47 that have not yet been incorporated in earlier models.

Collagen protein force field parameterization

Collagen is the sole protein that features hydroxyproline (HYP), a non-standard amino acid, resulting from the post-translational hydroxylation of proline. As it is rarely found, HYP is not parameterized in common biomolecular force fields, such as CHARMM48. However, a set for HYP has been developed49 based on quantum mechanical models and have been subsequently used to derive the atomistic parameters that best match the quantum-mechanics calculation, with particular focus on the correct modelling of the pucker of the HYP ring. These parameters have been then successfully incorporated into an extended CHARMM force field as used for collagen modelling in a series of previous works30,50,51. This model is also used here.

Crystal geometry and HAP force field parameterization

We generated the hexagonal HAP crystal unit cell by using Material Studio 4.4 (Accelrys, Inc.), which features the following lattice parameters: a=9.4214 Å, b=9.4214 Å, c=6.8814 Å, α=90°, β=90° and γ=120°. On the basis of this unit cell (with 44 atoms per unit cell), HAP crystals of varying size are generated for the purpose of this study, using the in silico mineralization method described below. Force fields for biomaterial simulations like CHARMM48 do not include parameters for mineral crystal, such as HAP. Therefore, to model protein–HAP composites, we extend the CHARMM force field and use bond, angle and dihedral parameters following those reported earlier52, which are based on both quantum mechanical calculations and empirical data. For non-bonded parameters, we adopt a Born–Mayer–Huggins model. For non-bonded terms, we use data from Bhowmik et al.53 in which the authors fitted the Born–Mayer–Huggins potential52 with a simpler Lennard–Jones potential.

Mineralized collagen microfibril model: in silico mineralization

The collagen microfibril model used here is based on Protein Data Bank entry 3HR2 that was used as a basis for the non-mineralized collagen fibril model reported in 30. For further information, we refer the reader to a previous paper with the details on the development of the collagen fibril atomistic model30 and to the work that describes the details of the X-ray crystallography54. The mineralized models are built starting from the non-mineralized model by filling the model’s periodic box of HAP unit cells. Then, the HAP unit cells that have at least one atom overlapping with the atoms of collagen are removed. In this way, the HAP is left only in void space of the collagen fibril box. The overlapping condition is controlled by a chosen cutoff distance between HAP atoms and collagen atoms. By selecting the appropriate cutoff distance, it is possible to add more or less HAP unit cells and, thus, reach varied degrees of mineralization. This procedure is implemented through a tcl script used in the Visual Molecular Dynamics programme55. The choice of physiologically relevant amounts of HAP is chosen based on the literature2,56, which show that the content of the mineral phase is 60–70% w/w of bone tissue. In addition to the crystallites associated with the collagen fibrils (intrafibrillar mineralization), a significant amount of mineral is believed to be located in the spaces between the fibrils (extrafibrillar mineralization). It has been reported that the extrafibrillar mineral fraction can account to as much as 75% of the total mineral in bone tissues57,58,59. Therefore, in our models, which account only for interfibrillar mineralization, we consider degrees of mineralization of up to 40%, which allows us to study moderately to highly mineralized fibrils. The orientation and surface of HAP crystals also has a significant effect in controlling the collagen–mineral interactions. Recent studies using density functional theory shows that the mineral crystal surface (that is, either a OH- or a Ca-terminated surface) has a significant role in controlling the interaction between collagen and HAP60. In our model, the mineral crystals are oriented such that orientation of the OH groups of the crystal is parallel to the collagen fibril axis as shown in Fig. 1b. This orientation of crystals and the growth of the crystals along the fibril axis typically observed in experiments56 are, hence, consistent with our mineralized collagen microfibril model. The mineral distribution inside the unit cell and its densities are shown in Fig. 2a. To enable a direct comparison with experimental results for mineral distribution, we utilize the density distribution for non-mineralized and mineralized collagen as reported in Nudelman et al.26, and plot the difference of the two data sets in comparison with our 40% mineral-density case (Fig. 2c). It is important to note that the in vivo mineralization process is a very complex phenomenon, in which immature amorphous HAP clusters nucleate close to the C-terminal end (which is at the gap-to-overlap transition region), and then grow by further mineral deposition. Furthermore, the process is believed to be assisted by acidic non-collagenous proteins. This process is too complex to be reliably replicated with direct atomistic simulations. Hence, our approach is to mimic the mineralization process by filling the void space within the collagen unit cell, in the process described above. Future works could be devoted to improve the in silico model of the mineralization process, to mimic more closely the in vivo mature mineralization state.

Fibril equilibration

Molecular dynamics simulations are performed using the LAMMPS code61 and the modified CHARMM force field described above. The constructed collagen–HAP model is first geometrically optimized through energy minimization, and then an NVT equilibration is performed for 2 ns. The unit cell that comprises collagen and HAP has triclinic symmetry as described in earlier work54. Rigid bonds are used to constrain covalent bond lengths, thus allowing for an integration time step of 2 fs. Non-bonding interactions are computed using a switching function between 0.8 and 1.0 nm for van der Waals interactions, whereas the Ewald summation method62 is applied to describe electrostatic interactions. The system temperature is maintained at 300 K (room temperature). We confirmed that the root mean square deviation of the mineralized collagen is stable, and that no major changes occur in the structure towards the end of equilibration.

In silico mechanical testing

To assess the mechanical properties of mineralized collagen fibrils, we perform stress-controlled (NP ij T) molecular dynamics simulations with increasing tensile stress applied along the x axis of the unit cell as depicted in Fig. 3a. The unit cell is under constant atmospheric pressure along other two axes (y and z). We use an NP ij T ensemble63 for loading the samples with different stress states (σ=−P ij ) ranging from: (i) atmospheric pressure, (ii) 20 MPa, (iii) 60 MPa and (iv) 100 MPa. We observe that samples reached equilibration under applied load at ~6 ns. The strain is computed as ε=(L−L o )/L o , where L o is the equilibrium length identified at atmospheric pressure. To ensure that equilibrium is achieved, we monitor the pressure at equilibrium, root mean square deviation, and confirm that the size of the simulation cell reaches a steady-state value. Using the fibril strain ε associated with each applied stress σ, we obtain the stress–strain behaviour for each case by plotting σ over ε. We use a second-order polynomial function to fit the stress–strain data, which includes all positive coefficients. The general form of the equation is σ=a 1 ε2+a 2 ε. For both non-mineralized and mineralized data sets, the coefficients a 1 and a 2 are positive so that the fitting procedure is consistent. The modulus is computed from the first derivative of a polynomial function that is fitted to the stress–strain data. The polynomial functions that are used to fit the stress–strain data are as follows: (a) for the 0% case, 2874.9x2+240.61x; (b) for the 20% case, 8101.8x2+787.89x; and (c) for the 40% case, 13651x2+1094.2x (with appropriate dimensions to express stress in MPa).The computational time requirement for equilibration is 0.2 ns per week using 24 processors for the 40% mineral-density case.

Visualization

We use Visual Molecular Dynamics55 to visualize the snapshots of simulation and compute hydrogen bonds. We count the hydrogen bonds within a cutoff distance of 3.5 Å and an angle range of 30°. MATLAB (MathWorks, Inc.) is used to plot the stress distribution inside the collagen microfibril model.

Strain analysis

We find the distance between the centre of mass of glycine residues in each chain of the collagen triple helix, and then utilize the distance between the centre of mass to compute the total deformation in collagen microfibril. This approach of computing the strain in a collagen molecule has been discussed in detail in earlier work64, and the readers are referred to this article for further information. Once the deformation is obtained for the whole microfibril, we compute the deformation of the gap and overlap region.

Theoretical models

We compare the effective modulus predicted by atomistic simulations with analytical models by Halpin–Tsai36,38 and Gao et al.37, and also with the experimental values reported in Hang and Barber31. According to the Halpin–Tsai model, the longitudinal modulus of a unidirectional plane, parallel platelet-reinforced composite is given by

where E m is the matrix modulus (here we use 0.5 GPa as the value obtained from our non-mineralized sample) and Φ is the mineral volume fraction. The constants A and B are given by

where ρ is the aspect ratio of the mineral (for biological materials, such as bone ρ~30) and E p is the modulus of the mineral platelet (100 GPa). According to Gao et al.37, the effective modulus for a nanocomposite is given by

where Φ is the mineral volume fraction, E mg is the elastic modulus of mineral (100 GPa), ρ is the aspect ratio of the mineral (for bone ρ~30) and G p =0.03 GPa is the shear modulus of collagen. The comparison of the analytical models, simulation and experiments are shown in Fig. 7.

Force field files, scripts and atomic structure files

The force-field files for LAMMPS, structure files in PDB format and all scripts used in this study are available from the corresponding author upon request.