A Functional Monomer Is Not Enough: Principal Component Analysis of the Influence of Template Complexation in Pre-Polymerization Mixtures on Imprinted Polymer Recognition and Morphology

In this report, principal component analysis (PCA) has been used to explore the influence of template complexation in the pre-polymerization phase on template molecularly imprinted polymer (MIP) recognition and polymer morphology. A series of 16 bupivacaine MIPs were studied. The ethylene glycol dimethacrylate (EGDMA)-crosslinked polymers had either methacrylic acid (MAA) or methyl methacrylate (MMA) as the functional monomer, and the stoichiometry between template, functional monomer and crosslinker was varied. The polymers were characterized using radioligand equilibrium binding experiments, gas sorption measurements, swelling studies and data extracted from molecular dynamics (MD) simulations of all-component pre-polymerization mixtures. The molar fraction of the functional monomer in the MAA-polymers contributed to describing both the binding, surface area and pore volume. Interestingly, weak positive correlations between the swelling behavior and the rebinding characteristics of the MAA-MIPs were exposed. Polymers prepared with MMA as a functional monomer and a polymer prepared with only EGDMA were found to share the same characteristics, such as poor rebinding capacities, as well as similar surface area and pore volume, independent of the molar fraction MMA used in synthesis. The use of PCA for interpreting relationships between MD-derived descriptions of events in the pre-polymerization mixture, recognition properties and morphologies of the corresponding polymers illustrates the potential of PCA as a tool for better understanding these complex materials and for their rational design.

However, the relationships between the factors influencing the recognition qualities of MIPs are difficult to study, since evaluation with conventional univariate methods is limited. This limitation can be overcome by the use of multivariate chemometric methods, in which mathematical and statistical methods are applied to chemical data. The methods allow for the simultaneous assessment of several factors, as opposed to the more common single variable studies of univariate methods. One multivariate method is principal component analysis (PCA), a data analysis technique for the reduction of multi-dimensional datasets to lower dimensions [30]. A PCA can be used to discover patterns in complex data matrices and can reveal relationships between measurements and variables, as well as between the variables themselves, which otherwise are likely to remain undetected [31].
We have recently reported on the influence of polymer composition on morphology and template recognition [27]. The polymer systems used in the present study were copolymers of the crosslinking agent ethylene glycol dimethacrylate (EGDMA) prepared with either methacrylic acid (MAA) or methyl methacrylate (MMA) as the functional monomer and imprinted with the local anesthetic, bupivacaine. In total, 16 MIP and corresponding non-imprinted reference polymer (REF) systems were included in the analysis. The polymers were prepared with different molar ratios of the respective monomers, though the stoichiometries of the MAA and MMA polymer series were equivalent (Tables 1 and 2). The choice of MMA as a functional monomer was based upon an interest in observing the impact of removing the capacity of the functional monomer to donate hydrogen bonds. PCA was applied to data extracted from MD simulations of the respective all-component prepolymerization mixtures, morphology characterization (surface area and porosity) and recognition behavior of the corresponding bulk polymers.  The analysis revealed that the molar fraction functional monomer in the MAA-polymers contributed to describing both the binding and morphology features. Moreover, competition in the hydrogen bond contact to the template between the functional monomer and the crosslinker in the pre-polymerization phase was revealed in both polymer series. Interestingly, weak correlations between the swelling behavior and the rebinding characteristics of the MAA-MIPs were observed. Furthermore, the results highlight the potential of using MMA as a non-crosslinking analogue of EGDMA for future studies on the impact of the degree of cross-linking on the behavior of EGDMA-based polymers.

Results and Discussion
In this study, the influence of template-monomer complexation in the pre-polymerization phase on polymer-template recognition and morphology was investigated through the use of PCA. The polymer systems studied were bupivacaine-imprinted EGDMA-copolymers containing different amounts of either MAA or MMA as the functional monomer. The analyzed data were obtained from equilibrium binding experiments (Tables S1 and S2), gas sorption measurements (Tables S3 and S4), swelling studies (Tables S5 and S6)  The variables describing the data extracted from MD simulation trajectories were measures of the hydrogen bond contact between the molecules present expressed as the percentage of the total simulation time (Tables S11 and S12): the average hydrogen bond contact between bupivacaine and the respective functional monomer MAA or MMA (%Bup-funcM), the average hydrogen bond contact between bupivacaine and the crosslinker EGDMA (%Bup-EGDMA), bupivacaine self-association (%Bup-Bup), contact between the amide proton functionality (HAB) of bupivacaine and the carbonyl oxygens (OAD) of either MMA or MAA (%HAB-OAD), as well as the contact between the amide proton of bupivacaine and a carbonyl oxygen (OAF) of EGDMA (%HAB-OAF). Both the total average hydrogen bond contact between the different molecules and the contact between the specific interaction points on the individual molecules were employed as variables in the PCA, since the total average hydrogen bond contact is the sum of the specific interactions between the individual molecules. From Tables S11 and S12, it can be seen that the hydrogen bond contact between the specific interaction points on the individual molecules contributes differently to the total average contact between the molecule species, in particular to the total average contact between the template bupivacaine and the two different functional monomers, MAA and MMA. In a previous report, it was suggested that the strength of the hydrogen bond contact between the amide proton functionality of bupivacaine and the carbonyl oxygen of MAA was steering both binding capacity and the specific binding in a series of MAA-polymers [27]. See Figure 2 for atom designations.
In PCA, multivariate data are projected onto a low-dimensional plane, where, e.g., groups of observations and trends can be revealed. The different principal components (PCs) are vectors in the n-dimensional variable space computed in such a way as to best approximate the data in the least squares sense. The first principal component (PC1) passes through the average point (the origin) in the data space pointing in the direction of maximum variance in the variable dataset. The second PC (PC2) is orthogonal to PC1, passes also through the origin and points in the direction of the second maximum variance in the dataset. A plot of the residual variance is shown in Figure S1. The two PCs define a two-dimensional plane onto which each data point is projected with a new coordinate value. The projections of the observed data (measured values of a sample) are called scores, and the projections of the variables are called loadings. Accordingly, the loading plot ( Figure 1A) reveals correlations between the variables, and the score plot ( Figure 1B) describes correlations between the observed data. Variables with high loadings (far from the origin) describe most of the variation in the data. The variables that are located closely together are highly correlated [30,31]. If they share the same sign, the correlation is positive, and if they have opposite signs, the correlation is negative. The score plot describes the measured properties of the samples (polymers), where grouping of the scores indicates that the samples share the same characteristics [30,31]. The morphological characteristics of the polymers (surface area and pore volume) together with the binding properties, the molar fraction of MAA and the hydrogen bond contact between functional monomer and bupivacaine in the pre-polymerization phase account for most of the variation in the data underlying PC1 (55%). The evident separation of the loadings corresponding to the binding properties and those corresponding to the morphological features in PC1 demonstrates the importance of taking morphological aspects into account when characterizing the binding performance of MIPs.
The loading describing the averaged bupivacaine-functional monomer hydrogen bond interactions (%Bup-funcM) is located close to that describing the molar fraction MAA in the polymers (MfracMAA), suggesting strong positive correlation between these factors. This can be explained by the fact that the averaged hydrogen bond contact between MAA and bupivacaine is higher than that between MMA and bupivacaine due to the additional interaction points provided by the acidic proton at MAA's carboxyl functionality (Figure 2). Further, it can be seen from the loading plot ( Figure 1A) that the molar fraction MAA in the polymers is positively correlated to the binding of bupivacaine to both the MIPs and the REFs in accordance with the results presented previously [27], implying higher binding with increased molar fraction MAA. The molar fraction MMA, in contrast, is negatively correlated to the variables describing binding and positively correlated to surface area and pore volume. However, the loading of MfracMMA on PC1 is lower compared to the loading of MfracMAA, indicating that the molar fraction MMA in the polymers does not contribute as much as the molar fraction MAA to the variation in the dataset. It is worth noting that the descriptors of the surface area for both the MIPs and the REFs derived with the BET and the Langmuir methods contribute equally to the variation in the data. Moreover, weak positive correlations between the swelling behavior and the rebinding characteristics of the MAA-MIPs were revealed, as the corresponding variables are located in the same quadrant of the loading plot. PC2 denotes 16% of the variation in the data that can predominantly be ascribed to the interactions between functional and crosslinking monomers and the template in the pre-polymerization phase. These interactions arise from hydrogen bond contacts between the carbonyl oxygens of MAA, MMA or EGDMA and the amide proton of bupivacaine (Figure 1, Tables S11 and S12). The results clearly demonstrate the presence of competition for hydrogen bonding to the template between the carbonyl interaction points of the respective functional monomer, MAA or MMA, and EGDMA, since the descriptors of these interactions are positioned in diagonally located quadrants of the plot, thus being negatively correlated. This provides a unique illustration of the diversity of the interactions between pre-polymerization mixture components and the template. Interestingly, the binding of the template to the MIPs (B/TMIP) and the specific binding (B/TSpec) were separated in PC2 from the binding to the non-imprinted reference polymers (B/TREF).
In the score plot, three groupings could be identified ( Figure 1B). The groupings of the scores are related to which functional monomer was used during synthesis and also to its molar fraction. Interestingly, the scores corresponding to the MMA polymers (Polymers 10-15, red scores) and those corresponding to the polymer not containing any functional monomer (Polymer 0, black scores) can be considered as one group, which demonstrates that the MMA polymers and the polymer prepared without the functional monomer share similar properties. These polymers are correlated to the descriptors of morphology where the spread of their scores in PC1 is of minor magnitude. In contrast, the spread of the scores corresponding to these polymers in PC2 is more pronounced. This can be ascribed to the increasing molar fraction MMA in the polymers. These results indicate that the increase of the molar fraction MMA in the polymers affects the morphological features (surface area and pore volume) only to a minor extent. Further, the similarity of the properties of the MMA polymers and the polymer prepared without a functional monomer suggests that MMA may be considered as a non-crosslinking analogue of EGDMA, of interest for studies of the influence of the degree of crosslinking in these polymers, a topic of on-going work.
In the case of the MAA polymers ( Figure 1B, blue and green scores), the scores were spread in PC1 according to the increasing molar fraction MAA in the polymers, the positively correlated rebinding behavior to both the MIPs and the REFs and the positively correlated hydrogen bond contact between template and functional monomer. Polymer 9, containing the highest molar fraction MAA within the studied polymer series, displayed the highest score values in PC1. This polymer has previously been shown to display the highest binding capacity and the most frequent hydrogen bond contact between template and functional monomer in the pre-polymerization phase [27]. The distribution of the scores in relation to the location of the loadings suggests that binding to both the MIPs and the REFs, as well as the specific binding, increases with increasing molar fraction MAA in the polymers, in accordance with results presented previously [27].
In summary, the polymers containing MMA as a functional monomer and the one prepared without functional monomer share the same properties. The rebinding capacity of these polymers can be interpreted as equally low, and the morphological characteristics, such as the surface area and pore volume, were comparable and rather independent from the molar fraction MMA used in the synthesis. The molar fraction of functional monomer in the MAA-polymers, in contrast, contributed to describing both binding and morphology features. Interestingly, the analysis exposed competition in the hydrogen bond contact to the template between the functional monomer and the crosslinker in both the simulated MMA-and MAA pre-polymerization mixtures.

Data Analysis
The PCA was performed using the software package The Unscrambler v.10.2 (Camo Software AS, Oslo, Norway). A total of 2,964 data points was used in the study. For the binding study 9 replicates were performed; for the morphology experiments 1 batch of bulk polymer was used for each composition; and the MD data were extracted from molecular trajectories obtained from five simulations on each system. Prior to data analysis, the variables were mean centered and then scaled using the standard deviation in order to provide equal weighting to the variables.

Molecular Dynamics (MD) Simulations
All-atom MD simulations of the polymer systems were performed as described previously [25][26][27] using the AMBER (v.10.0 UCSF, San Francisco, CA, USA) platform of programs [49,50]. The simulated pre-polymerization mixtures differed in composition through variations of the molar ratio of the monomers employed. The molar fraction functional monomer (either MMA or MAA) was increased over the range of simulated molar ratios, while the other components where held effectively constant. For detailed information regarding the compositions of the simulated all-component pre-polymerization mixtures, as well as for equilibration and production run data, see Tables S7-S10. All systems were simulated in quintuplicate, each covering 10 ns of recorded trajectory data for each mixture (totally 50 ns). Final equilibrated trajectories were analyzed using the PTRAJ module implemented in AmberTools (v.1.3 UCSF, San Francisco, CA, USA). Hydrogen bond interactions were extracted from the trajectories using a cut-off distance and angle of 3.0 Å and 120, respectively. The structures of the molecular species and the designations of atoms potentially participating in hydrogen bond interactions are presented in Figure 2.
For the average hydrogen bond occupancies for all analyzed atom-pairs in each system and the corresponding averaged lifetimes, see Tables S11 and S12.

Polymer Synthesis
A series of bupivacaine-imprinted polymers containing either MMA or MAA as the functional monomer were synthesized and treated as described previously [27,51]. Corresponding non-imprinted reference (REF) polymers were prepared following the same procedure. The compositions of the polymers (Tables 1 and 2) and their stoichiometries were equivalent to those of the pre-polymerization mixtures simulated by MD (Tables S7 and S8).

Equilibrium Binding Study
Binding to both the MIP and the REF polymers was studied at a polymer concentration of 0.75 mg/mL. The radioligand, [ 3 H]-(R,S)-bupivacaine (15 pmol/mL), and polymer suspended in toluene were mixed in a total volume of 1 mL. The samples were incubated on a rocking table at 293 K for 3 h. Thereafter, they were centrifuged (7,000× g for 5 min), and the supernatant (600 μL) was mixed with the scintillation cocktail (2 mL, Beckman Ready Safe, Beckman Caulter, Bromma, Sweden). The activity was then measured for 2 min by scintillation counting (Packard Tri-Cab 2100TR liquid scintillation counter, PerkinElmer, Zaventem, Belgium). Control samples containing no polymer were also prepared and treated identically. All samples were analyzed in triplicate (Tables S1 and S2).

Swelling Studies
A volume of 1 mL dry polymer (V0) was measured in a glass cylinder, and after the addition of an excess of toluene, the cylinder was sealed. The volume of the swollen polymer (VSW) was recorded after 24 h. The swelling ratio (SW) was calculated according to Equation (1) (Tables S5 and S6).

Examination of Gas Accessible Surface Areas and Porosities
Polymer surface area and porosity were examined with nitrogen sorption studies by the Brunauer, Emmet and Teller [52] (BET), the Langmuir [53] and the Barrett, Joyner and Halenda [54] (BJH) methods (Tables S3 and S4). Prior to the measurements, samples were degassed at 323 K for 24 h to remove adsorbed gases and moisture. BET-and Langmuir surface areas were calculated from the adsorption data using 0.162 nm 2 as the molecular cross-sectional area for adsorbed nitrogen molecules. The BJH method was applied to calculate the pore volumes of the polymers from the desorption branch of the isotherms. All measurements were performed on an ASAP 2004 instrument (Micromeritics, Norcross, GA, USA) at 77 K.

Conclusions
In this study, a PCA was performed on data extracted from all-component MD simulation trajectories of a series of molecularly imprinted pre-polymerization mixtures together with data describing the surface characteristics and rebinding behaviors of analogous synthesized bulk polymers. Unique insights concerning the influence of template-monomer complexation on recognition were gained.
The polymers prepared with MMA as a functional monomer were, through the analysis, revealed to share the same properties as the polymer prepared without any functional monomer, which highlights the potential of the use of MMA as a non-crosslinking analogue of EGDMA. The rebinding capacities of these polymers were equally low, and the morphological characteristics, such as surface area and pore volume comparable and rather independent of the molar fraction MMA used in synthesis.
The molar fraction functional monomer in the MAA-polymers, on the other hand, contributed to describing both the binding, surface area and pore volume. Interestingly, weak positive correlations between the swelling behavior and the rebinding characteristics of the MAA-MIPs were revealed.
Finally, the results presented in this report demonstrate the potential of multivariate methods for analyzing the characteristics of complex materials, e.g., MIPs, where unique insights of importance for the further development of these materials may be gained