On Estimation of Adhesive Strength of Implants Bioactive Coating with Titanium by Density Functional Theory and Molecular Dynamics Simulations

One of the ways to improve and accelerate osseointegration of a surgical implant with bone is application of biocompatible coatings, in particular, hydroxyapatite (HAp). Since the cases of delamination of the coating take place in dental practice, it is very important to estimate the adhesive strength of HAp with the implant. A measure of the coating-to-substrate bond strength is the energy of this bond. In this research, quantum chemistry is used to calculate the binding energy of functional groups (anions) of hydroxyapatite and titanium 2+, which is a standard implant material. First, using Density Functional Theory with Becke three-parameter Lee-Yang-Parr hybrid exchange-correlation functional, the lowest potential energy surface is calculated. Then, by ab initio molecular dynamics, the reaction path, the reaction products, frequencies of oscillations, the activation energies and binding energies between various combinations of component anions of HAp and Ti(II) are calculated.


Introduction
The key requirement of dental implantation is strong osseointegration (fusion of an implant with a bone - Figure  1) 1,2 . To achieve this, implants are often covered with bioactive coatings 3,4,5 . One of the most popular materials is hydroxyapatite (HAp), with the molecular formula Ca 10 (PO 4 ) 6 (OH) 2 ( Figure 2) 6,7,8 . The composition of HAp is similar to the mineral constituent of bone. The nanocrystalline structure of HAp gives a microrelief of the surface, favorable for osseointegration [9][10][11][12][13] .
However, such coatings showed a number of drawbacks: cases of peeling of the coating from the titanium base 15,16 , a moderate rate of osseointegration, and exposure to the environment. An important task therefore, among others, is to study the bond strength of the HAp coating with titanium 17,18 . The strength characteristics of the bond-substrate interaction is the energy of this bond 19 . The aim of the present work is to determine the binding energies between hydroxyapatite functional groups (anions) and Ti(II) (titanium cation with plus 2 charge) -the standard material for implants. This is accomplished with computational methods. In the future, these components will be used to calculate the total binding energy of the unit cell HAp and the cation Ti(II).

General computational method
In computational chemistry, the Density Functional Theory (DFT) is used to describe the state of many-particle systems and to determine their geometric and chemical properties 20,21 . In this paper it is used to determine the ground state energies of polyatomic complexes in the Ti (II) -hydroxyapatite system. The ultimate goal is the theoretical calculation of the binding energy of the HAp coating and Ti(II).
Using the formulas obtained with the help of DFT, we calculate the ground-state energies in the computational chemistry software suite, Gaussian 09 (Revision C.01 was used 22 ). The suite allows us to perform geometric optimizations of structures of interest, which result in obtaining incredibly specific nuclei coordinates in three-dimensional space. These geometry optimizations give the necessary positions of nuclei that deliver a global minimum on the potential energy surface. The program calculates and interprets the partial differential equations that deliver the solution with minimum energy to the Schrodinger equation. All calculations for the polyatomic complexes in the Ti (II) -hydroxyapatite system were performed using the DFT functional B3LYP (Becke three-parameter Lee-Yang-Parr functional) in Gaussian 09 system.

The time independent Schrödinger equation
The time independent (general) Schrödinger equation describes the interactions of N particles. It states that when the Hamiltonian, Ĥ, acts on a wave function, Ψ, it is proportional to the stationary state of the same wave function with a proportionality constant E. For our HAp calculations, we determine the lowest possible energy that the systems can have: (1) where E is the total energy for the system, and Ĥ is a Hamilton operator for a molecular system containing N nuclei and electrons: (2) where T is a kinetic energy and V is the potential energy, or where the two first terms are kinetic energies for electrons and nuclei, and the last three terms are potential energies of nuclei-electron attractions, nuclei-nuclei repulsions and electron-electron repulsions correspondingly.

Density functional theory and Born-Oppenheimer approximation
In Born-Oppenheimer (BO) approximation, which is an important aspect of Density Functional Theory (DFT), we assume that the kinetic energy of electrons is much greater than the kinetic energy of nuclei, and because of that the nuclei are believed to be fixed 21 . Their kinetic energy is zero, and the potential energy due to nucleus-nucleus repulsion is a constant. From this concept, for each arrangement of nuclei, the ground-state energy, E, (or full energy of a manyparticle system) can be found on the base of DFT. Thus, at each moment in time the system is settled at the ground state energy, or at absolute zero, and therefore the velocity of nuclei are close to zero and the kinetic energy of nuclei is negligible small ( T nuc = 0). Then the expression (3) becomes (4) or is called the Electronic Hamiltonian. It is easy to see that, in this case, in the Schrödinger equation (5) the eigenvalues E eff will give the ground state energy of the nuclei, or ab-initio lowest energy E min , and that is a unique function of R v , the coordinates of the nuclei. The ab initio lowest total energies can be evaluated for many configurations, which form the minimum energy surface, or ab-initio potential-energy surface (see Figure 3), for a given reaction 23-28 . According to DFT, the reactants and products are both at the positions of the absolute minima on the potential energy surface. From energy principals, the path between reactants and products goes in the direction of minimum spending energy, which is viewed as a path through a saddle point between two global minima on the potential energy surface (see Figure 3). The difference between the saddle point energy and the reactant energy will give the activation energy barrier, E act , or the energy that is needed to be overcome to form a bond, and the difference in reactants energy and the product energy, ΔG, will give the binding energy, or the energy stored in the product compared with the reactants. In this paper, the subject of our interest is the binding energy only. In the future however, we will be analyzing our results for the activation energies as well, which are related with , , ... , , ...
the manufacturing technologies and processes for spattering nano-coatings on the implant.

Density functional theory basics (electronic optimization)
and taking into account the normalization of wave functions, this expression can be rewritten (6) which means the minimization of quadratic functional with respect to Ψ elec will give us the true energy minimum, E eff . In Density Functional Theory Energy ρ (r), the electron density is introduced as (7) Then the quadratic functional taking into account expression for the Hamiltonian (4) can be written as (8) Hohenberg-Kohn Theorems 21 : is a unique functional of ρ (r). HK2 Functional minimized with respect to ρ (r) will give true ground state energy (which means E min ), or E[ ρ ] reaches its minimum value for the real ground state density, ρ (r).
In the expression (7) Ψ i is the Kohn-Sham orbital satisfying the following orbital equation (9) with V eff as the Kohn-Sham effective potential that depends in turn on the electron density. The expression (8) for total energy can be expressed as (10) where (11) and (12) E xc is the exchange-correlation functional and (13)

Calculation method
This paper represents the first step in solving the general problem of determining the binding energy of HAp with titanium. We determine the binding energies between individual fragments (anions) of hydroxyapatite and titanium. In general, a Ti 2+ atom was paired with various anions from HAp (see Figure 2). In Gaussian, it requires the initial "guess" , set up of positions of the atoms. The more complex the final structure is, the more important the initial guess. This is due to the fact that as the geometry gets more complicated, the charge distribution among the atoms in the structure plays an important role in the structure of the final optimized geometry. The charge distribution may cause atoms in the final geometry to repel and attract one another, resulting in particular geometries (e.g. angles and bond lengths). First, by using the B3LYP hybrid-exchange-correlation functional, we calculated ground state energies for different configurations of hydroxyapatite, namely for OH -, PO 4 3-(which are anions), and then for Ti 2+ (which we took with charge 2+ to have attraction between ions) ( Figure 2). Then, we calculated various combinations of Ti 2+ , OH -, PO 4 3-(that is when the initial guesses become more important). Once these energies were found, we calculated the binding energies (BE) by subtracting the ground state energies of products from the ground state energies of reactants.
To achieve convergence during calculations, we used fine integration grid and basis set 3-21G during optimization. Also, frequency calculations were performed before geometric optimization. If frequencies are negative, the initial configurations are not correct and ground state energy is not found and the optimized position of nuclei is not at global minimum. If that is the case, we need to figure out the new initial guess configurations and repeat the optimization process from the beginning. For all the combinations considered, the main energy levels and binding energies are calculated, with particular emphasis on the equilibrium bond lengths and the angles of Ti-O bonds. In the next stage of work, this data will be used in calculating the total binding energy of the anion Ti (II) and the whole unit cell HAp. Depending of the complexity of the final configuration, the calculation time varied from several minutes to several hours.

Results of calculations
The results are shown in Fig. 4-8. For example, let us briefly summarize the general results of one of the calculations. According to calculations of the complex [TiOH]+ ( Figure  4), the charge of the Ti (II) ion decreases. In terms of the Density Functional Theory (DFT), this means that Ti (II) experiences an increase in the electron density due to its attraction with negatively charged hydroxide ions, OH -. The length of the reference bond is about 1.8 Å. For the vibrational frequency, which indicates as interaction between Ti (II) and the oxygen atom, the calculation gives a reasonable value of about 970 cm -1 . The Hessian index contains all positive eigenvalues of a constant-force matrix and corresponds to real vibrational frequencies of bonds within molecules, or in this case, interacting ions. In the first structure, the Titanium (II) ion experiences a decrease of almost half of its initial charge and demonstrates an expected charge distribution. The initial charge of Ti(II) is 2+, however, after its reaction with the negatively charged hydroxide (the oxygen atom in particular), the Titanium ion loses approximately .876 charge (Table 1 and Figure 4). This information provides several implications about the function of titanium in the [TiOH] + structure and its centralized relationship with the oxygen atom. The orbital space around titanium is now host to more electrons, which suggests an increase in electron density around the titanium atom 29 . The [TiOH] + structure is a useful reference as each configuration grows in complexity, especially to make inferences about the maximum amount of phosphate ions that can form a bond with a Ti (II) ion and the overall prediction of the binding energy for HAp. This structure demonstrates the lowest binding energy with Ti (II) at .734 atomic units (Table 1).

Titanium (II) and two hydroxide ions: Ti(OH) 2
A neutrally charged compound, Ti(OH) 2 introduces electrons from another oxygen atom to the orbital space  (Table 1 and Figure 5). As mentioned earlier, the bond length between oxygen and titanium in [TiOH] + is approximated at 1.8 Ǻ. The bond length between oxygen and titanium in Ti(OH) 2 is approximated at 1.73 Ǻ, which suggests a stronger bonding strength between titanium and oxygen in the structure. This is additionally supported by a stronger binding energy than [TiOH] + at 1.28 atomic units. The binding energy of the structure nearly doubles and this can primarily be attributed to the strong electronic differences between the oxygen and titanium atoms.

Titanium (II) and the phosphate ion: [TiPO 4 ] -
In this structure, the titanium atom loses a significant amount of charge due to electrons from three oxygen atoms occupying its orbital space. The oxygen atoms lose much of their electrons and the titanium atom results in a charge of 0.124 during the interaction ( Figure 6). The effect that this change in electron density has on the bond lengths between the oxygen atoms, titanium, and phosphorus is such: three of the oxygen atoms that contain most of the negative charge in the phosphate ion are pulled towards the titanium atom and, in consequence, the remaining oxygen atom develops a stronger bond to the phosphorus atom in [TiPO 4 ] 2-. As a reference for later binding energies between phosphate, titanium, and hydroxide, it is important to note that the binding energy for [TiPO4] 2is calculated at approximately 1.88 atomic units.

Titanium (II), hydroxide, and phosphate: [Ti(OH)PO 4 ] 2-
In this configuration, this is the first instance where the titanium atom will have a negative charge (Figure 7). The trend in electron density of titanium is such that titanium consistently gains electron density at a much more dramatic  rate than the phosphorus atom. The phosphorus atom is also consistently gaining electron density from the oxygen atoms, but at much smaller amounts than the titanium atom.

Titanium (II), two hydroxides, and one phosphate: [Ti(OH) 2 PO 4 ] 3-
In the final and most complex structure, the oxygen atoms are contributing even more electron density to the orbital space of the primarily titanium. This results in a centrally negative structure, with all of the positive charge belonging to the hydrogen and phosphorus atoms (Figure 8). Most of the positive charge belongs to phosphorus with a charge of 1.01. The charge of phosphorus in the [TiPO 4 ]was 1.19, which means phosphorus will experience a maximum charge loss of only 0.180. This is considerably smaller than the maximum charge loss of titanium in the [TiOH] + structure, with a maximum charge decrease of approximately 2.00. All of this seems to suggest very favorable trends and interactions with the titanium (II) ion and the component ions of Hydroxyapatite. This particular structure may also demonstrate the maximum amount of phosphate that can form a bond with Ti (II), due to an insufficient positive charge and energy of Ti (II).

Titanium(II) and the hydroxide ion: [TiOH] +
Once again [TiOH] + is a convenient initial structure to act as a reference for the following structures.
[TiOH] + is a linear structure with a bond length of 1.8 Ǻ between the oxygen and the titanium atoms ( Figure 4). The repulsion between the positively charged atoms titanium and hydrogen are likely the cause for this linearity. The calculated frequency interaction between titanium and oxygen is approximately 970 cm -1 , a reasonable vibration 30 .

Titanium (II) and two hydroxide ions: Ti(OH) 2
In the Ti(OH) 2 structure, higher frequencies are calculated between the titanium and oxygen atoms. The geometry of the structure changes to that of being slightly larger than the trigonal planar geometry ( Figure 5). The geometry is slightly larger due to the oxygen-oxygen repulsion of Ti(OH) 2 . Additionally, it is important to recognize that the bond length between the oxygen atoms and the titanium atom has decreased to 1.73 Ǻ when compared to the [TiOH] + structure, whose Titanium-oxygen bond length is 1.80 Ǻ. A larger interaction exists between oxygen and titanium in the Ti(OH) 2 structure, which suggests a favorably higher binding energy to the titanium atom.

Titanium (II) and the phosphate ion: [TiPO 4 ] -
[TiPO 4 ] -( Figure 6) is useful as a reference for the other phosphate configurations in terms of their geometries. Initially, the bond lengths between titanium and oxygen are 1.94 Ǻ per bond. Additionally, the bond angles between titanium and oxygen are 74 degrees per bond. There exists one oxygen atom that is not paired to the titanium atom and instead develops a stronger and shorter bond to the phosphorus atom after reaction ( Figure 6).

Titanium (II), hydroxide, and phosphate: [Ti(OH)PO 4 ] 2-
This structure is a special geometry of interest, as it exhibits a quasi-chair structure (Figure 7). This is a favorable geometry due to a more mechanically reliable structure 31 . When compared to the standard vibrational behavior of phosphate, the increase in frequency between the titanium and oxygen interaction suggests a higher energy in the structure. Additionally, the bond lengths between titanium and oxygen have shortened, also suggesting a more favorable binding energy.

PO 4 ] 3-
Now let us consider the case where another hydroxide reacts with [Ti(OH)PO 4 ] 2-. A tetrahedral geometry is formed and the titanium atom continues to pull oxygen atoms away from the phosphorus atom. However, since the titanium atom also has electrons from two hydroxide ions occupying its orbital space, the bond lengths between the titanium and oxygen atoms are slightly larger than that of the [Ti(OH) PO 4 ] 2structure.

Trends of electron density
Regarding the information from Table 1, the change in charge of Ti(II) appears to occur at a constant rate. As the structure gains more oxygen atoms, the electron density increases much more significantly than that of the phosphorous atom's charge. This seems to suggest that the most influential interaction regarding the binding energy of Ti(II) and HAp is the interaction between titanium and oxygen. As the electron density of titanium increases at a constant rate for each structure, there is a direct relationship with other aspects of the structures such as stronger bonds, higher interaction frequencies, and more favorable geometries.

Geometries and binding energies for each structure.
Additionally, the geometries become more favorable as the complexity of the structures grow. With the [Ti(OH) PO 4 ] 2- (Figure 7) the quasi-chair geometry contains a fourmember ring with low bond angles and high energy. There is an open space on top of the structure near the hydroxide, which would be an area available and energetically favorable for additional interactions. This most likely explains why the [Ti(OH) 2 PO 4 ] 3structure (refer to Figure 8) contains the final hydroxide on the opposite side of the titanium atom, as the titanium atom contains a more positive charge than either of the two oxygen atoms at the top of the structure. This is the closest spot the hydroxide can bond to the lower energy part of [Ti(OH) 2 PO 4 ] 3-(the top space that contains the open oxygen and hydroxide). The [Ti(OH) 2 PO 4 ] 3structure displays the most complex geometry as well as the lowest ground state energy, which according to DFT suggests the highest binding energy compared to the other structures. ] 3-(Ti -larger grey atom, O -red atom, P -orange atom, H -smaller grey atom) (a) longitudinal oscillation centralized between Ti and P (b) transverse oscillation between Ti and Phosphate (c) longitudinal oscillation centralized between Ti and P with strong Hydrogen-Hydrogen repulsion -boat structure, which is the most stable configuration, since Titanium and P interact because of Ti having now a negative charge.

There seems to be a constant relationship between vibrational frequencies, bond lengths, geometries, and electron densities with the overall binding energy for Ti(II) for each compound.
As stated in 5.1, as the electron density of titanium increases with each structure, the bonding strength increases between titanium and its neighbor oxygen atoms. Titanium is predominantly influential with its drastic change in charge compared to the phosphorous atom, which has very little change in charge per structure. As the structures increase in complexity, the titanium atom's charge becomes increasingly more negative, which results in a more centrally negative structure and a lower ground state energy. Since the orbital space of titanium is increasingly occupied with electrons from neighbor oxygen atoms, there is a limited number of phosphate ions that can form a bond with titanium. This suggests that the phosphate-titanium (II) interactions may act as a parameter for this model, which will constrain its scope.

What can the conclusions suggest about the overall strength of HAp and what is the "right" bonding strength for the interaction with bone, so that it doesn't interfere with necessary bodily healing functions.
In biological and medical practices, it is typical to consider if the binding energy of chemicals that attach to compounds and tissues found in the human body will interfere with necessary biological functions of the body 32,33 . Determining an accurate maximization optimization of the binding energy of HAp to Ti(II) will be useful for practicalities such as ensuring that the bonding strength of HAp and Ti(II) doesn't interfere with the natural healing process. However, considering the grafting and healing capabilities that the HAp has with bone, it may be of best interest to maximize the binding energy of HAp with Ti(II) for medical procedures.

Acknowledgements
The work was done with partial support of RFBR grants No. 17-08-01579 and No. 17-08-01312 and of state grant #AAAA-A17-117021310386-3. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562. We thank consultants for their assistance with porting, optimization and visualization, which was made possible through the XSEDE Extended Collaborative Support Service (ECSS) program 34,35 .