Substitution of Strontium, Beryllium, and Magnesium into the Matrix of Hydroxyapatite: A Computer Modeling Study

R. Tony Billingsley

Senior Chemistry Research

May 16, 1997

Dr. Karen Stevens

Abstract

A scaled-down version of hydroxyapatite, the unit cell of bone, was used for geometry optimization and frequency calculations. The calcium atoms that naturally occur in bone were then replaced with strontium, beryllium, and magnesium and the geometry and frequency calculations were then repeated. Errors occurred which prevented the beryllium and magnesium species from being analyzed completely, however, the vibrational modes of the strontium species were compared to those of the calcium species. From this data it was discovered that strontium can weaken the hydroxyapatite matrix through resonance vibrations when it is present in addition to calcium.

Theory

The unit cell of bone, hydroxyapatite, is composed of PO₄, OH, and Ca (figure 1). The hydroxyapatite matrix is composed of repeating hexagonal columnar crystals (figures 2 and 3). These crystals pack together three dimensionally as shown in figures 2 and 3. Because calcium, beryllium, magnesium and strontium are all in group 2A of the periodic table, they share certain physical properties.

The group 2A elements, commonly called the alkaline earth metals, all have ns² valence electron configurations. For this reason, they are all able to be substituted into hydroxyapatite in a fashion that is similar to that of calcium. The group 2A elements, like other groups in the periodic table, have atomic radii that increase with increasing atomic number. Figure 4 shows the pattern of increasing atomic radius for the group 2A elements. Although the radii do increase, it is interesting to note that the difference between the radii of Ca and Sr is the smallest out of all of the elements studied. Table 1 shows the percent increase in the radii of the group 2A elements.

Strontium does not occur as the free element. Strontium is softer than calcium and decomposes water more vigorously. Strontium-90 has a half-life of 28 years. It is a product of nuclear fallout and presents major health problems. Because calcium and strontium share the properties of group 2A elements, strontium can replace calcium in hydroxyapatite, and hence in bone, without much difficulty. This substitution, however simple, should lead to some differences in the properties of the bone. Because strontium is much larger than calcium it might seem logical that its incorporation into the matrix could cause a bulge that would lead to a looser fit of the strontium-containing cell with surrounding hydroxyapatite cells. This could lead to an increased amount of flexibility in that section of the bone as well as changes in the dipole moment and other physical characteristics of hydroxyapatite.

However, since hydroxyapatite is such a large molecule and has a large number of ionic bonds and associations (as compared to covalent bonds), it is difficult to model using Spartanä . For this reason, the theory of the model was extrapolated to a scaled down version that included two PO₄ groups with three cations between them (figure 4). Calculations done with the different metals would be expected to give results that could be used to understand the changes in hydroxyapatite as a whole.

In order to determine whether the addition of strontium or other metals does change the properties of bone, vibrational calculations were carried out. Because the vibrations of a molecule depend on the molecule’s rigidity, molecules that have larger ranges of vibrations can be understood to be more flexible (less rigid). This increased flexibility of the molecules could be considered to correlate to the softness of the molecule as a whole. Thus, vibrational analysis can be used to analyze the properties of bone.

Procedure

The apparatus used for this research project consisted of an Alltek Pentiumä 200 computer running Windows95ä and PC Spartanä , a Macintosh 7200/75 computer running Spartan Plusä , and a Macintosh 4400/200 computer also running Spartan Plusä .

Spartanä offers five levels of calculational accuracy: SYBYL, AM1, PM3, 3-21G ^(*), and 6-31G^*. Each of these types of calculations is good for certain types of molecules. The explanations of each of these types are directly from both the Spartanä manual (2) and from the Spartanä on-line help program. SYBYL is a molecular modeling technique that is very simple, and hence is applicable to a wide variety of compounds. With this method, the total energy of a molecule is described in terms of a sum of contributions arising from distortions from "ideal" bond distances ("stretch contributions"), bond angles ("bend contributions"), and dihedral angles ("torsion contributions"), together with contributions due to "non-bonded" interactions (2). This can be understood by looking at equation 1 below. In this equation, the first three summations are over all "bonds," all "bond angles," and all "dihedral angles," respectively, while the last summation is over all pairs of atoms which are not bonded.

Equation 1 - The Molecular Mechanics Model Framework

Stretch and bend terms are usually given as simple quadratic ("Hook’s Law") forms as can be seen in equations 2 and 3 below. In these equations, r_i^equilibrium and a _i^equilibrium are idealized bond lengths and bond angles respectively. In addition, k_i^stretch and k_i^bend are so-called stretch and bend "force constants," respectively.

Equations 2 and 3 - Stretch and Bend terms in simple quadratic form

Other factors, such as torsion, van der Waals interactions, and electrostatic terms are also taken into consideration, but will not be discussed in this paper. Molecular mechanics models do incorporate a range of factors in their calculations, but they are relatively simplistic compared to other models. Because of this simplicity SYBYL calculations can lead to poor results when high precision is required.

The ab-initio methods offered in Spartanä are 3-21G^(*) and 6-31G*. These models use an approximation of the Schrö dinger equation. The Schrö dinger equation () has been solved for the hydrogen atom, but exact solutions of many-electron atoms, much less whole molecules, are beyond reach. Therefore, to go from the full Schrö dinger equation to a practical approximation, the following approximations are made.

There needs to be a separation of nuclear and electron motions. This is accomplished by implementing the Born Oppenheimer approximation which assumes that the nuclei are stationary from the point of view of the electrons. Using this approximation eliminates the nuclear kinetic energy term in the Hamiltonian and leads to a constant nuclear-nuclear potential energy term. Thus eliminating the mass dependence in the Schrö dinger equation.

The electron motions must also be separated. This is done by implementing the Hartree-Fock approximation. This is used to represent the many-electron wavefunction as a sum of products of one-electron wavefunctions, the spatial parts of which are termed "molecular orbitals."

Finally, the individual molecular orbitals must be represented in terms of liner combinations of atom-centered basis functions, or atomic orbitals. This is termed the LCAO approximation and helps to reduce the problem of finding the best functional form to the much simpler problem of finding the best set of linear coefficients.

3-21G^(*) calculations are generally successful in accounting for equilibrium structures of organic molecules. They are also successful when incorporating second row and heavier main-group elements into molecules. 3-21G^(*) calculations are the simplest method of choice for widespread application while maintaining a high level of accuracy. They are usually superior for performing calculations on all systems except for those requiring a distinction between cis and trans.

AM1 and PM3 are semi-empirical methods of calculation. Semi-empirical models follow from ab-initio models except that they use a single approximation, the NDDO (neglect of diatomic differential overlap) approximation. This completely eliminates overlap of atomic basis functions of different atoms. However, it leads directly to a reduction in computation effort, thus saving computer time. Semi-empirical calculations are generally good for calculations of equilibrium geometry. Despite an increased accuracy over SYBYL, AM1 and PM3 calculations are not as reliable as even the simplest ab-initio 3-21G^(*) calculation. They can have bondlength and bond angle errors that are twice as large as those from 3-21G^(*) calculations (2). Also, they are not as successful for calculations on molecules with second row and heavier main-group elements than they are for calculations with first row atoms. AM1 and PM3 calculations also perform somewhat less well for ions and free radicals than for conventional molecules. This is most likely due to sparse representation of these types of molecules in the sets used for parameterization, rather than any limitations of the model itself. AM1 and PM3 calculations are also considered to be generally unreliable in descriptions of hydrogen-bonded systems.

For each of the molecules studied, a model was built using the Spartan Plus (or PC Spartan) builder program. This molecule was then used for calculations of geometry optimization and frequency. The geometry optimization was always carried out at the SYBYL level. Frequency calculations were conducted at the lowest available level of accuracy. For the (PO₄)₂Be₃ and (PO₄)₂Mg₃ this was the PM3 level, however, for (PO₄)₂Ca₃ and (PO₄)₂Sr₃ the frequency calculations were carried out at the 3-21G* level.

The numerous problems encountered (discussed in detail in the following section) each called for some variation in the procedure. Those variations will not be discussed here because this section is meant to describe only the general procedure used. When a solution was found, it is described in the section on problems.

Problems Experienced

As is often the case with research, a number of problems were experienced during the course of this study. For this project the problems were divided into three main areas: problems with Spartanä , problems with hardware, and problems with scheduling and resource availability. Each of these difficulties played a major role in inhibiting the productivity of this project, and often the presence of one of these problems not only led to, but compounded the effect of other problems. The details of the difficulties experienced, as well as their solutions, are described here.

Spartanä , regardless of whether it was running on a PC or a Macintosh, was the cause of a number of problems throughout the sequence of this project. Obviously this entire research project could not have been done without the program, but its limitations were definitely apparent. The most common error received while performing calculations with Spartanä was "SCF Failed to Converge." The Spartanä manual is not exactly detailed concerning how to solve problems such as this one. The details on convergence say only that you can turn off pseudodiagonalization (in semi-empirical methods) in response to convergence difficulties (2); It does not, however, define pseudodiagonalization. It also says that the SCF density can be set to F by using the "convergence = F" keyword, but again it does not define what this (2). Some more detailed investigative work reveals that increasing the ‘x’ in the "maxcycle = x" keyword can help with convergence difficulties. This was tried with the value of ‘x’ being increased from 200 to 2,000,000 cycles with no success. Finally, the manual indicates that checking the "converge" box in the calculation dialogue can help with convergence difficulties, but may add time to calculations. Using this method, a successful frequency calculation was finally completed for the (PO₄)₂Ca₃ molecule. Although this worked for the (PO₄)₂Ca₃ species, the same problems were encountered when calculating the frequency of (PO₄)₂Sr₃ and (PO₄)₂Mg₃. For (PO₄)₂Sr₃, the problem was finally solved by running the calculation with the "converge" box checked. Although this extended the calculation to 99 hours 46 minutes, it solved the problem. The problem with (PO₄)₂Mg₃ was never solved, and as a result only partial frequency data was obtained. The other result of this was that the vibrations could not be viewed using Spartanä , so a visual comparison of the vibrations of (PO₄)₂Mg₃ with the other molecules studied could not be performed.

This was not, however, the only error received from Spartanä . There were a number of other errors, most of which are insignificant. However, when performing the frequency calculation for (PO₄)₂Be₃, the error "constraint involves stationary atoms" was encountered. Spartanä allows constraints to be placed on specific atoms when the calculation requires a fixed distance between atoms. In this case not only was that not required, but it was not wanted because the spatial arrangement of the atoms in the molecule is crucial to the molecule’s vibrational flexibility. After making sure that no constraints had been placed on any of the atoms the calculation was attempted a second time with the same results. This problem was also never solved, and the Spartanä manual mentions nothing helpful on the subject. The result of this problem is the same as that for (PO₄)₂Mg₃, only partial frequency data was obtained and visual vibrational comparisons could not be performed.

The final area of Spartanä problems was the inability of the Spartanä program to perform calculations on molecules converted from the Spartan Plusä that contained inorganic atoms. Originally a new Pentiumä computer was purchased for this research, however when it was found that calculations involving inorganic atoms had to be done using Spartan Plusä , a Macintosh system had to be used instead. The reason for this is that Spartan Plusä has not yet been written for the PC. This problem with the conversion of molecules between Spartan Plusä and Spartanä actually led to most of the hardware and resource availability problems described in the remainder of this section.

Hardware problems were actually minimal in this study. Because a Macintosh system had to be used, the calculations were initially carried out using the Macintosh 7200/75 computer in the science auditorium. However, as the problems with scheduling (described in the next paragraph) escalated, a Macintosh 4400/200 was used instead. The availability of this machine was quite fortunate because it allowed calculations to be carried out on a computer that was not being used by anyone else at the time. Unfortunately, however, none of the calculations involving heavy atoms (such as strontium) could be carried out on this computer because it had only 16 megabytes of physical RAM. Spartan Plusä kept reporting that the calculation could not proceed because the Spartanä engine had "run out of memory."

Finally, scheduling and resource availability problems comprised a large portion of the difficulties encountered while performing this study. The use of the Macintosh 7200/75 in the science auditorium was greatly appreciated, however that computer is available for use by any class that meets in that room. One class in particular met there every Tuesday and Thursday and used the computer every time the class met. Because of this, approximately twenty to thirty calculations that were in progress were stopped before they had finished. Since most of the calculations done took more than forty-eight hours to complete, most of the calculations had to be done between Thursday afternoon and Tuesday morning.

Data and Calculations

For all of the calculations performed in this study, Spartanä or Spartan Plusä was used. SYBYL was used as the level of accuracy for all geometry optimizations and either the PM3 or 3-21G* methods were used for frequency calculations. The calculations lasted anywhere from ten seconds to nearly 100 hours depending on the molecule, the type of calculation, and the level of accuracy. Tables one and two list the vibrational frequencies obtained for (PO₄)₂Ca₃ and (PO₄)₂Sr₃. Paired with each of the frequencies is a description of the molecular vibrations. The descriptions can be understood in the following manner. Sync/Async describes whether the proceeding action happened synchronously or asynchronously between the top and bottom halves of the molecule. Symm/Asymm describes whether the proceeding action happened at the same time or at opposite times around the species designated (i.e., symm stretch of O-Ca-O means that both Ca-O bonds elongate and withdraw simultaneously). A stretch is the elongation and shrinking of a bond. A wag is the movement of an atom radially back and forth with a fixed bond distance. Finally, in the table describing (PO₄)₂Ca₃, the asymmetric O-Ca-O group (asymm O-Ca-O) refers to the side of the (PO₄)₂Ca₃ molecule that is asymmetric to the rest of the molecule (see figure 6).

Table four shows a list of all of the vibrational frequencies obtained for all four molecules studied. For (PO₄)₂Ca₃ and (PO₄)₂Sr₃ (the two molecules for which the vibrations could be visualized in Spartanä ) the colored boxes indicate vibrational modes that correspond between the two molecules.

Table 4 – Vibrational Frequencies for each of the tested molecules

Table 3 shows that (PO₄)₂Sr₃ has a number of repeated vibrational frequencies whereas Table 2 shows that (PO₄)₂Ca₃ does not. The double frequencies can be attributed to vibrations that differ by the axis of vibration. Generally, one of the vibrations will vibrate along the axis of one of the bonds while the other will vibrate between bond axes. This can be seen in figure 5. The reason that this is not seen in the (PO₄)₂Ca₃ species is because Spartanä calculated it’s geometry as asymmetric (figure 6).

Figure 6 - Side and Top Views of (PO₄)₂Ca₃ - Showing Asymmetry

Figure 7 - Side and Top Views of (PO₄)₂Sr₃ - Showing Symmetry

Figure 8 - Side and Top Views of (PO₄)₂Be₃ - Showing Asymmetry

Figure 9 - Side and Top Views of (PO₄)₂Mg₃ - Showing Asymmetry

Error Analysis

There are two main sources of error that occurred in this study. The first source of error is in the Spartanä calculations of the geometry. As can be seen in figures 6, 7, and 9, the geometry of three of the species studied is asymmetrical. This most likely is the result of the repeated stopping and restarting of the calculations. This error actually led to the second source of error, vibrational comparisons. Although the greatest effort was made to match up the vibrational motifs of (PO₄)₂Ca₃ and (PO₄)₂Sr₃, the asymmetry of the (PO₄)₂Ca₃ species made the task very difficult. If there had been a symmetry to the (PO₄)₂Ca₃ species, it is likely that there would have been more vibrational modes that could have been correlated between the two molecules. The same stopping and repeated restarting of the calculations may also have led to the errors associated with the vibrations of the (PO₄)₂Be₃ and (PO₄)₂Mg₃ species. As mentioned already, the calculations for the frequency of these two species produced errors which prevented the vibrations of these two molecules from being viewed using the Spartanä program. Because of this, there could be no reliable comparisons made to the other two species. Although these errors were significant, the study was completed despite them. With more time and computer availability it is probable that these problems could be overcome through additional studies.

Results

Looking at tables 2-4 it can be seen that the vibrations of (PO₄)₂Ca₃ always occur at a higher vibrational frequency that the corresponding vibrations in (PO₄)₂Sr₃. This can be understood by looking at the reduced mass for both species. The equation for the reduced mass of a bond is shown below in equation 4. To compare the (PO₄)₂Ca₃ and (PO₄)₂Sr₃ species the reduced mass of a Ca-O and Sr-O bond (respectively) was calculated. The results are shown in table 7.

Equation 4 - The reduced mass of a bond

If the force constant k is assumed to be a constant in equation 5, and w is understood to be equal to (û )*(2ð c) (where û is the wavenumber of corresponding vibrations), then it can be seen that the vibrations of (PO₄)₂Sr₃ should occur at a lower wavenumber than those of (PO₄)₂Ca₃. Indeed this is what is seen in table 4.

Equation 5 - Relating frequency to reduced mass and force constant

To determine what the value of k is for both species, equation 5 was solved for k as shown in equation 6.

Equation 6 - solving for the force constant k

Solving for k for both (PO₄)₂Ca₃ and (PO₄)₂Sr₃ yields the values shown in table 7. Although these values are not identical, they are comparable. One factor which could have affected this value is the asymmetry of the (PO₄)₂Ca₃ species. The asymmetry of one of the sides definitely had some effect on the average force constant of the Ca-O bonds, that effect could have been to increase the value of k. It is important is to keep in mind that these values for k are similar, and can therefore help to make sense of the fact that the vibrational frequencies of (PO₄)₂Sr₃ were consistently lower than those of (PO₄)₂Ca₃.

At times, the order of the vibrational modes changes, for instance in (PO₄)₂Ca₃ the asynchronous wag of the P=O groups accompanied with an asynchronous, asymmetric stretch of two O-Ca-O groups occurs at a lower vibrational frequency than the synchronous, symmetric stretch of all of the O-Ca-O groups. In (PO₄)₂Sr₃ however, the order is reversed. This happened with a number of vibrational frequencies (as can be seen in table 4) and could very well be due to the fact that the (PO₄)₂Ca₃ molecule was calculated with an asymmetrical geometry. This asymmetry (shown in figure 6) made it difficult to compare vibrational motifs with (PO₄)₂Sr₃, which was symmetric (figure 7). The asymmetry of the (PO₄)₂Ca₃, (PO₄)₂Be₃, and (PO₄)₂Mg₃ species (figures 6, 8, and 9) is believed to have resulted from the repeated stopping and restarting of the calculations that occurred as a result of the computer being used by a class. The damage to the calculations was much more extensive in the (PO₄)₂Be₃, and (PO₄)₂Mg₃ species because Spartanä was not able even to show their vibrational frequencies. The (PO₄)₂Mg₃ species shows the most asymmetry of any of the molecules studied. Because of these factors the (PO₄)₂Be₃, and (PO₄)₂Mg₃ species will not be considered for the results of this study. It is, however, probable that these two species would follow the same pattern as the (PO₄)₂Ca₃ and (PO₄)₂Sr₃ species, with the vibrational frequency of corresponding vibrational modes increasing above those of (PO₄)₂Ca₃ for the (PO₄)₂Mg₃ species and increasing even more for (PO₄)₂Be₃.

So, how does all of this data fit into the understanding of what atomic substitutions in hydroxyapatite do to the stability of bone? To understand this the vibrational frequencies of PO₄ (table 5) must be compared with the data for (PO₄)₂Ca₃ and (PO₄)₂Sr₃. The stability of the hydroxyapatite matrix, like any matrix, depends partially on the vibrations of each of the involved species. If resonance frequencies exist within the matrix, it is generally less stable. Because this model has been scaled down for this study, the vibrations of (PO₄)₂Ca₃ and (PO₄)₂Sr₃ were compared to those of PO₄ in order to get a general understanding of the stability of each relative to PO₄. Looking at these vibrations, there are no exact matches between multiples of PO₄ vibrations and (PO₄)₂Ca₃ or (PO₄)₂Sr₃ vibrations. So it would appear that the hydroxyapatite matrix would be approximately equally stable regardless of whether (PO₄)₂Ca₃ or (PO₄)₂Sr₃ were involved. However, people who are exposed to radioactive fallout, including fallout involving strontium, never have all of the calcium in their bones replaced with strontium. Therefore it is important to look at the vibrational interference between (PO₄)₂Ca₃ and (PO₄)₂Sr₃. Looking at table 4 it can be seen that there is definite resonance between these two species. The (PO₄)₂Sr₃ frequencies of 241.69 and 262.63 (cm^-1) are almost exactly half of the (PO₄)₂Ca₃ frequencies of 483.93 and 524.16 (cm^-1). Similarly, the (PO₄)₂Ca₃ frequencies of 121.35, 298.60, and 324.16 (cm^-1) are almost exactly half of the (PO₄)₂Sr₃ frequencies of 241.69, 600.35, and 646.12 (cm^-1). A comparison of these vibrations and their descriptions can be found in table 6. Although a few of the resonant vibrations do not involve the same vibrational motifs, two of them do. The synchronous, symmetric stretching of all of the O-Ca-O and the synchronous, asymmetric stretching of the O-Sr-O groups, occurring at 298.60 and 600.35 cm^-1 respectively are in resonance. Also, the synchronous, asymmetric stretching of all O-Ca-O groups (483.93 cm^-1) and the synchronous, symmetric stretching of all O-Sr-O groups (241.69 cm^-1) are also in resonance. In both cases these are major vibrational modes. The stretching of these groups can be seen in a number of other modes as well. Therefore, interference by resonance of this type of vibration could lead to a decreased stability of the hydroxyapatite matrix when both calcium and strontium are present.

Conclusions

Because it is chemically similar to calcium, strontium can be substituted into the matrix of bone relatively easily. In fact, Sr⁸⁹ is often intentionally placed into the bones of elderly people to relieve some of the bone pain associated with cancer (3). Strontium-89 is an injectable compound that the body mistakes for calcium and absorbs into the bone. Doctors believe the radiation decreases the swelling around the bone, caused by the cancer, therefore relieving the pain. At M.D. Anderson Cancer Center in Houston, about 400 patients have been injected with Strontium-89. Seventy-five percent have had partial or complete relief from their pain (4). The results of this study, however, would seem to warrant additional research concerning the long-term effects of combining strontium and calcium in bone. Certainly it would be worth investigating how much damage to the bone matrix this strontium treatment can cause. The results of this study suggest that, although strontium would most likely not cause damage to a fully strontium-saturated hydroxyapatite matrix, it has the potential to cause resonance-related damage to a hydroxyapatite matrix when both calcium and strontium are present. However, most of the radioactive strontium found in man and animals is ingested either by eating contaminated plants or drinking contaminated milk (5). So if a person lives near an area where strontium-90 is in high concentrations it can be difficult to avoid. Currently, there is no known way to remove strontium from bone without virtually demineralizing the entire bone (5). Unfortunately, this study will not help those who already have strontium-90 in their bones, but it does shed some light on one way that the strontium may be causing damage. Hopefully this knowledge will be helpful to those people who are looking for ways to fight strontium poisoning.

Works Cited

1 - Zumdahl, Steven S. Chemistry. D.C. Heath and Company: 1986.

2 - Wavefunction, Inc. MacSpartan Plus. Tutorial and User’s Guide. Wavefunction, Inc.: 1996.

3 – Conquest, Spring 1995 Volume 9 Issue 4. (Internet Source).

4 - Ivanhoe Broadcast News, Inc. Bone Pain Relief. Television News Service/Medical Breakthroughs. Ivanhoe, Broadcast News, In. 1997 (Internet source).

5 - McLean, Franklin. Radiation, Isotopes, and Bone. New York, Academic Press: 1964.

General Reference

Grosch, Daniel E, and Larry E. Hopwood. Biological Effects of Radiations. New York, Academic Press: 1979.