In silico study on probing atomistic insight into structural stability and tensile properties of Fe-doped hydroxyapatite single crystals

doi:10.21203/rs.3.rs-1989927/v1

Download PDF

Article

In silico study on probing atomistic insight into structural stability and tensile properties of Fe-doped hydroxyapatite single crystals

https://doi.org/10.21203/rs.3.rs-1989927/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Hydroxyapatite (HA, Ca₁₀PO₄(OH)₂) is a widely explored material in the experimental domain of biomaterials science, because of its resemblance with natural bone minerals. Specifically, in the bioceramic community, HA doped with multivalent cations (e.g., Mg²⁺, Fe²⁺, Sr²⁺, etc.) has been extensively investigated in the last few decades. Experimental research largely established the critical role of dopant content on the changes in mechanical and biocompatibility properties. The plethora of experimental measurements of mechanical response on doped HA is based on compression or indentation testing of polycrystalline materials. Such measurements, as well as computational predictions of mechanical properties on single crystalline (doped) HA are scarce. On that premise, the present study aims to build atomistic models of Fe²⁺-doped HA, a model system, with varying Fe content (10, 20, 30, and 40 mol%) and to explore their uniaxial tensile response by means of molecular dynamics (MD) simulation, together with the calculation of IR spectrum. In the equilibrated unit cell structures, Ca(1) sites were found to be energetically favourable for Fe²⁺ substitution. The local distribution of Fe²⁺ ions significantly affects the atomic partial charge distribution and chemical symmetry surrounding the functional groups. These signatures are reflected in the significant decrease in the intensity of IR peaks found in the Fe-doped HA, together with peak splitting because of the symmetry change in the crystal structure. Another important objective of this work is to computationally predict the mechanical response of doped HA in their single crystal format. An interesting observation is that the elastic anisotropy of undoped HA was not compromised with Fe-doping. Tensile strength (TS) is systematically reduced in doped HA with Fe²⁺ dopant content and a decrease in TS with temperature can be attributed to the increased thermal agitation of atoms at elevated temperatures. The physics of the tensile response was rationalized in terms of the strain dependent changes in covalent/ionic framework (Ca-P distance, P-O bond strain, O-P-O angular strain, O-H bond distance). Further, the dynamic changes in covalent bond network were energetically analyzed by calculating the changes in O-H and P-O bond vibrational energy. Summarizing, the current work develops our foundational understanding of the atomistic phenomena involved in the phase stability and tensile response of Fe-doped HA single crystals.

Hydroxyapatite

iron-doping

MD simulation

Mechanical response

IR spectrum

Hydroxyapatite (HA,Ca₁₀PO₄(OH)₂) is the principal inorganic phase in natural bone tissue ^1,2. It is interspersed in the collagen matrix to form a composite structure, and such composite architecture allows the natural bone to withstand the loads experienced by the surrounding tissue ³. In view of its resemblance with the inorganic composition of the natural bone, HA is widely used to develop synthetic grafts and implants, aimed toward hard tissue regeneration application. ^4–7. HA can be artificially synthesized with different Ca/P ratios. The porous scaffolds of HA exhibit moderate strength and better integration with the surrounding tissue (osseointegration) ^8,9. For example, the implant surfaces coated with porous hydroxyapatite are reported to have better osseointegration ¹⁰.

While pure HA has been used extensively for implants and coatings, many studies have explored the effects of various dopants on the structure and properties of Hydroxyapatite. Traditionally, the properties of undoped HA were controlled by factors, such as grain size, porosity, surface roughness, etc. ¹¹. Doping provides the opportunity to modify the chemical and physical properties of HA. Various ions, such as Ca²⁺, PO₄³⁻, and OH⁻ ions in the apatite crystal structure can be substituted by cations (Fe²⁺/ Fe³⁺, Sr²⁺, Mg²⁺, Zn²⁺, Ag²⁺), and anions (Cl⁻, F⁻, CO₃⁻ ) ^12–17. Some dopants, like Ag⁺ are used to inhibit bacterial activity ¹⁷. Overall, the incorporation of these dopants, alone or in combination has been shown to improve bioactivity and osseointegration ^18,19. For example, Fe-doped HA, in particular, appears to improve osseointegration, bone regeneration, and hemocompatibility ^12,20,21. The magnetic behavior of iron-doped HA nanoparticles has been utilized in the field of targeted nanomedicine and MRI imaging ^22,23.

For load bearing applications, the bulk modulus of an implant plays a crucial role in its successful integration with surrounding bone tissue. Stress shielding occurs in a situation, wherein the modulus of the implant material is greater than the surrounding tissue ^24,25. The implant with the higher modulus tends to transfer much less load to the surrounding tissue. This leads to resorption of the bone tissue around the implant, eventually leading to implant failure. While doping undoubtedly improves biocompatibility and specific mechanical properties ^26,27, it is important to study the extent to which doping may improve these properties. While mechanical testing of the polycrystalline doped HA-based bioceramics has been accomplished in many lab-scale studies, the in silico prediction of the mechanical properties of the single crystalline doped HA can have many advantages, and such advantages are relatively unexplored ^28,29.

Molecular dynamics is an in-silico simulation method for studying the temporal evolution of atoms in molecular systems. It employs Newton’s equations of motion to predict the trajectory of various atoms in a system, and allows one to determine the interactions, that atoms have with each other in a given system ³⁰. MD simulations have been conducted to study the mechanical properties of HA crystals in a few studies ^31–33. For example, Qin et al. found that HA crystals with sizes up to 2 nm in a collagen fibril matrix can increase mechanical strength ³¹. Using the polymer consistent forcefield (PCFF) interatomic potentials, they studied the mechanical behavior of monoclinic HA crystals under different strain rates ³². Their study indicated anisotropic mechanical response, with the Z direction being the weakest. In compression, the application of strain in the Z direction resulted in better fracture resistance than that in X or Y direction ³². The temperature was not found to significantly influence the mechanical behavior of the HA crystal. In a different study, Synder et al. studied nanoscale mechanical deformations in hexagonal hydroxyapatite using MD simulations ³³. While using a consistent valence force field (CVFF), they too found that the fracture strength under tension was weakest along Z axis and strongest under compression. While comparing the mechanical properties of systems containing 44000 atoms and 5500 atoms, they found very little variation in properties.

While the mechanical properties of pure HA have been studied from a molecular dynamics perspective, few studies have explored the mechanical properties of doped HA.^34,35 Particularly, there are no molecular dynamics studies exploring the effect of iron doping in HA, and its mechanical response at different temperatures. As HA experiences different ranges of temperature, during manufacturing, sterilization, and implantation, it becomes important to study the temperature-dependent and compositional-dependent mechanical behaviour of Fe-doped HA. Against the above backdrop, this study aims to perform a molecular dynamics study on iron-doped HA, specifically to find the changes in the FTIR spectra and the mechanical properties in comparison to pure HA. Also, brittle materials, like HA are more susceptible to failure under tension and exhibit much better properties under compression. We present here the results obtained under tensile loading of undoped /doped HA single crystals, when strained along X, Y, or Z axes.

The unit cell of HA consists of 44 atoms with hexagonal symmetry. The lattice parameters are described as a = 9.417 Å, b = 9.417 Å, and c = 6.875 Å, with interior angles α = β = 90° and γ = 120° ^36,37. In the present study. HA unit cell was modelled using BIOVIA Materials Studio 6.0.³⁸ Fig. 1 shows a ball and stick representation of the structure. A supercell consisting of 1180 atoms was created for the study of the mechanical response. The structure was then geometrically optimized using the smart algorithm, implemented in Materials Studio package.³⁹ The optimized structure was then equilibrated first at a desired temperature for 100ps. This was followed by a second phase of equilibration at 1 atmospheric pressure for 100 picoseconds. Berendsen thermostat and barostat were employed for temperature and pressure coupling, respectively.⁴⁰ The cell time constant was kept at 1.0 ps and the decay constant was kept at 0.1 ps.

2.1 in silico tensile response analysis

The equilibrated structure was used in LAMMPS to evaluate the uniaxial tensile response, in silico ⁴¹. The PCFF-interface forcefield was used to conduct the simulations ⁴². Earlier studies have reported that the ‘Interface Forcefield’ provides reasonable accuracy in determining bulk modulus ^32,33. The tensile loading simulations were performed using a model of 1188 atoms. The periodic boundary conditions were applied in all 3 directions. The cutoff distance for Lennard-Jones and Coulombic pair potentials was set to 9.5 Å, similar to the study by Snyder et al. ³³. Ewald summation was used in all simulations, to ensure the acceptable convergence of the system ⁴³. On MD simulation platform, the system was uniaxially and virtually loaded in the X, Y, and Z directions at 10 K, 300 K, and 400 K to study the tensile properties. The rationale for the selection of different temperatures can be explained as follows. Since we have used different force field equations for undoped and Fe-doped HA, we wanted to initially assess the reliability of our results against the published results with undoped HA at 10 K.³² During synthesis, HA or doped HA experiences sub-zero temperatures during freeze-drying. For bioactivity assessment, HA/doped HA is tested for biomineralization at 300 K. Similarly, prior to biocompatibility testing, HA-based materials undergo steam autoclave-based sterilization at 400K. Therefore, it was important to understand the impact of such variation in temperature on interatomic bond configuration changes as well as the mechanical response under tension, which primarily causes the failure of brittle ceramics like HA at much lower stress or strain values, compared to that under compression. Also, a few published papers reported the compression failure of undoped HA single crystals, using MD simulations and therefore, the present in silico study is expected to add new results on HA system. ³³

In each case, the system was first equilibrated at the respective temperature for 100 picoseconds before the load application. A strain rate of 1×10¹⁰ /s was applied in each case, to study the stress response of the structure. This was then used to derive the tensile properties of the system.

2.2. in silico structural analysis

The IR spectra analysis was carried out using LAMMPS and available python code ^41,44. LAMMPS was first used to calculate the net dipole moment of the system. A python script was then used to calculate the dipole-dipole autocorrelation function and to perform the Fourier transform to obtain the IR line shape. The IR spectra were then calculated by multiplying the lineshape with an electromagnetic field-dependent prefactor.

The Fourier Transform was used to calculate the IR lineshape, which is given by the following equations:

$$I\left(\omega \right)= \frac{1}{2\pi } {\int }_{-\infty }^{\infty }{e}^{-i\left(\omega -\left(\omega \right)\right)t}\varphi \left(t\right)dt$$

...

Where, $\varphi \left(t\right)=⟨{\overrightarrow{\mu }}_{01}\left(0\right).{\overrightarrow{\mu }}_{01}\left(t\right){e}^{i{\int }_{0}^{t}\delta \omega \left(\tau \right)d\tau }⟩{e}^{-\left|t\right|/2{T}_{1}}$ ...(2)

Here, $\delta \omega \left(t\right)=\omega \left(t\right)-⟨\omega ⟩$, ${\overrightarrow{\mu }}_{01}\left(t\right)={\overrightarrow{\mu }}_{01}\left(t\right){\overrightarrow{\text{e}}}_{OH}\left(t\right)$ is the $0 \to 1$ (fundamental) transition dipole vector for the OH oscillator directed along the bond, and T₁ is the vibrational energy relaxation time.

In order to study the Fe doped hydroxyapatite structure, the Ca²⁺ ion was replaced with Fe²⁺ ion in either Ca(1) position or in Ca(2) position in the HA unit cell. Energy minimization was performed on both the structures to acquire the most stable structures together with an insight into preferable substitution sites of Fe²⁺ ions in apatite structure. The stable structures thus obtained were subjected to mechanical testing and spectral analysis, following the above-mentioned protocols. Throughout the text, xFeHA designates x mol% Fe-doped HA and the sample designation are mentioned in Table 1.

Table 1

Lattice Parameters and interior angles for single snit sell of Fe-doped structure
Structure	Designation	a (Å)	b (Å)	c (Å)	α (°)	β (°)	γ (°)
0 mol% Fe-HA	0FeHA	9.417	9.417	6.875	90	90	120
10 mol% Fe-HA	10FeHA	9.614	9.614	7.200	90	90	120
20 mol% Fe-HA	20FeHA	9.669	9.669	7.247	90	90	120
30 mol% Fe-HA	30FeHA	9.803	9.803	7.197	89.99	90	120
40 mol% Fe-HA	40FeHA	9.939	9.939	7.148	90	89.99	120

We shall analyze the results in two important aspects. First, the IR spectrum will be calculated on the basis of the structural parameters as well as the unit cell configuration for different composition of Fe-doped Hydroxyapatite. The second aim is to understand the anisotropy in the mechanical response as a function of Fe²⁺ content. Furthermore, the effects of temperature on mechanical properties will also be quantitatively assessed.

3.1 Insight into phase stability in Fe-doped HA single crystals

The initial structure of HA unit cell has been presented in Fig. 1. Primarily, the difference in energy states between the Fe²⁺ substituted structures at Ca(1) and Ca(2) positions were explored, and the final energy values, obtained after the minimization was depicted in Fig. 2. From Fig. 2, it can be easily deduced that the Ca(1) positions are slightly favourable for Fe²⁺ substitution, irrespective of the dopant content, probably because of the smaller volume of Ca(1) site ³⁴. Another important aspect is the stability of the apatite phase after Fe²⁺ doping. It can be clearly seen from Fig. 2 that, the minimum energy of the HA phase systematically decreased with an increment in Fe²⁺ content. In harmony with the experimental evidences, the substitution of Fe²⁺ ions negatively impacts the phase stability of HA⁴⁵. The energetically minimized structures (Fe²⁺ doped at Ca(1) positions) were shown in Fig. 3. All further calculations were carried out using the structures doped in Ca(1) positions, as shown in Fig. 3.

The substitution of Fe²⁺ results in the change in different crystallographic parameters and those are summarized in Table 1. A systematic increase in lattice parameters a and b with Fe²⁺ content has been observed, while a fluctuating trend was noticed along ‘c’-axis (Table 1). For instance, lattice parameter (a) increased from 9.417 Å in pure HA to 9.669 Å in 20FeHA (Table 1). A slight reduction (9.939 Å) was observed in 40FeHA (Table 1). For the lattice parameter (c), it first increased from 6.875 Å in pure HA to 7.247 Å in 20FeHA and again reduced to 7.148Å in 40FeHA (Table 1). It is interesting to note that with Fe²⁺ substitution, the ‘c’ axis changes most significantly compared to ‘a’ and ‘b’ (Table 1). For example, the change in a in 40FeHA was only 0.5% compared to pure HA, whereas the change in the c axis was found to be 3.9%. On the other hand, lattice angles remain almost constant with very marginal changes (Table 1). Hence, there are no major distortions in the lattice structure. This reconfirms the well-known ability of apatite structure to retain its crystal symmetry after doping ³⁴.

3.2 Spectroscopical analysis

The calculated IR spectra of undoped and different Fe²⁺ doped HA have been calculated from the dipole-dipole autocorrelation function and have been presented in Fig. 4. The peak positions of the calculated spectra of pure HA corroborate well with the experimental results ⁴⁶⁴⁷. The most prominent peaks in the calculated spectra in the region of 1090–1110 cm^− 1 correspond to the asymmetric stretching mode (v₃ ) of the phosphate group (Fig. 4).⁴⁷ Moreover, small peaks correspond to the symmetric bending mode (v₂ ) of the phosphate group, which can also be seen in the region of 497–533 cm^− 1(Fig. 4). Besides this, the peak at ~ 3600 cm^− 1 represents the hydroxyl stretching (Fig. 4).

The number and the intensity of the IR peaks were found to be affected by the Fe²⁺ doping inside the apatite structure. The most notable change appeared in the (v₃) vibrational mode of the phosphate group (Fig. 4). It can be easily seen that the intensity of the (v₃) peaks was greatly reduced with the Fe²⁺ doping and the lowest intensity was found for 10 and 30 mol % Fe-doped HA (Fig. 4). Together with that, prominent peak splitting was also recorded for 10 and 30 mol % Fe-doped HA.

The calculated change of the IR spectra can be explained in the qualitative manner. Among several factors, the number and intensity of the IR bands depend on the symmetry surrounding the functional group and on the partial charge distribution of the atoms. The presence of symmetry helps to reduce the number of peaks. With the substitution of Ca²⁺ ions by Fe²⁺ ions, it is observed that the symmetry around the phosphate groups along with the partial charge distribution significantly changed, as seen in Fig. 2 and Table 3. It is evident from Fig. 3 that, the symmetry surrounding the phosphate group was most affected in the case of 10 and 30 mol% of Fe²⁺ doping. This breakage of symmetry caused the reduction in the intensity together with the splitting of the v₃ vibrational peaks in these two Fe-doped HA.

3.3 Uniaxial Tensile Response of single crystals of Fe-doped HA

MD simulation was successfully implemented to predict the mechanical stress-strain response under tension. As mentioned earlier, we have conducted simulations at three different temperatures namely 10 K, 300 K, and 400K. We have first benchmarked the simulation results at 10 K with the results published in the literature by Snyder et al. and Ou et al. ^32,33. A comparison of the simulation results with the existing literature shows that the mechanical response of undoped HA along the X, Y, and Z-axis follows a similar pattern. We have considered a system consisting of 1188 atoms to conduct our simulation within a reasonable time period. However, simulations in the existing literature have used larger systems (44000 atoms) ³³. The usage of this smaller system results in significant undulations in the post-TS stress-strain response. In the current study, an important parameter viz. the tensile strength, which represents the maximum load-bearing capacity of the single crystals, has been extracted from the stress-strain response and has been analyzed with respect to both temperature and dopant concentration.

Figure 5 shows the stress-strain response of both doped and undoped HA at 10 K. The common characteristic is that the system undergoes a linear response, followed by a nonlinear response. Post yielding, there is a sharp drop in the stress values in almost all the cases, except when the system is strained along Z-axis (Fig. 5). One notable observation is the appearance of two prominent peaks in the stress-strain curve for Fe-doped HA at 10K, when loading direction was along the X axis (Fig. 5). This can be attributed to the breakage and rearrangement of the ions inside the material during the tensile test. The tensile strengths obtained from these computed stress-strain curves have been tabulated in Table 2. It can be noted that the maximum value of tensile strength (point of maximum stress) was always recorded when the load was applied along Y-direction (Table 2, Fig. 5). The tensile strength also remained lower in doped HA, compared to the pure one at 10K (Fig. 5, Table 2).

Table 2

Tensile strength of systems at different temperatures and dopant concentrations
Structure Designation	Tensile strength (Gpa)
	T = 10 K			T = 300 K			T = 400 K
	X-axis	Y-axis	Z-axis	X-axis	Y-Axis	Z-axis	X-axis	Y-axis	Z-axis
0FeHA	5.65	6.20	4.90	4.65	5.50	4.00	4.50	5.05	3.75
10FeHA	4.60	5.60	4.70	4.30	4.75	3.52	4.00	4.50	3.10
20FeHA	4.80	5.60	4.50	4.00	4.60	3.48	3.75	4.35	3.15
40FeHA	3.75	4.50	4.30	3.20	3.70	3.75	3.10	3.50	3.20

Similar stress-strain response was recorded at higher temperatures (Figs. 6 and 7). One notable difference is the presence of undulations in the stress-strain curve of pure HA with loading along the X-axis at 400K (Fig. 7). The values of tensile strength (TS) decreased systematically with an increase of temperature, independent of the loading direction and dopant content (Table 2). In most cases, the TS always remained lower in doped HA than in pure HA (Table 2).

Table 3

partial charges on different atoms in undoped and Fe-doped Ha, calculated using charge equilibration (Qeq) method.⁵⁶
Element Name	Partial changes (a.u)
Element Name	0FeHA	10FeHA	20FeHA	30FeHA	40FeHA
Ca	1.210–1.250	1.233–1.272	1.251–1.282	1.269-1.300	1.287
O(H) (Oxygen in OH⁻)	(-)0.837	(-)0.818- (-)0.816	(-)0.796- (-)0.797	(-)0.775- (-)0.776	(-)0.754
O(P) (Oxygen in PO₄³⁻)	(-)0.562- (-)0.590	(-)0.560- (-)0.601	(-)0.537- (-)0.581	(-)0.546- (-)0.573	(-)0.527- (-)575
P	0.480	0.494–0.499	0.513	0.527–0.533	0.546–0.547
H	0.230	0.236–0.243	0.242-255	0.254–0.261	0.267
Fe		0.735	0.750	0.752–0.773	0.775

In our efforts to develop a deeper insight into the physics of tensile response, we shall now analyze the atomistic insights and the rearrangement or restructuring of the covalent bonds in doped HA single crystals.

4.1 Atomistic insights into tensile response of single crystalline doped HA

In order to obtain atomistic insight into the stress-strain response, the distance between Ca²⁺/Fe²⁺ and phosphate ions has been calculated. For each case, the distance was found to be monotonously increasing with the applied strain, almost in a linear manner (Fig. 11, Fig. S4 and S5 in supporting information). This increment of the ionic distance is the major contributor to the mechanical response of doped HA. Moreover, it is clear from the above trends that the tensile response of undoped and Fe-doped HA strongly depends on the loading direction. Different atomic arrangement along the orthogonal crystallographic directions is the underlying reason behind the observed anisotropic nature of HA. Profoundly, different sets of atoms are displaced in the crystal during the mechanical response, depending on the loading directions. For instance, in the case of pure HA, OH⁻ ions are observed to move apart, during loading along the X-axis (Fig. 9). This indicates the effective reduction of the attractive interaction between OH⁻ ions and Ca²⁺ ions, which ultimately led to material failure. Similarly, the interactions among Ca²⁺ ions and PO₄³⁻ /OH groups were noted to be reduced, when the load was applied along the Y-direction (Fig. 10). As more ionic bonds were broken along Y-direction, the tensile strength was higher along the Y-axis (Table 2).

On the other hand, the mechanical response of pure HA along the Z-direction strongly depends on the ionic bonds between Ca²⁺ ions and PO₄³⁻ groups, which started to break during tensile testing (Fig. S2, S3 in Supporting information). The above-mentioned trend is observed independent of the temperature and the doping amount of HA.

Besides this, the redistribution of the atomistic charges in the presence of the foreign metal ions also alters the mechanical behaviour, which has been substantiated further in the present work. When Fe²⁺ ions were introduced in the apatite structure, the partial charge distribution of the constituent atoms changed significantly (Table 3). Partial charges on Ca²⁺, H⁺ and P⁵⁺ ions generally increased with an increase in Fe^+ 2 content, whereas the same on O(H) generally decreased (Table 3). On the other side, partial charges on O(P) atoms fluctuate as a function of dopant content (Table 3). The observed decrease in partial charges on O(H) resulted in the change in the attractive electrostatic interactions among Ca²⁺ and O(H) atoms, while the increased charges on Ca and H atoms enhanced the repulsive coulombic forces (Table 3). The effects of these changes got reflected in the decreasing nature of the tensile strength with dopant content (Table 2). It is also evident from Table 3 that the Fe²⁺ ions carry lesser partial charges than Ca²⁺ ions. Hence, the presence of Fe²⁺ ions gave rise to a weaker electrostatic attractive interaction among Fe²⁺- PO₄³⁻ or Fe²⁺- OH⁻ ions compared to that among the Ca²⁺- PO₄³⁻ or Ca²⁺- OH⁻ groups, present in pure HA. This also promotes the recorded decline in the tensile strength in doped HA.

As seen in Figs. 5, 6, and 7 and also in Table 2, the temperature also negatively impacts the tensile strength. With an increase in the temperature, the atomic/molecular vibration increases because of the enhanced thermal movement. This, in turn, affects the tensile strength in a negative manner. Other than the above-mentioned bond breakage at TS point, the deformations of several other bonds are also responsible for the linear response of HA under tensile loading (see Fig. 9, 10 and Fig. S1-S3 in Supporting information). The strength of these bonds also changed with doping content and temperature. The synchronized influence of these changes gave rise to the non-monotonous behaviour of the tensile strength.

In order to quantitatively analyze the differences in the mechanical response of doped HA in different loading directions, we have calculated the anisotropic variation in the tensile strength for undoped and Fe-doped HA, at different temperatures. In analyzing the anisotropic strength behaviour, Fig. 8 shows the ratios of tensile strength plotted as σ_z/σ_x and σ_z/σ_y. A non-monotonous trend was observed for tensile strength (Fig. 8). A general observation is that the anisotropic nature of the tensile strength slightly decreased in the doped HA with an increase in temperature (Fig. 8). Like before, all these observations can be attributed to the altered partial charge distribution among the atoms in the presence of a foreign metal ion (Fe).

The anisotropic nature of the bone tissue is reported in the literature ⁴⁸. This arises from the preferential orientation of the collagen fibers and mineral crystals (calcium phosphate like HA) along the bone growth direction ⁴⁹. The anisotropic property of natural bone helps to maximize the stiffness and strength along the major load bearing axis, while keeping the bone mass minimum ⁵⁰. Hence, an ideal bone implant material should exhibit elastic anisotropy. HA is known to exhibit anisotropic mechanical response [49], and our study establishes that the anisotropic nature of HA is retained after Fe²⁺ doping. Moreover, the elastic modulus of Fe²⁺-doped HA (data not shown) is close to that of cortical bone (7-30GPa) ⁵¹, which establishes its efficacy as bone implant material.

4.2 Rearrangement of covalent bonds during tensile loading

To study the effects of the tensile force on the covalent bonds, bond strain and the change in vibrational bond energy have been calculated.

First, we analyze the changes in the characteristic properties of covalent bonds, (P-O) due to the application of the tensile force at both 10 K and 300 K. We have computed the bond strain and the effective change in its vibrational bond energy as a function of applied strain. These quantities are calculated using the following expressions,

Δr = r-r₀ (1)

ΔE = 0.5×k(r² - r₀² ) (2)

where k is the force constant and r₀ is the equilibrium bond length.

The average bond strain of the P-O bonds, present in pure HA at 10K has been presented in Fig. 12. The P-O bonds were found to be contracted at the beginning of the application of the tensile force and then stretched with the increment of the applied force (Fig. 12). After the TS point, the bonds again contracted and remained in that state for the rest of the simulation time (Fig. 12). In corroboration with the stress-strain curve (Figs. 5–7), the above-mentioned trend is most prominent when the force was applied along X and Y-axis. On the other hand, the O-P-O angle slightly decreased with the application of the tensile force, and then it maintained a steady state from 0.05% strain onward (Fig. 13).

The changes in the P-O vibrational bond energy at 10K took place in complete harmony with the bond strain (Fig. 14). For X and Y-loading directions, the change was maximum initially and then reached minima near the TS point. Beyond that, the vibrational bond energy was again moved further away from its value at the equilibrium position, r₀ (Fig. 14). The bond energy under stress was recorded to be closer to the equilibrium bond energy value near TS point, when the force was applied along Y-axis, compared to that for other loading directions (Fig. 14). In agreement with the stress-strain curve (Figs. 5–7), the observed pattern of the bond energy change under stress was significantly different for Z- axis (Fig. 14). The anisotropic nature of the HA crystal is the underlying reason behind the differences recorded in the change in bond energy under stress. Because of the different atomic arrangements along the various crystallographic directions, the P-O bonds were experiencing significantly divergent resultant forces, when exposed to the same external stress along the orthogonal loading direction. Such phenomenon results in distinct anisotropic behaviour, as presented in Fig. 14.

When the temperature was elevated to 300 K, the general trend of both the bond strain and change in bond vibrational energy for pure HA remains identical with additional oscillations (Fig. S6 in Supporting Information). This is attributed to the thermal motion at high temperature. At both 10K and 300K, the behaviour of P-O bond strain/bond energy of Fe-doped HA under external tensile stress remains mostly identical to that of pure HA at 10K and 300K, respectively (Fig. S7 in Supporting Information). The only significant change observed in P-O bond strain was for 40FeHA at 10K with loading direction along the Z axis (Fig. 12). Throughout the tensile loading, the P-O bonds were found to be expanded, probably because of the high Fe²⁺ content (Fig. S7 in Supporting Information). Similar behaviour of P-O bonds for 40FeHA was absent at 300K due to the thermal agitation (Fig. S7 in Supporting Information).

The signature of the bond deformation with increased temperature can be found in the O-P-O angular strain in terms of oscillatory nature as well (Fig. S8 in Supporting Information). At 10K, the angular strain became much smaller for Fe-doped HA, compared to pure HA (Fig. 13), whereas similar angular stain was recorded for each material at 300K (Fig. S8 in Supporting Information).

As far as the O-H bond strain of the pure HA is concerned, it exhibited oscillatory behaviour under the external stress at 10K (Fig. 15), in corroboration with the change of O-H bond vibrational energy (Fig. 16). It is seen from the obtained atomistic trajectory during the tensile loading that the OH⁻ groups underwent orientational changes under the tensile stress. The oscillatory bond strain is a signature of the force experienced by the atoms of the hydroxyl groups during the orientational change. The frequency of oscillation increased in the presence of Fe²⁺ ions at 10K, because of the altered interactions among the hydroxyl groups and metal ions (Fig. 15). For pure and Fe-doped HA, the amplitude of oscillation increased together with a decrease (reduction) in the oscillation frequency as well as with the increment of temperature, because of thermal motion (Fig. 15).

On the other hand, the oscillatory nature of the change in the O-H bond vibrational energy is more prominent at 300K for 0FeHA and 20FeHA, because of the greater change in its bond length at higher temperatures (Fig. S9 and S10 in Supporting Information). For 40FeHA, the oscillatory nature is visible in both 10 K and 300 K, although the amplitude is larger at 300K, due to greater (larger) thermal vibration (Fig. S10 in Supporting Information). Summarizing the results in this section, we can conclude that although the application of external tensile stress results in the change in both covalent bond length (Δr) and its corresponding vibrational energy (ΔE) in the apatite structure, such changes do not appear to result in the bond breakage. Moreover, the presence of dopant ions and the change in temperature also significantly influence the behaviour of the covalent bonds in apatite structure.

Most importantly, our study has provided an atomistic comprehension to the tensile response of single crystalline HA. On the experimental aspects, Gupta and co-workers performed a series of synchrotron-based studies to quantitatively analyze the deformation mechanism of natural bone, which is a composite of collagen and HA^52–54. Since the elastic stiffness or strength properties largely depend on HA content in natural bone, we believe that the computational output from studies like the present study can be used further to validate the synchrotron-based experimental results. Also, synchrotron-based in situ study on single crystalline HA can be performed in future.

At the closure, our study shed light to the fracture mechanism of the single crystalline doped HA using a model system, which will be helpful to understand the planned experimental results. Similar in silico studies can be used to generate significant data for biomaterialomics approach.⁵⁵

In the current study, we have successfully implemented MD simulation platform to understand the structural changes in Fe-doped HA and to establish an atomistic level understanding of the dopant-dependent changes in ionic arrangement together with its impact on the mechanical response. Following conclusions can be drawn from the present work.

Between two crystallographically non-equivalent Ca sites, Ca(1) sites are slightly preferable over Ca(2) sites for Fe²⁺ substitution inside the apatite structure. The lattice parameters significantly change along the ‘c’ axis after doping.
The presence of Fe²⁺ ions inside the doped HA lattice restricts the chemical symmetry around the functional groups and this corroborates well with the changes in the number and intensity of the calculated IR bands.
The tensile strength of HA decreases with Fe²⁺ doping because of the changes in the atomic charge distribution inside the crystal structure. The elastic modulus of Fe-doped HA is comparable with the natural cortical bone. An increase in temperature induces higher bond vibration and thermal movement, which significantly decreases the tensile strength of HA.
The tensile deformation of Fe-doped HA at the atomistic level has been quantitatively analyzed using the dynamic changes in the major physical parameters of the covalent/ionic bond framework, including Ca-P distance, P-O bond strain, O-P-O angular strain, O-H bond distance with strain. Also, the dynamic changes in vibrational bond energy for O-H and P-O bonds are evaluated at various temperature.
No significant breakage of the covalent bonds in the apatite structure is recorded during the tensile test. O-H bond strain exhibited oscillatory nature under tensile load, in a manner dependent on Fe²⁺dopant content.

Taken together, the current study presents a molecular-level insight into the altered structure and mechanical behaviour and establishes the efficacy of MD simulation platform to gain deeper atomistic insights into the mechanical response of single crystalline doped HA.

Acknowledgements

One of the authors, BB acknowledges the financial support provided by the Abdul Kalam National Innovation Fellowship of Indian National Academy of Engineering, and ongoing research projects, “Investigation on Machine learning approaches for optimization of process parameter for additive manufacturing and accelerated design of patient specific hip implants, based on Finite Element Analysis’ (CRG/2020/001145). SN acknowledges Prof. S. Yashonath, Solid State and Structural chemistry Unit, IISc, Bangalore for his useful suggestions.

Data Availability

The data that support the findings of this study are submitted to the materials Cloud archive (https://archive.materialscloud.org/) with DOI: https://doi.org/10.24435/materialscloud:13-sa. (Under Review).

Rey, C., Combes, C., Drouet, C. & Glimcher, M. J. Bone mineral: update on chemical composition and structure. Osteoporos. Int. 20, 1013–1021 (2009).
Hueb De Menezes Oliveira, M. A. et al. Microstructure and mineral composition of dental enamel of permanent and deciduous teeth. Microsc. Res. Tech. 73, 572–577 (2010).
Fratzl, P., Gupta, H. S., Paschalis, E. P. & Roschger, P. Structure and mechanical quality of the collagen–mineral nano-composite in bone. J. Mater. Chem. 14, 2115–2123 (2004).
Kuriakose, T. A. et al. Synthesis of stoichiometric nano crystalline hydroxyapatite by ethanol-based sol–gel technique at low temperature. J. Cryst. Growth 263, 517–523 (2004).
Lin, K., Chang, J., Cheng, R. & Ruan, M. Hydrothermal microemulsion synthesis of stoichiometric single crystal hydroxyapatite nanorods with mono-dispersion and narrow-size distribution. Mater. Lett. 61, 1683–1687 (2007).
Gao, C. et al. Current progress in bioactive ceramic scaffolds for bone repair and regeneration. Int. J. Mol. Sci. 15, 4714–4732 (2014).
Szcześ, A., Hołysz, L. & Chibowski, E. Synthesis of hydroxyapatite for biomedical applications. Adv. Colloid Interface Sci. 249, 321–330 (2017).
Sopyan, I., Mel, M., Ramesh, S. & Khalid, K. A. Porous hydroxyapatite for artificial bone applications. Sci. Technol. Adv. Mater. 8, 116–123 (2007).
Tripathi, G. & Basu, B. A porous hydroxyapatite scaffold for bone tissue engineering: Physico-mechanical and biological evaluations. Ceram. Int. 38, 341–349 (2012).
Xue, W. et al. Preparation and cell–materials interactions of plasma sprayed strontium-containing hydroxyapatite coating. Surf. Coatings Technol. 201, 4685–4693 (2007).
Maté Sánchez de Val, J. E. et al. Influence of hydroxyapatite granule size, porosity, and crystallinity on tissue reaction in vivo. Part A: synthesis, characterization of the materials, and SEM analysis. Clin. Oral Implants Res. 27, 1331–1338 (2016).
Mercado, D. F. et al. Paramagnetic iron-doped hydroxyapatite nanoparticles with improved metal sorption properties. A bioorganic substrates-mediated synthesis. ACS Appl. Mater. Interfaces 6, 3937–3946 (2014).
Capuccini, C. et al. Interaction of Sr-doped hydroxyapatite nanocrystals with osteoclast and osteoblast-like cells. J. Biomed. Mater. Res. Part A 89A, 594–600 (2009).
Stanić, V. et al. Synthesis, characterization and antimicrobial activity of copper and zinc-doped hydroxyapatite nanopowders. Appl. Surf. Sci. 256, 6083–6089 (2010).
Rulis, P., Ouyang, L. & Ching, W. Y. Electronic structure and bonding in calcium apatite crystals: Hydroxyapatite, fluorapatite, chlorapatite, and bromapatite. Phys. Rev. B 70, 155104 (2004).
Tampieri, A., Celotti, G. C., Landi, E. & Sandri, M. Magnesium doped hydroxyapatite: synthesis and characterization. in Key Engineering Materials 264, 2051–2054 (Trans Tech Publ, 2004).
Riaz, M. et al. Synthesis of monophasic Ag doped hydroxyapatite and evaluation of antibacterial activity. Mater. Sci. Eng. C 90, 308–313 (2018).
Vranceanu, D. M. et al. Osseointegration of sputtered SiC-added hydroxyapatite for orthopaedic applications. Ceram. Int. 42, 10085–10093 (2016).
Wang, J.-Y. et al. Flame-sprayed strontium- and magnesium-doped hydroxyapatite on titanium implants for osseointegration enhancement. Surf. Coatings Technol. 386, 125452 (2020).
Gloria, A. et al. Magnetic poly(ε-caprolactone)/iron-doped hydroxyapatite nanocomposite substrates for advanced bone tissue engineering. J. R. Soc. Interface 10, 20120833 (2013).
Sarath Chandra, V. et al. Blood compatibility of iron-doped nanosize hydroxyapatite and its drug release. ACS Appl. Mater. Interfaces 4, 1200–1210 (2012).
Panseri, S. et al. Magnetic labelling of mesenchymal stem cells with iron-doped hydroxyapatite nanoparticles as tool for cell therapy. J. Biomed. Nanotechnol. 12, 909–921 (2016).
Ribeiro, T. P., Monteiro, F. J. & Laranjeira, M. S. PEGylation of iron doped hydroxyapatite nanoparticles for increased applicability as MRI contrast agents and as drug vehicles: A study on thrombogenicity, cytocompatibility and drug loading. Eur. Polym. J. 137, 109934 (2020).
Huiskes, R., Weinans, H. & van Rietbergen, B. The relationship between stress shielding and bone resorption around total hip stems and the effects of flexible materials. Clin. Orthop. Relat. Res. 124–134 (1992).
Nagels, J., Stokdijk, M. & Rozing, P. M. Stress shielding and bone resorption in shoulder arthroplasty. J. Shoulder Elb. Surg. 12, 35–39 (2003).
Mansour, S. F., El-dek, S. I., Dorozhkin, S. V & Ahmed, M. K. Physico-mechanical properties of Mg and Ag doped hydroxyapatite/chitosan biocomposites. New J. Chem. 41, 13773–13783 (2017).
Vladescu, A. et al. Mechanical properties and biocompatibility of the sputtered Ti doped hydroxyapatite. J. Mech. Behav. Biomed. Mater. 63, 314–325 (2016).
Butler, K. T., Frost, J. M., Skelton, J. M., Svane, K. L. & Walsh, A. Computational materials design of crystalline solids. Chem. Soc. Rev. 45, 6138–6146 (2016).
Lau, D., Jian, W., Yu, Z. & Hui, D. Nano-engineering of construction materials using molecular dynamics simulations: Prospects and challenges. Compos. Part B Eng. 143, 282–291 (2018).
Andersen, H. C. Molecular dynamics simulations at constant pressure and/or temperature. J. Chem. Phys. 72, 2384–2393 (1980).
Qin, Z., Gautieri, A., Nair, A. K., Inbar, H. & Buehler, M. J. Thickness of Hydroxyapatite Nanocrystal Controls Mechanical Properties of the Collagen–Hydroxyapatite Interface. Langmuir 28, 1982–1992 (2012).
Ou, X. & Han, Q. Molecular dynamics simulations of the mechanical properties of monoclinic hydroxyapatite. J. Mol. Model. 20, 1–12 (2014).
Snyder, A. D. & Salehinia, I. Study of nanoscale deformation mechanisms in bulk hexagonal hydroxyapatite under uniaxial loading using molecular dynamics. J. Mech. Behav. Biomed. Mater. 110, 103894 (2020).
Basu, S. & Basu, B. Doped biphasic calcium phosphate: synthesis and structure. J. Asian Ceram. Soc. 7, 265–283 (2019).
Ullah, I. et al. Simultaneous co-substitution of Sr2+/Fe3 + in hydroxyapatite nanoparticles for potential biomedical applications. Ceram. Int. 44, 21338–21348 (2018).
Hughes, J. M. & Rakovan, J. The crystal structure of apatite, Ca5(PO4)3(F,OH,Cl). Phosphates Geochemical, Geobiol. Mater. Importance 48, 1–12 (2019).
Kay, M. I., Young, R. A. & Posner, A. S. Crystal structure of hydroxyapatite. Nature 204, 1050 (1964).
BIOVIA, D. S. Materials Studio, version 6.0. San Diego, CA, USA (2012).
Sharma, S. et al. Overview of BIOVIA materials studio, LAMMPS, and GROMACS. Mol. Dyn. Simul. Nanocomposites using BIOVIA Mater. Stud. Lammps Gromacs, Amsterdam, Netherlands Elsevier 39–100 (2019).
Berendsen, H. J. C., Postma, J. P. M., Van Gunsteren, W. F., Dinola, A. & Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81, 3684–3690 (1984).
Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 117, 1–19 (1995).
Heinz, H., Lin, T.-J., Kishore Mishra, R. & Emami, F. S. Thermodynamically Consistent Force Fields for the Assembly of Inorganic, Organic, and Biological Nanostructures: The INTERFACE Force Field. Langmuir 29, 1754–1765 (2013).
Ewald, P. P. Die Berechnung optischer und elektrostatischer Gitterpotentiale. Ann. Phys. 369, 253–287 (1921).
Braun, E. Open Source Code: Calculating an IR Spectra from a LAMMPS Simulation. (2016). doi:10.5281/zenodo.154672
Basu, S. & Basu, B. Unravelling Doped Biphasic Calcium Phosphate: Synthesis to Application. ACS Appl. Bio Mater. 2, 5263–5297 (2019).
Basu, S., Ghosh, A., Barui, A. & Basu, B. (Fe/Sr) Codoped Biphasic Calcium Phosphate with Tailored Osteoblast Cell Functionality. ACS Biomater. Sci. Eng. 4, 857–871 (2018).
Rehman, I. & Bonfield, Wjj. Characterization of hydroxyapatite and carbonated apatite by photo acoustic FTIR spectroscopy. J. Mater. Sci. Mater. Med. 8, 1–4 (1997).
Bonfield, W. & Grynpas, M. D. Anisotropy of the Young’s modulus of bone. Nature 270, 453–454 (1977).
Novitskaya, E. et al. Anisotropy in the compressive mechanical properties of bovine cortical bone and the mineral and protein constituents. Acta Biomater. 7, 3170–3177 (2011).
Kang, J. et al. Anisotropy characteristics of microstructures for bone substitutes and porous implants with application of additive manufacturing in orthopaedic. Mater. Des. 191, 108608 (2020).
Basu, B. Biomaterials for Musculoskeletal Regeneration. (Springer, 2017). doi:10.1007/978-981-10-3059-8
Hang, F., Gupta, H. S. & Barber, A. H. Nanointerfacial strength between non-collagenous protein and collagen fibrils in antler bone. J. R. Soc. Interface 11, 20130993 (2014).
Gupta, H. S. et al. Intrafibrillar plasticity through mineral/collagen sliding is the dominant mechanism for the extreme toughness of antler bone. J. Mech. Behav. Biomed. Mater. 28, 366–382 (2013).
Karunaratne, A. et al. Symmetrically reduced stiffness and increased extensibility in compression and tension at the mineralized fibrillar level in rachitic bone. Bone 52, 689–698 (2013).
Basu, B., Gowtham, N. H., Xiao, Y., Kalidindi, S. R. & Leong, K. W. Biomaterialomics: Data science-driven pathways to develop fourth-generation biomaterials. Acta Biomater. 143, 1–25 (2022).
Rappe, A. K. & Goddard III, W. A. Charge equilibration for molecular dynamics simulations. J. Phys. Chem. 95, 3358–3363 (1991).

No competing interests reported.

supportinginfo210822.docx

Download PDF

Editorial decision: Major revision
12 Oct, 2022
Reviews received at journal
29 Sep, 2022
Reviews received at journal
07 Sep, 2022
Reviewers agreed at journal
07 Sep, 2022
Reviewers agreed at journal
01 Sep, 2022
Reviewers invited by journal
01 Sep, 2022
Editor assigned by journal
01 Sep, 2022
Editor invited by journal
01 Sep, 2022
Submission checks completed at journal
01 Sep, 2022
First submitted to journal
23 Aug, 2022

You are reading this latest preprint version

In silico study on probing atomistic insight into structural stability and tensile properties of Fe-doped hydroxyapatite single crystals

Status:

Version 1

Abstract

Figures

1. Introduction

2. Methods

2.1 in silico tensile response analysis

2.2. in silico structural analysis

3. Results

3.1 Insight into phase stability in Fe-doped HA single crystals

3.2 Spectroscopical analysis

3.3 Uniaxial Tensile Response of single crystals of Fe-doped HA

4. Discussion

4.1 Atomistic insights into tensile response of single crystalline doped HA

4.2 Rearrangement of covalent bonds during tensile loading

5. Conclusions

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1