Systematic Design of Molecularly Imprinted Polymers for Triclosan Using Design of Experiments and Molecular Dynamics Simulations

1. Introduction

The extensive application of antimicrobial compounds such as triclosan (TCS) in various fields has increased its ubiquitous occurrence in the environment as well as in biological organisms, thus increasing concern about the adverse effects that TCS can exert on ecosystems due to its toxicity, bioaccumulation capability, and possible negative impacts [1]. Triclosan is a hydrophobic and chlorinated aromatic compound bearing a phenol group, allowing it to undergo a variety of interactions, ranging from hydrogen bonding to hydrophobic and π-π interactions. The physicochemical characteristics of triclosan make it a very difficult compound to separate selectively from water, and its handling poses many challenges [2,3]. Within the available strategies for the selective recognition and elimination of target compounds, MIPs have been proven to be one of the most interesting materials due to their unique property of generating specific binding sites through polymerization processes with the incorporation of a template molecule. Such materials present enhanced selectivity, chemical robustness, and recyclability when compared to traditional sorbents such as activated carbon or other types of polymers [4]. In addition, MIPs showed superior performance compared to traditional sorbents in TCS-related studies [5].

In this context, methacrylic acid (MAA) is among the most popular functional monomers employed in MIP production. This is attributed to the fact that it creates interactions, which allow creating highly specific materials. Nevertheless, MAA systems are characterized by high density and low diffusion of polymers, limiting their performance in terms of adsorption kinetics and desorption. As a solution to the problem, hydrophilic co-monomers have been added to increase the polarity of a polymer, its swelling and accessibility of target molecules. For instance, 2-hydroxyethyl methacrylate (HEMA) is a good candidate for this purpose. Moreover, another important factor influencing MIP performance concerns the influence of the porogenic solvent used during MIP synthesis on its functionality, as well as pore structure and particle morphology [6,7]. Porogen-assisted polymerization has already been reported to impact adsorption kinetics and adsorption equilibrium via altering the internal structure of the polymer matrix. In triclosan adsorption systems, fast equilibrium processes and good adsorption capacity are commonly achieved with MIPs with developed porosity and high surface area.

While some previous research focused on one of the influencing factors, such as specific monomers or porosity, a comprehensive characterisation of the main contributing factors to develop efficient MIPs for triclosan is still missing. Most existing MIP systems were obtained based on trial-and-error procedures, which do not provide full understanding about the mutual impact between monomer composition and functional characteristics of the resulting MIP. Currently there is no study to systematically optimize all three main parameters simultaneously, i.e, functional monomer interactions (MAA), hydrophilicity (HEMA) and the effect of the porogen solvent (ACN). In this study, a design of experiments (DOE) method utilizing the Box-Behnken design was applied to explore the impact of the contents of MAA, HEMA, and porogens on the performance of triclosan imprinted polymers. The correlation among the compositions, morphologies, and performances of MIPs was explored by combining the synthetic process, structural analysis, and the adsorption process. The findings reveal that the optimal formula can provide a high imprinting factor (around 3.0) and release efficiency (around 65%-70% after 24 hours). This paper presents a systematic strategy for creating MIPs with adjustable adsorption and desorption capabilities, which will facilitate their applications in selective adsorption and controlled release of hydrophobic bioactive molecules like triclosan.

2. Materials and Methods

2.1. Materials

Triclosan (TCS) was used as the template. Methacrylic acid (MAA) and 2-hydroxyethyl methacrylate (HEMA) were employed as functional and hydrophilic comonomers, respectively. Trimethylolpropane triacrylate (TMPTA) served as the crosslinking agent, and azobisisobutyronitrile (AIBN) was used as the thermal initiator. Acetonitrile (ACN) was used as the porogenic solvent. All reagents were of analytical grade, purchased from Sigma and used without further purification.

2.2. Experimental Design (Doe)

To effectively explore the influence of significant formulation parameters on the behavior of the molecularly imprinted polymers (MIPs), a response surface method (RSM) was used. Three independent factors were chosen because of their importance in the efficiency of imprinting and structural characteristics of the produced polymers: proportion of the functional monomer (MAA, X₁), proportion of the hydrophilic comonomer (HEMA, X₂), and porogenic solvent volume (ACN, X₃). The content of the crosslinker (TMPTA) remained unchanged in all experimental trials. The chosen independent variables were varied using three-level factorial analysis, and the experimental design matrix was built following the principles of the Box-Behnken technique. In such a manner, the influence of independent variables and their interactions on the response could be efficiently explored. As the main response variable, the imprinting factor (IF) was used [8].

2.3. Synthesis of Molecularly Imprinted Polymers

The preparation of MIPs was performed using the bulk free radical polymerization technique. Triclosan was first dissolved in the acetonitrile solvent in the appropriate amount followed by the addition of MAA and HEMA according to the DOE matrix. The cross-linker TMPTA was then added to the above mixture, followed by stirring until it was homogenized in the system. Initiation of the reaction was carried out using AIBN in the amount of 0.030 g and nitrogen gas was passed in the mixture for 10 min to get rid of oxygen. The bulk polymerization and curing was done at 60^oC for 24 h and 70^oC for 2 h respectively. The triclosan was extracted using a mixture of methanol and acetic acid in a ratio of 9:1, v/v. Extraction was stopped when there was no trace of triclosan found in the solvent by UV-Vis spectroscopy.

2.4. Characterization

FTIR analysis was done to determine the presence of the polymer structure and removal of the template. The spectra were taken by observing the peaks related to the functional groups (C=O, O-H) and disappearance of the templates. Scanning electron microscope (SEM) was employed for the analysis of the surface morphology and particle structure. The samples were gold-coated prior to taking pictures. The particle size and morphological details were observed based on the SEM images. The triclosan concentration was measured by UV–Vis spectroscopy. Calibration curves were drawn for triclosan concentrations. A linear regression model was established: Y=0.0154X+0.193 (R^2=0.9914) where X=concentration and Y=absorbance.

2.5. Adsorption and Release Experiments

The batch adsorption experiment was done by adding a known weight of polymer to a solution containing triclosan with a known concentration. After equilibration, the polymer was then separated through centrifugation, after which UV-Vis spectroscopy was done for analysis of supernatant. Calculation of the adsorption capacity ($q_e$) is represented by:

$$q_e={\displaystyle \frac{(C_0-C_e)V}{m}}$$

Removal efficiency (%) was calculated as:

$$\%Removal={\displaystyle \frac{C_0-C_e}{C_0}}\times 100$$

The imprinting factor (IF) was calculated as:

$$IF={\displaystyle \frac{q_{e,MIP}}{q_{e,NIP}}}$$

For release experiments, triclosan loaded MIPs was suspended in water, followed by sampling at specific times using UV-Vis spectroscopy for analysis. Cumulative release (%) was computed from the amount of triclosan released as compared to the total amount loaded [9,10].

2.6. Molecular Dynamics Simulations

A pre-polymerization molecular dynamics model was performed containing triclosan as template, methacrylic acid (MAA) as functional monomer, 2-hydroxyethyl methacrylate (HEMA) as comonomer, and acetonitrile as porogenic solvent. The simulated mixture contained 1 triclosan molecule, 7 MAA molecules, 2 HEMA molecules, and 100 acetonitrile molecules. The ratios were picked to represent condition 1 with 0.1 g triclosan, 0.2 g MAA and 0.1 g HEMA. Small-molecule structures were generated from SMILES strings using RDKit, with explicit hydrogens added prior to 3D coordinate generation. The SMILES used were triclosan, Clc2cc(Cl)ccc2Oc1ccc(Cl)cc1O; methacrylic acid, CC(C(O)=O)=C; 2-hydroxyethyl methacrylate, O=C(OCCO)\C(=C)C; and acetonitrile, CC#N. Initial 3D conformers were generated with RDKit and were subsequently optimized with the MMFF force field [11,12]. Partial charges were assigned with the AM1-BCC method using the OpenFF toolkit [13,14,15,16].

All molecules were described with the Open Force Field small-molecule force field openff-2.2.1.offxml [15,16,17]. The full simulation system was assembled in OpenMM via the OpenFF toolkit using the parameterized component library [15,16,17,18,19]. Acetonitrile was modeled explicitly and parameterized in the same manner as the other small molecules. The initial configuration was generated in a periodic cubic box of side length 5.0 nm by random insertion of 1 triclosan, 7 MAA, 2 HEMA, and 100 acetonitrile molecules, using a minimum placement padding of 0.16 nm between molecular bounding spheres.

Molecular dynamics simulations were performed in OpenMM 8.5 [18,19]. A Langevin middle integrator was used with a time step of 1.0 fs, a friction coefficient of 1.0 ps^-1, and a production temperature of 300 K. Pressure was maintained at 1.0 bar using a Monte Carlo barostat with volume-move attempts every 25 steps. Periodic boundary conditions were applied in all three dimensions. The initial system was energy-minimized for 2000 iterations, followed by a 5 ps equilibration at 50 K, and a 100 ps equilibration at 300 K. The production trajectory was propagated for 30 ns. Coordinates were saved every 100 ps.

Because this starting configuration was intentionally loose, the system was allowed to relax under NPT conditions at 300 K and 1 bar. The box contracted rapidly during the early stage of the simulation and then fluctuated around a stable liquid-like density. Across the 30 ns production trajectory, the average density was 0.83 g mL^-1. The average density over the first 1 ns of the production run was 0.826 g mL^-1, and the final density was 0.824 g mL^-1, indicating that no substantial drift occurred during the simulation. The final cubic box length was 2.197 nm.

Trajectory analysis focused on triclosan-centered interaction patterns rather than exact structural microstates. All saved frames were processed using a custom Python analysis workflow based on MDTraj [20]. The trajectory was first imaged under periodic boundary conditions and aligned on triclosan to remove rigid-body translation and rotation before interaction analysis. Two triclosan interaction regions were defined around its two oxygen-centered sites: a phenolic oxygen site (O1) represented by atoms O1x, C1x, C2x, and C6x, and a bridging ether oxygen site (O2) represented by atoms O2x, C7x, C8x, and C12x. For each frame, binary interaction descriptors were assigned according to whether at least one MAA or HEMA molecule formed a contact with either triclosan site and whether hydrogen-bond geometries were present between triclosan and the monomer species.

Contact patterns were defined using site-specific heavy-atom proximity between triclosan and MAA or HEMA. Hydrogen-bond patterns were identified geometrically by evaluating donor-hydrogen-acceptor triplets involving triclosan and monomer oxygen atoms. Each frame was therefore assigned to an interaction-pattern family defined by the set of contact and hydrogen-bond events present in that frame. For each interaction pattern, the number of assigned frames and its fraction of the total trajectory were calculated. Representative structures were selected by constructing a triclosan-centered local feature vector for each frame and identifying the frame nearest the centroid of all frames belonging to the same interaction pattern. The molecular dynamics simulations were visually inspected using ChimeraX [21,22].

3. Results and Discussion

The IR spectra of the resultant MIP indicate the presence of all the characteristic peaks of an acrylate-based polymer network. The broad peak in the range of 3400 cm⁻¹ indicates the presence of the O-H group, which is due to the formation of carboxylic groups; they play a vital role in forming bonds with the template. The peaks located in the ranges of 2950-2850 cm⁻¹ are attributed to the C-H stretching vibration in the polymer backbone. The peak observed at 1720-1730 cm⁻¹ indicates the presence of the C=O group (Figure 1) [23].

The adsorption capacity of the synthesized MIPs and their non-imprinted analogs towards triclosan were determined using batch rebinding experiments. It is obvious from the results obtained that MIPs possess greater adsorption capacity than NIPs. This proves that the formation of selective cavities for triclosan molecules took place. The selectivity of MIPs towards triclosan may be explained by the formation of complementary cavities with respect to triclosan molecules which result from the polymerization process. In this regard, NIPs show significantly lower binding capacity due to absence of specific cavities as there are no cavities for specific interaction with target molecules. Thus, nonspecific adsorption due to van der Waals forces occur in this case. On the molecular level, the binding process depends on the formation of complexes between triclosan and methacrylic acid due to bonding between triclosan and carboxyl groups. Besides, the addition of 2-hydroxyethyl methacrylate improves the wettability of polymers, thus increasing the adsorption capacity of the resulting material (Figure 2) [24].

Optimized MIP surface morphologies were studied using SEM analysis. Microphotograph in Figure 3 depicts a heterogeneous, rugged, and very aggregated morphology of optimized MIP, which is typical for bulk MIP systems. The particles form nodular domains that are connected into porous clusters. Using SEM magnification scale, it was determined that the size of primary particles was within the range of about 0.4-1.2 µm, and an average primary particle size was ~0.8 µm, whereas larger agglomerates had sizes within the range of 5-25 µm. The existence of submicron-sized primary particles means that the MIP polymerization network was comprised of domains rather than large compact particles, which is common for conventional bulk-type MIP systems. Based on the observation of this surface morphology, it can be stated that this optimized MIP exhibits a porous surface structure, which probably arises from a high content of porogens (ACN). The optimized formulation, corresponding to MAA ≈ 0.50 g, HEMA ≈ 0.35 g, and ACN ≈ 14–16 mL, exhibited a predicted imprinting factor (IF) of approximately 2.8–2.9, consistent with the response surface analysis. The relatively high IF can be directly correlated with the observed morphology. The submicron particle size and porous structure reduce diffusion limitations and allow triclosan molecules to access imprinted cavities more efficiently. The improved accessibility in the present system enables effective utilization of binding sites, resulting in enhanced selectivity and higher IF values. At the same time, the presence of well-defined cavities formed through imprinting ensures that nonspecific adsorption remains limited, maintaining selectivity [25].

The release behavior of triclosan from the optimized molecularly imprinted polymer (MIP) was evaluated over a 48 h period, and the fractional release (M_t/M_∞) is presented in Figure 4. The experimental data show a characteristic two-stage release profile, consisting of an initial moderate release followed by a slower, sustained release phase. In the early stage (0–8 h), approximately 45–50% of triclosan was released, which can be attributed to the diffusion of weakly bound or surface-accessible molecules. This is followed by a slower release phase, reaching approximately 75–80% release at 48 h, indicating the presence of bound triclosan within imprinted cavities. Moderate initial release followed by sustained diffusion-controlled release suggests that the system effectively retains triclosan within the polymer network, consistent with template–monomer interactions.

To elucidate the release mechanism, the experimental data were fitted to four commonly used kinetic models: zero-order, first-order, Higuchi, and Korsmeyer–Peppas models. The fitting results and corresponding correlation coefficients ($R^2$) are summarized below:

Table 1. Release Model.

Model	Equation	R2	Interpretation
Korsmeyer-Peppas	Q= K.tn	0.9754	Best fit; describes complex diffusion.
Higuchi	Q= K.t0.5	0.925	High fit; indicates matrix-based diffusion.
First-order	Ln(1-Q) = -K.t	0.8441	Moderate fit; release depends on concentration.
Zero-order	Q=K.t	0.2864	Poor fit; release is not at a constant rate.

Table 1. Release Model.

Model	Equation	R2	Interpretation
Korsmeyer-Peppas	Q= K.tn	0.9754	Best fit; describes complex diffusion.
Higuchi	Q= K.t0.5	0.925	High fit; indicates matrix-based diffusion.
First-order	Ln(1-Q) = -K.t	0.8441	Moderate fit; release depends on concentration.
Zero-order	Q=K.t	0.2864	Poor fit; release is not at a constant rate.

Korsmeyer-Peppas model (R2=0.9754) showed the best correlation with experimental results, which means that it describes the release process most accurately among other models. This model is used when release processes are not ideal and several mechanisms work together. High correlation coefficient of Korsmeyer-Peppas model proves that a combination mechanism controls the release from the studied matrix. This mechanism consists of the following steps: Desorption from the imprinted sites; interaction of triclosan molecules with functional groups (MAA and HEMA) through interactions; breaking of these interactions takes time, hence, slowing down the second stage of release. Diffusion within the polymer matrix: porous nature of the polymer due to presence of porogen (ACN) causes gradual diffusion of triclosan molecules. This can be proved by a good fitting of Higuchi model (R2=0.9250), which corresponds to diffusion-controlled release from the matrix. Also, a relatively good fit of the first order model (R2=0.8441) shows that concentration gradient effects take place during the process. Zero-order model (R2=0.2864) was not correlated with data well, which means that a constant rate release mechanism did not occur in the system [26].

Finally, we perform molecular dynamics simulations to understand the most significant interaction patterns between MAA, HEMA and triclosan. Molecular dynamics simulations can be used to investigate MIP behavior by modeling the dynamic template-monomer-solvent interactions that govern pre-polymerization complex formation, binding-site development, and selective rebinding [27,28,29,30]. Analysis of the trajectory showed that triclosan interacted more frequently with MAA than with HEMA. It should be noted that the simulated mixture contained seven MAA molecules per triclosan, so only a subset of MAA molecules could interact directly with the template at any given time, consistent with a saturable local binding environment around the available triclosan interaction sites. The most populated single state was with no specific interactions between triclosan and the monomers (36% of frames), indicating that the local environment remained dynamic and that persistent binding was not maintained throughout the whole simulation. The most common MAA pattern was simultaneous MAA contact with both the phenolic hydroxyl region and the diaryl ether oxygen region of triclosan (15.3% of frames; Figure 5a). HEMA also interacted substantially with triclosan, with the most common HEMA pattern being simultaneous HEMA contact with both the phenolic hydroxyl region and the diaryl ether oxygen region (10.3% of frames; Figure 5a). Overall, MAA-only patterns accounted for 31% of frames, compared with 15.7% for HEMA-only patterns, while mixed MAA+HEMA environments accounted for 17.3% of frames. These results indicate that MAA remains the dominant recognition partner for triclosan in this pre-polymerization mixture, but that recognition is governed mainly by dynamic contact-based association around the phenolic and diaryl ether regions of triclosan rather than by persistent specific bonding.

4. Conclusion

In this present work, we employed DOE and MD simulations for the preparation and improvement of triclosan MIPs. Based on the Box-Behnken design method, the effects of MAA, HEMA, and ACN on imprinting were efficiently assessed, and an optimal formula was identified which resulted in a relatively high imprinting factor of 2.8–3.0. As seen from the results, the optimal MIPs had a higher adsorption capacity and selectivity than NIPs because of the creation of recognition sites. The porous morphology could improve the mass transfer rate and access to the binding sites. The release behavior illustrated the controlled diffusion process. MD simulations offered further insights into the imprinting process at the molecular level, demonstrating that MAA was essential in triclosan recognition via its interaction with the template molecules. The correspondence between the optimal formula obtained via DOE and MD indicated the advantages of using a combination of DOE and MD. In summary, this study contributes towards a comprehensive approach to develop MIPs with adjustable absorption and desorption characteristics. The current methodology may further be applied in the development of MIPs for other target molecules.