Multiscale computational insights into 5-fluorouracil delivery via zeolite imidazole frameworks (ZIFs)

Abstract

Zeolitic Imidazolate Frameworks (ZIFs) are considered as potential nanocarriers in biomedical applications such as storage and transportation of drugs, due to their low toxicity, high internal load and controlled release. In this work, the adsorption of the anticancer drug 5-Fluorouracil (5-FU) in selected ZIFs is studied by employing a multiscale computational scheme that includes semi-empirical computational techniques, Grand Canonical Monte Carlo (GCMC) and molecular dynamics simulations. Our investigation is based on ZIF-8 which is characterized by pH-sensitive controlled drug release. In order to improve the capacity and the drug interaction with the framework we expanded the parent ZIF-8, by replacing each imidazole linker with 3-(1H-pyrrol-3-yl)-pyridine. 5-FU shows enhanced binding affinity (34 kcal mol⁻¹) in modified ZIF with respect to the parent ZIF-8 (12 kcal mol⁻¹) (PM7). Moreover, GCMC simulations were employed to determine the loading of 5-FU in both ZIF compounds under different thermodynamic conditions and over a wide range of pressures, where it is revealed that the loading for the modified ZIF (2311 mg g⁻¹) is four times higher than the loading for ZIF-8 (560 mg g⁻¹). Moreover, molecular dynamics simulations offer a detailed comparison between the two ZIF materials in aqueous systems and provide deeper insight into their drug-material interactions.

---

Introduction

Cancer is a complex disease characterized by uncontrolled cell proliferation, resistance to apoptosis, and invasive growth. Its development is influenced by genetic, environmental, and lifestyle-related factors such as exposure to carcinogens, poor diet, and obesity.^1,2 Chemotherapy remains a primary treatment strategy, especially for metastatic cancers, although its efficacy is often limited by off-target toxicity and systemic side effects. 5-Fluorouracil (5-FU), a widely used anticancer agent commonly known as Adrucil, exerts its effect by disrupting DNA synthesis but is associated with adverse effects such as inflammation and hair loss.^3,4 Due to its clinical relevance and small molecular size, 5-FU serves as a useful model compound in computational drug delivery studies.^5–8

To improve therapeutic selectivity and reduce side effects, nanoparticle-based carriers are increasingly employed in chemotherapy to encapsulate, protect, and control the release of anticancer drugs, such as 5-FU.^9,10 Key attributes for these nanoparticles include non-toxicity, stability, biocompatibility, and controlled drug release. Nanoparticles achieve targeted delivery via passive accumulation in tumor tissue or active targeting through specific tumor cell interactions.^{11,12,13–20,21–23} Despite challenges such as low drug-loading capacity and insufficient release control, Zeolitic Imidazolate Frameworks (ZIFs) offer promising advancements in nanoparticle-based drug delivery, potentially overcoming previous limitations. They are a subclass of MOFs^24,25 with zeolite-like topologies, consisting of tetrahedral transition metal ions interconnected by organic linkers rendering them promising materials studied for applications such as gas capture, drug storage and transport, heterogeneous catalysis, and chemical sensing.^26–32 ZIF-8, for instance, demonstrates potential in drug delivery systems due to its high stability in water and rapid decomposition in acidic environments, facilitating pH-sensitive drug release. This property is particularly beneficial for targeting cancerous tissues, which have a more acidic pH compared to normal tissues, minimizing undesired drug release during transport and enhancing targeted delivery.³³ ZIF-8's ability to retain 5-FU under neutral pH and release it rapidly under acidic conditions has been experimentally validated,³² with over half of the drug released within 12 hours at pH 5.0. Our study along with ZIF-8, also investigates a modified ZIF with a larger pore size than ZIF-8. By replacing ZIF-8's imidazole ring linker with 3-(1H-pyrrol-3-yl)-pyridine, the new linker forms a larger unit cell with more available pore volume, increasing drug capacity as indicated by our simulations when compared with other porous materials. The corresponding linkers are highlighted in red in Fig. 1. Calculations using PoreBlazer³⁴ reveal the unique characteristics of these materials, emphasizing their potential in improving drug interaction and expanding the capabilities of ZIF-8 for selective drug delivery, as shown in the Table 1.

Fig. 1 Schematic representation of the structure of ZIF-8 (left) and the modified ZIF (right). The corresponding linker is highlighted in red.

Table 1 Structural properties of ZIFs in study

Structural properties	ZIF-8^35	Modified ZIF
a PLD: pore limiting diameter. b LCD: largest cavity diameter.
PLD^a (Å)	2.93	13.05
LCD^b (Å)	11.36	19.82
Pore volume (Å3)	21 139	58 950

---

Computational chemistry applies computer science to chemical problems, using theoretical models to predict molecular structures and properties. Since solving quantum n-particle problems exactly is impractical, it complements experimental data and can even predict phenomena difficult to observe directly. This field is essential for designing drugs and materials, utilizing methods that vary in accuracy and computational demand based on system size. Approaches include ab initio methods for high accuracy but high cost, empirical methods for faster, approximate results, and semi-empirical methods that balance precision with efficiency. Method choice depends on the system's complexity and resource availability. For this study, a multiscale computational approach^36,37 consisting of three stages has been developed to investigate various aspects of drug-porous material systems, shown in Fig. 2. Initially, a semi-empirical method is used to examine the molecular interactions between 5-FU and the corresponding ZIFs by considering judiciously selected finite size models. In the next stage, Grand Canonical Monte Carlo classical simulations are performed focusing on the analysis of drug uptakes within the pores of the two ZIFs, in order to correlate how porosity features in ZIF materials, such as pore size and available free volume, influence drug loading. Finally, molecular dynamics simulations are employed to model the behavior of drug-loaded porous materials in an aqueous environment, mimicking the conditions found in body circulation, thus providing a more realistic perspective that sheds light on the kinetics of the system.

Fig. 2 The multi-scale approach of our drug delivery system: from quantum mechanics to molecular dynamics.

For the quantum chemical calculations, Gaussian 16,³⁸ a widely used software package for computational chemistry, was employed. Molecular interactions between 5-FU and ZIFs were calculated with the PM7 (Parametric Method 7).³⁹ Semi-empirical methods, like PM7, rely on empirical approaches where the correlation between two electrons in the Hamiltonian is not explicitly considered, using instead a simplified form of the Hamiltonian operator that integrates parameters derived from experimental data. While providing lower accuracy than ab initio methods, semi-empirical methods like PM7 are computationally faster and more suitable for handling larger systems of 50–100 atoms like in this case. Moreover, PM7 is suitable for the accurate prediction of the interactions of small molecules of this magnitude compared to more time-consuming QM methods. Grand Canonical Monte Carlo (GCMC) simulations operate under constant chemical potential, temperature, and volume (μVT ensemble), allowing for the exchange of particles between the simulated system and a reservoir. This dynamic approach enables the exploration of systems in equilibrium with their environment, providing valuable insights into phase transitions, adsorption phenomena, and other thermodynamic properties. For this study, RASPA⁴⁰ software was employed to perform the simulations, with van der Waals interaction parameters for the atoms of each ZIF framework obtained from the Universal Force Field (UFF).^41–44 The 5-FU drug molecule was parametrized using the TraPPE force field.⁴⁵ Lorentz–Berthelot mixing rules describing van der Waals interactions between the MOF framework and the drug molecule were used.⁴⁶ Electrostatic interactions were considered, with point charges assigned to the drug atoms, obtained from quantum mechanics calculations using the PBE level of theory and the def2-TZVP basis set.⁴³ Partial charges for both ZIFs were calculated for the crystallographic geometry in vacuum via a single-point quantum mechanics periodic simulation using the PBE method and the DZVP-MOLOPT-SR-GTH basis set for Zn and DZVP-MOLOPT-GTH for all other atoms, using the CP2K code.⁴⁷ This system was studied under several pressures at 298 K with 50 000 equilibration cycles and 500 000 production cycles. There are four types of moves that are possible to occur in each step with equal probability, insertion of a drug molecule at a random place, deletion of a molecule, translation of an existing molecule, rotation or reinsertion. MD simulations offer insights into the thermodynamics, kinetics, and structural dynamics of molecules, guiding the design of novel materials, drugs, and catalysts. For this study, LAMMPS software⁴⁸ was used for MD simulations, with the parametrization of both ZIFs based on GAFF⁴⁹ or most bonded interactions and van der Waals interactions. The MCPB.py Python algorithm,⁵⁰ compatible with AmberTools⁵¹ is a tool that performs calculations for organometallic systems using B3LYP/6-31G* (by default) and it was employed for bond, angle and dihedral coefficients between Zn and neighbouring atoms that could not be defined accurately by the GAFF force field. The same set of partial charges used in the GCMC calculations was assigned to each atom type here. Tests and comparisons for ZIF-8, including RMSD graphs and optimization results, showed good agreement between the crystallographic and simulated structures of ZIF-8, validating the accuracy of the parametrization set. The average RMSD deviation of ZIF-8 atoms from the crystallographic data was less than 0.2 Å. And finally, to implement water solution, the SPCE water model was selected for use in the simulations due to its efficiency and agreement with experimental data, providing a time-efficient model suitable for large systems.⁵² To systematically investigate the performance of porous materials in chemotherapeutic drug delivery, this study attempts to implement a multiscale computational workflow comprising three distinct stages. First, semi-empirical quantum chemical calculations are employed to evaluate the interaction energies between 5-fluorouracil and the two different ZIF frameworks—ZIF-8 and its structurally modified analogue—providing insights into drug-framework affinity. Second, Grand Canonical Monte Carlo (GCMC) simulations quantify the maximum drug loading capacities under varying thermodynamic conditions, revealing the effect of structural modifications on uptake efficiency. Lastly, Molecular Dynamics (MD) simulations offer a dynamic perspective on the mobility, retention behavior, and selectivity of drug molecules within the frameworks. This integrative approach enables a comprehensive comparison of the two ZIF systems and facilitates the rational assessment of their suitability for targeted and sustained drug delivery. The results of each stage are presented and analyzed in detail in the following section.

Results and discussion

PM7 results

Two carriers, namely ZIF-8 and the extended ZIF, were investigated for their potential as storage and delivery platforms for the 5-FU drug. The primary objective of the calculations at this level of theory was to identify the most suitable candidate in terms of drug-host interaction and analyse the interaction between each ZIF and the drug. Initially, instead of examining the entire unit cell of each ZIF, a small cluster (of 60 and 88 atoms, respectively) consisting of the metal ion and the coordinated four linkers (representing a pocket of the whole cell) was isolated. The coordinates of the nearest neighbor atoms of the metal corners were fixed in order to preserve the coordination of the metal and the angles at which the linkers are connected. Several calculations with different initial conformations were performed on these small structures, along with the drug 5-FU, to determine the spatial arrangement that yielded the minimum total energy. This corresponds to the configuration with the strongest interactions between the ZIF and the drug (as depicted in Fig. 3). Calculations were conducted considering various initial positions and angles to explore different configurations and identify the optimal structure that maximizes the interaction between the ZIF and the drug. The interaction energy was calculated from the formula:

ΔE = E_(ZIF-drug) − E_(ZIF) − E_(drug)

Fig. 3 The interaction energy between the expanded ZIF and the drug (right) is much higher than the ZIF-8 (left) interaction.

The interaction energy between ZIF-8 and 5-fluorouracil (5-FU) was calculated to be −12 kcal mol⁻¹ using the PM7 semi-empirical method. In contrast, the modified ZIF framework exhibited a significantly stronger interaction with 5-FU, with an interaction energy of −34 kcal mol⁻¹, highlighting its enhanced affinity toward the drug molecule. This threefold increase in binding strength can be primarily attributed to the structural features of the modified linker. The expanded linker architecture provides additional interaction sites and facilitates more extensive orbital overlap with 5-FU. Furthermore, the modified framework adopts a spatial conformation that more effectively encapsulates the drug molecule, thereby promoting stronger host–guest interactions. These findings underscore the first indication of the potential of the modified ZIF as a superior platform for drug encapsulation and controlled release applications. Also, the average optimal distance between the drug molecule and the metal of both ZIF compounds is around 6.0 Å, (Fig. 3). This value is also validated by molecular dynamics simulations described below. Additional PM7 method calculations were conducted to investigate the interactions between ZIF fragments and a water molecule, as well as between two drug molecules. To start with, calculations to analyse the interactions between a water molecule and the ZIF fragments were performed. The interaction energy for the ZIF-8 fragment with water was found to be −10 kcal mol⁻¹, while the modified ZIF fragment exhibited an interaction energy of −11 kcal mol⁻¹. This slight difference of 1 kcal mol⁻¹ lies within the accuracy threshold of semi-empirical simulations, rendering it insignificant. Additionally, the weak interaction energies suggest that the water molecule does not strongly bind to the porous material during drug delivery, implying that water molecules do not occupy binding positions instead of drugs within the pores. The interaction energy between two drug molecules was determined to be −9 kcal mol⁻¹, indicating a weak interaction, too. This suggests that there is no significant agglomeration or formation of strong bonds when drug molecules are in close proximity.

GCMC results

The subsequent phase of our study focuses on predicting the uptake of the 5-FU drug within the ZIF pores and determining how the drug molecules interact with the porous structure of the ZIF, providing crucial insight into the capacity and efficiency of the ZIF materials as drug carriers. For these simulations, the unit cell of each ZIF compound was used and gradually filled up randomly with the drug molecule according to the stochastic Grand Canonical Monte Carlo algorithm. We obtained the uptake isotherms for fugacities up to 1 Pascal and at an ambient temperature of 298 K using GCMC simulations until the system was well equilibrated. In chemical thermodynamics, the concept of fugacity represents the effective partial pressure of a real gas, replacing the mechanical pressure for more accurate computation of the chemical equilibrium constant. Essentially, fugacity equates to the pressure of an ideal gas under the same conditions, facilitating a more precise understanding of gas behavior in non-ideal scenarios.

The GCMC simulations were performed in vacuum and no solvent was considered. The drug loading in ZIFs obtained from GCMC simulations are presented in Fig. 4. For ZIF-8, the average maximum loading capacity was calculated to be 560 mg 5-FU per g ZIF-8, equivalent to approximately 13 molecules of 5-fluorouracil per unit cell. In contrast, the modified ZIF exhibited a considerably higher loading capacity of 2310 mg 5-FU per g modified ZIF under the same conditions, translating to about 37 drug molecules per unit cell. This fourfold increase in uptake capacity underscores the superior ability of the modified ZIF to accommodate the physical interaction with the drug molecules, likely due to both enhanced pore volume and the stronger molecular interactions, as corroborated by our earlier interaction energy calculations. Also, we should emphasize that there is physical and not chemical interaction between the ZIFs and the drug, there is no new bonds being formed during the adsorption process. Drug adsorption onto the modified ZIF occurs in two distinct stages before reaching full capacity, attributable to two main factors. In the initial stage, drug molecules preferentially occupy optimal positions near the metal corners of the pores, where interactions are most favorable. As fugacity increases, these preferred sites become saturated, leading to further drug uptake within the remaining pore space in a second stage. Another explanation is based on the internal pore geometry of ZIFs, which contain two distinct pore sizes. These results suggest that the smaller pores are the first to fill since they provide a more confined space with multiple interaction points. Once they are saturated, the larger cavities are subsequently occupied as the applied external pressure increases. These results are also contrasted with the drug uptake capacities of other porous materials such as ZIF-8, MOF-123, and MUT, computed under the same conditions.^53–55 Our findings suggest that the modified ZIF not only meets but surpasses the drug loading capacities of these well-studied materials. This indicates that the linker replacement introduced to the ZIF structure significantly enhances its capability to adsorb the drug molecules. In terms of volumetric uptake, a detailed comparison at the low loading limit reveals that the modified ZIF performs comparably to the parent ZIF-8. However, while ZIF-8 reaches its maximum volumetric uptake at 100 cm³ g⁻¹ framework, the modified ZIF achieves a significantly higher maximum value of 400 cm³ g⁻¹ framework at a higher fugacity. This substantial difference further underscores the enhanced capacity of the modified ZIF to accommodate larger quantities of the drug. The observed superior performance of the modified ZIF can be attributed to several key factors. The binding energy of 5-FU associated with the global minimum plays a crucial role as an indicator of the material's physical adsorption capabilities. In GCMC simulations, higher interaction energies are directly correlated with increased adsorption capacities, highlighting the importance of binding interactions in achieving high drug loading in a physical way. Overall, these computational results collectively demonstrate that the modified ZIF offers significantly improved drug loading and retention capabilities compared to both its parent ZIF-8 and other porous materials, positioning it as a highly effective platform for drug delivery applications.

Fig. 4 The uptake comparison of 5-FU for several porous materials. MUT (green line) and MOF-123 (red line) correspond to data derived from bibliography. The dashed blue line is the experimental uptake of 5-FU for ZIF-8. The uptake for ZIF-8 (blue line) and modified ZIF (orange line) was calculated by performing GCMC simulations.⁵³

Molecular dynamics results

The concluding phase of our multiscale approach involves conducting molecular dynamics simulations to validate preceding findings and analyze the interactions within the system comprising our ZIF compounds and 5-Fluoruracil. To construct a realistic representation closely resembling the human body environment, we developed a simulation model comprising of three main components: the grid of the porous material, multiple 5-FU drug molecules, and water molecules as the solvent. The simulation box was constructed with the use of Packmol⁵⁶ and Moltemplate.⁵⁷ The simulation box is periodic in all dimensions but in one dimension there is water at either side of the ZIF structure in order to simulate the aqueous environment and allow the drugs to have the option to enter or escape the pores in case the ZIF-drug interactions are not strong enough. For the ZIF-8 system depicted in Fig. 5, a scaffold of 2 × 2 × 4 cells was utilized, along with 2400 water molecules to achieve a density of approximately 1.00 gr cm⁻³, and varying quantities of drug molecules depending on the specific aim of the simulation. The drug molecules were randomly preloaded inside the ZIF pores.

Fig. 5 ZIF-8 system in NPT ensemble, black: ZIF-8, red: 5-FU, cyan: water.

Similarly, the system containing the modified ZIF featured a 2 × 2 × 3 cell structure and an equivalent amount of water, but due to the larger internal volume of the modified ZIF, a higher drug uptake was anticipated, necessitating the loading of up to 35 drug molecules per pore until maximum capacity was reached. The maximum amount of drug molecules corresponds to the uptake obtained from the GCMC simulations and they were inserted in the simulation box randomly in this stage. Both systems underwent a standardized simulation protocol after all molecules have been initially placed at their starting positions, starting with a minimization phase to allow molecules in close proximity with each other to relax. Subsequently, the annealing process commenced, gradually heating the system from 300 K to 1000 K, followed by a cooling phase back to 310 K at a rate of 1 fs per step. Each of these stages comprised 1 000 000 steps in the NVT ensemble and all ZIF coordinates were kept rigid because in high temperatures the structure would not remain intact and would not retain the crystallographic cell dimension. This stage is necessary to remove inhomogeneities or vacuum regions that may have occurred due to the random insertion of drug and water molecules. Also, starting from different initial configurations was examined on several occasions and resulted in similar production phases. This was validated by comparing thermodynamic properties such as density, pressure and RMSD analysis which showed that after the annealing process, all initial systems end up in different configurations, thermodynamically similar to each other. The system then entered the production phase, simulated under NPT conditions at 310 K and 1 atm for 10 nanoseconds, where every molecule in the system is free to move with no rigid constrains, with the results analysed from the final 5 nanoseconds of the simulation, ensuring equilibrated system parameters such as temperature, pressure, volume and density. The simulation time of 10 ns is adequate for the two different ZIF systems to show off their differences in drug delivery which are apparent in the analysis that was performed.

The primary outcomes gleaned from such simulations, crucial for understanding material interactions, are the Radial Distribution Function (RDF) graphs. Fig. 6 illustrates the radial distribution of various atoms from the 5-FU drug molecule, including N, O, and F, in proximity to the Zn atoms located at the metal corners of ZIF-8. The black line represents the average radial distribution of each atom of 5-FU around the Zn atoms. A notable peak at approximately 6 Å signifies the average distance of interaction between the drug and the metal corner within ZIF-8 pores. This observation validates the semi-empirical calculations conducted on the metal corner, confirming the most favorable drug-material interactions depicted in Fig. 2. Subsequent less prominent peaks denote additional favorable interactions, indicative of certain drug molecules occupying energetically less favorable positions within the pores due to their high abundance. Furthermore, the RDF reveals two distinct peaks concerning the radial distribution of oxygen atoms of 5-FU around Zn, corresponding to the presence of two oxygen atoms connected to the ring of 5-FU and their approximate 4 Å separation, aligning with the distances of the respective peaks.

Fig. 6 Radial distribution function of Zn atoms of ZIF-8 metal corners with atoms in 5-FU molecules.

The investigation into the behavior of systems loaded with varying amounts of drug molecules in ZIF-8 is elucidated by the RDFs depicted in Fig. 7. A comprehensive spectrum analysis of this material was conducted through simulations involving five different loadings of 5-FU, ranging from a single drug molecule to an excess exceeding the maximum capacity (as determined by GCMC simulations). These RDFs reveal a discernible trend: as the drug loading increases, the pores become densely packed with molecules, resulting in a compacted system that restricts the drug molecules from occupying solely the most favorable positions and interacting optimally with ZIF-8. Conversely, with lower drug quantities, each pore becomes less densely packed, allowing drug molecules to concentrate primarily at the most favorable distance of 6 Å relatively to the Zn metal corners, consistent with earlier findings. Notably, in instances where only one drug molecule occupies the pore, the two distinct peaks representing different interaction energy states become more evident.

Fig. 7 Radial distribution function of Zn atoms of ZIF-8 metal corner with 5-FU molecules on several drug loadings.

The comparative analysis between two systems featuring ZIF-8 and modified ZIF, each loaded with an equivalent quantity of drugs, is presented herein. Fig. 8 depicts the radial distribution function (RDF) of the Zn atom of the metal corner for both materials in conjunction with all 5-FU molecules on average. Notably, a discernible peak is observed in both cases at approximately 6 Å, indicative of the lowest energy state, corroborating earlier calculated interactions. However, a notable distinction arises in the intensity of ZIF-8's RDF, which markedly surpasses that of modified ZIF. This disparity can be rationalized by the fact that within ZIF-8, the quantity of drug molecules occupies the available pore space to its maximal capacity. In contrast, in the context of modified ZIF, the pores exhibit a less densely packed configuration, thereby affording a greater tolerance for the incorporation of additional 5-FU molecules. In the context of the modified ZIF pores, a substantial presence of water molecules persists within each pore, thereby occupying space that could otherwise accommodate additional 5-FU molecules. Conversely, in ZIF-8 pores, the presence of water molecules is significantly diminished or absent entirely following a 10 ns interval. This phenomenon aligns with the principle that in densely packed systems, only molecules that exhibiting stronger interactions endure, thereby rationalizing the observed outcome.

Fig. 8 Radial distribution function of Zn atoms (of both ZIFs) with 5-FU molecules in a system containing 200 drug molecules in total.

In the final step of the study, we conducted a mean square displacement (MSD) analysis to elucidate the diffusion dynamics of drug molecules within the system. It became evident that the drug molecules exhibited movement throughout the bulk of the water, a phenomenon commonly referred to as diffusion, which occurs naturally in equilibrium fluids. However, the scenario becomes more complex when drug molecules are entrapped within interactions with the porous material. In dense fluid environments, the motion of individual molecules is characterized by random collisions with neighboring molecules, resulting in a trajectory akin to a random walk. Albert Einstein famously analyzed such paths in his study of Brownian motion, demonstrating that the mean square displacement of particles following a random walk is proportional to the elapsed time. Our MSD plots (Fig. 9) reveal that in a simulation box containing pure water, 5FU molecules are able to move freely in a random manner, leading to a diffusion coefficient of the same order of magnitude as that typically observed for molecules dissolved in aqueous solution and of course being higher than systems containing ZIF materials, where the freedom of movement is hindered. The high diffusion coefficient of 1.5 × 10⁻⁹ m² s⁻¹ (derived from the red line of the figure) in this scenario indicates minimal resistance to the movement of drug molecules, allowing them to traverse the simulation box unhindered. This value is close to the self-diffusion coefficient of neat water that is 2.3 × 10⁻⁹ m² s⁻¹ at 298 K.⁵⁸ When comparing two modified ZIF systems containing 200 and 300 molecules of 5-FU (blue and green line respectively), respectively, it is expected that the more densely packed system exhibits a lower diffusion coefficient. This restriction results from the drug molecules frequently colliding with the internal surfaces of the pores, hindering their free movement. It becomes evident that the more densely packed system exhibits a lower diffusion coefficient (5.1 × 10⁻¹⁰ m² s⁻¹ for 200 5-FU in modZIF (orange) and 1.8 × 10⁻¹⁰ m² s⁻¹ for 300 5-FU in modZIF (green)). Furthermore, the diffusion phenomenon appears comparable between the ZIF and modified ZIF systems, when they are loaded with the same number of molecules. In the case of ZIF-8, the low diffusion coefficient can be attributed to the small pore apertures hindering drug molecule escape, leading to a low diffusion coefficient at 8.0 × 10⁻¹⁰ m² s⁻¹, thereby facilitating drug delivery until the structure is dissolved by the low pH environment of cancerous tissue. The similarity in diffusion coefficients between the modified ZIF and ZIF-8 systems indicates that while modified ZIFs offer larger cavities and pore sizes facilitating increased drug uptake and stronger interactions, the drugs remain bound inside the material and are less likely to escape once engaged in interactions with the framework. Overall, the MSD analysis highlights the interplay between molecular mobility and material structure. In bulk water, drug molecules diffuse freely, whereas in the constrained environments of ZIF-8 and modified ZIF, their movement is significantly hindered. This restriction is crucial for the controlled release of drugs in delivery applications, ensuring that the drug remains encapsulated until it reaches the desired site of action. The findings from our MSD analysis provide a deeper understanding of the diffusion dynamics in these systems, informing the design and optimization of ZIF-based drug delivery platforms.

Fig. 9 Mean square displacement plots and diffusion coefficients for 5-FU molecules in ZIF pores and in pure water.

Conclusions

This study comprehensively demonstrates the promise of ZIFs, particularly ZIF-8 and its structurally modified analogue, as nanocarriers for controlled chemotherapeutic drug delivery. Through a combination of multiscale computational techniques, we identify the modified ZIF as a markedly superior candidate for the encapsulation and sustained release of 5-FU, a widely used anticancer agent. Semi-empirical quantum chemical calculations revealed that the interaction energy of −34 kcal mol⁻¹ (PM7) between 5-FU and the metal node in the modified ZIF reached nearly three times stronger than the −12 kcal mol⁻¹ observed in pristine ZIF-8. This enhanced interaction is attributed to the structural modifications introduced in the framework, including the incorporation of an expanded linker, which increases the pore volume and facilitates greater orbital overlap. These features enable the formation of additional interaction points between the framework and the drug, promoting tighter binding and reducing the likelihood of premature drug desorption. Grand Canonical Monte Carlo (GCMC) simulations provided further validation of these findings, showing a fourfold increase in the maximum 5-FU uptake in the modified ZIF compared to ZIF-8 under equilibrium conditions. This considerable improvement in adsorption capacity is critical for increasing therapeutic loading efficiency and minimizing dosage frequency. Molecular Dynamics (MD) simulations offered molecular-level insights into the stability and selectivity of drug loading. The results confirmed that 5-FU remains preferentially bound within the framework's pores rather than interacting with surrounding water molecules, highlighting the selectivity of the porous matrix. Importantly, mean square displacement (MSD) analysis revealed a significant reduction in 5-FU mobility within the pores of the modified ZIF—approximately an order of magnitude lower than in bulk aqueous solution—emphasizing the framework's ability to retain the drug effectively. Such confinement is essential in maintaining local drug concentrations at therapeutic levels while minimizing systemic toxicity. Collectively, these results not only underscore the potential of the modified ZIF as a highly efficient platform for chemotherapeutic delivery but also establish a methodological framework for the rational design and computational evaluation of porous drug carriers. Future experimental studies are encouraged to validate these computational predictions and explore the broader applicability of the proposed structural modifications in clinical drug delivery systems. These results highlight the promise of ZIFs as effective nanocarriers for chemotherapeutic agents and provide a strong foundation for future experimental validation and optimization. Continued exploration of linker functionalization and framework modification may pave the way for the rational design of advanced nanoporous drug carriers. This approach can be readily extended to other drug molecules and porous materials, offering valuable pathways for tailored, targeted delivery strategies in oncology. Such developments could significantly enhance the efficacy of chemotherapeutic delivery while reducing toxicity to healthy tissues. Moreover, this approach can be extended to other clinically relevant chemotherapeutic drugs, supporting the development of targeted and controlled delivery platforms across various cancer types. Tailoring porous materials to the physicochemical properties of specific drugs may improve therapeutic outcomes and broaden the applicability of chemotherapy. Finally, the multiscale computational methodology employed in this study represents a robust framework for the systematic investigation of porous materials. As demonstrated here, it is particularly effective for evaluating newly developed materials alongside established benchmarks, providing valuable insights for material screening and design in drug delivery applications.

Author contributions

Michail Vlachos: writing – original draft, writing – review & editing, conceptualization, methodology, investigation, formal analysis, visualization. Giorgio Turtù: methodology, investigation, formal analysis. Marco Severi: investigation, formal analysis. Emmanuel Tylianakis: writing – review & editing, supervision, investigation. Emmanuel Klontzas: writing – review & editing, investigation, resources. Francesco Zerbetto: writing – review & editing, supervision, conceptualization, methodology. George Froudakis: writing – review & editing, supervision, conceptualization, methodology, funding aquisition.

Conflicts of interest

There are no conflicts to declare.

Data availability

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

Acknowledgements

This work was supported through a studentship in the University of Crete, Greece and a traineeship in the University of Bologna, Italy.

References

1. D. S. Shewach and R. D. Kuchta, Introduction to cancer chemotherapeutics, Chem. Rev., 2009, 109, 2859–2861.
2. S. Mosleh-Shirazi, M. Abbasi, M. R. Moaddeli, A. Vaez, M. Shafiee and S. R. Kasaee, et al., Nanotechnology Advances in the Detection and Treatment of Cancer: An Overview, Nanotheranostics, 2022, 6(4), 400–423.
3. S. Ghafouri-Fard, A. Abak, F. Tondro Anamag, H. Shoorei, F. Fattahi and S. A. Javadinia, et al., 5-Fluorouracil: A Narrative Review on the Role of Regulatory Mechanisms in Driving Resistance to This Chemotherapeutic Agent, Front. Oncol., 2021, 11, 1–21.
4. P. Xie, J. L. Mo, J. H. Liu, X. Li, L. M. Tan and W. Zhang, et al., Pharmacogenomics of 5-fluorouracil in colorectal cancer: review and update, Cell. Oncol., 2020, 43(6), 989–1001.
5. D. B. Longley, D. P. Harkin and P. G. Johnston, 5-Fluorouracil: Mechanisms of action and clinical strategies, Nat. Rev. Cancer, 2003, 3(5), 330–338.
6. Y. Zhang, Y. Qu and Y. Z. Chen, Influence of 6-shogaol potentiated on 5-fluorouracil treatment of liver cancer by promoting apoptosis and cell cycle arrest by regulating AKT/mTOR/MRP1 signalling, Chin. J. Nat. Med., 2022, 20(5), 352–363,  DOI:10.1016/S1875-5364(22)60174-2.
7. P. Guo, C. Pi, S. Zhao, S. Fu, H. Yang and X. Zheng, et al., Oral co-delivery nanoemulsion of 5-fluorouracil and curcumin for synergistic effects against liver cancer, Expert Opin. Drug Delivery, 2020, 17(10), 1473–1484,  DOI:10.1080/17425247.2020.1796629.
8. Y. Zhihong, Y. Ye, A. Pejhan, A. H. Nasr, N. Nourbakhsh and R. Tayebee, A theoretical study on the pure and doped ZnO nanoclusters as effective nanobiosensors for 5-fluorouracil anticancer drug adsorption, Appl. Organomet. Chem., 2020, 34(4), 1–10.
9. S. A. AlQahtani, G. I. Harisa, A. H. Alomrani, F. K. Alanazi and M. M. Badran, Improved pharmacokinetic and biodistribution of 5-fluorouracil loaded biomimetic nanoerythrocytes decorated nanocarriers for liver cancer treatment, Colloids Surf., B, 2021, 197, 111380,  DOI:10.1016/j.colsurfb.2020.111380.
10. Z. Liu, Z. Yu, L. Zhang, R. Tayebee, S. Pourmojahed and E. E. Shahri, et al., Studying loading of 5-fluorouracil onto the surface of gold–silver nanocomposite utilizing DLS and UV–Vis spectroscopy, Chem. Pap., 2023, 77(10), 6135–6145,  DOI:10.1007/s11696-023-02927-9.
11. X. Sun, G. Huang, F. Shao, R. Tayebee, Z. Jalili and E. E. Shahri, Au-doped Fe3O4 nanoparticle as an effective nanocarrier for delivery of 5-fluorouracil anti-liver cancer drug, Inorg. Chem. Commun., 2024, 167, 112685,  DOI:10.1016/j.inoche.2024.112685.
12. J. Du, L. A. Lane and S. Nie, Stimuli-responsive nanoparticles for targeting the tumor microenvironment, J. Controlled Release, 2015, 219, 205–214.
13. C. V. Soares, A. A. Leitão and G. Maurin, Computational evaluation of the chemical warfare agents capture performances of robust MOFs, Microporous Mesoporous Mater., 2019, 280, 97–104.
14. P. Horcajada, C. Serre, G. Maurin, N. A. Ramsahye, F. Balas and M. Vallet-Regí, et al., Flexible porous metal-organic frameworks for a controlled drug delivery, J. Am. Chem. Soc., 2008, 130(21), 6774–6780.
15. M. Kotzabasaki, E. Tylianakis, E. Klontzas and G. E. Froudakis, OH-functionalization strategy in Metal-Organic Frameworks for drug delivery, Chem. Phys. Lett., 2017, 685, 114–118,  DOI:10.1016/j.cplett.2017.07.053.
16. S. Keskin and S. Kizilel, Biomedical Applications of Metal Organic Frameworks, Ind. Eng. Chem. Res., 2011, 50(4), 1799–1812.
17. P. Horcajada, C. Serre, M. Vallet-Regí, M. Sebban, F. Taulelle and G. Férey, Metal-organic frameworks as efficient materials for drug delivery, Angew. Chem., Int. Ed., 2006, 45(36), 5974–5978.
18. D. Giliopoulos, A. Zamboulis, D. Giannakoudakis and K. Triantafyllidis, Polymer/Metal Organic Framework (MOF) Nanocomposites for Biomedical Applications, Molecules, 2020, 25(1), 185.
19. M. C. Bernini, D. Fairén-Jimenez, M. Pasinetti, A. J. Ramirez-Pastor and R. Q. Snurr, Screening of Bio-Compatible Metal-Organic Frameworks as Potential Drug Carriers using Monte Carlo Simulations, 2014. Available from: https://people.ds.cam.ac.uk/df334.
20. M. Kotzabasaki and G. E. Froudakis, Review of computer simulations on anti-cancer drug delivery in MOFs, Inorganic Chemistry Frontiers. Royal Society of Chemistry, 2018, vol. 5, pp. 1255–1272.
21. M. Vassaki, K. E. Papathanasiou, C. Hadjicharalambous, D. Chandrinou, P. Turhanen and D. Choquesillo-Lazarte, et al., Self-sacrificial MOFs for ultra-long controlled release of bisphosphonate anti-osteoporotic drugs, Chem. Commun., 2020, 56(38), 5166–5169.
22. Y. Chen, X. Zhang, M. R. Mian, F. A. Son, K. Zhang and R. Cao, et al., Structural Diversity of Zirconium Metal-Organic Frameworks and Effect on Adsorption of Toxic Chemicals, J. Am. Chem. Soc., 2020, 142(51), 21428–21438.
23. C. Lin, H. He, Y. Zhang, M. Xu, F. Tian and L. Li, et al., Acetaldehyde-modified-cystine functionalized Zr-MOFs for pH/GSH dual-responsive drug delivery and selective visualization of GSH in living cells, RSC Adv., 2020, 10(6), 3084–3091.
24. C. G. Livas, E. Tylianakis and G. E. Froudakis, Enhancing NO Uptake in Metal-Organic Frameworks via Linker Functionalization. A Multi-Scale Theoretical Study, Chem, 2022, 4(4), 1300–1311.
25. D. Raptis, C. Livas, G. Stavroglou, R. M. Giappa, E. Tylianakis and T. Stergiannakos, et al., Surface Modification Strategy for Enhanced NO2 Capture in Metal–Organic Frameworks, Molecules, 2022, 27(11), 3448.
26. J. P. Dürholt, G. Fraux, F. X. Coudert and R. Schmid, Ab Initio Derived Force Fields for Zeolitic Imidazolate Frameworks: MOF-FF for ZIFs, J. Chem. Theory Comput., 2019, 15(4), 2420–2432.
27. K. Zhang, R. P. Lively, M. E. Dose, A. J. Brown, C. Zhang and J. Chung, et al., Alcohol and water adsorption in zeolitic imidazolate frameworks, Chem. Commun., 2013, 49(31), 3245–3247.
28. R. J. Verploegh, A. Kulkarni, S. E. Boulfelfel, J. C. Haydak, D. Tang and D. S. Sholl, Screening Diffusion of Small Molecules in Flexible Zeolitic Imidazolate Frameworks Using a DFT-Parameterized Force Field, J. Phys. Chem. C, 2019, 123(14), 9153–9167.
29. Y. G. Proenza and R. L. Longo, Simulation of the Adsorption and Release of Large Drugs by ZIF-8, J. Chem. Inf. Model., 2020, 60(2), 644–652.
30. M. K. Rana, G. B. Suffritti, P. Demontis and M. Masia, Simulation study of CO2 adsorption properties in small Zeolite Imidazolate Frameworks, Chem. Phys. Lett., 2013, 580, 99–102.
31. Q. Wang, Y. Sun, S. Li, P. Zhang and Q. Yao, Synthesis and modification of ZIF-8 and its application in drug delivery and tumor therapy, RSC Advances. Royal Society of Chemistry, 2020, vol. 10, pp. 37600–37620.
32. C. Y. Sun, C. Qin, X. L. Wang, G. S. Yang, K. Z. Shao and Y. Q. Lan, et al., Zeolitic imidazolate framework-8 as efficient pH-sensitive drug delivery vehicle, Dalton Trans., 2012, 41(23), 6906–6909.
33. X. Y. Sun, H. J. Zhang, X. Y. Zhao, Q. Sun, Y. Y. Wang and E. Q. Gao, Dual functions of pH-sensitive cation Zr-MOF for 5-Fu: large drug-loading capacity and high-sensitivity fluorescence detection, Dalton Trans., 2021, 50(30), 10524–10532.
34. L. Sarkisov, R. Bueno-Perez, M. Sutharson and D. Fairen-Jimenez, Materials Informatics with PoreBlazer v4.0 and the CSD MOF Database, Chem. Mater., 2020, 32(23), 9849–9867.
35. M. Bergaoui, M. Khalfaoui, A. Awadallah-F and S. Al-Muhtaseb, A review of the features and applications of ZIF-8 and its derivatives for separating CO2 and isomers of C3- and C4- hydrocarbons, J. Nat. Gas Sci. Eng., 2021, 96, 104289,  DOI:10.1016/j.jngse.2021.104289.
36. M. Kotzabasaki, I. Galdadas, E. Tylianakis, E. Klontzas, Z. Cournia and G. E. Froudakis, Multiscale simulations reveal IRMOF-74-III as a potent drug carrier for gemcitabine delivery, J. Mater. Chem. B, 2017, 5(18), 3277–3282.
37. C. G. Livas, D. Raptis, E. Tylianakis and G. E. Froudakis, Multiscale Theoretical Study of Sulfur Dioxide (SO2) Adsorption in Metal–Organic Frameworks, Molecules, 2023, 28(7), 3122.
38. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. a. Robb, J. R. Cheeseman, et al., G16_C01. Gaussian 16, Revision C.01, 2016, Gaussian, Inc., Wallin.
39. J. J. P. Stewart, Optimization of parameters for semiempirical methods VI: more modifications to the NDDO approximations and re-optimization of parameters, J. Mol. Model, 2013, 19(1), 1–32.
40. D. Dubbeldam, S. Calero, D. E. Ellis and R. Q. Snurr, RASPA: molecular simulation software for adsorption and diffusion in flexible nanoporous materials, Mol. Simul., 2016, 42(2), 81–101.
41. A. K. Rappe, C. J. Casewit, K. S. Colwell and W. A. Goddard, UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations, J. Am. Chem. Soc., 1992, 114(25), 10024–10035.
42. D. E. Coupry, M. A. Addicoat and T. Heine, Extension of the Universal Force Field for Metal-Organic Frameworks, J. Chem. Theory Comput., 2016, 12(10), 5215–5225.
43. M. K. Rana, F. G. Pazzona, G. B. Suffritti, P. Demontis and M. Masia, Estimation of partial charges in small zeolite imidazolate frameworks from density functional theory calculations, J. Chem. Theory Comput., 2011, 7(6), 1575–1582.
44. X. Wu, J. Huang, W. Cai and M. Jaroniec, Force field for ZIF-8 flexible frameworks: atomistic simulation of adsorption, diffusion of pure gases as CH₄, H₂, CO₂ and N₂, RSC Adv., 2014, 4, 16503–16511.
45. B. L. Eggimann, A. J. Sunnarborg, H. D. Stern, A. P. Bliss and J. I. Siepmann, An online parameter and property database for the TraPPE force field, Mol. Simul., 2014, 40(1–3), 101–105.
46. S. L. Mayo, B. D. Olafson and W. A. Goddard, DREIDING: a generic force field for molecular simulations, J. Phys. Chem., 1990, 94(26), 8897–8909.
47. T. D. Kühne, M. Iannuzzi, M. Del Ben, V. V. Rybkin, P. Seewald and F. Stein, et al., CP2K: An electronic structure and molecular dynamics software package -Quickstep: Efficient and accurate electronic structure calculations, J. Chem. Phys., 2020, 152, 194103.
48. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown and P. S. Crozier, et al., LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comput. Phys. Commun., 2022, 271, 108171,  DOI:10.1016/j.cpc.2021.108171.
49. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, Development and testing of a general amber force field, J. Comput. Chem., 2004, 25(9), 1157–1174.
50. P. Li and K. Merz, MCPB.py: A Python Based Metal Center Parameter Builder, J. Chem. Inf. Model., 2016, 56(4), 599–604.
51. D. A. Case, H. M. Aktulga, K. Belfon, D. S. Cerutti, G. A. Cisneros and V. W. D. Cruzeiro, et al., AmberTools, J. Chem. Inf. Model., 2023, 63(20), 6183–6191.
52. P. Mark and L. Nilsson, Structure and dynamics of the TIP3P, SPC, and SPC/E water models at 298 K, J. Phys. Chem. A, 2001, 105(43), 9954–9960.
53. R. K. Alavijeh, K. Akhbari, E. Tylianakis, G. E. Froudakis and J. M. White, Two-Fold Homointerpenetrated Metal-Organic Framework with the Potential for Anticancer Drug Loading Using Computational Simulations, Cryst. Growth Des., 2021, 21(11), 6402–6410.
54. M. Gomar and S. Yeganegi, Adsorption of 5-Fluorouracil and Thioguanine drugs into ZIF-1, ZIF-3 and ZIF-6 by simulation methods, Mater. Sci. Eng., C, 2019, 97, 461–466.
55. M. Gomar and S. Yeganegi, Adsorption of 5-fluorouracil, hydroxyurea and mercaptopurine drugs on zeolitic imidazolate frameworks (ZIF-7, ZIF-8 and ZIF-9), Microporous Mesoporous Mater., 2017, 252, 167–172.
56. L. Martinez, R. Andrade, E. G. Birgin and J. M. Martínez, PACKMOL: A package for building initial configurations for molecular dynamics simulations, J. Comput. Chem., 2009, 30(13), 2157–2164.
57. A. I. Jewett, D. Stelter, J. Lambert, S. M. Saladi, O. M. Roscioni and M. Ricci, et al., Moltemplate: A Tool for Coarse-Grained Modeling of Complex Biological Matter and Soft Condensed Matter Physics, J. Mol. Biol., 2021, 433(11), 166841,  DOI:10.1016/j.jmb.2021.166841.
58. S. H. Young and H. L. Cherl, Self-diffusion coefficient of water in tofu determined by pulsed field gradient nuclear magnetic resonance, J. Agric. Food Chem., 2006, 54(1), 219–223.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5pm00058k

---