Molecular modeling studies of the binding characteristics of phosphates to sevelamer hydrochloride--assessing a novel technique to reduce phosphates contamination.

Recent advances in the use of polymeric materials for the remediation of phosphate from aqueous systems, including biological environments, have prompted the use of computational techniques to determine which methods are feasible to model such complex systems and to model the mechanism of binding action. In particular, Sevelamer Hydrochloride (Renagel®) is used as our model polymer. A relatively simple system is constructed with a dimer of Sevelamer Hydrochloride and four phosphate ions used for capture. This work reports on molecular dynamics and Monte Carlo simulations used to determine average structure, and points of intermolecular interaction, and calculate changes in volume after the successful capture of phosphate in our model. Our resulting volume changes are of the order of 20–25%, in comparison to experimental swelling measurements on similar systems, which have an average swelling of 50–60% in aqueous media.


Introduction
Environmental concerns over phosphate contamination in aqueous systems have spawned the need for its targeted removal from lakes and ponds. Eutrophication of such water systems is a direct effect of phosphorus and other excessive elemental concentrations [1][2][3][4]. These compounds are found in industrial and domestic wastewater. Elements such as nitrogen, which is known to contribute to eutrophication, have been examined for remediation by use of nitrification processes.
However, phosphorus containing compounds have proven to be much more difficult to handle in a practical fashion [2][3][4]. One method which employs the use of a condensed phosphate template resin has been shown to absorb phosphate. This system was produced via an emulsion polymerization of oleylamine, oleyl alcohol, and divinylbenzene to form a surface imprinted binding site for the capture of phosphates [4]. This current work reports on a more simplified system which uses Sevelamer Hydrochloride to sequester phosphate ions.
The similarities between the functional groups in the structure of Sevelamer Hydrochloride and the oleylamine polymerized system allow for direct comparisons between phosphate binding structure and electrostatic/van der Waals interactions between the captured phosphate and amine groups. Use of Monte Carlo simulations to calculate molecular volumes provide a means for comparison of the in silico results with experimental results. This work reports on the culmination of two computational methods to reproduce the swelling in Renagel ® [5][6][7][8] and determine the binding interactions, and structure involved in phosphate capture. This simple model can then be extrapolated to more sophisticated systems for use in predicting polymer systems to be used in wastewater treatment. Figure 1 shows the schematic of Sevelamer Hydrochloride. The NH 2 groups are protonated by the HCl producing a set of positively charged cages to capture the phosphate anions. In addition, hydroxyl groups can function as supplemental sites for binding through hydrogen bond interactions. Using our simple model of 1.5 repeat units shown in Figure 1, molecular dynamics [9,10], and Monte Carlo methods [9,11], insight into the interactions and bonding structure of this system is defined. In particular, Monte Carlo simulations were used to conduct volume calculations to be compared with experimental swelling measurements by taking the difference between the calculated volumes of our model and the same model with four captured phosphate ions in the system.  Figure 1 was taken from the Renagel® website: www.renagel.com.

Experimental Methods (Swelling Measurements)
Particles for swelling studies are obtained by dissolving a Renagel® tablet (800 mg) in 100 mL of deionized water for 1 hour. The solid is filtered, collected, and dried overnight under pressure. Once dried, particles are spread over double-sided tape, and mounted in a flow cell. Solution flow through the cell is set to 100 mL/hr, which is regulated via a syringe pump. Particle size change is analyzed after flowing solutions of .1 M HCl through the cell for 1 hour. Swelling is monitored through an optical microscope at a 10x magnification, and measurements are taken using Zeiss imaging software. Readings are taken in triplicate, the third set using Image J software, and are found to be within ± 5 % of each other.

Molecular Dynamics
Vacuum System Sybyl 6.9 Software [12] was used to build the inhibitor in vacuum and perform minimization and molecular dynamics simulations on SGI Octane 2.
Energy minimization was performed with Tripos [12] force-field and Gasteiger-Huckel [13] charges with a non-bond cutoff distance of 8 Å using the conjugate gradient method until the root mean square deviation (RMS) < 0.05 kcal/mole/Å (where Å is angstrom). An NVT ensemble [9] molecular dynamic simulation was conducted in vacuum with Gasteiger-Huckel [6] charges assigned at 300 K for 400 ps with a time step of 1 fs. Every 200 fs a snapshot was collected.

Solvated System
The dielectric constant of the system was changed to a constant of 78 debye to mimic a solvated (H 2 O) system. The same non-bonded cutoff distance of 8Å was used in the solvated system, along with the same conjugate gradient method and convergence criteria. Again, Gasteiger-Huckel charges were used on the atoms of the system. Using this prescription the large dielectric constant diminishes the coulombic interactions of the charged atoms. This allows for the study of the van der Waals interactions absent strong coulombic effects. Therefore the effect of van der Waals attraction on phosphate capture could be determined.

Monte Carlo
The procedure used to determine the volume of each molecule was a simple Monte Carlo procedure. In broad terms, the Monte Carlo method was broken into four steps. First, a spherical shell was formed around the molecule.
Secondly, several random values were generated as points in space, bound by the dimensions of the sphere. This allowed all points generated to fall within the sphere; however, they might not have fallen within the molecule. Thirdly, an algorithm was used to assess if the values fell within the molecule. Lastly, the volume of the sphere was multiplied by the ratio of points falling within the molecule to the total number of points. The product was assessed to be the volume of the molecule.

Determining the Shell Volume
The spatial coordinates of the molecular system are exported from the molecular modeling program in the cssr format (Cambridge Structure Search and Retrieval). The cssr format is composed of a numbering column, an atom identification column, spatial coordinates, a connectivity matrix, each atom's potential charge and a few other pieces of chemical information.
The spatial coordinates are given in a Cartesian format. The center of the molecule (x',y',z') was determined by averaging the values in the coordinate space. The molecule was then translated by spatial subtraction of each value by the center. This resulted in the center of the molecule having the coordinate of (0,0,0).
The coordinate space is subsequently transformed into spherical coordinates. The largest value of r is determined. A value of 2 Å is added to the larges value of r. This value is termed r shell . From this, the volume of the shell is determined.

Generating Random Numbers
The random numbers used for this simulation are obtained via the stock random number generator of Matlab ® . For each point three numbers are generated. The first two numbers are bound by (-π, π), and the angles are given in radians. The last number is bound by the radius of the shell.

Determining the Ratio
To increase the efficiency of the assessment, the solution technique is parsed into three steps. The first step evaluates if the generated radius is greater than any radius existing on the molecule. If so, the value is assessed to be outside of the molecule and added only to the total count. If the value is less than all values within the molecule, it is counted in the molecule. For the other values, the proximity to an angle is assessed. If 1 or more points have a radius greater than the said angle, and is within 0.07 radians (~2 degrees) of any of those points having a larger radius, then it was assessed within the molecule. This algorithm was run several times such that the value of ratio as shown in Equation 2 became sufficiently stable over many runs.
The reportable final volume of the molecule is assessed by taking the average of several runs. This is performed by utilizing a nested for-loop. The first forloop is used to determine a ratio. The second for-loop is used to determine the stability of the Monte Carlo method in predicting the volume of the molecule. If the predicted value has a covariance less than 3%, then we conclude that the prediction method is valid. Figure 2 depicts experimental swelling results are shown below. These results showed swelling ranging between 10-200%, with the average being approximately 50-60%. As mentioned above, these results were determined by light microscopy using a light microscope at 1300 x10 magnification using Zeiss Magnification Software. The results presented at this point will compare the experimental swelling measurements to that of the calculated values from the simulations. The molecular dynamic simulations are used to construct a view of the average structure of our model system in the presence of four triphosphate ions. The average structure is then used in Monte Carlo simulation to compute the volume of the model system. The percent change of volume (% swelling) is calculated by comparing the volume of model system absent phosphates to that of the same system with the four triphosphate ions. To begin with we present the results of the molecular dynamic system in vacuum. Figure 3 shows the plot of the total energy of the system vs. the time for the NVT ensemblen [9] simulation in vacuum with a constant dielectric of 1 debye. As depicted in Figure 3, the simulation is stable for 1 ns. A space filled representation of the model system is shown in figure 4. It is clear from the representation that the phosphate ions appear trapped in NH 2 and NH 3 cages in the center of the system (3 out of the 4 phosphates) with the outer phosphate ion seemingly plugged into an amine socket on the surface of the model system. This observation leads us to believe that van der Waals interactions are the major contributors to ion capture in this model system. In order to further support this idea the interatomic distances were plotted between various nitrogen and phosphorus atoms. The following discussion shows data depicting the dynamic motion of between the nitrogen groups and phosphate ions being captured. Figure 5 shows the interatomic distances between a nitrogen atom labeled Nitrogen 126 and phosphorus atom labeled Phosphorous 237.

Results
The atom labels represent the atom type and atom number for the entire system. From the plot it is evident that upon the resolution of the reorganization energy the system stabilizes and maintains a bond distance between 2.7 and 4 Å. The plot also shows an adjustment of the system at around 480 ps to allow for the Nitrogen -Phosphorous separation to decrease below 3.5 Å, and eventually the system adjusts back to a larger average separation around 3.5 Å. This would allow for hydrogen bonding between the phosphate ions and the hydrogen atoms on the amine groups and/or van der Waals interactions ions.

Figure 5:
Interatomic distance between Nitrogen (126) and Phosphorous (237) in the model system Renagel ® and four triphosphate ions. The plot shows an adjustment of the system at about 480 ps to allow for the N-P separation to decrease below 3.5 Å and eventually the system adjusts back to a larger average separation around 3.5 Å. Figure 6 below is a similar type of plot with the interatomic distances of Nitrogen 204 and Phosphorous 237 being depicted. In this plot the Nitrogen group is bonded to the end of the polymer allowing for more freedom of motion. This group moves freely in a wagging motion. Therefore, the capture phosphate is held steady while the amine is allowed to wag in and out of the proximity of the phosphate, resulting in a periodic motion depicted in Figure 6. The periodic interatomic separation between the phosphate ion and the amine nitrogen is due to the wagging motion of the amine on the polymeric tail, and varies between a separation of ~5 Å for 200 ps and then returns to a separation of below 4.3 Å for 400 ps with only a brief change in that approximate distance during that period.    Figure 7 above shows the plot of the interatomic distance between Nitrogen (89) and Phosphorous (237). This figure displays the reorganization and stabilization that followed during the simulation. After the reorganization of the polymer the separation was maintained at about 3.5-4 Å. Nitrogen (89) and Phosphorous (237) show typical behavior of amines nitrogen groups and phosphate ions in the central cages of the polymers.
After measuring the distances between a set of interior (captured) phosphate ions the distances remain constant after the stabilization of the simulation. The plot shown in Figure 8 shows the stable separation between Phosphorous 237 and Phosphorous 216 of about 6-7 Å. In fact this type of stabilized separation is evident between all of the phosphorus atoms found in phosphates captured in the center of the polymer. Finally, we discuss the calculated volumes of the model system with four phosphates and the same system without the phosphates. In comparing the two aforementioned system volumes we are able to estimate the percent swelling. Using Monte Carlo simulations we were able to determine the volumes of these systems. Below we show the results of the Monte Carlo simulations using structures from the molecular dynamic simulations. The results of Monte Carlo simulations used to calculate volume of 733 and 933 cubic Å for the systems containing no phosphate ions and for the system containing four phosphate ions, respectively, are shown below in the plots.  With these calculated volumes shown in Figure 9 and 10, the percent swelling can be determined and compared to experimental data. The percent swelling was calculated to be 25% which is below the average percent swelling of 50-60% found in the experimental results.
Tables 1 and 2 summarize the interatomic distances of several of the Nitrogen-Phosphate pairs. One of the possible explanations for the discrepancy between our results and the experimental averages is the lack of solvent found in the experimental measurements. We used a constant dielectric of 78 debye to mimic a solvated system. We produced similar molecular dynamic simulations for several hundred picoseconds that showed internal and external capture of 4 phosphate ions by our model polymer system. The energy plots for these systems are not shown here, but a similar table (Table 2) to the one shown below is presented below for our solvated simulation. Both tables show the mean standard deviation, high, and low values of various inter-atomic separations are reported from the results of a NVT ensemble simulation in vacuum for the four phosphate system. Captured phosphates in the inner portion of the polymer, exhibit Nitrogen-Phosphorous separations on the order of 3-4Å with very small standard deviations. This suggests little or no movement of the central phosphate (captured). In the solvated system the system goes through considerable reorganization before being stabilized. Due to the large dielectric constant in the simulated solvated system, see the equation below, the coulombic interactions are decreased tremendously.  Table 1 the standard deviation of the Nitrogen -Phosphorous separations that are on the centrally captured phosphate amines is very low where as in Table  2 where the solvation is mimicked by the change in dielectric from 1 to 78 debye the deviations are much higher. This is due to the fact that the capture of phosphate is attributed to van der Waals interactions. In fact, the polymer must undergo structural rearrangement to properly incorporate the phosphates in an amine cage. Table 2 summarizes such results.
The hydrogen bonding found with chloride ions and hydrogen atoms in the NVT ensemble simulation with dielectric of 1 (vacuum) is not present in the solvated system (dielectric 78). This is evident when all of the chloride ion escape the system and must be removed from the system entirely in order to produce a stable simulation. The chloride ions rely on electrostatic interactions due to their smaller size as opposed to the larger phosphate ions which experience significant van der Waals interactions.
The calculated average volume from the molecular dynamic simulation of the solvated four phosphate systems produced a percent swelling of 140% which is good agreement with some of the extreme values of the experimental results.

Discussion
The results presented in this work conclude that a simple model of the Renagel ® compound Sevelamer hydrochloride can reproduce % swelling results. Our results derived from molecular dynamic and Monte Carlo simulations showed an average value of 25% in vacuum and over 140% in H 2 O mimicked (dielectric 78) simulations. These results are relatively consistent with experimental results of between 10-200%, with the average being approximately 50-60%. Our results also bring out the dominant nature of van der Waals bonding in this system by effectively eliminating electrostatic interactions by increasing the dielectric to 78 debye.
In this way our results are able to show the capture of the four phosphates with minimal electrostatic interactions. The internal capture of phosphate is apparently achieved through formation of amine cages around the phosphate in the center of the polymer while the external capture of phosphate was accomplished by interactions about one of the P-O bonds with the amines on the exterior of the polymer.