Cancer remains hard to treat, partially due to the non-specificity of chemotherapeutics. Metal–organic frameworks (MOFs) are promising carriers for targeted chemotherapy, yet, to date, there have been few detailed studies to systematically enhance drug loading while maintaining controlled release. In this work, we investigate which molecular simulation methods best capture the experimental uptake and release of cisplatin from UiO-66 and UiO-66(NH2). We then screen a series of biocompatible, pH-sensitive zeolitic imidazolate frameworks (ZIFs) for their ability to retain cisplatin in healthy parts of the patient and release it in the vicinity of a tumor. Pure-component GCMC simulations show that the maximum cisplatin loading depends on the pore volume. To achieve this maximum loading in the presence of water, either the pore size needs to be large enough to occupy both cisplatin and its solvation shell or the MOF–cisplatin interaction must be more favorable than the cisplatin–shell interaction. Both solvated and non-solvated simulations show that cisplatin release rates can be controlled by either decreasing the pore limiting diameters or by manipulating framework–cisplatin interaction energies to create strong, dispersed adsorption sites. The latter method is preferable if cisplatin loading is performed from solution into a pre-synthesized framework as weak interaction energies and small pore window diameters will hinder cisplatin uptake. Here, ZIF-82 is most promising. If it is possible to load cisplatin during crystallization, ZIF-11 would outcompete the other MOFs screened as cisplatin cannot pass through its pore windows; therefore, release rates would be purely driven by the pH triggered framework degradation.

In 2018, the World Health Organization (WHO) reported that cancer is responsible for 17% of deaths globally,1 and the number of cases are expected to rise by 2030.2 Chemotherapy is one of the primary cancer treatment methods, during which cytotoxic drugs, such as cisplatin [cis-diaminedichloroplatinum (II) (CPT)], are usually administered intravenously. Once incorporated into a cell, cisplatin activates and is able to form covalent bonds with nucleotides.3 This damages DNA, preventing cell mitosis and activating apoptosis (i.e., the controlled self-destruction of cells). Unfortunately, intravenous administration is non-specific, and therefore, cisplatin (among other cytotoxic drugs) can have adverse impacts on normal healthy cells. This causes the harmful side effects of chemotherapy, such as hemorrhages, fatigue, damage to the nervous system, kidneys, bladder, and so on. Furthermore, indirect administration reduces the dose reaching cancerous cells, increasing the chance of the tumor developing CPT resistance. To combat these issues, a non-toxic carrier can be used to deliver cisplatin directly to the tumor.

FIG. 1.

Structure of (a) ZIF-8 and (b) ZIF-11. View down the crystallographic c-axis. Hydrogens are omitted for clarity (dark gray = carbon, blue = nitrogen, light gray = zinc).

FIG. 1.

Structure of (a) ZIF-8 and (b) ZIF-11. View down the crystallographic c-axis. Hydrogens are omitted for clarity (dark gray = carbon, blue = nitrogen, light gray = zinc).

Close modal

Existing drug delivery carriers can be subcategorized as organic and inorganic. Organic carriers are typically biocompatible, and they include materials such as polymeric structures,4 liposomes,5 and nanocapsules.6 However, they usually have low loading efficiency and slow drug release. Inorganic carriers include materials such as ferromagnetic,7 gold,8 and mesoporous silica nanoparticles,9 the properties of which can be easily tailored for targeted drug delivery and high loadings. However, they tend to be inert and non-biocompatible.10 Metal–organic frameworks (MOFs) can offer the best of both worlds. MOFs consist of metal nodes coordinated by organic ligands forming highly porous materials with large specific internal surface areas and easily tunable compositions and functionalities. Weak coordination bonds between the primary building units allow some MOFs to deconstruct inside the body, and their large pore volumes and internal surface areas enable high drug loading capacities.11 

Various experimental groups have successfully loaded cisplatin into metal–organic frameworks, which enabled them to achieve controlled release rates. For example, Taylor-Pashow et al. successfully loaded CPT into MIL-101(Fe), and to reduce the release rate, they coated the MOF nanoparticles with sodium metasilicate, which slowed down the deconstruction of the framework.12 Several groups have focused on using members of the UiO series of MOFs [Zr6O4(OH)4 nodes] because the strong Zr–O bonds and high connectivity of the nodes increases their stability in aqueous environments. For example, He et al. incorporated cisplatin into UiO nanoparticles via conjugation (i.e., the covalent binding of a cisplatin prodrug to functional groups in the framework), and functionalized the surface for the attachment of small interfering RNAs to achieve CPT/siRNA co-delivery.13 Lin et al. incorporated missing ligand defects in UiO-66 to increase the specific cisplatin uptake compared to the defect free structure.14 Mocniak et al. compared encapsulation and conjugation as methods to uptake a cisplatin pro-drug into UiO-66 and the amine functionalized variant UiO-66(NH2). They achieved a higher uptake and a slower release rate using conjugation as opposed to encapsulation, since the conjugation method involves the formation of pro-drug framework covalent bonds. Compared to encapsulated cisplatin that can passively diffuse through the framework, conjugated cisplatin can only move through the framework, provided that the peptide bonds are broken via hydrolysis.15 

While a slow release rate is beneficial in healthy parts of the body, when in the vicinity of the tumor, it could be problematic, particularly if the carrier is not anchored in its vicinity. One way to speed up delivery is to induce framework degradation. Some MOFs undergo structural changes in response to variations in external stimuli, for example, changes in pH can hydrolyze metal–ligand bonds.16 For the application of chemotherapy, this is advantageous as the extracellular pH (pHe) of tumor tissues are often acidic (pHe = 6.5–6.9) due to their enhanced lactate secretion from anaerobic glycolysis.17 On the other hand, the pHe surrounding normal healthy cells is heavily regulated at pHe = 7.2–7.5.18 If the carrier MOF can retain cisplatin while intact in normal areas of the body and rapidly release it in the vicinity of a tumor as the framework degrades, this could result in effective targeted drug delivery. This pH sensitive controlled release has been demonstrated experimentally for ZIF-8 for a variety of anti-cancer drugs, e.g., for combined cisplatin and oleanolic acid,19 5-FU,20 DOX,21 and benzonidazole.22 

In this work, we first assess a variety of simulation methods for their suitability to provide information about cisplatin uptake and release by comparing the simulation results to available experimental results of cisplatin uptake and release rates for UiO-66 and UiO-66(NH2).15 We then computationally screen a series of biocompatible, pH sensitive zeolitic imidazolate frameworks (ZIFs—8/11/68/70/78/79/82) that are stable at pH > 7 and that will deconstruct in more acidic environments16 and highlight property–performance relationships with respect to the ability of each ZIF to uptake and retain cisplatin.

Out of the ZIFs screened, ZIF-8 [Zn2(mIm)12] was studied because it is one of the most widely investigated MOFs for drug delivery applications as highlighted in the review by Wang et al.23 Toxicology reports on ZIF-8 show that below a threshold concentration, ZIF-8 is biocompatible and will not cause cell toxicity. Above the threshold concentration, ZIF-8’s toxicity mechanism is caused by the Zn2+ ions as the framework deconstructs, as opposed to the organic ligand that has a significantly lower IC50 number (this is the case for a wide variety of frameworks highlighted in the work by Tamames-Taber et al.24).25 ZIF-8 has a sodalite (SOD) topology that consists of uniform 11.6 Å pores constructed of 24 vertices (Zn nodes) each coordinated by four methyl imidazole ligands (mIm) forming eight 6 MR [(MR) membered rings] and six 4 MR/cavity. The pore window diameter of ZIF-8 (∼3.4 Å) can change depending on interactions between the framework atoms and guest molecules.26 This window diameter is comparable to the size of cisplatin (∼5 Å); thus, CPT diffusion will be substantially hindered by gate-opening energetic barriers. (In comparison, ZIFs 68–82 have substantially larger pore windows). ZIF-11 consists of 14.6 Å cavities constructed of benzene imidazole ligands (bIm) that coordinate 48 Zn nodes into a RHO topology. This forms pores that consist of 6 octahedral, 8 hexagonal, and 12 square windows and octagonal channels that connect the main cavities.27 It was selected for screening as the largest pore-window size in ZIF-11 is ∼3 Å, and similar to ZIF-8, the bIm ligands can rotate, providing another example of diffusion hinderance due to framework displacement. The structures of ZIF-8 and ZIF-11 are shown in Fig. 1.

ZIFs-68/70/78/79/82 were selected for screening because they have the same topology, therefore allowing us to compare the influence of functionalization on cisplatin uptake and retention. They consist of three cage types: (i) hpr, (ii) gme, and (iii) kno. The cages consist of Zn nodes, each coordinated to four imidazole ligands (two nIm and two substitute-Im, as shown in Fig. 2). Each structure consists of channels formed by the larger kno cages (24 Zn connecting of three 4 MR, three 8 MR, and two 12 MR) down the c-axis. These channels are surrounded by four parallel rows of gme cavities (18 Zn connecting nine 4 MR, two 6 MR, and three 8 MR), which are connected by hpr cages (12 Zn connecting six 4 MR and two 6 MR).28 

Several research groups have used molecular simulations to screen MOFs for drug delivery. For example, Bernini et al. used grand-canonical Monte Carlo (GCMC) simulations to determine the adsorption behavior of ibuprofen in a selection of non-biocompatible MOFs, for which they had experimental results to validate their simulations (MIL-53, MIL-100, and MIL-101). They then used their simulations to screen biocompatible MOFs, highlighting thermodynamic details behind their performance as drug delivery carriers, based on their maximum loading capacity and adsorption energies.29 Similarly, Erucar and Keskin performed Monte Carlo simulations of adsorption of MTX and 5-FU, two chemotherapy drugs, in a selection of MOFs and compared their simulation results to experimental data where available. They also used molecular dynamics (MD) simulations to study the rate of drug molecule diffusion.30 

To the best of our knowledge, the majority of theoretical research studies on drug delivery in MOFs use pure-component grand-canonical Monte Carlo (GCMC) or equilibrium molecular dynamics (MD) simulations.31–34 However, it is unclear if these are enough to assess MOFs for drug delivery. Here, we present one of the first MOF–drug theoretical research papers in which a wide range of equilibrium and non-equilibrium simulation techniques (GCMC, equilibrium and non-equilibrium MD, umbrella sampling, ab initio optimizations and static energies, geometric analysis, and simulation of polyhedra software GASP) with and without water are tested for their ability to capture experimental cisplatin uptake and release rates. Making use of experimental data available for cisplatin release, we initially use a wide range of simulation techniques on UiO-66 and UiO-66(NH2) to see how well each method compares to the available experimental data.15 Based on this assessment, we then identify techniques most suitable for screening MOFs for cisplatin delivery. We then use our validated methodology to screen pH sensitive ZIFs to highlight key structural parameters to look for in cisplatin drug delivery carriers.

FIG. 2.

Top: ligands used in ZIFs-68/70/78/79/82.28 Bottom: the structure of ZIF-70 as a representative example of the topologies for the isostructural ZIFs-68/70/78/79/82 (view down the crystallographic c-axis; pink = hydrogen, dark gray = carbon, blue = nitrogen, red = oxygen, light gray = zinc).

FIG. 2.

Top: ligands used in ZIFs-68/70/78/79/82.28 Bottom: the structure of ZIF-70 as a representative example of the topologies for the isostructural ZIFs-68/70/78/79/82 (view down the crystallographic c-axis; pink = hydrogen, dark gray = carbon, blue = nitrogen, red = oxygen, light gray = zinc).

Close modal
FIG. 3.

(a) Structure of UiO-66 [Zr6O4(OH)4BDC6]. Representation of the accessible pore volume for a 2 Å probe (pink) generated using Mercury Software.55 Color representations are Zr = turquois, C = gray, O = red, and N = blue. Hydrogens are omitted for clarity. (b) Schematic of models used in alchemical simulations [a = one cisplatin molecule in bulk water, b = cisplatin and water molecules in UiO-66/UiO-66(NH2)]. ∆Ga→b and ∆Gb→a are the changes in Gibbs free energies when moving one cisplatin molecule from a → b and b → a. At equilibrium, ∆Ga→b = ∆Gb→a. (c) Hess cycle used to determine the equilibration loading of cisplatin and water in UiO-66 and UiO-66(NH2). λ = 0 and λ = 1 refer to the beginning and end states of cisplatin during decoupling.

FIG. 3.

(a) Structure of UiO-66 [Zr6O4(OH)4BDC6]. Representation of the accessible pore volume for a 2 Å probe (pink) generated using Mercury Software.55 Color representations are Zr = turquois, C = gray, O = red, and N = blue. Hydrogens are omitted for clarity. (b) Schematic of models used in alchemical simulations [a = one cisplatin molecule in bulk water, b = cisplatin and water molecules in UiO-66/UiO-66(NH2)]. ∆Ga→b and ∆Gb→a are the changes in Gibbs free energies when moving one cisplatin molecule from a → b and b → a. At equilibrium, ∆Ga→b = ∆Gb→a. (c) Hess cycle used to determine the equilibration loading of cisplatin and water in UiO-66 and UiO-66(NH2). λ = 0 and λ = 1 refer to the beginning and end states of cisplatin during decoupling.

Close modal

In our classical simulations, cisplatin was modeled using the forcefield parameters developed by Yesylevskyy et al.35 Since this forcefield was parameterized with the TIP3P model for water, the same water model was also used in this work.36 For the frameworks, the structures were taken from the Cambridge Structural Database and periodic boundary conditions were used in all directions to model infinite crystals. Non-bonded Lennard-Jones and bonded parameters for the frameworks were obtained from the Universal forcefield.37 Electrostatic interactions were calculated using Mulliken charges, determined using ab initio static energy calculations on fragments of the bulk (as described in the supplementary material). While Mulliken charges are known to be dependent on the molecular geometry and the basis set,38 we will show later by comparing results from classical simulations with ab initio simulations that our choice of combining UFF with Mulliken charges is adequate to compare different MOF structures.

Grand-canonical Monte Carlo (GCMC) simulations at 310 K were used to calculate the maximum cisplatin loading (nGCMC) and the MOF–CPT interaction energies (UMOF-CPT), which represents the affinity of cisplatin to the framework hence provides a semi-quantitative measure of release rates at body temperature. The GCMC simulations were carried out inside the static frameworks using the multipurpose simulation code (MuSiC).39 MOF–CPT and CPT–CPT Coulombic interactions were calculated using Ewald summation40 and the Wolf summation method, respectively.40 All non-Coulombic interactions were calculated using the Lennard-Jones 12–6 potential. 10 × 106 iterations of insertion, deletion, translation, and rotation moves were used for each pressure point, and the first 40% of microstates were ignored to ensure that the microstates sampled are at equilibrium.

Most of our methods (with the exception of GCMC, GASP and ab initio simulations) involve molecular dynamics (MD). When screening all frameworks, MD simulations were used to pre-equilibrate all structures and characterize the energy barriers of intracell and intercell cisplatin movement. For the UiO-66 variants, for which experimental data are available, MD simulations were also used to calculate cisplatin diffusivity with and without water to determine the influence of water during encapsulation and to compare forcefield based and ab initio based MOF–CPT interaction energies.

All molecular dynamics (MD) simulations were carried out using the Groningen Machine for Chemical Simulations (GROMACS, 2019.2) code41–46 keeping the MOF frameworks flexible. The MOF topologies were generated using OBGMX, which models all interactions using the Universal forcefield.37,47 No additional constraints were imposed on the framework. Long-ranged electrostatic interactions were calculated using the particle mesh Ewald algorithm.48,49 Long range dispersion corrections were applied to the energy and pressure terms. Prior to all MD simulations, the structures were pre-equilibrated using a steepest decent energy minimization followed by a 100 ps NVT equilibration. Separate temperature coupling baths were used for cisplatin and the framework, along with a modified Berendsen coupling method.50 Depending on the ensemble used during production simulations, a further 100 ps NPT equilibration was carried out using a Parrinello–Rahman barostat. During the production simulations, temperature and (depending on the ensemble) pressure coupling was achieved using a Nosé–Hoover thermostat and a Parrinello–Rahman barostat.51,52

Since GCMC simulations cannot easily be used to determine the impact of water on cisplatin loading from solution, we used Bennett’s Acceptance Ratio (BAR) simulations to model cisplatin + water uptake in UiO-66/UiO-66(NH2). In the BAR simulations, a single cisplatin molecule was initially placed in either the octahedral or tetrahedral pores of UiO-66 and UiO-66(NH2) [Fig. 3(a)]. We then randomly inserted NSOL water molecules (NSOL being the number of water molecules in one of a selection of simulations, over which the water content inside the pores is varied) and equilibrated the system at 298 K (representing encapsulation temperature). The NPT ensemble was used in these simulations so that the framework could expand if necessary, at higher water loadings. At equilibrium, there is no net transfer of cisplatin between system (a) (cisplatin in bulk water) and system (b) (MOF + cisplatin + water). Hence, the net change in Gibbs free energy when moving a cisplatin molecule between (a) and (b) is nil, as shown conceptually in Fig. 3(b). By gradually decoupling a single cisplatin molecule from systems (a) and (b), it is possible to use the BAR method to calculate the free energy change associated with cisplatin uptake from solution using the Hess cycle in Fig. 3(c). To decouple cisplatin (corresponding to turning off all non-bonded interactions between cisplatin and the rest of the system), a series of configurations were generated between the fully coupled and decoupled states. The BAR algorithm53 was then used to calculate the change in Gibbs free energy between the two states: ΔG=01dλdHλdλ (where H is the parameterized Hamiltonian and λ refers to the position along an integration path between states 0 and 1). In this work, state 0 corresponds to a cisplatin molecule that fully interacts with the system and state 1 corresponds to no interactions between cisplatin and the system.54 To prevent molecules of opposite charge becoming trapped, the Coulombic interactions were decoupled before the LJ interactions. In total, 20 configurations were used to decouple the non-bonded interactions with varying λ between 0 and 1 in intervals of 0.1 to ensure each stage of decoupling shares part of the same phase space.

Diffusion is a direct measure of cisplatin release rates. We calculated cisplatin diffusivity based on the gradient of mean-square displacement plots of cisplatin in solvated and non-solvated UiO-66 and UiO-66(NH2) at elevated temperatures of 398 K (to capture interpore movement within reasonable simulation times). This allowed us to determine the influence of water on release rates. We extended these simulations to include wider-pore variants of UiO-66 (UiO-67 and UiO-68) to check whether the influence of water is the same when diffusion energy barriers are dominated by non-bonded as opposed to bonded interaction energies.

Helmholtz total free energy barriers and bonded/non-bonded potential energy barriers associated with cisplatin diffusion are an indirect measure of cisplatin release rates, and they were calculated using umbrella sampling and steered-MD (SMD) simulations, respectively. Since the bulk crystal would confine the volume of a unit cell, these simulations were conducted using the NVT ensemble. A single cisplatin molecule was placed in the center of the largest pore of the crystal structure, and the system was relaxed using a steepest descents energy minimization and 100 ps NVT pre-equilibration. The temperature was regulated using a modified Berendsen coupling method.50 Cisplatin was then dragged through the path of least resistance (assumed to be the 1D line in which cisplatin passes through the fewest and largest pore windows) using the GROMACs pull code.56 In this work, “the paths of least resistance” are 1D routes that solely pass through the 6 MRs of ZIF-8, 8 MRs of ZIF-11, or through the kno channels of ZIFs 68, 70, 78, 79, and 82. The position of cisplatin relative to the path was defined relative to its center of mass. The simulation time was adjusted to ensure that cisplatin moves from one pore to the same position in a neighboring pore. For all simulations, a pull rate of 0.01 nm/ps was used with a time step of 0.001 ps. A spring constant of 1000 kJ/mol nm2 was used for all ZIFs with the exception of ZIF-11 for which a tighter constraint potential was needed to adequately sample cisplatin movement at the window. Sampling windows for the umbrella simulations were collected at least every 0.1 nm along the simulation path, though more sampling windows were taken as necessary in the high energy regions of the pore windows. The sampling configurations were not pre-equilibrated but each was subjected to 10 ns of MD simulations in which cisplatin is constrained by the harmonic potential.57 The weighted histogram method was used to obtain the free energy profiles along the pull direction.58 At least six pull and umbrella sampling repeats were carried out for each MOF [more repeats were used in systems where the potential mean force (PMF) curves were less consistent]. During both the pull and umbrella simulations, separate temperature coupling of CPT and the framework was achieved using a Nosé–Hoover thermostat.51,52 Both the time-dependent potential mean force (PMF) curves and associated energy histograms are shown in Sec. S6 in the supplementary material.

The intrinsic flexibility of the MOF frameworks were determined using the GASP (geometric analysis and simulation of polyhedra) software.59 GASP was used in combination with the aforementioned SMD simulations to determine the potential energy barrier associated with cisplatin moving through the frameworks, to identify the bonded and non-bonded contributions to this energy barrier, and to determine whether cisplatin can pass through the pore windows without severely distorting the framework covalent bonding. In this approach, the local atomic geometry of the input structure is constrained by a series of superimposed templates while the simulation cell parameters are varied. This identifies a range of cell parameters, termed a “flexibility window,” within which the global geometry can vary without alterations to the local geometry. The variants produced using GASP were then analyzed using PoreBlazer60 and ab initio static energy calculations (detailed information on the ab initio and GASP simulations can be found Secs. S1 and S2 in the supplementary material).

FIG. 4.

Equilibrium MD simulations at 398 K in UiO-66, UiO-67, and UiO-68 [with BDC and BDC(NH2) ligands]. (a) Cisplatin diffusivity in solvated MOFs (1CPT/o-pore) compared to non-solvated MOFs (1CPT/unit cell). (b) Cisplatin diffusivity in solvated frameworks as a function of ∆UMOF-CPT,SMD (@298 K) calculated in non-solvated frameworks using SMD simulations.

FIG. 4.

Equilibrium MD simulations at 398 K in UiO-66, UiO-67, and UiO-68 [with BDC and BDC(NH2) ligands]. (a) Cisplatin diffusivity in solvated MOFs (1CPT/o-pore) compared to non-solvated MOFs (1CPT/unit cell). (b) Cisplatin diffusivity in solvated frameworks as a function of ∆UMOF-CPT,SMD (@298 K) calculated in non-solvated frameworks using SMD simulations.

Close modal

To understand which of our simulation methods are appropriate for screening pH-sensitive MOFs for cisplatin delivery, we tested them on UiO-66 and UiO-66(NH2), for which experimental data are available.15 To date, most of MOF–drug simulation papers utilize pure component GCMC simulations to determine the drug molecule uptake (nGCMC).31 The maximum uptake of cisplatin in UiO-66 and UiO-66(NH2) using pure-component GCMC simulations is 0.67:1 and 0.50:1 CPT:Zr, respectively (which is in line with their pore volumes of 0.46 to 0.42 cm3/g60,61). The GCMC loading in UiO-66(NH2) matches the experimental loading in UiO-66(NH2) via conjugation (i.e., diffusion and covalent bonding of functionalized, pre-activated cisplatin molecules into the MOF from solution15). As an efficient, widely used31 preliminary screening tool, we therefore used GCMC simulations to obtain a maximum loading, as our first screening criteria to get information about the relative ranking of different materials.

However, it is important to note that the cisplatin loadings obtained using pure component GCMC simulations are the maximum possible loading, and whether the maximum loading is achieved experimentally will strongly depend on the procedure used. For example, experimental loading via encapsulation (i.e., diffusion of cisplatin from solution into the framework) yields uptakes of 0.06:1 CPT:Zr in both UiO-66 and UiO-66(NH2)15 (which is tenfold lower than nGCMC). Encapsulation cannot achieve the maximum loading because (as shown by Bennett’s acceptance ratio simulations in Sec. S3 of the supplementary material) cisplatin favors its solvation shell to an empty framework. As water enters the octahedral pores, cisplatin + solvation shell + MOF interactions are more favorable than cisplatin in bulk water. However, the tetrahedral pores are too small to occupy cisplatin and its solvation shell; hence, at equilibrium the theoretical cisplatin loading in UiO-66 and UiO-66(NH2) is 1 CPT per octahedral pore (0.17:1 CPT:Zr), which is significantly closer to the encapsulation loading 0.06:1 CPT:Zr15 than nGCMC. Although our BAR simulations of cisplatin adsorption in the presence of water provide thermodynamic insight into the experiments, we did not account for water when screening pH sensitive ZIFs because, as a preliminary screening tool, BAR is too computationally demanding. Furthermore, GCMC simulations of water and cisplatin mixtures are inefficient because packing inside the pores would lead to high rejection probabilities. To account for the lack of water, we regard wider cavities [>8 Å, the size of tetrahedral pores in UiO-66/UiO-66(NH2)] to be beneficial for uptake via encapsulation.

Experimental cisplatin release rates are fourfold higher in UiO-66 compared to UiO-66(NH2), likely due to enhanced amine–cisplatin interaction energies.15 While we cannot not directly compute the release rate easily, there are a number of other computational indicators of release rates, including (a) cisplatin diffusivity through the framework (DCPT), (b) cisplatin diffusion energy barriers, and (c) MOF–CPT interaction energies (UMOF-CPT). To decide which indicators to use for screening pH sensitive ZIFs, we calculated each indicator for cisplatin release from UiO-66 and UiO-66(NH2) and compared the results to experimental release rates.

Cisplatin diffusivity is a direct measure of release rates and was calculated using equilibrium MD simulations (both with and without water) at experimental temperatures. As shown in Sec. S4 in the supplementary material, regardless of the water loading, cisplatin diffusivity is slightly lower in UiO-66(NH2) compared to UiO-66 due to amine–water–cisplatin and amine–cisplatin hydrogen bonding. However, limiting sub-diffusion was witnessed in both frameworks due to high energy barriers associated with cisplatin movement through the pore windows. In theory, diffusivities calculated from equilibrium MD provide a direct measure of release rates; however, we deemed them inappropriate for screening ZIFs because interpore movement could not be captured within reasonable simulation times.

The necessity to include water when calculating the energy barrier for diffusion therefore depends on whether water alters trends in cisplatin diffusivity. To determine the impact of water on cisplatin diffusivity, we performed equilibrium MD simulations with water (@1CPT/octahedral pore and the optimal water loading described in Sec. S4 of the supplementary material) and without water (@1CPT/unit cell to prevent cisplatin molecules grouping together) at the elevated temperature of 398 K to capture interpore cisplatin diffusion in UiO-66 and UiO-66(NH2). As cisplatin could potentially be carried away in a “water channel,” if water disrupts MOF–CPT interactions in an effectively uninhibited space, we also extended these simulations to include isoreticular variants of UiO-66 with larger pores [UiO-67, UiO-67(NH2), UiO-68, and UiO-68(NH2)62]. As shown in Fig. 4(a), while the presence of water considerably slows down diffusion, the ranking of the diffusivities remains the same regardless of the presence of water. Furthermore, the non-bonded MOF–CPT diffusion energy barriers of cisplatin movement through each non-solvated UiO framework using SMD simulations (ΔUMOF-CPT,SMD, which is described in detail below) is strongly correlated with the cisplatin diffusivities calculated in solvated frameworks, as seen in Fig. 4(b). Therefore, when screening pH sensitive ZIFs, we compared diffusion energy barriers in non-solvated frameworks as semi-quantitative measurements of release rates.

TABLE I.

Summative ranking of screening criteria used to test pH sensitive ZIFs as drug delivery carriers. The overall MOF rank will be calculated as score × ranking position, where a higher rank corresponds to a better performing material.

Screening criteriaDescription/SimulationRanking scoreReason
nGCMC Maximum cisplatin loading/pure-CPT GCMC Adequately captures the maximum potential uptake. 
UMOF-CPT MOF–CPT interaction energy/pure-CPT GCMC CPT affinity to the framework is important for release/uptake to/from solution. Equilibrium depends on CPT affinity to the framework vs water. Values are subject to CPT and water loading. 
ΔAumbrella Helmholtz total free energy barrier/umbrella sampling Paints the complete picture (bonded + non-bonded free energies), but it cannot reproduce experimental release rates without high levels of uncertainty. 
ΔUbonded Bonded potential energy barrier/SMD or GASP + ab initio a Comparison of SMD to the GASP flexibility window provides an immediate indication as to whether CPT can move through the framework without distortions. 
ΔUMOF-CPT Non-bonded potential energy barrier/SMD Adequate indicator for capturing experimental release rates. 
Screening criteriaDescription/SimulationRanking scoreReason
nGCMC Maximum cisplatin loading/pure-CPT GCMC Adequately captures the maximum potential uptake. 
UMOF-CPT MOF–CPT interaction energy/pure-CPT GCMC CPT affinity to the framework is important for release/uptake to/from solution. Equilibrium depends on CPT affinity to the framework vs water. Values are subject to CPT and water loading. 
ΔAumbrella Helmholtz total free energy barrier/umbrella sampling Paints the complete picture (bonded + non-bonded free energies), but it cannot reproduce experimental release rates without high levels of uncertainty. 
ΔUbonded Bonded potential energy barrier/SMD or GASP + ab initio a Comparison of SMD to the GASP flexibility window provides an immediate indication as to whether CPT can move through the framework without distortions. 
ΔUMOF-CPT Non-bonded potential energy barrier/SMD Adequate indicator for capturing experimental release rates. 
a

Only used as a screening criterion if the uptake method is from solution into a pre-synthesized framework. If ΔUbonded results in breaking framework bonds, we immediately remove this framework from our selection of potential candidates.

To determine the energy barrier that a CPT molecule is experiencing when moving through the framework, we calculated the Helmholtz total free energy barrier consisting of both bonded (resulting from the framework) and non-bonded contributions (resulting from the CPT/MOF interactions). For this, we used umbrella sampling simulations in which a single cisplatin molecule was moved along the path of least resistance through otherwise empty UiO-66/UiO-66(NH2). As shown in detail in Sec. S6 in the supplementary material, the calculated Helmholtz total free energy barriers are ΔAumbrella = 22 kJ/mol for UiO-66 and UiO-66(NH2). These energy barriers are dominated by bonded energies in the framework and compare reasonably well to the π-flip energy barrier of empty UiO-66 calculated by Kolokolov et al. who used nuclear magnetic resonance experiments (∼30 kJ/mol63). The Helmholtz total free energy barrier was therefore used as a screening criterion.

As an alternative to the Helmholtz total free energy barrier, we also calculated the individual bonded and non-bonded contributions to the total potential energy barrier separately to see which term correlates with experimentally seen release rates. For this purpose, we performed SMD simulations [again pulling a single cisplatin molecule through the “path of least resistance” in non-solvated UiO-66/UiO-66(NH2)]. The SMD bonded potential energy barrier was higher in UiO-66 (ΔUbonded,SMD = 406 kJ/mol) compared to UiO-66(NH2) (ΔUbonded,SMD = 184 kJ/mol). We then compared the energy barrier corresponding to maximum framework deformation from SMD simulations to the bonded potential energy barrier, ΔUbonded,GASP, calculated as a function of framework flexibility using a combination of GASP (to determine structures within the flexibility window) and static energy calculation using ab initio simulations. As the values are very similar to the results from the SMD simulations [ΔUbonded,GASP = 443 kJ/mol for UiO-66 and ΔUbonded,GASP = 157 kJ/mol for UiO-66 (NH2)], it showed that our choice of forcefield and charges is adequate. Furthermore, the extent of deformation from SMD simulations was found to be within the flexibility window of both UiO-66 and UiO-66(NH2), meaning cisplatin can move through the pore windows without breaking framework bonds. The higher bonded potential energy barrier for UiO-66 agrees with trends observed by Pakhira et al. who argued that the energy barrier to rotate polar ligands is lower due to electrostatic repulsion.64 The combination of SMD, ab initio, and GASP simulations provides information whether cisplatin can pass through the framework without the need to break framework bonds. This was used as a screening test on pH sensitive ZIFs that show particularly high energy barriers due to small pore window sizes.

FIG. 5.

Maximum loading of cisplatin, nGCMC, as a function of (a) accessible pore volume and (b) MOF–CPT interaction energy at maximum loading, UMOF-CPT. Symbols: UiO-66 (×), UiO-66(NH2) (△), ZIF-8 (□), ZIF-11 (○), ZIF-68 (+), ZIF-70 (∗), ZIF-78 (◇), ZIF-79 (▬), and ZIF-82 ().

FIG. 5.

Maximum loading of cisplatin, nGCMC, as a function of (a) accessible pore volume and (b) MOF–CPT interaction energy at maximum loading, UMOF-CPT. Symbols: UiO-66 (×), UiO-66(NH2) (△), ZIF-8 (□), ZIF-11 (○), ZIF-68 (+), ZIF-70 (∗), ZIF-78 (◇), ZIF-79 (▬), and ZIF-82 ().

Close modal

While the bonded potential energy barriers can tell us if it is feasible for cisplatin to move through a framework, they do not compare well to trends observed experimentally for the release rates (where UiO-66 shows a higher release rate) most likely due to the fact that the framework undoubtedly deforms when solvated. Instead, the MOF–CPT non-bonded potential energy barriers provide a more meaningful, indirect measure of release rates. Using SMD simulations, the non-bonded potential energy barriers were calculated to be ΔUMOF-CPT,SMD = 55 kJ/mol for UiO-66 and ΔUMOF-CPT,SMD = 67 kJ/mol for UiO-66(NH2), respectively. Stronger interaction energies between the polar arms of cisplatin and the amine groups in UiO-66(NH2) result in higher non-bonded MOF–CPT energy barriers compared to UiO-66. Unlike the bonded contribution to the total energy barrier, this is in agreement with the higher release rates from UiO-66.15 Therefore, the non-bonded potential energy barrier was included as a screening criterion. Table I summarizes our choice of screening criteria to select pH sensitive ZIFs for cisplatin drug delivery. Other factors that influenced our decisions when comparing ZIFs included the pore volume (uptake via encapsulation will require the pore volume to occupy cisplatin and its solvation shell), the limiting diameter along the path of least resistance (if it is too small, both release and uptake will be hindered), and the presence of electronegative groups (hence the ability to form MOF–cisplatin or MOF–cisplatin–water hydrogen bonds and the possibility to incorporate cisplatin experimentally via conjugation).

We used GCMC simulations to calculate the maximum loading (nGCMC) and MOF–CPT energy (UMOF-CPT) associated with cisplatin uptake in the pH sensitive ZIFs. The maximum uptake should be as high as possible to increase the dosage to the tumor for a given amount of MOF. The adsorption energy acts as a driving force for encapsulation, and a high absolute value improves retention and prevents premature release of cytotoxic drugs in healthy areas of the patient. Figure 5(a) shows that the maximum loading is dominated by the accessible pore volume since the more cisplatin molecules can pack into larger void spaces (for calculated and experimental structural properties, see Sec. S7 in the supplementary material). ZIF-70 has the highest uptake because its small nIm + Im ligands result in larger pore volume. Figure 5(b) shows that the adsorption energy decreases with maximum cisplatin loading because a lower density of framework atoms results in larger pore space and weaker interaction energies. For example, ZIF-78 has the strongest average adsorption energy (hence highest affinity for cisplatin) and lowest uptake attributed to its polar ligands in combination with small pore sizes. However, since UMOF-CPT decreases with higher cisplatin loading and a lower density of framework atoms, there is a clear trade-off between nGCMC and UMOF-CPT.

Adsorption isotherms and energy histograms (i.e., the fraction of GCMC configurations that correspond to a particular adsorption energy) provide further insight into the adsorption mechanism. Here, it is specifically of interest whether the cisplatin loading is associated with strong MOF–CPT energies or a pore-filling mechanism. The latter is less useful for drug delivery since the encapsulation relies on cisplatin energetically favoring the pore space as opposed to solution.

Cisplatin adsorption in ZIF-8 exhibits a type 1 isotherm [Fig. 6(a)] due to a pore filling mechanism with weak MOF–CPT interaction energies of UMOF-CPT = −36 kJ/mol in its hydrophobic 11 Å pore cavities [Fig. 6(b)]. In ZIF-11, MOF–CPT interaction energies are comparatively strong at the cavity edges where larger bIm ligands form rings [UMOF-CPT < −90 kJ/mol Fig. 6(b) and Sec. S8 in the supplementary material]. These high energy sites are the first to become occupied at low pressures and followed pore filling, resulting in the steps in Fig. 6(a) and dual peaks in Fig. 6(b).

FIG. 6.

(a) Cisplatin adsorption isotherms (error bars show the standard deviation from three repeat simulations) and (b) UMOF-CPT histograms at 310 K (i.e., the statistical distribution of configurations displaying a particular MOF–CPT interaction energy within 1 kJ/mol bin sizes). Open triangles/dashed line = ZIF-8. Filled triangles/solid line = ZIF-11.

FIG. 6.

(a) Cisplatin adsorption isotherms (error bars show the standard deviation from three repeat simulations) and (b) UMOF-CPT histograms at 310 K (i.e., the statistical distribution of configurations displaying a particular MOF–CPT interaction energy within 1 kJ/mol bin sizes). Open triangles/dashed line = ZIF-8. Filled triangles/solid line = ZIF-11.

Close modal

ZIFs-68/70/78/79/82 consist of large kno cages surrounded by channels of smaller gme-hpr cages that are parallel to the c-axis (diameter ≈ 3.6 Å65). The functional groups of the ligands range from polar (ZIF-78 and ZIF-82) to non-polar (ZIF-68, 70, and 70) pointing into the kno cages, resulting in a variation in pore aperture from 3.8–13.1 Å and pore diameter from 7.1–15.9 Å.28 Adsorption isotherms and interaction energies [Figs. 7(a) and 7(b) and also Sec. S7 in the supplementary material] show a small step in the uptake at low pressures as cisplatin occupies the high energy adsorption sites (typically UMOF-CPT = −150 kJ/mol in the smaller gme pores—except in ZIF-78 where it initially fills the high energy sites next to nbIm ligands in its kno pores). At higher pressures, there is a larger step in the isotherms due to adsorption in the low energy sites (UMOF-CPT < −100 kJ/mol typically in the larger kno cages). This “two step” behavior is similar to that seen by Van der Perre who performed adsorption experiments of polar guest molecules in ZIF-68.65 

FIG. 7.

(a) Adsorption isotherms (error bars show the standard deviation from three repeat simulations) and (b) UMOF-CPT histograms in ZIF-68, ZIF-70, ZIF-78, ZIF-79, and ZIF-82 at 310 K (i.e., the statistical distribution of configurations displaying a particular MOF–CPT interaction energy within 1 kJ/mol bin sizes). Corresponding energy maps and energies as a function of pressure are shown in the supplementary material.

FIG. 7.

(a) Adsorption isotherms (error bars show the standard deviation from three repeat simulations) and (b) UMOF-CPT histograms in ZIF-68, ZIF-70, ZIF-78, ZIF-79, and ZIF-82 at 310 K (i.e., the statistical distribution of configurations displaying a particular MOF–CPT interaction energy within 1 kJ/mol bin sizes). Corresponding energy maps and energies as a function of pressure are shown in the supplementary material.

Close modal

The benefits of large cisplatin uptakes are forfeit if the payload is released prematurely in the healthy tissues of chemotherapy patients. While MOF–CPT interaction energies, UMOF-CPT, are indicative of cisplatin’s affinity to the framework over solution, they cannot be directly related to release rates. As discussed above, better criteria are the energy barriers associated with moving a single cisplatin molecule down the path of least resistance (i.e., a 1D line in which cisplatin moves through the widest pores/channels with the fewest pore windows in each ZIF).

As a first step, we compared the extent of framework deformation (from SMD simulations, when cisplatin moves through the pore-limiting diameter) to the framework’s flexibility window (from GASP simulations) for each ZIF to determine whether cisplatin can fit through the narrowest parts of the path of least resistance without severely distorting the framework’s bonding. As shown in Sec. S5 in the supplementary material, interpore movement of cisplatin in ZIF-11 requires the pore windows to open beyond the extent of its flexibility window (i.e., pore-hopping is likely to distort the localized bonding). Therefore, in the case where cisplatin uptake is from solution into a pre-synthesized framework, we remove the ZIF-11 framework from our list of potential drug delivery carriers as encapsulation would be limited (that said, if cisplatin could be captured in the pores of ZIF-11 as the crystal grows, it would be an ideal cisplatin drug delivery carrier). All other MOFs allow interpore cisplatin movement without severely distorting localized bonding in the framework.

Figure 8(a) shows the Helmholtz total free energy barriers associated with cisplatin movement along the path of least resistance in each framework as a function of limiting pore diameter. The Helmholtz total free energy barriers (∆Aumbrella) increase rapidly as the framework pore limiting diameter (PLD) decreases below 4 Å. This rapid increase is due to the bonded energy penalty associated with displacing framework atoms to allow interpore cisplatin movement. ZIFs 8 and 11 have the smallest PLD out of the frameworks screened (<3 Å) and hence the highest Helmholtz total free energy barrier and the best retention properties. The Helmholtz free energy barrier is also relatively constant in MOFs with PLD ≥ 4 Å at ∼20 kJ/mol. UiO-66 (which has been shown experimentally to be a promising drug delivery carrier since its rate of release is ideal for reducing tumor cell viability15) also has a Helmholtz total free energy barrier of ∆Aumbrella = 20 kJ/mol. The pore windows in UiO-66 are slightly smaller than the collision diameter of cisplatin (3.6 Å); therefore, the Helmholtz total free energy barrier in UiO-66 is governed by the bonded contribution to the energy barrier (i.e., the energy penalty of displacing framework atoms when cisplatin moves through its small pore windows). Since the kno channels are generally larger than the pore windows in UiO-66, it would be expected that they pose even lower Helmholtz total free energy barriers. However, this is not the case as the non-bonded MOF–CPT interaction energy barriers (∆UMOF-CPT) increase with PLD > 4 Å (as discussed below). In other words, larger pore windows have lower bonded energy barriers, but larger pores have higher non-bonded energy barriers, which result in similar overall energy barriers (∆Aumbrella). This is important as larger pores would enable enhanced encapsulation while maintaining similar levels of cisplatin retention as UiO-66.

FIG. 8.

Difference in the Helmholtz total free energy barrier, ∆Aumbrella, associated with cisplatin movement through the path of least resistance in the framework at 310 K as a function of pore limiting diameter (PLD). Symbols: UiO-66 (×), UiO-66(NH2) (△), ZIF-8 (□), ZIF-11 (○), ZIF-68 (+), ZIF-70 (∗), ZIF-78 (◇), ZIF-79 (▬), and ZIF-82 ().

FIG. 8.

Difference in the Helmholtz total free energy barrier, ∆Aumbrella, associated with cisplatin movement through the path of least resistance in the framework at 310 K as a function of pore limiting diameter (PLD). Symbols: UiO-66 (×), UiO-66(NH2) (△), ZIF-8 (□), ZIF-11 (○), ZIF-68 (+), ZIF-70 (∗), ZIF-78 (◇), ZIF-79 (▬), and ZIF-82 ().

Close modal
FIG. 9.

Maximum cisplatin uptake as a function of the non-bonded MOF-CPT potential energy barrier (ΔUMOF-CPT) (red: PLD > 4 Å, blue: PLD < 4 Å). Symbols: UiO-66 (×), UiO-66(NH2) (△), ZIF-8 (□), ZIF-11 (○), ZIF-68 (+), ZIF-70 (∗), ZIF-78 (◇), ZIF-79 (▬), and ZIF-82 ().

FIG. 9.

Maximum cisplatin uptake as a function of the non-bonded MOF-CPT potential energy barrier (ΔUMOF-CPT) (red: PLD > 4 Å, blue: PLD < 4 Å). Symbols: UiO-66 (×), UiO-66(NH2) (△), ZIF-8 (□), ZIF-11 (○), ZIF-68 (+), ZIF-70 (∗), ZIF-78 (◇), ZIF-79 (▬), and ZIF-82 ().

Close modal

Figure 9 shows the non-bonded MOF–CPT interaction potential energy barriers (∆UMOF-CPT), which provided a direct correlation to the release rates for UiO-66 and UiO-66(NH2). In general, the MOF–CPT interaction energies (UMOF-CPT) are enhanced in the presence of larger ligands (e.g., in the kno channels of ZIF-79) or ligands with polar functional groups (e.g., in the kno channels of ZIF-78 and ZIF-82, Fig. 7). However, to increase the potential energy barrier, there needs to be a large variation in UMOF-CPT i.e., large ΔUMOF-CPT, caused by strong adsorption sites spread out between weak adsorption sites. As seen in Sec. S9 in the supplementary material, UMOF-CPT is uniformly distributed in ZIF-78 and ZIF-79, and so despite favorable interaction energies, the non-bonded potential energy barriers (ΔUMOF-CPT) are low, meaning cisplatin can move relatively freely in these frameworks. In ZIF-82 however, there are strong interaction energies in the vicinity of the C≡N groups and weak interaction energies elsewhere, which enhances ΔUMOF-CPT. The phenomena of large steps between favorable adsorption sites (i.e., large variations in energy) increasing the energy barrier was also seen by Amirjalayer and Schmid who simulated benzene diffusion through variants of IRMOF-1.66 

Overall, we draw the following recommendations. Depending on the method of loading cisplatin into the MOF, a release rate hindered by the PLD may or may not be a drawback to using that MOF as a drug delivery carrier. For example, windows connecting the cavities in ZIF-8 (pore limiting diameter, PLD = 2.9 Å) and ZIF-11 (PLD = 2 Å) are very narrow compared to the diameter of CPT (d ≈ 5 Å) and even smaller than in UiO-66 (PLD = 3.6 Å) where encapsulation takes more than 48 h.15 Combined with weak UMOF-CPT, there is little driving force and a large free energy barrier (∆Aumbrella:ZIF-8 = 42 kJ/mol and ∆Aumbrella:ZIF-11 = 41 kJ/mol) for cisplatin to enter the cavities of these MOFs. We therefore conclude that these structures will not exhibit good cisplatin loading by encapsulation. However, if there were a method of growing the crystals around the drug molecules, ZIF-11 would make an excellent drug delivery carrier due to its retention properties. In situations where cisplatin uptake is hindered by encapsulation, it may be beneficial to use ZIFs with larger kno pores that can display a similar rate of retention (∆Aumbrella) as in UiO-66, which has been shown to have ideal release rates for killing cancer cells.15 To further enhance encapsulation and to prevent premature release rates caused by water displacement, ΔUMOF-CPT can be enhanced by the addition of polar groups—as in the case of ZIF-82. As a next step, the intrapore and interpore diffusivities would need to be calculated in the presence of water to fully determine the influence of the solvation shell prior to experimental screening. However, such calculations are computationally expensive and beyond the scope of this paper, which is to screen ZIFs and identify promising materials that can then be further investigated in more detail.

TABLE II.

Summary of cisplatin uptake and release properties in the screened MOFs. Note that the highest rank corresponds to the best performing MOF and that the score is calculated as rank × ranking score from Table I.

Screening criterianGCMC (NCPT3) rank (score)UMOF-CPT (kJ/mol) rank (score)ΔAumbrella (kJ/mol) rank (score)ΔUMOF-CPT (kJ/mol) rank (score)Total score
UiO66 0.0018 −45 22 55 32 
5 (15) 2 (4) 5 (5) 2 (8) 
UiO66(NH20.0014 −61 22 67 48 
3 (9) 7 (14) 5 (5) 5 (20) 
ZIF8 0.0025 −36 42 75 59 
8 (24) 1 (2) 9 (9) 6 (24) 
ZIF11a 0.0018 −90 41 120 69 
5 (15) 5 (10) 8 (8) 9 (36) 
ZIF68 0.0018 −150 23 55 41 
5 (15) 6 (12) 6 (6) 2 (8) 
ZIF70 0.0027 −150 20 80 64 
9 (27) 3 (6) 3 (3) 7 (28) 
ZIF78 0.0011 −150 30 55 36 
1 (3) 9 (18) 7 (7) 2 (8) 
ZIF79 0.0014 −150 16 57 43 
3 (9) 8 (16) 2 (2) 4 (16) 
ZIF82 0.0021 −51 13 90 62 
7 (21) 4 (8) 1 (1) 8 (32) 
Screening criterianGCMC (NCPT3) rank (score)UMOF-CPT (kJ/mol) rank (score)ΔAumbrella (kJ/mol) rank (score)ΔUMOF-CPT (kJ/mol) rank (score)Total score
UiO66 0.0018 −45 22 55 32 
5 (15) 2 (4) 5 (5) 2 (8) 
UiO66(NH20.0014 −61 22 67 48 
3 (9) 7 (14) 5 (5) 5 (20) 
ZIF8 0.0025 −36 42 75 59 
8 (24) 1 (2) 9 (9) 6 (24) 
ZIF11a 0.0018 −90 41 120 69 
5 (15) 5 (10) 8 (8) 9 (36) 
ZIF68 0.0018 −150 23 55 41 
5 (15) 6 (12) 6 (6) 2 (8) 
ZIF70 0.0027 −150 20 80 64 
9 (27) 3 (6) 3 (3) 7 (28) 
ZIF78 0.0011 −150 30 55 36 
1 (3) 9 (18) 7 (7) 2 (8) 
ZIF79 0.0014 −150 16 57 43 
3 (9) 8 (16) 2 (2) 4 (16) 
ZIF82 0.0021 −51 13 90 62 
7 (21) 4 (8) 1 (1) 8 (32) 
a

Rejected from screening process if uptake is from solution into pre-synthesized framework due to large bonded energy barrier and bond breaking.

Table II summarizes the findings for each framework screened and emphasizes the importance of including a variety of properties, even in an initial screening process, to identify promising materials. It should be noted that our total scores in Table II depend on the ranking scores defined in Table I, but nevertheless they show that it is not enough to screen based on a single criterion alone. As an example, many MOF-drug screening papers only focus on the drug molecule loading. In this work, ZIF-8 exhibits the second-highest loading attributed to its large pore volume. However, when we consider ZIF-8’s weak interaction energies with cisplatin (due to the hydrophobic nature of ZIF-8, cisplatin will favor bulk water), mediocre non-bonded cisplatin diffusion energy barrier (attributed to the lack of adsorption hotspots), and high bonded cisplatin diffusion energy barrier (caused by small pore windows which will hinder cisplatin uptake from solution), this popular framework becomes less promising for cisplatin drug delivery.

Overall, ZIFs-82, 70, and 11 are most promising as they exhibit good uptake (due to their large pore volumes) and large diffusion energy barriers. However, if encapsulation is carried out from solution into a pre-synthesized framework, we omit ZIF-11 from this list because cisplatin would not be able to enter the pores without breaking framework bonds. Compared to ZIF-11, the window diameters of the kno channels in ZIFs-82/70 are larger, and so there would be less diffusion limitation during encapsulation. Although the gme pore windows are smaller (3.6 Å) and will hinder cisplatin diffusion, strong MOF-CPT interactions in the polar cavities will act as a driving force to encapsulate cisplatin. This effect was witnessed by Van der Pierre et al., who occupied the gme pores of ZIF-68 with TMB molecules, which are larger than cisplatin and contain similar functional groups.65 Compared to ZIF-70, ZIF-82 also has polar groups in the kno channel creating strong, dispersed cisplatin adsorption hotspots, so at this stage, it would appear to be the most promising framework for enhanced cisplatin encapsulation.

In this work, we assessed a variety of simulation methods for their suitability to provide information about cisplatin uptake and release by comparing the simulation results to available experimental results of cisplatin (CPT) uptake and release rates for UiO-66 and UiO-66(NH2). For screening applications, it would be beneficial not to include water in the simulations to reduce the computational effort. Pure-component GCMC simulations adequately capture the maximum reported experimental loading of CPT. However, MD simulations in the presence of water reveal favorable interaction energies between cisplatin and its solvation shell, which might result in this theoretical loading not being achieved experimentally unless the pore size is large enough to occupy cisplatin with its solvation shell or the MOF–cisplatin interaction energies are more favorable than the cisplatin-shell energy. For the release of CPT, the diffusion barrier is of importance and our non-equilibrium MD simulations show that the diffusion energy barriers are the same regardless of water. We therefore conducted our screening simulations without water. Apart from the maximum uptake, we identified four other important properties to take into account when screening porous materials for drug delivery, all of which provide indirect information about the retention and release rate. The MOF–CPT adsorption energies describe the affinity of cisplatin to the framework as opposed to bulk solution (i.e., the likelihood of it being displaced by water). Bonded potential energy barriers can be used to determine whether cisplatin can physically pass through the MOF without breaking framework bonds (hence whether encapsulation from solution into a pre-synthesized framework is possible). Non-bonded energy barriers provide a measure for the distribution of adsorption hotspots and hence the difficulty for cisplatin to pass between them and move through the framework. Finally, the Helmholtz free energy barrier describes this overall contribution of bonded and non-bonded diffusion barriers.

We then screened a selection of biocompatible pH sensitive ZIFs for their ability to achieve a high maximum loading of cisplatin as well as a high energy barrier to increase the payload reaching the tumor in chemotherapy patients and decrease release rates in healthy parts of the body. The maximum uptake is predominately governed by the accessible pore volume. The pore-hop energy barrier (i.e., the energy barrier associated with cisplatin movement between pores, which is inversely proportional to the drug molecule release rates from the framework) can be enhanced by decreasing the pore window diameter at a compromise that this will reduce accessibility of cisplatin entering the framework during encapsulation. Instead, similar energy barriers can be achieved without compromising encapsulation rates in frameworks with large open pores but with enhanced MOF–cisplatin interaction energies.

Ultimately, this work has highlighted the structural properties to aim for when designing an ideal cisplatin drug delivery carrier. From the MOFs screened in this work, ZIF-11, ZIF-70, and ZIF-82 all exhibit good maximum uptakes. Out of these, the pore window diameter in ZIF-11 is the bottleneck for diffusion, thus reducing release rates, yet this will also hinder encapsulation i.e., the uptake of cisplatin into the pre-synthesized framework. In ZIF-70, there are fewer polar groups and therefore weaker interaction energies will hinder encapsulation while enabling premature release rates. ZIF-82 is the most promising material due to its large uptake, variation in interaction energies, and the addition of polar groups to anchor cisplatin and prevent premature drug release.

Details of the ab initio, GASP, BAR, equilibrium MD, umbrella sampling, and SMD simulations as well as MOF structural characteristics and potential energy maps can be found in the supplementary material.

This work was supported by funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant Agreement No 648283 “GROWMOF”). This research made use of the Balena High Performance Computing (HPC) Service at the University of Bath.

The data that support the findings of this study are available within the article and its supplementary material.

1.
See http://tiny.cc/rzwpjz for WHO, Cancer, 2020 (cited 2020 07/02/2020).
2.
L. B. P.
Boyle
, World Cancer Report No. 2008,
IARC Press
,
Lyon
,
2008
.
3.
M. J.
Lind
, “
Principles of cytotoxic chemotherapy
,”
Medicine
36
(
1
),
19
23
(
2008
).
4.
D.
Moreno
et al., “
Pharmacodynamics of cisplatin-loaded PLGA nanoparticles administered to tumor-bearing mice
,”
Eur. J. Pharm. Biopharm.
74
,
265
274
(
2009
).
5.
J. M.
Terwogt
et al., “
Phase I and pharmacokinetic study of SPI-77, a liposomal encapsulated dosage form of cisplatin
,”
Cancer Chemother. Pharmacol.
49
(
3
),
201
210
(
2002
).
6.
K. N. J.
Burger
et al., “
Nanocapsules: Lipid-coated aggregates of cisplatin with high cytotoxicity
,”
Nat. Med.
8
(
1
),
81
84
(
2002
).
7.
M.
Kettering
et al., “
Characterization of iron oxide nanoparticles adsorbed with cisplatin for biomedical applications
,”
Phys. Med. Biol.
54
(
17
),
5109
5121
(
2009
).
8.
L.
Ren
et al., “
Cisplatin-loaded Au-Au2S nanoparticles for potential cancer therapy: Cytotoxicity, in vitro carcinogenicity, and cellular uptake
,”
J. Biomed. Mater. Res., Part A
85
(
3
),
787
796
(
2008
).
9.
P.
Yang
,
S.
Gai
, and
J.
Lin
, “
Functionalized mesoporous silica materials for controlled drug delivery
,”
Chem. Soc. Rev.
41
(
9
),
3679
3698
(
2012
).
10.
X.
Duan
et al., “
Nanoparticle formulations of cisplatin for cancer therapy
,”
Wiley Interdiscip. Rev.: Nanomed. Nanobiotechnol.
8
(
5
),
776
791
(
2016
).
11.
O. M.
Yaghi
et al., “
Reticular synthesis and the design of new materials
,”
Nature
423
,
705
714
(
2003
).
12.
K. M. L.
Taylor-Pashow
et al., “
Postsynthetic modifications of iron-carboxylate nanoscale metal−organic frameworks for imaging and drug delivery
,”
J. Am. Chem. Soc.
131
(
40
),
14261
14263
(
2009
).
13.
C.
He
et al., “
Nanoscale metal–organic frameworks for the co-delivery of cisplatin and pooled siRNAs to enhance therapeutic efficacy in drug-resistant ovarian cancer cells
,”
J. Am. Chem. Soc.
136
(
14
),
5181
5184
(
2014
).
14.
S.-X.
Lin
et al., “
Effective loading of cisplatin into a nanoscale UiO-66 metal–organic framework with preformed defects
,”
Dalton Trans.
48
(
16
),
5308
5314
(
2019
).
15.
K. A.
Mocniak
et al., “
Incorporation of cisplatin into the metal–organic frameworks UiO66-NH2 and UiO66—Encapsulation vs. conjugation
,”
R. Soc. Chem. Adv.
5
(
102
),
83648
83656
(
2015
).
16.
A. J.
Howarth
et al., “
Chemical, thermal and mechanical stabilities of metal–organic frameworks
,”
Nat. Rev. Mater.
1
(
3
),
15018
(
2016
).
17.
Y.
Kato
et al., “
Acidic extracellular microenvironment and cancer
,”
Cancer Cell Int.
13
(
1
),
89
(
2013
).
18.
M.
Chen
et al., “
Extracellular pH is a biomarker enabling detection of breast cancer and liver cancer using CEST MRI
,”
Oncotarget
8
(
28
),
45759
45767
(
2017
).
19.
D.
Chen
et al., “
Cancer cell membrane-decorated zeolitic-imidazolate frameworks codelivering cisplatin and oleanolic acid induce apoptosis and reversed multidrug resistance on bladder carcinoma cells
,”
ACS Omega
5
(
2
),
995
1002
(
2020
).
20.
C.-Y.
Sun
et al., “
Zeolitic imidazolate framework-8 as efficient pH-sensitive drug delivery vehicle
,”
Dalton Trans.
41
(
23
),
6906
6909
(
2012
).
21.
H.
Ren
et al., “
Polyacrylic acid@zeolitic imidazolate framework-8 nanoparticles with ultrahigh drug loading capability for pH-sensitive drug release
,”
Chem. Commun.
50
(
8
),
1000
1002
(
2014
).
22.
L. R.
de Moura Ferraz
et al., “
ZIF-8 as a promising drug delivery system for benznidazole: Development, characterization, in vitro dialysis release and cytotoxicity
,”
Sci. Rep.
10
,
16815
(
2020
).
23.
Q.
Wang
et al., “
Synthesis and modification of ZIF-8 and its application in drug delivery and tumor therapy
,”
RSC Adv.
10
(
62
),
37600
37620
(
2020
).
24.
C.
Tamames-Tabar
et al., “
Cytotoxicity of nanoscaled metal–organic frameworks
,”
J. Mater. Chem. B
2
(
3
),
262
271
(
2014
).
25.
M.
Hoop
et al., “
Biocompatibility characteristics of the metal organic framework ZIF-8 for therapeutical applications
,”
Appl. Mater. Today
11
,
13
21
(
2018
).
26.
D.
Fairen-Jimenez
et al., “
Opening the gate: Framework flexibility in ZIF-8 explored by experiments and simulations
,”
J. Am. Chem. Soc.
133
(
23
),
8900
8902
(
2011
).
27.
K. S.
Park
et al., “
Exceptional chemical and thermal stability of zeolitic imidazolate frameworks
,”
Proc. Natl. Acad. Sci. U. S. A.
103
(
27
),
10186
10191
(
2006
).
28.
R.
Banerjee
et al., “
Control of pore size and functionality in isoreticular zeolitic imidazolate frameworks and their carbon dioxide selective capture properties
,”
J. Am. Chem. Soc.
131
(
11
),
3875
3877
(
2009
).
29.
M. C.
Bernini
, “
Screening of bio-compatible metal–organic frameworks as potential drug carriers using Monte Carlo simulations
,”
J. Mater. Chem. B
2
(
7
),
766
774
(
2014
).
30.
I.
Erucar
and
S.
Keskin
, “
Computational investigation of metal organic frameworks for storage and delivery of anticancer drugs
,”
J. Mater. Chem. B
5
(
35
),
7342
7351
(
2017
).
31.
M.
Kotzabasaki
and
G. E.
Froudakis
, “
Review of computer simulations on anti-cancer drug delivery in MOFs
,”
Inorg. Chem. Front.
5
(
6
),
1255
1272
(
2018
).
32.
J.-Q.
Liu
et al., “
A combined experimental and computational study of novel nanocage-based metal–organic frameworks for drug delivery
,”
Dalton Trans.
44
(
44
),
19370
19382
(
2015
).
33.
J.
Wang
et al., “
Combined experimental and theoretical insight into the drug delivery of nanoporous metal–organic frameworks
,”
R. Soc. Chem. Adv.
5
(
104
),
85606
85612
(
2015
).
34.
F.
Li
et al., “
Encapsulation of pharmaceutical ingredient linker in metal–organic framework: Combined experimental and theoretical insight into the drug delivery
,”
R. Soc. Chem. Adv.
6
(
53
),
47959
47965
(
2016
).
35.
S.
Yesylevskyy
et al., “
Empirical force field for cisplatin based on quantum dynamics data: Case study of new parameterization scheme for coordination compounds
,”
J. Mol. Model.
21
(
10
),
268
(
2015
).
36.
W. L.
Jorgensen
et al., “
Comparison of simple potential functions for simulating liquid water
,”
J. Chem. Phys.
79
(
2
),
926
935
(
1983
).
37.
A. K.
Rappe
et al., “
UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations
,”
J. Am. Chem. Soc.
114
(
25
),
10024
10035
(
1992
).
38.
S.
Hamad
et al., “
Atomic charges for modeling metal–organic frameworks: Why and how
,”
J. Solid State Chem.
223
,
144
151
(
2015
).
39.
A.
Gupta
et al., “
Object-oriented programming paradigms for molecular modeling
,”
Mol. Simul.
29
(
1
),
29
46
(
2003
).
40.
D.
Wolf
et al., “
Exact method for the simulation of Coulombic systems by spherically truncated, pairwise r−1 summation
,”
J. Chem. Phys.
110
(
17
),
8254
8282
(
1999
).
41.
H. J. C.
Berendsen
,
D.
van der Spoel
, and
R.
van Drunen
, “
GROMACS: A message-passing parallel molecular dynamics implementation
,”
Comput. Phys. Commun.
91
(
1
),
43
56
(
1995
).
42.
H.
Bekker
 et al, “
GROMACS: A parallel computer for molecular dynamics simulations
,” in
Physics Computing ’92
(
World Scientific Publishing
,
1993
), pp.
252
256
.
43.
E.
Lindahl
,
B.
Hess
, and
D.
van der Spoel
, “
GROMACS 3.0: A package for molecular simulation and trajectory analysis
,”
Mol. Model. Annu.
7
(
8
),
306
317
(
2001
).
44.
D.
Van Der Spoel
et al., “
GROMACS: Fast, flexible, and free
,”
J. Comput. Chem.
26
(
16
),
1701
1718
(
2005
).
45.
B.
Hess
et al., “
GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation
,”
J. Chem. Theory Comput.
4
(
3
),
435
447
(
2008
).
46.
S.
Pronk
et al., “
GROMACS 4.5: A high-throughput and highly parallel open source molecular simulation toolkit
,”
Bioinformatics
29
(
7
),
845
854
(
2013
).
47.
G.
Garberoglio
, “
OBGMX: A web-based generator of GROMACS topologies for molecular and periodic systems using the universal force field
,”
J. Comput. Chem.
33
(
27
),
2204
2208
(
2012
).
48.
T.
Darden
,
D.
York
, and
L.
Pedersen
, “
Particle mesh Ewald: An N · log(N) method for Ewald sums in large systems
,”
J. Chem. Phys.
98
(
12
),
10089
10092
(
1993
).
49.
U.
Essmann
et al., “
A smooth particle mesh Ewald method
,”
J. Chem. Phys.
103
(
19
),
8577
8593
(
1995
).
50.
G.
Bussi
,
D.
Donadio
, and
M.
Parrinello
, “
Canonical sampling through velocity rescaling
,”
J. Chem. Phys.
126
(
1
),
014101
(
2007
).
51.
S.
Nosé
, “
A unified formulation of the constant temperature molecular dynamics methods
,”
J. Chem. Phys.
81
(
1
),
511
519
(
1984
).
52.
W. G.
Hoover
, “
Canonical dynamics: Equilibrium phase-space distributions
,”
Phys. Rev. A
31
(
3
),
1695
1697
(
1985
).
53.
C. H.
Bennett
, “
Efficient estimation of free energy differences from Monte Carlo data
,”
J. Comput. Phys.
22
(
2
),
245
268
(
1976
).
54.
M. R.
Shirts
et al., “
Extremely precise free energy calculations of amino acid side chain analogs: Comparison of common molecular mechanics force fields for proteins
,”
J. Chem. Phys.
119
(
11
),
5740
5761
(
2003
).
55.
C. F.
Macrae
et al., “
Mercury 4.0: From visualization to analysis, design and prediction
,”
J. Appl. Crystallogr.
53
,
226
235
(
2020
).
56.
E.
Darve
,
D.
Rodríguez-Gómez
, and
A.
Pohorille
, “
Adaptive biasing force method for scalar and vector free energy calculations
,”
J. Chem. Phys.
128
(
14
),
144120
(
2008
).
57.
G. M.
Torrie
and
J. P.
Valleau
, “
Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling
,”
J. Comput. Phys.
23
(
2
),
187
199
(
1977
).
58.
S.
Kumar
et al., “
THE weighted histogram analysis method for free-energy calculations on biomolecules. I. The method
,”
J. Comput. Chem.
13
(
8
),
1011
1021
(
1992
).
59.
S. A.
Wells
and
A.
Sartbaeva
, “
GASP: Software for geometric simulations of flexibility in polyhedral and molecular framework structures
,”
Mol. Simul.
41
(
16–17
),
1409
1421
(
2015
).
60.
L.
Sarkisov
and
A.
Harrison
, “
Computational structure characterisation tools in application to ordered and disordered porous materials
,”
Mol. Simul.
37
(
15
),
1248
1257
(
2011
).
61.
W.
Zhang
et al., “
Cooperative effect of temperature and linker functionality on CO2 capture from industrial gas mixtures in metal–organic frameworks: A combined experimental and molecular simulation study
,”
Phys. Chem. Chem. Phys.
14
(
7
),
2317
2325
(
2012
).
62.
J. H.
Cavka
et al., “
A new zirconium inorganic building brick forming metal organic frameworks with exceptional stability
,”
J. Am. Chem. Soc.
130
(
42
),
13850
13851
(
2008
).
63.
D. I.
Kolokolov
et al., “
Probing the dynamics of the porous Zr terephthalate UiO-66 framework using 2H NMR and neutron scattering
,”
J. Phys. Chem. C
116
(
22
),
12131
12136
(
2012
).
64.
S.
Pakhira
, “
Rotational dynamics of the organic bridging linkers in metal–organic frameworks and their substituent effects on the rotational energy barrier
,”
R. Soc. Chem. Adv.
9
(
65
),
38137
38147
(
2019
).
65.
S.
Van der Perre
et al., “
Adsorptive characterization of the ZIF-68 metal-organic framework: A complex structure with amphiphilic properties
,”
Langmuir
30
(
28
),
8416
8424
(
2014
).
66.
S.
Amirjalayer
and
R.
Schmid
, “
Influence of pore dimension on the host–guest interaction in metal–organic frameworks
,”
J. Phys. Chem. C
120
(
48
),
27319
27327
(
2016
).

Supplementary Material