Evaluation of polymer–preservative interactions for preservation efficacy: molecular dynamics simulation and QSAR approaches†

Qisong Xu^a, Pui Shan Chow^a, Erte Xi^b, Randy Marsh^b, Shikar Gupta^c and Krishna M. Gupta*^a
aInstitute of Sustainability for Chemicals, Energy and Environment (ISCE2), Agency for Science, Technology and Research (A*STAR), 1 Pesek Road, Jurong Island, Singapore 627833, Republic of Singapore. E-mail: krishna_mohan_gupta@isce2.a-star.edu.sg
^bProctor & Gamble, Winton Hill Business Center, 6280 Center Hill Ave., Cincinnati, OH 45224, USA
^cProcter & Gamble International Operations SA SG Branch, Singapore 138547, Singapore

First published on 14th August 2024

---

---

1. Introduction

Cosmetic industries reached a global market value of US $379.7 billion in 2022 and this is expected to grow at a rate of 5.50% per year during 2022–2028.¹ However, as cosmetic product formulations contain mostly water and organic compounds, such formulations inevitably increase the risks of microbial contamination in aqueous media.² In addition, the toxins arising from microbial contamination may also cause chemical, biological and physical deterioration of cosmetic products, which may pose detrimental health risks to consumers.³ Therefore, the addition of appropriate preservatives not only ensures long shelf-life and the safety of cosmetic products but also increases product quality.⁴ Indeed, preservatives are frequently employed in cosmetic products,⁵ and in various pharmaceutical and food industries.^6,7 The development of an effective preservative system should also adhere to good manufacturing practices (GMPs) during and after manufacturing through primary and secondary preservation strategies.⁵ Therefore, the development of new preservative systems with high efficacy continues to attract considerable interest for various consumer products.

To date, the advancement of high efficacy preservative systems in consumer products still poses significant challenges and relies predominantly on experimental approaches.⁵ Typically, the performance of these systems depends on the characteristics of polymer fibers, preservatives, and their intricate interactions in water. For example, Dubrovski and Brezocnik demonstrated that the increase in volume porosity of hydrophilic viscose and hydrophobic polyester fibers resulted in higher vertical wicking rate for all fiber contents used.⁸ Using spun lacing, Jain et al. successfully fabricated various non-woven structures of fabric with significant differences in pore diameter and thickness by varying the waterjet pressure, delivery speed, web mass and composition.⁹ To guide the design of anti-fouling fibers, Kim and coworkers comprehensively investigated the bacterial adhesion on polystyrene and polylactic acid fibers, which was highly affected by various fiber characteristics, including fiber spatial distribution, wettability, pore volume and pore size.¹⁰ Salama et al. also highlighted the importance of the nature and combinations of polyester and viscose on the interactions between wet wipes and preservatives, thereby influencing the efficacy of the preservative system for wet wipes.¹¹ Different classes of preservatives have also been studied and developed. For example, natural preservatives from plant extracts and essential oils were screened for the formulation of cosmetic products to produce safe preservative systems with high efficacy.¹² In addition, multifunctional ingredients, such as chelating agents, were also explored in the formulation of preservative systems to disrupt bacterial growth and activity.^13,14 Furthermore, the impact of a series of phenolic preservatives in consumer products on human health and the environment was systematically evaluated using hazard assessment and computational modelling.¹⁵

To quantify the efficacy of the preservative system, an antimicrobial test was initially conducted by assessing the performance of antimicrobials to impede bacterial growth.¹⁶ Subsequently, the challenge test was internationally adopted to determine the stability and efficacy of the preservative system, where the products were inoculated with specified amounts of microorganisms at room temperature over 28 days.¹⁷ The decimal time reduction (D-value) method is a rapid technique to evaluate the preservative efficacy.^18,19 In particular, the D-value represents the inactivation rate of microorganisms (e.g., bacteria, yeast, and mold) after the inoculation of prepared preservative systems. Furthermore, the capacity test was evaluated to gauge the effectiveness of the preservative concentration in the preservative systems, and to measure the resulting bacterial activity in a range of formulations (e.g., creams, suspensions, and solutions).²⁰

Apart from the aforementioned experiments, several computational studies have also been widely employed as promising tools to investigate the molecular interactions between polymers and small molecules in various formulation systems.^21–23 For instance, in the context of the current topic using molecular dynamics (MD) simulations, Wang and coworkers modelled 13 different small molecules in amorphous polyethylene terephthalate (PET) to estimate their diffusion coefficients, which matched well with experimental measurements.²⁴ They further simulated the diffusion of three additives in PET to reveal key properties critical to a diffusion mechanism, namely, the interaction energy between additives and PET, fractional free volume, molecular shape and size, and self-diffusion of the polymer.²⁵ For polypropylene (PP), the effect of the types of amorphous homopolymers and copolymers on the diffusion behaviour of limonene in PP was investigated; larger free volume clusters in random copolymers enhance limonene transport.²⁶ In the temperature range of 293–343 K, five additives (BHT, UV-P, UV-327, UV-329 and UV-531) in isotactic PP were also systematically simulated to highlight key properties of materials contributing to molecular diffusion.²⁷ Furthermore, a two-phase MD simulation model was proposed to evaluate the diffusion of five additives from PP to food simulants (modelled as an ethanol/water mixture and pure isooctane, respectively).²⁸ In addition, Mazeau and Vergelati employed the Monte Carlo technique to model the adsorption configuration and energy of benzophenone on amorphous and crystalline cellulose.²⁹ By complementing molecular modelling and experiments, they further determined and compared the adsorption enthalpies of various aromatic molecules on cellulose.³⁰ Using density functional theory calculations, Todde et al. computed the adsorption interaction energies and configurations of five small molecules in various cellulose derivatives.³¹ However, to the best of our knowledge, previous molecular studies are not only largely limited to a small number of molecules in specified polymers, but also to simulations under dry conditions. As a result, these studies may not be able to capture the critical role of water in mediating interactions in realistic scenarios.

In this study, we report a comprehensive evaluation of polymer–preservative interactions in water for preservation efficacy using MD simulation and quantitative structure–activity relationship (QSAR) approaches. While past QSAR models of preservatives have been reported for their antibacterial activity and risk assessment,^32–34 the computational methods in this work are demonstrated for the molecular investigation of three polymers of distinct nature (PET, PP, and cellulose) and fifteen preservatives from three chemical classes (aliphatic, aromatic and organic acids).

2. Materials and methods

2.1. Force field validation of polymers

The molecular structures of individual polymer chains (PET, PP and cellulose) were first generated using a web-based simulation builder (CHARMM-GUI).³⁵ As shown in Fig. 1, for a given polymer, 18 chains in total were used to assemble the 3D atomistic model of an amorphous polymer. In each chain, different numbers of monomers were also defined to ensure that the total mass of a polymer model was approximately 87,000 g (see Table S1†). An MD simulation protocol was employed to assemble polymer chains to form a polymer model.³⁶ The simulation protocol consists of a 7-step compression and relaxation scheme (Table S2†) to obtain an equilibrated structure of the polymer. Three repeat simulation runs were also performed at 298 K for each polymer. The simulated polymer densities obtained from the simulation protocol were subsequently compared to the literature data. It is worth noting that the simulation protocol in this work is applicable to the construction of pure polymers, while simulation protocols have also been reported for polymer composites containing nanomaterials, including carbon and boron nitride nanotubes.^37,38

2.2. Adsorption model

To study the adsorption phenomena of a preservative at a polymer surface, a 2D polymer model with periodicities along the x- and y-axis was first constructed using the same 7-step simulation protocol described in section 2.1.³⁶ As opposed to the 3D models, the construction of 2D polymers with smooth surfaces was facilitated by two additional atomic walls in the z-direction during the simulation protocol. Fig. S1† shows the final MD snapshots of the 2D polymer models for adsorption.

Table 1 illustrates the list of fifteen preservatives investigated in this study. Broadly, these preservatives belong to three different chemical classes, namely, aliphatic, aromatic and organic acids. An adsorption model was constructed, as shown in Fig. 2, to simulate the dynamic adsorption behaviour of individual preservatives. First, an equilibrated aqueous-preservative model was constructed for all preservatives based on the experimental conditions (a diluted solution with each preservative conditioned at pH = 3.8). A total of 10 molecules (preservative and its associated ions) and an appropriate number of sodium ions were inserted into the water box. Table S3† lists the pK_a of preservatives, and the number of neutral and ionic forms of preservatives at pH 3.8 (calculated using the Henderson–Hasselbalch equation). Next, the water box was positioned at the two interfaces of a dry polymer. Finally, a 50-ns NPT run was performed to monitor the dynamic adsorption of the preservatives to the polymer. For each polymer, three repeated simulation runs were performed for reproducibility. During the NPT run, the trajectories of the preservative molecules and the associated ions were monitored over the 50 ns period. The trajectories of the preservatives along the direction normal to the polymer surface were tracked to obtain molecular insights into the interaction behaviour at the interface. To quantitatively analyse the adsorption of the preservatives, two interaction energies (IEs) were computed from the NPT, namely, (1) IEs between the preservatives and the polymers at the interface, and (2) IEs between the preservatives and water. For preservatives that were strongly adsorbed to polymers, the last 5 ns of the NPT run was used to calculate the IEs. On the other hand, for preservatives that were weakly adsorbed to polymers, the last 20 ns of the NPT run was used to compute the IEs; the longer duration for analysis provides more statistically significant IEs for the weakly adsorbed preservatives. Furthermore, the average IEs were computed from three repeated simulation runs.

---

2.3. Diffusion model

To calculate the interior diffusion coefficient of the preservatives in polymers, a selected number of preservatives were first artificially dragged from the bulk water phase to the interior of the polymers. As illustrated in Fig. 3, once the preservatives are present in the polymer's interior, a set of MD simulations (10 ns NPT followed by 20 ns NVT) was performed. The 10 ns NPT ensured that the preservatives were well-equilibrated in the polymers. Subsequently, the 20 ns NVT was performed to allow sufficiently long simulation times to obtain the interior diffusion coefficients of the preservatives and their IEs with the polymers. For each combination of preservatives and polymers, five repeated NVT runs were also performed for reproducibility. In each NVT run, the initial 10 ns was employed to obtain a more precise estimate of diffusion coefficient with less statistical error. On the other hand, for the calculation of IEs of the preservatives in the polymer interior, the last 5 ns was employed. Finally, the average internal diffusion coefficients and IEs for each combination of preservatives and polymers were computed from the five repeated NVT runs.

All the simulations were performed using GROMACS version 2018.³⁹ The CHARMM force field parameters were adopted to model all polymers, molecules and ions.⁴⁰ Energy minimization was conducted using the steepest descent method before performing the MD simulations. The initial velocities for MD simulations were generated using the Maxwell–Boltzmann distribution. For pure polymer simulations, the temperature was controlled by a Nosé–Hoover thermostat with a relaxation time of 0.1 ps. The pressure was controlled by a Berendsen barostat via semi-isotropic coupling. For adsorption and diffusion simulations, the temperature was controlled by the Nosé–Hoover thermostat whereas the pressure was controlled by a Parrinello–Rahman barostat via semi-isotropic coupling. The coulombic interactions were calculated using the particle-mesh Ewald summation method with a grid spacing of 1.2 Å. The LJ interactions were calculated using a cut-off of 12 Å. The equations of motion were integrated by the leap-frog algorithm with a time step of 1 fs. Overall, from the MD simulation of adsorption and diffusion, four quantitative target properties were calculated: (T1) IEs between the preservatives and the polymers at the interface; (T2) IEs between the preservatives and water; (T3) interior IEs between the preservatives and the polymers; and (T4) diffusion coefficients of the preservatives in polymers.

2.4. QSAR modelling

3. Results and discussion

3.1. Force field validation of polymers

Using PET polymer as an example, the impact of the number of polymer chains and the number of monomeric units in a chain on polymer packing was first evaluated. The variation in both numbers of polymer chains and monomeric units in a chain showed an insignificant effect on density of the PET polymer (see Table S6†). Therefore, for PP and cellulose, the same number of chains (18) was adopted while varying the number of monomeric units appropriately to construct polymer models of similar mass. Indeed, as tabulated in Table 2, the simulated densities for all three polymers agreed well with those reported in previous simulations and experiments, thereby validating the force fields for these polymers. Furthermore, the variation in densities is attributed to the difference in the molecular structure of these polymers. Particularly, the low density of PP reflects weak van der Waals interactions among its long hydrophobic aliphatic chains.⁴³ As an aromatic ester, PET possesses both polar and aromatic functional groups that enhance stronger intermolecular interactions among chains, which increase polymer density.⁴⁴ Owing to the large number of hydroxyl groups present in glucose units, the polymer chains in cellulose can form strong intermolecular hydrogen bonds, which further reduce the interchain distance and increases cellulose density.⁴⁵

---

3.2. Adsorption behaviour of the preservatives

From the MD simulations, distinctive adsorption behaviours were observed for various preservatives at the interface of each polymer. Notably, the unique adsorption of preservatives is a complex interplay of interactions among different types of polymers, classes of preservatives and water.

For PET, the preservatives were either adsorbed at the PET interface or existed in a dynamic state between the aqueous phase and the interface depending on their orientation and interaction strength with respect to the interface. As described in section 2.2, the aqueous preservative mixture contains both neutral and ionic forms of the preservatives; therefore, the adsorption behaviour of both forms was evaluated. Fig. 4 shows the representative snapshots of an adsorbed (benzoic acid) and a non-adsorbed (citric acid) neutral preservative on the PET surface as well as the trajectories of these preservatives during simulations. Using benzoic acid as the adsorbed example, most classes of preservatives were adsorbed to the PET surface as observed in Fig. 4a (green coloured molecules). As is evident from the molecular trajectories in Fig. 4b, the preservatives were indeed diffusing towards the PET surface. Due to the presence of both hydrophobic phenyl groups and hydrophilic ester groups in PET, these structural fragments facilitate favourable molecular interactions with most preservatives, thereby enabling the adsorption of preservatives to PET. However, there were a few preservatives that did not adsorb to the PET surface, as observed in Fig. 4c (green coloured molecules). Notably, these preservatives included some organic acids (e.g., citric acid, succinic acid and levulinic acid), which have several hydrophilic groups (e.g., carboxyl and hydroxyl) that interact more favourably with the bulk water phase, as illustrated in their molecular trajectories during simulation, as shown in Fig. 4d. Ionic forms of all classes of preservatives were in the dynamic state between the aqueous phase and the interface. Owing to the neutral PET surface, ionic forms of preservatives interacted weakly with the polymer. In addition, during the adsorption of all preservatives on PET, water molecules were observed to diffuse into the PET interior and interact with the ester groups of PET. To further demonstrate microscopic insights into the benzoic acid adsorption behaviour on PET, simulation snapshots near the interface at 2 ns and 4 ns were extracted, as shown in Fig. 5. The molecular configurations of the benzoic acid at the PET surface indicate that the phenyl ring of the benzoic acid (green coloured molecules) is closely oriented towards the phenyl rings of the PET. Such favourable π-stacking is a ubiquitous phenomenon in numerous chemical and biochemical systems that arises from the noncovalent interactions between aromatic rings.^51,52 Furthermore, our computational findings were consistent with observations from previous MD simulations, where π-stacking between aromatic carbon rings of polymers and the surfaces of nanomaterials (carbon and boron nitride nanotubes) was also significant.^53,54

For PP, all classes of preservatives were observed to adsorb to the polymer surface. Fig. 6a illustrates a representative snapshot of an adsorbed preservative (benzoic acid) on the PP surface. Furthermore, as visualized in Fig. 6b, the trajectories of benzoic acid demonstrate a favourable affinity and dynamic behaviour towards PP. This observation was attributed to the hydrophobic nature of PP that interacts weakly with water molecules, thereby exposing further PP surface for interactions with preservatives. Indeed, all studied classes of preservatives have nonpolar skeletal carbons, which have high affinity with hydrophobic PP. Interestingly, for succinic acid and citric acid, both preservatives were observed to adsorb to PP but not to PET. As illustrated in Fig. 7, there were favourable molecular configurations of these organic acids that enhanced their adsorption to the PP surface. Specifically, the carbon backbones of these organic acids were preferentially attached to the hydrophobic PP while their carboxyl groups were oriented towards the hydrophilic water phase. During the adsorption of all preservatives to PP, no water molecules were observed to diffuse into the hydrophobic interior of PP. Similar to PET, the ionic forms of all preservatives were in the dynamic state between the aqueous phase and the interface.

For cellulose, none of the preservatives were adsorbed to the polymer surface. Using benzoic acid as an example, Fig. 8a demonstrates the representative snapshot of the non-adsorbing preservative on the cellulose surface. As illustrated in the trajectories of Fig. 8b, benzoic acid was in the dynamic state between the cellulose surface and the water phase; the molecules did not adsorb to cellulose but remained in the water during the simulation. This is because water molecules have stronger interactions with hydrophilic cellulose (−9.50 kJ mol⁻¹) than with hydrophobic PP (−2.19 kJ mol⁻¹). Indeed, unlike PET and PP, cellulose has a higher number of hydroxyl groups in its structure, resulting in a highly hydrophilic surface that interacts favourably with water.⁵⁵ These water molecules tend to adsorb to cellulose, thereby reducing the adsorption of the preservatives to cellulose. Furthermore, during the adsorption processes of all preservatives for cellulose, some water molecules were observed to diffuse into cellulose and interact with its hydroxyl groups.

3.3. Preservative–polymer and preservative–water interaction energies

To quantitatively describe the difference in the adsorption behaviour of preservatives at the polymer surface, two types of IEs at the interface were computed: (1) IEs between the preservatives and the polymers, and (2) IEs between the preservatives and water. It is also worthwhile noting that the IEs were computed based on a single preservative solute, which provides a straightforward comparison of molecular interactions among the different types of polymers and preservatives.

In Fig. 9, the IEs between the preservatives and the polymers at the interface were compared for various polymers and preservatives. Generally, for all preservative classes, the magnitude of these IE values (negative means favourable interactions) indicated more favourable interactions between the preservatives and the polymers, starting from cellulose to PP, followed by PET. As the IEs between the preservatives and the polymers were the lowest for cellulose, the preservatives were less likely to be adsorbed onto cellulose. Furthermore, the most favourable IEs were observed for PET as this polymer possesses both hydrophobic phenyl and hydrophilic ester groups that interact with both hydrophobic fragments and hydrophilic hydroxyl groups of preservatives. In terms of the preservatives, for aliphatic and aromatic classes (Fig. 9a and b), the IEs generally increased with the molecular weight of the preservatives; larger molecules have higher surface areas to enhance molecular interactions with polymers. In Fig. 9c, while all simple organic acids (levulinic acid, succinic acid and citric acid) have the lowest IEs with cellulose, it is worthwhile noting that levulinic acid and succinic acid interact more favourably with PP, as compared to citric acid. As previously illustrated in Fig. 7, unlike citric acid having larger numbers of hydroxyl groups, the carbon backbones of succinic acid were more preferentially adsorbed to the PP surface to enhance the interactions between succinic acid and PP. Furthermore, compared to the simple organic acids, aromatic organic acids (benzoic acid and p-anisic acid) have the highest IEs with PET as these preservatives contain both phenyl and carboxylic groups, which interact more favourably with PET.

Fig. 10 shows the IEs between the preservatives and water at the interface for various polymers and preservatives. Generally, the magnitude of these numerically negative IE values indicated more favourable interactions between the preservatives and water, starting from PET to PP, followed by cellulose. Such an observation is opposite to the trend of IEs between the preservatives and the polymers as shown in Fig. 9. Specifically, if a preservative interacts strongly with a polymer, the corresponding IEs between the preservatives and water would be lower, thereby leading to adsorption to the polymer. Furthermore, owing to the strong cellulose–water affinity, as discussed previously, the preservatives were less likely adsorbed on cellulose but more readily dissolved in water. Therefore, the IEs between the preservatives and water for cellulose were the highest for all preservatives. In terms of the preservatives, for aliphatic and aromatic classes (Fig. 10a and b), the IEs between the preservatives and water generally increased with the molecular weight of the preservatives. This could be attributed to the presence of more hydroxyl groups and larger surface areas in the preservatives, which have high affinity for water. Furthermore, aliphatic preservatives tend to have higher IEs with water than aromatics; the phenyl rings in aromatic preservatives interact weakly with water. Similar observations were also evident for organic acids, as shown in Fig. 10c; the IEs between simple organic acids and water increased from levulinic acid to succinic acid and citric acid, which correlated with the higher number of hydroxyl groups present in the preservatives.

3.4. Diffusion behaviour of preservatives in polymers

As illustrated in Fig. 11a, the diffusion coefficients of preservatives in polymers generally increased in the order of cellulose < PET < PP. As compared in Table 2, polymer density follows the opposite trend, i.e. increasing in the order of PP < PET < cellulose. Due to the decrease in polymer density, polymer porosity is increased, which leads to an increase in the diffusion coefficients of the preservatives. Specifically, PP has the lowest polymer density of 0.85 g cm⁻³ as its hydrophobic polymer chains were only associated through weak van der Waals interactions. This resulted in a more porous PP that facilitates high diffusion for all preservatives. On the other hand, cellulose has hydrophilic hydroxyl groups that enable strong hydrogen-bonding among its polymer chains, thereby reducing the voids in cellulose (density of 1.43 g cm⁻³) and hindering the diffusion of preservatives.

Furthermore, the variations in the diffusion coefficients of the preservatives can also be attributed to the difference in the interior IEs between the preservatives and the polymers. As compared in Fig. 11b, these IE values generally suggested more favourable interior interactions between the preservatives and the polymers, starting from PP to PET, followed by cellulose. The opposite trend was observed for the diffusion coefficients of preservatives in polymers. Specifically, the stronger a preservative interacts with a polymer, the lower are its mobility and diffusion coefficient values within that polymer. For instance, cellulose not only has high polymer density that reduces voids for diffusion but also has a large number of hydroxyl groups that interact favourably with the polar functional groups of preservatives, thereby reducing their diffusion coefficients.

In addition, the molecular weight, size, and shape of preservatives can also influence their diffusion coefficients in polymers. In Fig. 11c, the diffusion coefficients decrease with increasing molecular weights of the preservatives. For example, the diffusion coefficients of the preservatives in PET decrease in the order from benzyl alcohol (MW 108.14 g mol⁻¹) to phenethyl alcohol (MW 122.17 g mol⁻¹) and phenyl propanol (MW 136.19 g mol⁻¹) even though the respective interior IEs between these three aromatic preservatives and PET were similar. Furthermore, it has been observed that a preservative with a smaller size leads to higher diffusion. For example, the molecular size of succinic acid is smaller than that of levulinic acid (Fig. S2a†); hence, the diffusion coefficient of succinic acid (D_{succinic acid} = 3.44 × 10⁻⁷cm s⁻²) in PET is higher than that of levulinic acid (D_{levulinic acid} = 1.68 × 10⁻⁷ cm s⁻²). Lastly, it was observed that the linear shape of preservatives contributes to higher diffusion coefficients compared to those of bulky-shaped preservatives. For example (Fig. S2b†), linear-shaped octanediol shows a higher diffusion coefficient than bulky-shaped phenoxyethanol in PET, despite octanediol (MW of 146.23) having a higher molecular weight than phenoxyethanol (MW of 138.16). Overall, the diffusion coefficients of preservatives in different polymers depend on various factors, including (1) the nature of the polymer, (2) interior IEs between preservatives and polymers, and (3) the structural properties of preservatives.

3.5. Development of QSAR models

As shown in Fig. S3–S6,† the sensitivity analysis of each target property in PET, PP, and cellulose polymers was comprehensively performed. Generally, as the number of molecular descriptors increases, the R² and RMSE values of the explored QSAR models increase and decrease, respectively. Furthermore, as the total number of descriptors increases beyond two, R² reaches at least 0.9, together with a lower RMSE value. Therefore, two unique descriptors were deemed sufficient to quantitatively predict four target properties in all polymers. These two unique descriptors identified for each target property in polymers are subsequently tabulated in Table 3. Apart from distinctive descriptors that were unique to a type of polymer, some common descriptors were also found to be critical to the prediction of a target property. Most of the descriptors that appear in the developed QSAR models belong to spatial descriptors. For example, when predicting the (T1) IEs between the preservatives and the polymers at the interface, “Jurs_PPSA_3” spatial descriptors (descriptors based on partial charges mapped on a surface area)⁵⁶ were found to be important to both PET and cellulose while “Molecular_Volume”, another 3D spatial descriptor that defines the molecular volume inside the contact surface, was found to be important to both PP and cellulose, which suggests that size, shape and partial charge on the preservatives are critical for the interaction between preservatives and polymers. To estimate the (T2) IEs between preservatives and water, “Molecular_PolarSurfaceArea” and “Jurs_RPCS” were found to be important to all three polymers, which suggests that the polarity on the preservative molecules is a key factor for the interaction between the preservatives and water. For the (T3) interior IEs between the preservatives and the polymers, “Jurs_PPSA_3” was found to be important to all three polymers. To compute the (T4) diffusion coefficients of the preservatives in the polymers, other spatial descriptors of shadow indices (geometric descriptors that help to characterize the shape of the molecules), such as “Shadow_XYfrac”, “Shadow_XZfrac” and “Shadow_Ylength” were identified as key descriptors. The presence of shadow index descriptors in the QSAR models clearly suggests that diffusion coefficients of preservatives in polymers depend on the shape of the preservatives.⁵⁷

---

Finally, the parity plots for all four target properties in various polymers are compared in Fig. 12. Using the respective two unique descriptors given in Table 3, the predicted IEs and diffusion coefficients of preservatives from QSAR models agreed well with the actual values calculated from MD simulations. All QSAR models provide high prediction power (R² > 0.9) and low prediction error (low RMSE) for all IEs and diffusion coefficients.

4. Conclusions

By complementing MD simulation and QSAR modelling, we have fundamentally elucidated the key microscopic insights and material properties that influence polymer–preservative interactions for preservation efficacy. In particular, the molecular interactions between three polymers (PET, PP, and cellulose) and fifteen preservatives from three chemical classes (aliphatic, aromatic, and organic acids) were investigated. The simulated densities of the three polymers were found to agree well with the literature, which validated the force field of the polymers. During adsorption in PET, most preservatives (except organic acids such as citric acid, succinic acid and levulinic acid) were observed to interact favourably with the phenyl groups and the ester groups in PET. For PP, all classes of preservatives were adsorbed to the polymer surface. Interestingly, the carbon backbones of organic acids (citric acid and succinic acid) were adsorbed to hydrophobic PP while its carboxyl groups interacted with water at the same time. On the other hand, no preservatives were adsorbed to the cellulose surface as water molecules interacted more with the cellulose surface and displaced preservatives to the water phase. To quantitatively verify the adsorption behaviour, the IEs between preservatives and polymers were calculated; for all preservatives, IEs were found to generally increase as cellulose < PP < PET. In contrast, a corresponding opposite trend in IEs between preservatives and water is observed in the polymers (PET < PP < cellulose). Furthermore, the diffusion coefficients of preservatives were found to decrease with increasing polymer density from PP to PET and cellulose. The favourable interaction between the polar groups of preservatives and the hydroxyl groups of cellulose can also impede the diffusion of preservatives. In addition, linear-shaped preservatives with low molecular weight and size tend to exhibit higher diffusion coefficients in polymers. From the sensitivity analysis and QSAR modelling, key molecular descriptors of the preservatives were identified to predict their IEs and diffusion coefficients in the polymers. Overall, a holistic computational approach is reported herein from deriving atomic-resolution insights to identifying key materials properties and estimating the molecular interactions between polymers and preservatives in water. This complementary method can potentially provide a systematic approach for rationally guiding the selection of polymers and preservatives, and the design of high efficacy preservative systems for pharmaceutical, food and cosmetic products.

Nevertheless, we should also note the assumptions and limitations associated with the computational approaches in this work: (1) the current work has focused on specific industrially relevant polymers (PET, PP and cellulose) and preservatives, which may limit the generalizability of the findings to other polymers and preservatives. In the future, simulation studies for other broader classes of polymers and preservatives should also be required for more generalizable insights. (2) Additional experiments could be considered to validate the computational findings. The adsorption and diffusion simulations of polymers and preservatives are typically simulated on shorter time and length scales than experiments; therefore, interaction energies and diffusion coefficients from computations might not be straightforwardly compared with those obtained from experiments. Nevertheless, the simulation findings could still provide critical microscopic trends of the intricate types of molecular interactions among different classes of polymers and preservatives, and the diffusional behaviour of preservatives in polymers, which may be intractable in experiments. Despite these assumptions and limitations, this work demonstrates the integrated strategy of MD simulation and QSAR modelling to quantitatively and rapidly predict adsorption and diffusion of preservatives in polymers, and fundamentally unravel molecular insights of polymer–preservative interactions in water, which may streamline and accelerate the development of high efficacy preservative systems.

Data availability

The data supporting this article have been included as part of the ESI.†

Conflicts of interest

The authors declare no competing financial interest.

Acknowledgements

This research was supported jointly by the Agency for Science, Technology and Research (A*STAR) and Proctor & Gamble (P&G) under Enterprise Gap Funding (Project Number: I22D4AG003). This research was also supported by the A*STAR under its Career Development Fund (Project Number: C210812028). This work was supported by the A*STAR Computational Resource Centre through the use of its high-performance computing facilities. The computational work for this article was partially performed on resources of the National Supercomputing Centre (NSCC), Singapore (https://www.nscc.sg). Q. Xu thanks the National Research Foundation Singapore for funding through the SGUnited Jobs Initiative (P20J3d1014).

References

1. . Global Cosmetics Market: Industry Trends, Share, Size, Growth, Opportunity and Forecast 2023–2028. URL: https://www.researchandmarkets.com/categories/cosmetics.
2. J. J. Kabara and D. S. Orth, Principles for product preservation, in Preservative-Free and Self-Preserving Cosmetics and Drugs: Principles and Practices, Marcel Dekker, New York, 1997, pp. 1–14 Search PubMed .
3. E. Neza and M. Centini, Microbiologically Contaminated and Over-Preserved Cosmetic Products According Rapex 2008–2014, Cosmetics, 2016, 3(1), 3 CrossRef .
4. J. C. Toler, Preservative stability and preservative systems, Int. J. Cosmet. Sci., 1985, 7(4), 157–164 CrossRef .
5. N. Halla, I. P. Fernandes, S. A. Heleno, P. Costa, Z. Boucherit-Otmani, K. Boucherit, A. E. Rodrigues, I. C. F. R. Ferreira and M. F. Barreiro, Cosmetics Preservation: A Review on Present Strategies, Molecules, 2018, 23(7), 1571 CrossRef .
6. P. Rathee, R. Sehrawat, P. Rathee, A. Khatkar, E. K. Akkol, S. Khatkar, N. Redhu, G. Türkcanoğlu and E. Sobarzo-Sánchez, Polyphenols: Natural Preservatives with Promising Applications in Food, Cosmetics and Pharma Industries; Problems and Toxicity Associated with Synthetic Preservatives; Impact of Misleading Advertisements; Recent Trends in Preservation and Legislation, Materials, 2023, 16(13), 4793 CrossRef PubMed .
7. C. Novais, A. K. Molina, R. M. V. Abreu, C. Santo-Buelga, I. C. F. R. Ferreira, C. Pereira and L. Barros, Natural Food Colorants and Preservatives: A Review, a Demand, and a Challenge, J. Agric. Food Chem., 2022, 70(9), 2789–2805 CrossRef PubMed .
8. P. D. Dubrovski and M. Brezocnik, Porosity and nonwoven fabric vertical wicking rate, Fibers Polym., 2016, 17(5), 801–808 CrossRef .
9. R. K. Jain, S. K. Sinha and A. Das, Structural investigation of spunlace nonwoven, Res. J. Text. Appar., 2018, 22(3), 158–179 CrossRef CAS .
10. T. Hemmatian, H. Lee and J. Kim, Bacteria Adhesion of Textiles Influenced by Wettability and Pore Characteristics of Fibrous Substrates, Polymers, 2021, 13(2), 223 CrossRef CAS .
11. P. Salama, A. Gliksberg, M. Cohen, I. Tzafrir and N. Ziklo, Why Are Wet Wipes So Difficult to Preserve? Understanding the Intrinsic Causes, Cosmetics, 2021, 8(3), 73 CrossRef CAS .
12. A. Kerdudo, P. Burger, F. Merck, A. Dingas, Y. Rolland, T. Michel and X. Fernandez, Development of a natural ingredient – Natural preservative: A case study, C. R. Chim., 2016, 19(9), 1077–1089 CrossRef CAS .
13. A. Varvaresou, S. Papageorgiou, E. Tsirivas, E. Protopapa, H. Kintziou, V. Kefala and C. Demetzos, Self-preserving cosmetics, Int. J. Cosmet. Sci., 2009, 31(3), 163–175 CrossRef .
14. W. Siegert, Boosting the antimicrobial efficiency of multifunctional additives by chelating agents, Int. J. Appl. Sci., 2014, 140, 1–6 Search PubMed .
15. H. L. Buckley, W. M. Hart-Cooper, J. H. Kim, D. M. Faulkner, L. W. Cheng, K. L. Chan, C. D. Vulpe, W. J. Orts, S. E. Amrose and M. J. Mulvihill, Design and Testing of Safer, More Effective Preservatives for Consumer Products, ACS Sustainable Chem. Eng., 2017, 5(5), 4320–4331 CrossRef .
16. A. M. Cundell, S. Chatellier, P. Schumann and R. Lilischkis, Equivalence of quality control strains of microorganisms used in the compendial microbiological tests: are national culture collection strains identical?, PDA J. Pharm. Sci. Technol., 2010, 64(2), 137–155 Search PubMed .
17. C. Booth, Antimicrobial effectiveness testing validation strategies, J. GXP Compliance, 2014, 18, 1–12 Search PubMed .
18. D. S. Orth, Linear regression method for rapid determination of cosmetic preservative efficacy, J. Soc. Cosmet. Chem., 1979, 30, 321–332 CAS .
19. D. S. Orth and D. C. Enigl, Preservative efficacy testing by a rapid screening method for estimation of D-values, J. Soc. Cosmet. Chem., 1993, 44(6), 329–336 Search PubMed .
20. S. V. Sutton and P. A. Geis, Antimicrobial preservative efficacy and microbial content testing, in Cosmetic microbiology, CRC Press, 2020, pp. 67–93 Search PubMed .
21. D. M. Walden, Y. Bundey, A. Jagarapu, V. Antontsev, K. Chakravarty and J. Varshney, Molecular Simulation and Statistical Learning Methods toward Predicting Drug–Polymer Amorphous Solid Dispersion Miscibility, Stability, and Formulation Design, Molecules, 2021, 26(1), 182 CrossRef CAS PubMed .
22. K. M. Gupta, X. Chin and P. Kanaujia, Molecular Interactions between APIs and Enteric Polymeric Excipients in Solid Dispersion: Insights from Molecular Simulations and Experiments, Pharmaceutics, 2023, 15(4), 1164 CrossRef CAS PubMed .
23. K. M. Gupta, S. Das, A. B. H. Wong and P. S. Chow, Formulation and Skin Permeation of Active-Loaded Lipid Nanoparticles: Evaluation and Screening by Synergizing Molecular Dynamics Simulations and Experiments, Langmuir, 2023, 39(1), 308–319 CrossRef CAS .
24. Z.-W. Wang, P.-L. Wang and C.-Y. Hu, Molecular dynamics simulation on diffusion of 13 kinds of small molecules in polyethylene terephthalate, Packag. Technol. Sci., 2010, 23(8), 457–469 CrossRef CAS .
25. B. Li, Z.-W. Wang, Q.-B. Lin and C.-Y. Hu, Molecular dynamics simulation of three plastic additives’ diffusion in polyethylene terephthalate, Food Addit. Contam.: Part A, 2017, 34(6), 1086–1099 CrossRef CAS .
26. Z.-W. Wang, P.-L. Wang and C.-Y. Hu, Investigation in Influence of Types of Polypropylene Material on Diffusion by Using Molecular Dynamics Simulation, Packag. Technol. Sci., 2012, 25(6), 329–339 CrossRef CAS .
27. Z.-W. Wang, B. Li, Q.-B. Lin and C.-Y. Hu, Molecular Dynamics Simulation on Diffusion of Five Kinds of Chemical Additives in Polypropylene, Packag. Technol. Sci., 2018, 31(5), 277–295 CrossRef CAS .
28. Z.-W. Wang, B. Li, Q.-B. Lin and C.-Y. Hu, Two-phase molecular dynamics model to simulate the migration of additives from polypropylene material to food, Int. J. Heat Mass Transfer, 2018, 122, 694–706 CrossRef CAS .
29. K. Mazeau and C. Vergelati, Atomistic Modeling of the Adsorption of Benzophenone onto Cellulosic Surfaces, Langmuir, 2002, 18(5), 1919–1927 CrossRef .
30. D. Da Silva Perez, R. Ruggiero, L. C. Morais, A. E. H. Machado and K. Mazeau, Theoretical and Experimental Studies on the Adsorption of Aromatic Compounds onto Cellulose, Langmuir, 2004, 20(8), 3151–3158 CrossRef .
31. G. Todde, S. K. Jha, G. Subramanian and M. K. Shukla, Adsorption of TNT, DNAN, NTO, FOX7, and NQ onto cellulose, chitin, and cellulose triacetate. Insights from Density Functional Theory calculations, Surf. Sci., 2018, 668, 54–60 CrossRef .
32. X. Shu-Wei, S. Wei, Y. Liang-Min and H. Zhe, QSAR Studies on the Antibacterial Activity of Some Substituted 3-Isothiazolinones, Acta Chim. Sin., 2007, 65(23), 2707–2714 Search PubMed .
33. T. Hisaki, M. Aiba née Kaneko, M. Yamaguchi, H. Sasa and H. Kouzuki, Development of QSAR models using artificial neural network analysis for risk assessment of repeated-dose, reproductive, and developmental toxicities of cosmetic ingredients, J. Toxicol. Sci., 2015, 40(2), 163–180 CrossRef CAS PubMed .
34. M. Bhardwaj, N. Masand, J. Sahoo and V. M. Patil, Risk Assessment of Cosmetic Preservatives Using QSAR: QSAR of Cosmetic Preservatives, Int. J. Quant. Struct. Relatsh., 2020, 5(1), 44–62 CAS .
35. S. Jo, T. Kim, V. G. Iyer and W. Im, CHARMM-GUI: A web-based graphical user interface for CHARMM, J. Comput. Chem., 2008, 29(11), 1859–1865 CrossRef CAS PubMed .
36. Q. Shi, K. Zhang, R. Lu and J. Jiang, Water desalination and biofuel dehydration through a thin membrane of polymer of intrinsic microporosity: Atomistic simulation study, J. Membr. Sci., 2018, 545, 49–56 CrossRef CAS .
37. P. P. Deshpande, M. S. Radue, P. Gaikwad, S. Bamane, S. U. Patil, W. A. Pisani and G. M. Odegard, Prediction of the Interfacial Properties of High-Performance Polymers and Flattened CNT-Reinforced Composites Using Molecular Dynamics, Langmuir, 2021, 37(39), 11526–11534 CrossRef CAS .
38. G. Sachdeva, S. U. Patil, S. S. Bamane, P. P. Deshpande, W. A. Pisani, G. M. Odegard and R. Pandey, Mechanical response of polymer/BN composites investigated by molecular dynamics method, J. Mater. Res., 2022, 37(24), 4533–4543 CrossRef CAS .
39. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers, SoftwareX, 2015, 1–2, 19–25 CrossRef .
40. R. B. Best, X. Zhu, J. Shim, P. E. M. Lopes, J. Mittal, M. Feig and A. D. MacKerell Jr., Optimization of the Additive CHARMM All-Atom Protein Force Field Targeting Improved Sampling of the Backbone ϕ, ψ and Side-Chain χ¹ and χ² Dihedral Angles, J. Chem. Theory Comput., 2012, 8(9), 3257–3273 CrossRef CAS PubMed .
41. F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss and V. Dubourg, Scikit-learn: Machine learning in Python, J. Mach. Learn. Res., 2011, 12, 2825–2830 Search PubMed .
42. T. Head, G. L. MechCoder and I. J. V. v. Shcherbatyi, scikit-optimize/scikit-optimize: v0. 5.2, 2018, vol. 5 Search PubMed .
43. V. Busico and R. Cipullo, Microstructure of polypropylene, Prog. Polym. Sci., 2001, 26(3), 443–533 CrossRef CAS .
44. A. B. Raheem, Z. Z. Noor, A. Hassan, M. K. Abd Hamid, S. A. Samsudin and A. H. Sabeen, Current developments in chemical recycling of post-consumer polyethylene terephthalate wastes for new materials production: A review, J. Cleaner Prod., 2019, 225, 1052–1064 CrossRef CAS .
45. B. Medronho and B. Lindman, Brief overview on cellulose dissolution/regeneration interactions and mechanisms, Adv. Colloid Interface Sci., 2015, 222, 502–508 CrossRef CAS .
46. M. Sangkhawasi, T. Remsungnen, A. S. Vangnai, R. P. Poo-arporn and T. Rungrotmongkol, All-Atom Molecular Dynamics Simulations on a Single Chain of PET and PEV Polymers, Polymers, 2022, 14(6), 1161 CrossRef CAS .
47. S.-P. Ju, C.-C. Chen, T.-J. Huang, C.-H. Liao, H.-L. Chen, Y.-C. Chuang, Y.-C. Wu and H.-T. Chen, Investigation of the structural and mechanical properties of polypropylene-based carbon fiber nanocomposites by experimental measurement and molecular dynamics simulation, Comput. Mater. Sci., 2016, 115, 1–10 CrossRef CAS .
48. M. Fukuda and S. Kuwajima, Molecular dynamics simulation of water diffusion in atactic and amorphous isotactic polypropylene, J. Chem. Phys., 1998, 108(7), 3001–3009 CrossRef CAS .
49. H. F. Mark, Encyclopedia of Polymer Science and Technology, John Wiley & Sons, Inc., New York, NY, USA, 2nd edn, 1982, pp. 35–68 Search PubMed .
50. Z. Ren, R. Guo, X. Zhou, H. Bi, X. Jia, M. Xu, J. Wang, L. Cai and Z. Huang, Effect of amorphous cellulose on the deformation behavior of cellulose composites: molecular dynamics simulation, RSC Adv., 2021, 11(33), 19967–19977 RSC .
51. S. E. Wheeler and J. W. G. Bloom, Toward a More Complete Understanding of Noncovalent Interactions Involving Aromatic Rings, J. Phys. Chem. A, 2014, 118(32), 6133–6147 CrossRef CAS .
52. C. R. Martinez and B. L. Iverson, Rethinking the term “pi-stacking”, Chem. Sci., 2012, 3(7), 2191–2201 RSC .
53. B. Damirchi, M. Radue, K. Kanhaiya, H. Heinz, G. M. Odegard and A. C. T. van Duin, ReaxFF Reactive Force Field Study of Polymerization of a Polymer Matrix in a Carbon Nanotube-Composite System, J. Phys. Chem. C, 2020, 124(37), 20488–20497 CrossRef .
54. S. S. Bamane, M. B. Jakubinek, K. Kanhaiya, B. Ashrafi, H. Heinz and G. M. Odegard, Boron Nitride Nanotubes: Force Field Parameterization, Epoxy Interactions, and Comparison with Carbon Nanotubes for High-Performance Composite Materials, ACS Appl. Nano Mater., 2023, 6(5), 3513–3524 CrossRef .
55. A. Etale, A. J. Onyianta, S. R. Turner and S. J. Eichhorn, Cellulose: A Review of Water Interactions, Applications in Composites, and Water Treatment, Chem. Rev., 2023, 123(5), 2016–2048 CrossRef PubMed .
56. D. T. Stanton and P. C. Jurs, Development and use of charged partial surface area structural descriptors in computer-assisted quantitative structure-property relationship studies, Anal. Chem., 1990, 62(21), 2323–2329 CrossRef CAS .
57. R. H. Rohrbaugh and P. C. Jurs, Descriptions of molecular shape applied in studies of structure/activity and structure/property relationships, Anal. Chim. Acta, 1987, 199, 99–109 CrossRef CAS .

---

Footnote

† Electronic supplementary information (ESI) available: Molecular simulation; QSAR modelling; simulation results; QSAR model results. See DOI: https://doi.org/10.1039/d4nr02162b

This journal is © The Royal Society of Chemistry 2024