Molecular Simulations of Recognitive Behavior of Molecularly

Simulation of the imprinted network formation was done using the all-atom kinetic gelation ... For a more comprehensive list of citations to this arti...
1 downloads 0 Views 290KB Size
Download PDF
6084

Ind. Eng. Chem. Res. 2007, 46, 6084-6091

Molecular Simulations of Recognitive Behavior of Molecularly Imprinted Intelligent Polymeric Networks David B. Henthorn† and Nicholas A. Peppas*,‡ Department of Chemical and Biological Engineering, UniVersity of Missouri-Rolla, 143 Schrenk Hall, Rolla, Missouri 65409, and Departments of Chemical and Biomedical Engineering and DiVision of Pharmaceutics, UniVersity of Texas at Austin, 1 UniVersity Station, C0400, Austin, Texas 78712

A method simulating the formation of densely cross-linked polymeric networks was developed that incorporates both intramolecular as well as intermolecular interactions and the subsequent effects they have on the end network structure. The all-atom nature of the model allows for the simulation of network formation in a variety of conditions including differing solvent qualities, presence of inert species, as well as nonlocal effects such as polymerization in the presence of a template molecule. We employed an all-atom kinetic gelation technique that utilized an off-lattice approach that tracked the position and interaction of all atoms throughout the simulation. This model was then used to study the formation of polymeric networks capable of recognizing and binding a specific molecule out of a host of competing species. Simulation of the imprinted network formation was done using the all-atom kinetic gelation framework, which helped identify the interactions central to recognition. These results were verified by comparison with previous experimental results. 1. Introduction One of the most regarded characteristics of biological molecules is the ability to recognize and bind specific molecules out of a host of competing species. This ability is key to the body’s signaling, regulation, catalysis, foreign body defense, and transport. As knowledge is furthered regarding how molecules such as enzymes and antibodies are able to recognize molecules, there is growing interest in creating synthetic materials able to mimic their functions. These biomimetic materials would see use in a host of applications such as separation media, sensors, regulators, and even potentially synthetic enzymes. In our group, hydrophilic polymeric networks have been developed1-5 that are capable of binding small molecules, much as natural proteins and enzymes do.6 In this work, we seek to investigate a modeling technique which will allow researchers to study or even design synthetic recognitive materials. Several synthesis strategies have emerged for the fabrication of synthetic materials with recognitive abilities, including designed peptide sequences,7-9 randomly folded heteropolymers,10-14 and molecular imprinting.15-19 The production of designed peptide sequences is hindered by two immense challenges: evaluation of the structures (an inverse protein folding problem) will require exorbitant amounts of computational resources, and synthesis techniques currently limit the number of residues allowed in a peptidic chain.10 As a result, strategies that rely on the production of random heteropolymers have emerged due to their relative ease of synthesis. Binding sites are built into these materials during formation through the use of a template molecule. Pande and co-workers11-14 modeled systems where random monomers were allowed to interact with a template molecule. After a relaxation period, the lowest energy conformation was frozen by “polymerizing” the monomers. Their work showed the imprinting success to be highly de* To whom correspondence should be addressed. Tel.: 512-4716644. Fax: 512-471-8227. E-mail: [email protected]. † University of Missouri-Rolla. ‡ University of Texas at Austin.

pendent on the reaction temperature and strength of interaction. At temperatures where the average thermal energy was greater than the strength of interaction, greatly reduced recognition was observed.13 Experimentally, the idea of molecular imprinting has been used to produce molecularly imprinted polymers (MIP) that exhibit molecular recognition from amino acids,20,21 glucose,1,2,22 and even protein epitopes.23-25 In molecular imprinting, a template molecule is placed in a prepolymerization mixture and allowed to interact. Polymerization of this mixture in the presence of a cross-linking agent results in production of a crosslinked three-dimensional structure with stabilized binding sites. The resulting MIP is then washed before use to remove the template molecule in order to open the binding sites. Cross-linking density is a crucial parameter in the MIP preparation. Yu and Mosbach21 studied its effect on the recognition and separation of with the use of an enantiomeric pair. As cross-linking density decreased, the ability for the material to discern between the two molecules decreased. This was in line with earlier work done by Wulff et al.26 using covalently imprinted MIPs. This need for densely cross-linked materials and the importance of intermolecular interactions makes modeling of MIP formation difficult using traditional methods. The various expressions for rates of initiation, polymerization, and termination may be integrated simultaneously using various numerical techniques. Nonidealities that occur in the formation of densely cross-linked materials may be added to these models, such as varying rate constants and the effects of autoacceleration.27,28 However, information on the network structure itself is still not available. In an effort to develop a model that would accurately model densely cross-linked materials, a random walk computer simulation was first proposed by Manneville and de Seze29 in 1981 and later demonstrated by Boots et al.30 in 1985 and Simon et al.31 in 1989. Initial models were limited in that molecules were fixed on a lattice and were not allowed to diffuse. In addition, all initiator molecules were activated simultaneously at the beginning of the reaction. Boots and Dotson32 later added

10.1021/ie061369l CCC: $37.00 © 2007 American Chemical Society Published on Web 05/08/2007

Ind. Eng. Chem. Res., Vol. 46, No. 19, 2007 6085 Table 1. Terms of the Intra- and Intermolecular Potential Energy Function Used in the All-Atom Kinetic Gelation Simulation interaction type

mathematical description

Lennard-Jones

VLJ )

coulombic

Vc ) qiqj/(4π0rij)

bond distance bond angle dihedral angle improper angle

Vbond ) Kr(r - req)2 Vθ ) Kθ(θ - θeq)2 Vdihedral ) Kχ(1 + cos(nχ - δ)) Vimproper ) Kimp(φ - φeq)2

4[(σ/r)12

-

(σ/r)6]

description used to model van der Waals attractions between atoms with strength of interaction, , and atomic radius, σ models ionic interactions between ions i and j with respective charges qi and qj interacting over a distance of rij through a medium of permittivity 0 harmonic potential for fluctuations about the equilibrium bond length, req, with force constant Kr deviations from the equilibrium bond angle, θeq, with rotational force constant Kθ calculation of potential of rotation around angle χ, with phase δ, and torsional force constant Kχ harmonic potential for deviations of the improper angle φ perturbed from the equilibrium angle φeq with force constant Kimp

polymer diffusion to the model. Bowman and Peppas33,34 substantially furthered the kinetic gelation model when they included the ability for monomers to occupy more than one lattice site, allowed monomers to be of different molecular types, and allowed for more realistic initiation involving exponential decay of the initiator molecules. Ward and Peppas35 further improved the kinetic gelation model to include the following: simulation of monomers of varying sizes including the simulation of macromonomers, reaction in the presence of inert species such as solvent molecules or dissolved drug, and simulation of living free radical polymerizations initiated with a class of molecules known as iniferters36 that are able to form reversible chain transfer bonds. Also, the cubic lattice was extended to a face-centered lattice that incorporated 12 nearest neighbors. The bulk of all kinetic gelation simulations are done with noninteracting spheres representing each monomer segment. With the increase in computational power, researchers have moved kinetic gelation simulations off-lattice and incorporated intermolecular interactions (Table 1).37,38 Addition of the Lennard-Jones potential allows for simulation of dispersion forces between various reacting and inert species. In developing a model for the simulation of recognitive, cross-linked polymeric networks, one essential feature will be the use of atomic-scale simulations. Since molecular-level phenomena govern the final structure and recognitive ability of the network, atomic-level interactions are added to the kinetic gelation model to handle both inter- and intramolecular interactions. 2. Methods 2.1. Force Field Development for Methacrylate and Acrylate Monomers. In order to simulate the polymerization of the various monomers at an atomic level, it was first necessary to develop a force field suitable for the methacrylate and acrylate class of monomers. Other force fields, such as CHARMM39-41 or AMBER,42 have been developed primarily to handle proteins and DNA. Few other classes of molecules have conjugated functional groups similar to those seen in the methacrylates and acrylates; specifically, the close proximity of the carbonyl carbon to the double bond. Figure 1 shows the representative methacrylate molecule, 2-hydroxyethyl methacrylate (HEMA), and the conjugation of the carbonyl carbon with the double bond. Because of the unique chemical environment present in these molecules, it was necessary to determine the force field parameters. Derivation of the force field parameters followed a procedure similar to that prescribed by McKerrell and co-workers41,43,44 in their development of the CHARMM force field. To do this, quantum mechanical simulations were performed to quantify the parameters from first principles. First, the molecular geometry of the monomer was optimized using the semiempirical AM1 method to find a lowest energy conformation. This conformation was then used as the starting point for more

accurate calculations, namely geometry optimization followed by frequency analysis (to ensure a true minimum) with the B3LYP DFT method and the 6-311G+(d,p) basis set. Partial charges were calculated with the natural bond order analysis method (NBO)45 because of its good performance when dealing with diffuse and polarized basis sets. The harmonic terms (deviations from equilibrium bond lengths, bond angles, and improper angles) were parametrized by perturbing the value from the equilibrium one and fitting via a least-squares procedure. Dihedral values were determined by rotating about the central atoms and fitting the corresponding energies to a five term cosine series. Lennard-Jones radii and interaction energies, the most difficult terms to accurately quantify, were adopted from similar atomic species parametrized in the most recent version of CHARMM.41 All electronic structure calculations were carried out with the Gaussian 9846 suite of programs. Cosine series fitting was done in Matlab (Version 6.0). 2.2. Kinetic Gelation Simulation Development. A computational approach has been utilized to understand the interaction between template and monomer by Piletsky et al.47 In the work, the researchers used computational chemistry tools to study interaction of the desired template with a library of potential monomers. The binding energy was then determined and the best monomer chosen for experimental validation. In this work, we aim to study the interaction of template with monomer, but we also include the polymerization process to ultimately study the final polymer network. In order to simulate network formation, the kinetic gelation method, a model based on percolation theory, was used. This model, originally developed by Manneville and de Seze29 and further formalized for simulation of multifunctional monomers by Bowman and Peppas,34 was extended in several ways. First, previous simulations grouped the entire repeating unit of a polymer chain together as a single bead, which then was allowed to take on only discrete positions in space by conducting the simulation on a lattice. While recent simulations35 have extended the lattice from a simple cubic one to a more flexible hexagonal one, this model still restricts the polymer repeating unit to possibly aphysical conformations. Our simulation was therefore moved to an off-lattice type, in which any possible conformation was possible and the polymer chains were therefore allowed to sample the full simulation volume. In addition, the repeating unit based “bead” approach was replaced by the full atomic structure of the repeating unit. This allows for interactions between atoms to occur and for the flexibility of less-energetic conformations to be screened. Solvent molecules were also handled explicitly. In the case of water, Jorgensen’s TIP3P48 water model was used while other models for dimethylsulfoxide49,50 were also utilized. In the previous collection of kinetic gelation simulations, interactions between polymer chains, monomers, solvent molecules, and other inert species were neglected. However, in order to capture the effects that these interactions have on the network

6086

Ind. Eng. Chem. Res., Vol. 46, No. 19, 2007

Figure 1. Structure of the monomer 2-hydroxyethylmethacrylate (HEMA).

Figure 3. Thermodynamic cycle used in the free-energy perturbation calculations. The difference in the free energy associated with binding species S2 instead of species S1 is given by ∆∆G1f2 ) ∆G2 - ∆G1. However, calculation of ∆G1 and ∆G2 is difficult computationally. Since G is a state function, an alternative thermodynamic pathway may be chosen. The work required to change species S1 into S2 is calculated, which while not physically realizable can be calculated by simulation, and ∆∆G1f2 ) ∆G4 - ∆G3.

Figure 4. Atomic definitions for parametrized HEMA molecule.

Figure 2. Kinetic gelation reaction scheme. After an initial equilibration, the reaction is started. At each reaction step, radicals are created, allowed to propagate, or terminated. This is followed by a diffusion and equilibration period. Reaction continues until a final conversion is achieved at which time the final network structure is given.

structure over the course of the reaction, a molecular mechanics style force field was included in the model. The force field developed for our monomers was inserted into the kinetic gelation model. The outline of the kinetic gelation simulation methodology is shown in Figure 2. Initially, a randomly created simulation volume containing monomers, solvent, initiator, and any nonreacting species was created. This mixture was then allowed to equilibrate prior to any reaction. In this time, the molecules diffuse and interact at the given reaction conditions. After this initial equilibration, the reaction steps begin. Reaction steps include initiation, propagation, and termination as shown previously for free-radical polymerizations in Figure 1. After each reaction step, the molecules are then again allowed to diffuse and relax. The reaction-equilibration cycle continues until either a set conversion is met or all propagating radicals are quenched or no longer active. The final network structure is then saved for postprocessing, which involves calculation of network structure quantities, and, in the case of recognitive networks, a free energy analysis of binding. While previous work has focused on the use of Monte Carlo methods to allow for the pseudodiffusion of all the species, it became necessary to switch to a molecular dynamics approach due to the all-atom nature of this approach. In the Monte Carlo methodology, molecules are moved at random, and an acceptance value is calculated according to the energetic penalty associated with the move. With systems such as extended, crosslinked polymer networks, the acceptance rate of these moves becomes exceedingly small due to the large energetic penalty

incurred. While work has been done by various researchers to simulate cross-linked polymer networks using Monte Carlo methods,51-53 acceptance rates are still low, especially for an all-atom simulation. Therefore, our work has focused on using molecular dynamics simulations to handle relaxation and diffusion as it alleviates the need for complex algorithms to improve Monte Carlo trial move acceptance. Before any equilibration or diffusion step, a conjugate gradient minimization was first conducted with 5000 steps to remove any atomic overlap that may have resulted during the construction of the simulation volume. This step is important since any overlap has the potential to cause an instability in the molecular dynamics steps that follow. After the minimization, a standard molecular dynamics simulation was used to allow the molecules to diffuse and interact. Specifically, either constant NVT or constant NPT ensembles were utilized, and a cubic simulation cell with periodic boundary conditions was employed. The equations of motions were integrated with a velocity Verlet method utilizing a time step of 1 femtosecond (fs). The initial simulation employed 40000 time steps for a total simulation time of 40 picoseconds (ps). Subsequent equilibrations between reaction events used 5000 time steps. Full electrostatic interactions were calculated via a particle mesh Ewald algorithm that allows for a reduction of the computational time from order N2 to order N log N.54 Temperature control was handled with velocity rescaling every 100 fs to the target temperature. For simulation in the NPT ensemble, pressure was controlled via coupling with a Berendsen pressure bath. All MD simulations were done using the NAMD molecular dynamics program55,56 because of its excellent scalability on parallel computers.57 In addition, NAMD was linked with the FFTW58 library to provide fast routines for the particle mesh Ewald algorithm. Kinetic gelation was handled with the developed program, “reaction”, written in C++. MD simulations were started from within reaction and their corresponding results read into memory after the completion of the simulation run. New MD runs were formulated after all reaction events that occurred during a given reaction time step were completed. For calculation of initiation

Ind. Eng. Chem. Res., Vol. 46, No. 19, 2007 6087 Table 2. Atomic Charges and CHARMM Definitions for the HEMA Molecule atom name

CHARMM type

atomic charge (units of e)

C1 H11 H12 H13 C2 C3 H31 H32 C4 O1 O2 C5 H51 H52 C6 H61 H62 O3 HO1

CT3 HA HA HA CE0 CE2 HE2 HE2 C O OS CT2 HB HB CT2 HB HB OH1 H

-0.576 0.154 0.154 0.154 0.478 -0.526 0.144 0.144 -0.016 -0.284 -0.073 -0.173 0.174 0.174 -0.146 0.149 0.149 -0.330 0.250

Table 3. Equilibrium Bond Lengths for the HEMA Monomer, as Determined by the DFT Method bond pair

equilibrium bond length (Å)

C1-C2 C2-C3 C2-C4 C4-O1 C4-O2 O2-C5 C5-C6 C6-O3 O3-HO1 (H11,H12,H13)-C1 (H31,H32)-C3 (H51,H52)-C5 (H61,H62)-C6

1.5208 1.3534 1.5062 1.2351 1.3752 1.4683 1.5399 1.4558 0.9792 1.1032 1.0952 1.1012 1.1063 Figure 5. Prereaction interactions seen between HEMA molecules and a single glucose molecule. Results from MD simulation (top) suggest a competitive interaction between water molecules (red and white triangular molecules) and HEMA molecules for the hydroxyl groups of the glucose molecule (center). All molecules within a 5 Å sphere of the glucose molecule are shown here. Three distinct interaction sites (red lines) are seen (bottom) between the monomer molecules and the template. All are of a hydrogenbonding type.

Table 4. Equilibrium Bond Angles for the HEMA Monomer, as Determined through DFT Method bond angle triplet

equilibrium bond angle (deg)

C1-C2-C3 H1*-C1-H1* H1*-C1-C2 C1-C2-C4 C2-C3-H3* C2-C4-O1 O1-C4-O2 C4-O2-C5 O2-C5-H5* C5-C6-O3 H5*-C5-C6 C5-C6-H6* H6*-C6-O3 C6-O3-HO1 H3*-C3-H3* H3*-C3-H3*

123.90 108.16 110.75 114.70 121.25 123.26 123.32 116.41 106.96 105.41 109.84 110.05 111.16 107.87 117.46 121.1

rates, initiator molecules were broken into radicals corresponding to eq 1

( τt )

I ) I0 exp -

(1)

where I is the initiator concentration, I0 is that of the initial conditions, t is the time step of the reaction simulation, and τ is the initiator decay constant. During the propagation steps, all possible reactive molecules within a reaction radius of 7 Å were considered for the reaction. In the absence of any reactivity data, a monomer was selected at random from those within the reaction radius. However, if specific monomer reactivity ratios,

defined as the probability of one monomer being added over another, were known then the reaction was biased accordingly. Termination was done solely by combination; two propagating radicals were required to diffuse within reaction distance and react to annihilate their respective radicals. Radicals were marked as trapped if no reaction had occurred at that center within 20 reaction time steps. 2.3. Free Energy Analysis of Binding for Recognitive Networks. Ultimately, understanding the binding event by a recognitive network involves quantifying the free energy change that is associated with a molecule displacing solvent molecules from the binding area and interacting with the functional groups present on the binding site. Two types of simulation techniques exist that permit the calculation of this free energy change: Gibbs ensemble Monte Carlo (GEMC) and free energy perturbation techniques. In the GEMC simulation, two simulation volumes are used to represent the two distinct phases.59 Temperature, pressure, and the total number of particles are held constant during the simulation, but particles may be transferred from one simulation box to the other, and by doing so, the chemical potential between the simulation volumes equilibrates. In order to simulate the binding of an analyte by a recognitive network, a GEMC simulation may be conducted with one

6088

Ind. Eng. Chem. Res., Vol. 46, No. 19, 2007

Figure 6. Cross-linked recognitive network structure with bound glucose (blue molecules).

simulation volume containing the network while the other is used as a solvent-analyte bath that imposes a chemical potential on the gel. However, like most Monte Carlo simulations of large or complex molecules, the particle transfer acceptance rate and, therefore, the efficiency of the simulation are expected to be low. Free energy perturbation (FEP)60,61 offers an alternative to GEMC and may be conducted within a molecular dynamics simulation. Absolute values of free energy changes are difficult to calculate; however, with this technique differences in ∆G between two states (∆∆G) are readily calculated (Figure 3). FEP is used frequently to study binding of small molecules by proteins with good results. In particular, rational drug design has seen success with FEP.62,63 First, a binding site is identified in a protein, and a series of possible ligands are created. Next, the possible ligands are inserted into the binding site, and ∆∆G is calculated for the series. When rank-ordered, the species match well with experimental trends.64 3. Results and Discussion Our work has focused on the development of the off-lattice kinetic gelation model that incorporates an all-atom potential to model realistic interactions between molecules. In addition to this, because of the large increase in computational resources needed by these techniques over past kinetic gelation simulations, emphasis was also placed on developing computational methods that scale well with size or can be run on parallel computers when possible. 3.1. Methacrylate/Acrylate Force Field Development. The representative molecule, 2-hydroxyethyl methacrylate (HEMA), was fully parametrized at the B3LYP/6-311G+(d,p)//B3LYP/ 6-311G+(d,p) level. The equilibrium geometry was then recorded and used to create a molecular topology file for later

use in the kinetic gelation simulation. All structural properties (bond and bond angle force constants, dihedral rotation penalties, improper energy constants, and atomic partial charges) were placed in a CHARMM compatible parameter file to be read in by the MD simulation engine. Parameter values are given in Tables 2-4 and reference the atomic configuration shown in Figure 4. 3.2. Kinetic Gelation Simulation and Correlations with Experiments. The all-atom kinetic gelation was used to simulate the polymerization reaction. Simulations were conducted on a parallel Linux cluster setup utilizing eight AMD Athlon MP 2000+ processors interconnected with 100BaseT Ethernet. In order to test the applicability of the kinetic gelation model to predict molecular arrangement and ordering, the first run of the kinetic gelation program, reaction, was done on a recently formulated and studied recognitive polymer network of Oral and Peppas65,66 designed to bind glucose. These materials were designed such that at equilibrium there was a 54% increase in glucose uptake when compared to control (nonimprinted) materials. In addition, the glucose-imprinted material showed excellent specificity: galactose, a structurally similar sugar, showed a 34% decrease in uptake with respect to the control material. According to her optimization, Oral found an optimal formulation to be a 1:1 ratio of HEMA molecules to glucose molecules in the reaction mixture. These materials are also densely cross-linked at 80 mol %. For the simulation of a glucose recognitive polymer based on the hydrophilic HEMA monomer, a reaction volume was created with the following: 300 difunctional ethylene glycol dimethacrylate (EGDMA) molecules, 800 water molecules, 160 D-glucopyranose molecules, 160 HEMA molecules, and 20 initiator molecules for a total of 35820 atoms. At an equilibration temperature of 300 K, the prereaction mixture was observed,

Ind. Eng. Chem. Res., Vol. 46, No. 19, 2007 6089

ability to bind specific molecules while rejecting others. In the case of the glucose-binding network, Oral showed these gels exhibit a large preference for glucose over its carbon-4 epimer, galactose. Between these two molecules, only the hydroxyl at the carbon-4 position has been rotated, as is shown in Figure 8. Results from our kinetic gelation model predict that three interaction sites are important to the binding of glucose by the HEMA moiety, including the hydroxyl group of carbon-4 in the down position. Therefore, of interest was the calculation of the change in free energy of binding between glucose and galactose, ∆∆Gglucosefgalactose. Using the free energy perturbation method, ∆∆Gglucosefgalactose was calculated to be +1.8 kcal/mol, which overpredicts galactose binding by 16% when compared to the values of Oral and Peppas. 4. Conclusions Figure 7. ATR-FTIR spectrum of glucose imprinted recognitive network from the Ph.D. thesis of Oral4 showing the shift of the carbonyl peak caused when binding occurs.

Figure 8. Structure of R-D-glucose (left) and R-D-galactose (right) showing epimerization of the carbon-4 hydroxyl group. Both are shown in the pyranose ring form, which is dominant in solution over the open form.

and interactions between the HEMA molecules and the glucose molecules were noted. The glucose molecules were predominantly solvated by a halo of water molecules; however, as shown in Figure 5, three specific interaction sites between the HEMA monomers and the glucose template molecule were seen. Specifically, the tail hydroxyl functional group of the HEMA repeating unit is seen hydrogen bonding with the hydroxyl groups based off of carbons 2 and 4 of D-glucopyranose. In addition, the other specific interaction seen is the hydrogen bonding between the oxygen of the HEMA carbonyl group and the hydroxyl group set on carbon 6 of the glucose molecule. The reaction of this mixture was then simulated in the kinetic gelation framework. Total run time was 5 days on eight processors of the Linux cluster. In addition, a control simulation was also run whereby the glucose template molecules were absent. Figure 6 shows the network structure formed for the recognitive network with bound glucose molecules. Analysis of the network-analyte system shows that both specific and nonspecific binding takes place. In the case of the specific binding, the same three interactions seen in the prereaction system were seen after the full network was formed. In order to validate this data, spectroscopic data from Oral and Peppas65 were used. She used attenuated total reflectanceFTIR (ATR-FTIR) to study any possible sites of interactions between the HEMA based network and the bound glucose molecules. In their studies, they showed that a significant shift of the carbonyl peak occurred during binding (Figure 7). The carbonyl group absorbance of an imprinted gel shifted from 1715 to 1725 cm-1 upon introduction and binding of glucose. In addition, a shoulder at 1715 cm-1 is present in the bound network, demonstrating that both interacting and noninteracting carbonyl groups are present. 3.3. Free Energy Analysis of Glucose Binding. Of great interest in the study of recognitive polymeric networks is their

An all-atom kinetic gelation simulation technique was developed to study free-radical polymerizations that are not easily described by conventional techniques. The method was successfully applied in the study of recognitive polymeric network formation in which interactions between the reacting monomers and small molecules acting as templates dictate the final network structure. With this technique the specific atomic interactions that cause the final network conformation were first seen and later validated with spectroscopic methods. It was also possible to study the specificity of binding using the free energy perturbation technique to quantify binding free energy differences between conformers. As computational resources continue to improve, all-atom models such as these could prove useful in elucidating the structure, function, and response of imprinted polymeric networks. In addition, knowledge of templatemonomer and template-polymer interactions may aid in the creation of in silico methods of network design and optimization. Further studies may include varying reactivity parameters for trapped or pendant groups, and the inclusion of “living” free polymerization techniques as illustrated in the work of Vaughan et al.67 Acknowledgment This work was supported in part by a grant from the National Science Foundation. This work is dedicated to Professor Kenneth A. Smith of MIT on the occasion of his 70th birthday. Ken Smith’s visionary approach to the mathematical analysis of chemical engineering problems had a major impact in the education of the senior author (NAP), who served as a postdoctoral fellow with him and Clark K. Colton at MIT from 1974 to 1976. Literature Cited (1) Byrne, M. E.; Park, K.; Peppas, N. A. Molecular Imprinting within Hydrogels. AdV. Drug DeliVery ReV. 2002, 54, 149-155. (2) Byrne, M. E.; Oral, E.; Hilt, J. Z.; Peppas, N. A. Networks for Recognition of Biomolecules: Molecular Imprinting and Micropatterning Poly(Ethylene Glycol)-Containing Films. Polym. AdV. Technol. 2002, 13, 798-810. (3) Byrne, M. E.; Park, K.; Peppas, N. A. Biomimetic Materials for Selective Recognition of Biomolecules. In Biological and Biomimetic MaterialssProperties to Function; McKittrick, J., Aizenberg, J., Kittrick, J. M. M., Orme, C. A., Vekilov, P., Eds.; Materials Research Society (MRS): Pittsburgh, PA, 2002; Vol. 724, pp 193-199. (4) Oral, E. Recognitive Polymer Networks for Therapeutic and Diagnostic Devices. In Ph.D. Thesis, School of Chemical Engineering, Purdue University, West Lafayette, IN, 2002. (5) Oral, E.; Peppas, N. A. Responsive and Recognitive Hydrogels Using Star Polymers. J. Biomed. Mater. Res., Part A 2004, 68, 439-508.

6090

Ind. Eng. Chem. Res., Vol. 46, No. 19, 2007

(6) Fersht, A. Structure and Mechanism in Protein Science: A Guide to Enzyme Catalysis and Protein Folding; W.H. Freeman and Co.: New York, 1999. (7) Kirshenbaum, K.; Zuckermann, R. N.; Dill, K. A. Designing Polymers That Mimic Biomolecules. Curr. Opin. Struct. Biol. 1999, 9, 530560. (8) Koehl, P.; Levitt, M. De Novo Protein Design. I. In Search of Stability and Specificity. J. Mol. Biol. 1999, 293, 1161-1181. (9) Koehl, P.; Levitt, M. De Novo Protein Design. II. Plasticity in Sequence Space. J. Mol. Biol. 1999, 293, 1183-1193. (10) Shakhnovich, E. I.; Gutin, A. M. Engineering of Stable and FastFolding Sequences of Model Proteins. Proc. Natl. Acad. Sci. U.S.A. 1993, 90, 7195-7199. (11) Pande, V. S.; Grosberg, A. Y.; Tanaka, T. Folding Thermodynamics and Kinetics of Imprinted Renaturable Heteropolymers. J. Chem. Phys. 1994, 101, 8246-8257. (12) Pande, V. S.; Grosberg, A. Y.; Tanaka, T. Thermodynamic Procedure to Synthesize Heteropolymers That Can Renature to Recognize a Given Target Molecule. Proc. Natl. Acad. Sci. U.S.A. 1994, 91, 1297612979. (13) Pande, V. S.; Grosberg, A. Y.; Tanaka, T. Phase Diagram of Heteropolymers with an Imprinted Conformation. Macromolecules 1995, 28, 2218-2227. (14) Pande, V. S.; Grosberg, A. Y.; Tanaka, T. How to Create Polymers with Protein-Like Capabilities: A Theoretical Suggestion. Physica D 1997, 107, 316-321. (15) Mosbach, K.; Ramstrom, O. The Emerging Technique of Molecular Imprinting and Its Future Impact on Biotechnology. Biotechnology 1996, 14, 163-170. (16) Mosbach, K.; Haupt, K.; Liu, X.-C.; Cormack, P. A. G.; Ramstrom, O. Molecular Imprinting: Status Artis et Quo Vadere. In Molecular and Ionic Recogniton with Imprinted Polymers; Bartsch, R. A., Maeda, M., Eds.; American Chemical Society: Washington, DC, 1998; pp 345-356. (17) Mosbach, K. Toward the Next Generation of Molecular Imprinting with Emphasis on the Formation, by Direct Molding, of Compounds with Biological Activity (Biomimetics). Anal. Chim. Acta 2001, 435, 3-8. (18) Whitcombe, M. J.; Vulfson, E. N. Imprinted Polymers. AdV. Mater. 2001, 13, 467-480. (19) Komiyama, M.; Takeuchi, T.; Mukawa, T.; Asanuma, H. Molecular Imprinting from Fundamentals to Applications; Wiley-VCH: Weinheim, Germany, 2003. (20) Andersson, L.; Sellergren, B.; Mosbach, K. Imprinting of Amino Acid Derivatives in Macroporous Polymers. Tetrahedron Lett. 1984, 25, 5211-5214. (21) Yu, C.; Mosbach, K. Influence of Mobile Phase Composition and Cross-Linking Density on the Enantiomeric Recognition Properties of Molecularly Imprinted Polymers. J. Chromatogr., A 2000, 888, 63-72. (22) Oral, E.; Peppas, N. A. Dynamic Studies of Molecular Imprinting Polymerizations. Polymer 2004, 45, 6163-6173. (23) Kempe, M.; Mosbach, K. Separation of Amino Acids, Peptides and Proteins on Molecularly Imprinted Stationary Phases. J. Chromatogr., A 1995, 691, 317-323. (24) Cederfur, J.; Pei, Y.; Zihui, M.; Kempe, M. Synthesis and Screening of a Molecularly Imprinted Polymer Library Targeted for Penicillin G. J. Comb. Chem. 2003, 5, 67-72. (25) Henthorn, D. B.; Zheng, Y.; Peppas, N. A. Synthetic LigandReceptor Interactions in Delivery Systems. In Nanotechnology in Therapeutics: Current Technology and Applications; Peppas, N. A., Thomas, J. B., Hilt, Z., Eds.; Horizon Press: Hetherset, U.K., 2006. (26) Wulff, G.; Vietmeier, J.; Poll, H.-G. Enzyme-Analogue Built Polymers, 22: Influence of the Nature of the Crosslinking Agent on the Performance of Imprinted Polymers in Racemic Resolution. Makromol. Chem. 1987, 188, 731-740. (27) Mikos, A. G.; Takoudis, C. G.; Peppas, N. A. Kinetic Modeling of Copolymerization/Cross-Linking Reactions. Macromolecules 1986, 19, 2174-2182. (28) Godner, M.; Bowman, C. N. Modeling Primary Radical Termination and Its Effects on Autoacceleration in Photopolymerization Kinetics. Macromolecules 1999, 32, 6552-6559. (29) Manneville, P.; Seze, L. d. Percolation and Gelation by Additive Polymerization. In Numerical Methods in the Study of Critical Phenomena; Dora, J. D., Demongeot, J., Lacolle, B., Eds.; Springer: Berlin, 1981; pp 116-124. (30) Boots, H. M. J.; Kloosterboer, J. G.; Hei, G. M. M. v. d.; Pandey, R. B. Inhomogeneity During the Bulk Polymerization of Divinyl Compounds: Differential Scanning Calorimetry Experiments and Percolation Theory. Br. Polym. J. 1985, 17, 219-223. (31) Simon, G. P.; Allen, P. E. M.; Bennett, D. J.; Williams, D. R. G.; Williams, E. H. Nature of Residual Unsaturation During Cure of Dimethacry-

lates Examined by CPPEMAS 13C NMR and Simulation Using a Kinetic Gelation Model. Macromolecules 1989, 22, 3555-3561. (32) Boots, H.; Dotson, N. The Simulation of Free-Radical Crosslinking Polymerization: The Effect of Diffusion. Polym. Commun. 1988, 29, 346354. (33) Bowman, C. N.; Peppas, N. A. A Kinetic Gelation Method for the Simulation of Free-Radical Polymerizations. Chem. Eng. Sci. 1992, 47, 1343-1351. (34) Bowman, C. N.; Peppas, N. A. Initiation and Termination Mechanisms in Kinetic Gelation Simulations. J. Polym. Sci., Part A: Polym. Chem. 1991, 29, 1575-1583. (35) Ward, J. H.; Peppas, N. A. Kinetic Gelation Modeling of Controlled Radical Polymerization. Macromolecules 2000, 33, 5137-5142. (36) Otsu, T.; Yoshida, M. Role of Initiator-Transfer Agent-Terminator (Iniferter) in Radical Polymerizations: Polymer Design by Organic Disulfides and Iniferters. Makromol. Chem., Rapid Commun. 1982, 3, 133140. (37) Hutchison, J. B.; Anseth, K. S. Off-Lattice Approach to Simulate Radical Chain Polymerizations of Tetrafunctional Monomers. Macromol. Theory Simul. 2001, 10, 600-607. (38) Nosaka, M.; Takasu, M.; Katoh, K. Characterization of Gels by Monte Carlo Method Using a Model of Radical Polymerization with Cross Linkers. J. Chem. Phys. 2001, 115, 11333-11338. (39) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. CHARMM: A Program for Macromolecular Energy, Minimization, and Dynamics Calculations. J. Comput. Chem. 1983, 4, 187-217. (40) MacKerell, A. D.; Wiorkiewiczkuczera, J.; Karplus, M. An AllAtom Empirical Energy Function for the Simulation of Nucleic Acids. J. Am. Chem. Soc. 1995, 117, 11946-11975. (41) MacKerrell, A. D.; Bashford, J. D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; JosephMcCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, I. W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586-3616. (42) Pearlman, D.; Case, D. A.; Caldwell, J. W.; Ross, W. S.; Cheatham, I. T. E.; DeBolt, S.; Ferguson, D.; Seibel, G.; Kollman, P. AMBER, a Package of Computer Programs for Applying Molecular Mechanics, Normal Mode Analysis, Molecular Dynamics and Free Energy Calculations to Simulate the Structural and Energetic Properties of Molecules. Comput. Phys. Commun. 1995, 91, 1-41. (43) Yin, D.; MacKerrell, J. A. D. Combined ab Initio/Empirical Approach for Optimization of Lennard-Jones Parameters. J. Comput. Chem. 1998, 19, 334-348. (44) Chen, I. J.; Yin, D.; A. D. M., Jr. Combined ab Initio/Empirical Approach for Optimization of Lennard-Jones Parameters for Polar-Neutral Compounds. J. Comput. Chem. 2002, 23, 199-213. (45) Glendening, E. D.; Reed, A. E.; Carpenter, J. E.; Weinhold, F. NBO Version 3.1. (46) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A., Jr.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B. G.; Chen, W.; Wong, M. W.; Andres, J. L.; Head-Gordon, M.; Replogle, E. S.; Pople, J. A. Gaussian 98, revision x.x; Gaussian, Inc.: Pittsburgh, PA, 1998. (47) Piletsky, S. A.; Karim, K.; Piletska, E. V.; Day, C. J.; Freebairn, K. W.; Legge, C.; Turner, A. P. F. Recognition of Ephedrine Enantiomers by Molecularly Imprinted Polymers Designed Using a Computational Approach, Analyst 2001, 126, 1826-1830. (48) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Liquid Water. J. Chem. Phys. 1983, 79, 926-935. (49) Liu, H.; Mueller-Plathe, F.; Gunsteren, W. F. v. A Force Field for Liquid Dimethyl Sulfoxide and Physical Properties of Liquid Dimethyl Sulfoxide Calculated Using Molecular Dynamics Simulation. J. Am. Chem. Soc. 1995, 117, 4363-4366. (50) Luzar, A.; Chandler, D. Structure and Hydrogen Bonding of WaterDimethyl Sulfoxide Mixtures by Computer Simulations. J. Chem. Phys. 1993, 98, 8160-8173.

Ind. Eng. Chem. Res., Vol. 46, No. 19, 2007 6091 (51) Escobedo, F. A.; Pablo, J. J. D. Molecular Simulation of Polymeric Networks and Gels: Phase Behavior and Swelling. Phys. Rep. 1999, 318, 85-112. (52) Gromov, D. G.; Pablo, J. J. d. Phase Behavior of Polymer-Solvent Mixtures. Fluid Phase Equilib. 1998, 150-151, 657-665. (53) Meredith, J. C.; Sanchez, I. C.; Johnston, K. P.; de Pablo, J. J. Simulation of Structure and Interaction Forces for Surfaces Coated with Grafted Chains in a Compressible Solvent. J. Chem. Phys. 1998, 109, 64246434. (54) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N‚ Log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089-10092. (55) Nelson, M.; Humphrey, W.; Gursoy, A.; Dalke, A.; Kale´, L.; Skeel, R.; Schulten, K. NAMDsA Parallel, Object-Oriented Molecular Dynamics Program. Int. J. Supercomput. Appl. High Perform. Comput. 1996, 10, 251268. (56) Kale´, L.; Skeel, R.; Bhandarkar, M.; Brunner, R.; Gursoy, A.; Krawetz, N.; Phillips, J.; Shinozaki, A.; Varadarajan, K.; Schulten, K. NAMD2: Greater Scalability for Parallel Molecular Dynamics. J. Comput. Phys. 1999, 151, 283-312. (57) Phillips, J. C.; Zheng, G.; Kumar, S.; Kale´, L. V. NAMD: Biomolecular Simulation on Thousands of Processors. In Proceedings of the IEEE/ACM SC2002 Conference, Baltimore, MD, November 16-22, 2002; IEEE Computer Society Press: Washington, DC, 2002; pp 234236. (58) Frigo, M.; Johnson, S. G. FFTW: An Adaptive Software Architecture for the FFT. In Proceedings of the IEEE ICASSP Conference, Seattle, WA, 1998; pp 123-128. (59) Frenkel, D.; Smit, B. Understanding Molecular Simulation; Academic Press: San Diego, CA, 1996.

(60) Reddy, M. R.; Erion, M. D.; Agarwal, A. Free Energy Calculations: Use and Limitations in Predicting Ligand Binding Affinities. In ReViews in Computational Chemistry; Lipkowski, K. B., Boyd, D. B., Eds.; Wiley-VCH: New York, 2000; pp 134-137. (61) Jorgensen, W. L.; Briggs, J. M. A Priori pKa Calculations and the Hydration of Organic Anions. J. Am. Chem. Soc. 1989, 111, 4190-4197. (62) Rastelli, G.; Constantino, L.; Gamberini, M. C.; Corso, A. D.; Mura, U.; Petrash, J. M.; Ferrari, A. M.; Pacchioni, S. Binding of 1-Benzopyran4-one Derivatives to Aldose Reductase: A Free Energy Perturbation Study. Bioorg. Med. Chem. 2002, 10, 1427-1436. (63) Fox, T.; Scanlan, T.; Kollman, P. A. Ligand Binding in the Catalytic Antibody 17e8. A Free Energy Perturbation Calculation Study. J. Am. Chem. Soc. 1997, 119, 11571-11577. (64) Erion, M. D.; van Poelje, P. D.; Reddy, M. R. Computer-Assisted Scanning of Ligand Interactions: Analysis of the Fructose 1,6-Bisphosphatase-Amp Complex Using Free Energy Calculations. J. Am. Chem. Soc. 2000, 122, 6114-6115. (65) Oral, E.; Peppas, N. A. Responsive and Recognitive Hydrogels Using Star Polymers. J. Biomed. Mater. Res. 2004, 68A, 439-447. (66) Oral, E.; Peppas, N. A. Dynamic Studies of Molecular Imprinting Polymerizations. Polymer 2004, 45, 6163-6173. (67) Vaughan, A.; Sizemore, S. P.; Byrne, M. E. Enhancing Molecularly Imprinted Polymer Binding Properties Via Controlled/Living Radical Polymerization and Reaction Analysis. Polymer 2007, 48, 74-81.

ReceiVed for reView October 24, 2006 ReVised manuscript receiVed March 24, 2007 Accepted March 30, 2007 IE061369L