Tine
Curk‡
ab,
Jure
Dobnikar
*a and
Daan
Frenkel
b
aInternational Research Center for Soft Matter, Beijing University of Chemical Technology, Beijing, China. E-mail: jd489@cam.ac.uk
bDepartment of Chemistry, University of Cambridge, Cambridge, UK
First published on 17th September 2015
Molecular imprinting is the process whereby a polymer matrix is cross-linked in the presence of molecules with surface sites that can bind selectively to certain ligands on the polymer. The cross-linking process endows the polymer matrix with a chemical ‘memory’, such that the target molecules can subsequently be recognized by the matrix. We present a simple model that accounts for the key features of this molecular recognition. Using a combination of analytical calculations and Monte Carlo simulations, we show that the model can account for the binding of rigid particles to an imprinted polymer matrix with valence-limited interactions. We show how the binding multivalency and the polymer material properties affect the efficiency and selectivity of molecular imprinting. Our calculations allow us to formulate design criteria for optimal molecular imprinting.
MIPs have been used in applications such as solid-phase extraction,11 chiral separation,12 and catalysis.13 They can act as molecular sensors,8,10,14 and mimic antibodies or enzymes.8 They can selectively bind drugs,15–17 proteins,18 or even whole bacteria.19,20Fig. 1 shows a schematic representation of the imprinting and subsequent recognition process. The efficiency of the molecular recognition process depends on a number of parameters: (1) the initial ligand concentration c, (2) the template–ligand binding affinity KD, and (3) the stiffness of the polymer matrix kh.
Clearly, it is important to maximize the selectivity of MIPs, but in experiments MIPs are often optimized by trial and error. In fact, the theoretical picture is rather fragmented as existing theoretical models for MIPs do not consider the imprinting process as a whole, but rather tend to focus on individual steps in their mode of action.21–27 Moreover, the atomistic and coarse-grained simulations of molecular imprinting that have been reported28–33 focused mostly on specific MIPs and did not explore generic trends that would allow us to arrive at general design principles.
Here, we present a generic, coarse-grained statistical mechanics model that captures the key features of molecular recognition. The model provides an integrated description of the MIP formation process and of the subsequent binding of analytes. For the simplest case of divalent particles (i.e. particles with 2 receptors), we derive analytical expressions for the binding free energy and, from that, the adsorption isotherm of analytes on a MIP. For the general case of multivalent particles, we use Monte Carlo simulations to obtain the adsorption isotherms – representative snapshots from the simulations are shown in Fig. 3 below. Two important measures of the quality of a MIP are: (i) how much more efficient is binding of a given analyte to imprinted (MIP) than to a non-imprinted matrix (NIP), and (ii) how well can we separate two analytes that are only slightly different (e.g., of the same size and with the same number – but different spatial patterns – of the receptors). In order to address these questions, we evaluate standard measures of MIP specificity, such as the imprinting factor (IF) and the separation factor (SF), as a function of matrix and analyte properties and identify the optimal range of the control parameters such as the template and ligand concentrations and polymer matrix stiffness. Our work provides insight into the generic features of MIP operation and leads to a set of simple design principles.
In what follows we use the standard notation for the inverse temperature: β ≡ 1/kBT with kB the Boltzmann constant and T the absolute temperature. Hybridization free energies for ligand–receptor bonds can be deduced from experimental data for binary association in solution, which yield the dissociation constant KD = ρ0eβΔG, where ρ0 = 1 M is the standard concentration. We note that the ligand–receptor dissociation constant in the pre-polymerization solution D can in principle be different from the one in a formed MIP KD if the conditions such as solvent, pH, salt concentration, or temperature are different. Corrections for the fact that the ligands and receptors have an excluded volume and are tethered to a surface can be computed (see ref. 35 and 36).
Using standard chemical equilibrium theory we can compute the amount of ligands adsorbing to the templates in step A, which depends on the ligand dissociation constant under the conditions of the imprinting phase, D, on the concentrations of templates CT and on c, the original concentration of ligands in solution (see ESI†). The fractional occupancy of receptors fr (the probability that a given receptor is bound to a ligand) is
(1) |
In our model we must account for the effect of the deformability of the matrix and for the directionality of the ligand–template interaction. To this end, we model MIPs as soft, deformable matrices that contain cavities imprinted by specific ligands (Fig. 1(C)). The ligands grafted to the matrix can fluctuate around their equilibrium positions due to the thermal fluctuations. We call the equilibrium positions of the imprinted ligands their “anchors”. In order to bind to receptors on the analyte, the ligands need to be displaced from their anchors, which increases the elastic free-energy of the matrix by an amount Uh. Following Rubinstein and Colby,37 we assume that the ligand–matrix interaction only depends on the distance between the receptor and the anchor and that it can be replaced by a harmonic spring:
(2) |
(3) |
We see that kh depends weakly on the template (cavity) size σ because the relative fluctuations of 2 ligands will decrease if the two ligands are closer than the mesh size . For flexible polymer gels Me ≈ kBT/3 (ref. 37) and large particles σ > the spring constant is determined simply by the cross-linking distance kh ≈ πkBT/2. Eqn (3) is important because, as we show below, it allows us to predict the effect of the stiffness of the matrix on the selectivity of MIPs. We will focus on the case where particles are rigid, however, the present model can be applied also to soft particles (such as proteins) where receptor positions on the particle itself are fluctuating (fluctuations characterised by spring constant kparth). In this case we could map (to the first order) a soft particle and a soft matrix to the current model of a hard particle and an “effective” softer matrix 1/keffh = 1/kh + 1/kparth.
In what follows we first focus on a single cavity system and calculate the free energy of binding an analyte. Afterwards we extend the picture to the whole polymer matrix and calculate corresponding binding affinities and adsorption isotherms.
The position rp and orientation Ωp of a rigid particle uniquely define the positions of all binding sites on its surface rbsi. When the particle approaches the cavity, bonds between the binding sites i and the ligands j can be formed – forming a bond results in a lowering of the free energy by an amount equal to the hybridization free energy ΔGij, but forming a bond also costs free energy because (in general) the ligand must be displaced from its equilibrium position rancj to rligj = rbsi: the corresponding free energy cost is given by (2). To compute the overall binding free energy, we need to compute the ratio of the partition function for the case where one or more ligands are bound to the particle, to the one for the case with no ligands bound.
This expression does not yet include the partition function of the unbound particles: this will later be accounted for through the chemical potential of the free particles. The partition function depends on the distribution of the binding sites on the particle rbs and on the arrangement of ligand anchors in the cavity ranc. To obtain the partition function we must sum over all particle positions and orientation and over all possible bonding arrangements:
(4) |
(5) |
(6) |
As a single bond cannot create a static stress in the matrix, q1 does not depend on the matrix stiffness kh. A similar result was reported by Tanaka et al. for the case of imprinted hydrogels.25 The simplest non-trivial term is the two-bond partition function. For any chosen combination of two binding sites and two ligands , there are two possible ways of forming two bonds. The total two-bond partition function in a system where there is more than one way to make two bonds is a sum of all possible two-bond pairs ij, i′j′,
(7) |
(8) |
Fig. 2 (a) Multivalency. The cavity occupancy fcav as a function of the imprinting mismatch b/a. The red solid line represents analytical calculations (4)–(9). The symbols depict the results of Monte Carlo simulations: blue circles for the divalent, green squares for the hexavalent case. Parameters: matrix stiffness , analyte concentration in solution , and the bond energies in the divalent/hexavalent case , . (b) Incomplete cavities. The specific binding free energy of a hexavalent analyte as a function of the number of ligands imprinted in the cavity. The results are for two values of (red: soft gel, green: stiff gel) and for two values of ligand binding strength ΔG* (circles: weaker bonds, squares: stronger bonds). |
Fig. 3 Simulation snapshots. (a) Imprinting phase: ligands (small balls, red when free and blue when bound) bind to the divalent templates (orange balls with cyan receptors). The ligand positions are subsequently frozen and templates extracted to form the cavities (MIP) shown in (b) and (c) as transparent blobs. (b) Re-binding of templates to MIP. (c) Binding of templates to the non-imprinted matrix (NIP). Snapshots show only a small part of the simulated system. (d) Adsorption isotherms. Average number of bound divalent (tetravalent in the inset) analytes depending on the concentration of analytes in solution . Solid lines depict the analytical prediction (12)–(14), and the symbols are for simulation results (black circles: re-binding of templates (b = a) to MIPs as shown in the snapshot (b)), green diamonds: binding of different analytes to MIPs, and red squares: binding of templates to a NIP – as shown in the snapshot (c). The parameters: , c* = 0.02, CT = 0.01, . The snapshots represent a configuration of bound analytes at with the system size V = (12.6σ)3. Inset in (d) shows isotherms of tetravalent particles binding at . |
In the limit of soft matrices kh ≪ kBT/b2, where thermal fluctuations are greater then the analyte size, the specificity towards analyte geometry is lost . This is the regime where only the number of ligands in the cavity (or the number of receptors on the analyte) remains important, not their spatial positions. In this limit the rotational degrees of freedom can be neglected, such a model of molecular imprinting has been considered in ref. 25 and 26.
If we assume that all binding cavities are identical and independent, the number of adsorbed analyte particles is determined by the binding free energy F and by the chemical potential μ of the analyte. We recover the simple Langmuir expression for the fraction of occupied cavities:
(9) |
MIPs that have been cross-linked in the presence of specific template particles will adsorb analytes that have a structure similar to the template: the greater the similarity between the template and the target particle, the larger the average occupancy fcav of the cavities imprinted by the template particles. In the case of particles with two binding sites, we can compute fcav analytically (eqn (4)–(9)). In Fig. 2(a) we compare these analytical calculations with Monte Carlo simulations. The good agreement between analytical results and simulations allows us to validate the Monte Carlo approach. For a particle with 6 binding sites, the analytical calculations are no longer tractable, but the MC simulations in Fig. 2 show that imprinting with a higher valency leads to a much stronger discrimination. In both cases, the cavities have a fixed ligand distribution that is imprinted by the template particle. The analytes are assumed to have the same geometry as the template but their sizes are rescaled by a factor b/a (the “imprinting mismatch”). Not surprisingly, the average occupancy fcav has a maximum at b/a ≃ 1. The difference between two and six binding sites shows up when we increase the mismatch: for the case of two binding sites, a mismatch of 30% decreases fcav only by a factor two. In contrast, for particles with six binding sites, a 30% mismatch leads to a decrease of fcav by more than an order of magnitude.
Finally, Fig. 2(b) illustrates how the binding depends on the number of ligands in the cavity. There are many practical examples where the number of ligands in the cavity is less than the number of receptors on the template. Obviously, we expect to observe ‘under-coordinated’ cavities in case the cavities are only partially formed. But even if cavities are initially well formed, they may become under-coordinated if the MIP is ground after imprinting (during the grinding, which is a common procedure to enhance the accessibility of cavities in stiff matrices,17 a fraction anchored ligands is likely to be detached from the cavities). Not surprisingly, we find that the binding strength increases with the number of ligands. The dependence is steepest for stiffer matrices (larger kh).
BAcav ≡ Ccav〈KcavA〉 | (10) |
(11) |
In the case of NIPs, there are no imprinted cavities and the average in (10) must be taken over a random distribution of ligands in the cross-linked matrix. The result derived in the ESI,† is
(12) |
In the divalent case we can evaluate the free energy - and thus the binding affinities - analytically (eqn (4)–(8) and (12)). Assuming that imprinted cavities are randomly distributed throughout the MIP the cross-cavity term becomes equal to BAnip2, however, single bond terms need to be subtracted to prevent double counting. If the template extraction process is efficient the cavity concentration is equal to the template concentration in the imprinting phase Ccav = CT. The MIP binding affinity becomes
(13) |
In the more general case κ > 2, BAmip can be computed by numerical simulations. From eqn (9) and (10) we observe that for low concentrations the binding affinity essentially determines the ratio of concentrations of bound analytes CB to analytes free in solution CA
(14) |
The binding affinities BAnip and BAmip can be used to evaluate two key quantities used to assess the performance of MIPs in the literature: the imprinting factor IF(a), which is the ratio of the binding affinity of a template a to a MIP and to a NIP, and the separation factor SF(b1, b2; a) measuring the ability of a MIP (imprinted by a template a) to distinguish between two different analytes b1 and b2:
(15) |
Fig. 4 MIP design principles. (a) Schematic phase diagram of MIP efficiency, summarizing the design principles. The ligand concentration for which MIPs are efficient lies between the two limits: (17) and (16). The dashed line drawn somewhat arbitrarily at separates the regions of bond selectivity and geometrical recognition. In (b) and (c) the imprinting factor calculated for divalent templates is shown as a function of and c* at a stoichiometric ratio of ligands–receptors (b), and as a function of and c* at fixed matrix stiffness (c). We have assumed an equal binding strength in the imprinting and binding stage . |
(16) |
(17) |
In some applications it is desirable to separate analytes according to receptor composition but not according to the geometry. For instance, if targeting a whole class of analytes with similar receptor composition but varying (or unknown) patterns (e.g., viruses with high mutation rates), geometrical selectivity is unwanted. The optimal regime of matrix stiffness for such applications is the regime of bond selectivity: . In this regime imprinting makes a difference, however, the substrate is too flexible to allow for efficient recognition of the geometrical patterns. Hence, MIPs can only discriminate analytes with different ligand compositions.
In order to test the predictions that follow from Fig. 4(a), we compared them to the analytical results for the divalent binding model (12) and (13). To this end, we evaluated the imprinting factor IF(a) for a broad range of ligand concentrations c*, binding affinities , and matrix elasticities . Fig. 4(b) shows the dependence of IF(a; a) on the matrix stiffness and on the ligand concentration at a fixed stoichiometric ratio of ligands and templates in the imprinting phase . Furthermore Fig. 4(c) shows the dependence of IF on the stoichiometry at fixed matrix stiffness, showing that the stoichiometric ratio is close to optimal, the optimal concentrations being and . We have assumed that the dissociation constants for individual bonds remain the same in the imprinting and binding phase .
In accordance with our tentative design rules (Fig. 4), we observe a sharp increase in the imprinting quality in the parameter region where MIPs should function optimally, i.e. for intermediate values of ligand concentrations and stiff matrices .
Within the same parameter space we also computed the separation factor SF(a, b; a) for two slightly different analytes (b1 = a and b2 = 1.1a) on MIPs imprinted by a template with a separation a between the two receptors (Fig. 5(a)). The analyte separation is effective for matrix stiffness . This range can be extended if we compare analytes that are less similar. Again, the region where analyte separation is most efficient roughly coincides with the onset of geometrical recognition in Fig. 4(a).
Fig. 5 (a) Separation Factor. The separation factor SF(b1, b2; a) for the imprinted template b1 = a and an analyte with slightly larger inter-receptor distance b2 = 1.1a as a function of the parameters and c*. (b) Enhancing the selectivity. The best separation is achieved when the imprinting is slightly mismatched relative to the analyte. The separation factor SF(b1, b2; aopt) at the optimal value of the imprinted distance aopt as a function of and c*. Parameters are the same as in Fig. 4(b). |
Somewhat surprisingly, we find that the optimal separation of two analytes, say b1 = σ and b2 = 1.1σ, is not best achieved by imprinting with a template identical to the first analyte (a = b1). Rather, the separation factor SF(b1, b2; a) can be maximized by designing the cavity with an optimal imprinted distance aopt < b1. This result can be intuitively understood by noting that the binding free energy is approximately a quadratic function of the mismatch close to the minimum. If the imprinting is slightly mismatched, the binding affinity of the chosen analyte is slightly smaller, but at the same time it increases relative to the binding affinity of the other particle resulting in better separation capacity of the MIP. Fig. 5(b) displays the separation factor SF(b1, b2; aopt) for the same analytes as in Fig. 5(a) but with the template size a optimized at each point in the parameter space (see ESI†). We can clearly see an increase in the separation capability, however, the qualitative features of the phase diagram in Fig. 4(a) remain valid.
The first key observation is that the quality of imprinting depends on the concentration of ligands and templates in the imprinting phase and on the binding strength between them: the optimal imprinting is achieved with a nearly stoichiometric ratio of ligands vs. template receptors, e.g. for tri-valent templates the stoichiometric ratio of templates:ligands should be 1:3. Additionally, initial ligand concentration should be of a similar value as the corresponding bond dissociation constant. This provides the optimal tradeoff between, on the one hand, efficient cavity formation, and on the other hand, selective re-binding of analytes to imprinted cavities.
The second key observation is that stiffer matrices are more selective – suggesting that it should be beneficial to make the gel as rigid as possible. Consequently, polymer gels are not a suitable matrix for efficient molecular imprinting of small molecules: when geometrical recognition is required, stiffer systems such as glassy polymers should be considered. However, the non-specific terms in the free energy – which are not explicitly considered here – will have an opposite effect upon increasing the stiffness of the matrix: slowing down of the kinetics and impeding the particles' access to the cavities. A commonly implemented solution for this problem is grinding the stiff gels in order to expose the imprinted ligands. In this case, much stiffer matrices can be used, however, the procedure inevitably reduces the imprinting quality by grinding-off a fraction of ligands in the cavities. In specific MIP systems, it is likely that the opposing thermodynamic and kinetic trends result in an application-specific optimum gel stiffness.
Moreover, when imprinting soft macromolecules such as proteins or biopolymers onto a hard matrix, the intrinsic softness of the template has a similar effect to a reduced matrix stiffness (within our coarse-grained model a soft template on a hard substrate resembles a hard template on a softer substrate). Extremely large values of the matrix stiffness in our model are therefore unlikely to be relevant for experimental realizations. A reasonable choice of the matrix stiffness to design applications seems to be between , i.e. fluctuations of the ligand position relative to the size of the template between 10% and 30%. In such a case, for imprinting to work, the strength of the individual ligand–receptor bonds needs to be strong enough. Divalent templates cannot be imprinted effectively unless the bond dissociation constant is . This changes considerably if the templates are multivalent: in the tetravalent case imprinting can be achieved with relatively weak bonds , while in the hexavalent case the bonds can be extremely weak , which is in the realm of hydrogen bonds in aqueous solutions. This suggests that it should be possible to efficiently imprint molecules, such as proteins or drugs, onto polymers in aqueous solutions.
Footnotes |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c5sm02144h |
‡ Current institution: University of Cambridge. |
This journal is © The Royal Society of Chemistry 2016 |