DFT-Based Theoretical Model for ... - ACS Publications

Novel pKa/DFT-Based Theoretical Model for Predicting the Drug Loading and Release of a pH-Responsive Drug Delivery System. Shaukatali N. Inamdar† ...
0 downloads 0 Views 3MB Size
Article pubs.acs.org/JPCC

Cite This: J. Phys. Chem. C 2018, 122, 12279−12290

Novel pKa/DFT-Based Theoretical Model for Predicting the Drug Loading and Release of a pH-Responsive Drug Delivery System Shaukatali N. Inamdar,† Khalid Ahmed,† Nashiour Rohman,† and Adam A. Skelton*,†,‡ †

Department of Pharmaceutical Sciences, University of KwaZulu-Natal, Durban 4000, South Africa Department of Chemistry, University of Liverpool, Liverpool L69 7ZD, U.K.



Downloaded via UNIV OF TOLEDO on June 15, 2018 at 14:20:52 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: Herein, we present a theoretical model that combines classical pKa theory with quantum mechanical calculations to predict the extent of interaction between acid-/base-dependent species over a full range of pH conditions. To demonstrate the theoretical model, we have predicted the drug loading and release of a pH-responsive drug delivery system consisting of sulfasalazine, an anionic antiinflammatory drug molecule, loaded onto the positively charged trimethylammonium (TA)-functionalized mesoporous silica nanoparticle surface. The model relies on the possible combinations of pH-dependent states of the surface (S) and drug (D) molecules as neutral (0) and deprotonated (−1) states, whose relative probabilities depend on their pKa value and the desired pH. The four possible combinations were identified as S0D0, S0D−1, S−1D0, and S−1D−1, and periodic density functional theory calculations were performed for systems comprising drug fragments adsorbed onto a model TA-functionalized quartz pH surface to calculate the pH-dependent interaction energy (EpH int ). The Eint value was attractive at the acidic environment of stomach (pH 2−5), where a drug is loaded and repulsive at neutral or slightly basic pH in the small intestine where it is released; the behavior is in accordance with the experimentally reported data. We also validated the experimental necessity of surface functionalization.



only in response to the aforementioned changes in pH,4−10 thereby going directly to where the drug is required. An advantage of such targeted drug delivery is that lower doses of the drugs are required, which reduces the side effects, whereas a controlled drug release offers the advantage of less-frequent administration and better compliance.28,29 Mesoporous silica nanoparticles (MSNs) are of interest as pH-responsive drug delivery systems because of their biocompatibility, low cytotoxicity, biodegradability, excretion, ordered and uniform size, high surface areas, excellent stability, and chemically modifiable surface properties.30−35 Two types of mesoporous silica materials, MCM-4136,37 and SBA-15,38,39 can act as reservoirs for drug molecules because of large surface areas and silanol groups (Si-OH) on their internal pore surfaces, which provide adsorption sites for H-bonding interactions. Although the silanol groups are neutral at highly acidic pH, a negative surface charge develops at neutral pH as a result of deprotonation of silanol groups.40−42 It has been suggested that the pKa values of these silanol groups are 4.5 and 8.5.40,41 As the interaction with an adsorbate such as a drug molecule should be different for a neutral silanol group versus a deprotonated silanol group, the acid−base behavior provides an

INTRODUCTION The pH of a solution, which is the measure of the concentration of protons, has a direct effect on any protonbearing atom in a molecule or material in contact with the solution, which will in turn be protonated or deprotonated.1,2 The extent to which an acidic or basic moiety will be protonated/deprotonated will depend on its acid−base equilibria at the pH in question.1,3 The protonation behavior of molecular entities has profound effects on the chemistry and interactions of many different systems in biology,1,3 nanotechnology,4−10 and coordination chemistry.11,12 pH can, therefore, affect the rate and yield of chemical reactions,13 complex interactions in biochemical processes,14 the efficacy of drugs,15 pharmacology, and the stability of coordination compounds.11,12 The pH in the human body varies from highly acidic to basic in the digestive system. The pH is 1.0−3.0 in the stomach, 7.4− 7.5 in the blood plasma, and 7.5−8.0 in the lower intestine. Furthermore, on the tissue or cellular level, there are small pH variations inside cancerous tumors (pH 6.5−6.8),15 bacterial infections (pH 7.4−7.8),16 and particular organelles such as lysosomes (pH ∼ 4−5).17,18 Nanotechnology applications are being developed to manipulate these variations in pH.4−10 Drug molecules can be loaded onto nanosized formulations such as liposomes,19 polymers,20,21 micelles,22,23 dendrimers,24 and nanoparticles of silica,25,26 carbon,27 or Fe3O4 and released © 2018 American Chemical Society

Received: March 23, 2018 Revised: May 8, 2018 Published: May 9, 2018 12279

DOI: 10.1021/acs.jpcc.8b02794 J. Phys. Chem. C 2018, 122, 12279−12290

Article

The Journal of Physical Chemistry C

the free energy required to protonate a titratable group and pHresponsive protein folding. The theoretical work that has been developed has shown to link up with computational simulations but is intractable to most computational chemists. We propose a method for calculating the extent of binding between two molecular species at any required pH; this method combines density functional theory (DFT) calculations with the theoretical fraction derived from the rearranged Henderson−Hasselbalch equation.52 If both species can exist in both neutral and deprotonated states, there will be a percentage of each combination of states; these are deprotonated species 1 with deprotonated species 2, deprotonated species 1 with neutral species 2, neutral species 1 with deprotonated species 2, and neutral species 1 with neutral species 2. As the proposed method uses an analytical equation to rescale the results of standard energy values, it can be used with any molecular modeling software that can calculate the energy or specific properties of two molecules. To thoroughly test and validate our method, we applied the model to predict the extent of loading of sulfasalazine drug with positively charged TA-functionalized MSN over a range of pH values and compared our prediction with experimental drug loading reported by Lee et al.43 As the TA group is a quaternary amine, it always has a positive (+1) charge; therefore, the first species is one of two fragments of sulfasalazine, N-(2-pyridyl) benzene sulfonamide or salicylic acid, each of which is acidic and can be neutral or deprotonated. The second species is the TA-functionalized or -nonfunctionalized (101) surface of quartz,55,56 which, due to the ordered silanol groups on its surface, is used as a model for the internal surface of MSNs. The silica surface is either completely protonated or has one deprotonated silanol group. We have run a series of periodic density functional theory calculations related to the pHdependent nature of silica (S) and the drug (D) molecule, giving rise to four combinations, S0D0, S0D−1, S−1D0, and S−1D−1. Our method then enables us to calculate the pHdependent binding energy (EpH int ), the interaction energy over a range of pH values with reference to the pKa value of the drug and carrier surface. Apart from acting as a proof of concept for our novel method that can predict binding of any molecular species over a range of pH values, we show that our model has a significant value for our chosen application to calculate the drug loading and release pH. Our proposed model will aid in the design of a suitable functionalized nanoparticle/drug pair and will help to predict the optimal pH with which the drug will be released. Our method can be used to guide the experimental process and provide a fundamental insight into the mechanism of pHresponsive drug loading and release.

opportunity to modify the interaction and, therefore, the amount of loading in response to pH. It should be noted that the hydrogen bonding of a drug molecule with the neutral silanol group needs to be enhanced for significant loading, and, therefore, additional organic function groups, such as amino and carboxylic acids, can be introduced onto the surface to increase the loading capacity.43 An example of such a system is the pH-responsive adsorption and desorption of anti-inflammatory drug sulfasalazine on a positively charged trimethylammonium (TA)-functionalized silica surface.43 Lee et al. reported that the amount of adsorbed sulfasalazine on MSNs functionalized with tetraalkylammonium (TA) was high at pH 2 and systematically decreased until pH 6, whereas it leveled off at basic pH 11.43 It was asserted that the drug could be loaded at acidic pH conditions of the stomach because of the attraction between the anionic drug and the positive charge of the functional group (FG) (see Figure 2A for the structure of sulfasalazine). Such binding allows the drug to remain trapped in the MSN-TA framework while passing through the acidic pH of the stomach and allowed to be released only at pH 7.4 of the small intestine. The authors postulate that the drug is released because of the electrostatic repulsion between the anionic drug molecule and the negatively charged silica surface arising from the deprotonation of the silanol group.43 Lee et al. also showed that for nonfunctionalized MSNs at pH 5.5 the binding was negligible and increased with increased concentration of the TA functional group.43 Molecular modeling is being widely used to predict the behavior and interactions among molecules.44 Quantum mechanical (QM) methods are used to provide the detailed interactions between a small subset of the real system or often two molecular species.45−48 Researchers interested in using QM methods to calculate interactions at different pH values usually simulate the different protonation states and assume that at that pH a particular protonation state is dominant.49 In reality, there are fractions of protonated, deprotonated, and neutral species in a thermodynamic ensemble of a large number of molecules. In classical molecular dynamics simulations, we can treat such situations by protonating or deprotonating the necessary percentage of molecules or surface species42,50 to adhere to a realistic system. Another widely used strategy is constant pH molecular dynamics (cpH-MD), where the protonation state is dynamically altered on the basis of the free energy of protonation.3,51 Despite the elegance of cpH-MD, these methods require extra parameters that need to be accurate to produce meaningful results; moreover, sophisticated software is required to implement cpH-MD and, therefore, it is not trivial to use for any chosen molecular modeling package. In QM calculations, due to small system sizes, the strategy of modeling the ensemble of protonation states is not possible, as generally a single interaction will be monitored. To overcome this limitation, finding a way to account for the thermodynamic fraction of protonated/deprotonated species to explain the pHresponsive behavior of adsorption of a single molecule is required. The rearranged Henderson−Hasselbalch equation provides the fraction of molecular species that is protonated or neutral.52 Much work on the development of pKa theory has been carried out with the focus on calculating the overall titration curve of proteins with a large number of interacting titratable groups.53 Furthermore, theoretical work was used to predict the distribution of protonation states in a protein54 to calculate



METHODS pKa/DFT-Based Theoretical Binding Model. For an acidic species at any given pH, the fractional concentration (F) of each of protonated (F−1) and neutral states (F0) depends on their pKa value and is given by the rearranged Henderson− Hasselbalch equation,52 as shown in eq 1. 1 ; 1 + 10 pKa − pH 1 1 + 10 pKa − pH

F −1(deprotonated) = F0(neutral) = 1 −

(1)

If we now consider two acidic species, a surface (S) and a drug molecule (D), each of which can be deprotonated (S−1 and 12280

DOI: 10.1021/acs.jpcc.8b02794 J. Phys. Chem. C 2018, 122, 12279−12290

Article

The Journal of Physical Chemistry C D−1) or neutral (S0 and D0), there will be four possible combinations, as shown in Figure 1. The combinations are

S

FS0D−1 = F0 × F −1 =

10 pKa − pH D

S

(1 + 10 pKa − pH)(1 + 10 pKa − pH) (3)

FS−1D−1 = F −1 × F −1 =

1 (1 + 10

pK aD − pH

S

)(1 + 10 pKa − pH) (4)

S

FS0D0 = F0 × F0 =

D

(10 pKa − pH)(10 pKa − pH) D

S

(1 + 10 pKa − pH)(1 + 10 pKa − pH)

(5)

It can be noticed that there are three different forms for the equations, where FS−1D0 and FS0D−1 differ by only the numerator. The two molecular species in each combination will interact with each other to different extents, referred to as the interaction energies (ES0D0, ES0D−1, ES−1D0, and ES−1D−1). At a specific pH, each combination will be present according to its fraction and, therefore, the pH-dependent interaction energy (EpH BE ) will be a superposition of all energies scaled by its fractions and is given by eq 6 pH E int = (ES0D0 × FS0D0) + (ES0D−1 × FS0D−1)

Figure 1. Schematic to highlight the different combinations of protonation states.

+ (ES−1D0 × FS−1D0) + (ES−1D−1 × FS−1D−1)

The full form of eq 6 is shown in the Supporting Information (eq S1). It should be noted that eqs 2−5 are set up to be independent of binding because the strength of the binding (interaction energies) is calculated later and included in eq 6. In fact, eqs 2−5 provide the probability that the two species would come into contact if there were no specific attractive interaction and if the drug and surface groups came into contact by a random diffusion process. Moreover, eqs 2−5 provide the probability of the four different species combinations during this random diffusion process. To elaborate further, we have provided a demonstration that our equations can be applied to an ensemble of molecules to provide the number of molecules in each combination (“statistics of binding” section and Figure S1). The demonstration highlights that the present model provides the number of surface states available for any number of molecules.

S0D0, where the neutral surface interacts with the neutral drug molecule, S0D−1, where the neutral surface interacts with the deprotonated drug molecule, S−1D0, where the deprotonated surface interacts with the neutral drug molecule, and S−1D−1, where the deprotonated surface interacts with the deprotonated drug molecule. Each combination will result in a different interaction, and the relative probability that a particular interaction will occur at a particular pH will be the fraction of each combination (FS−1D0 is the fraction of the S−1D0 combination and so on); this is given by the multiplication of the fractions of the individual contributions and is given in eqs 2−5. D

FS−1D0 = F −1 × F0 =

10 pKa − pH D

(6)

S

(1 + 10 pKa − pH)(1 + 10 pKa − pH) (2)

Figure 2. Acid−base reactions of the drug fragments N-(2-pyridyl) benzene sulfonamide (PBS) and salicylic acid (SA) (B, C) and the silica surface (D). (A) is the whole sulfasalazine molecule. 12281

DOI: 10.1021/acs.jpcc.8b02794 J. Phys. Chem. C 2018, 122, 12279−12290

Article

The Journal of Physical Chemistry C

Figure 3. Illustration of the dissociation of (A) PBS (pKa 6.89), (B) SA (pKa 2.97), (C) silica (pKa 4.5), and (D) silica (pKa 8.5) versus pH. The red line represents the fraction of the dissociated species (deprotonated (−1)), and the black line indicates the undissociated (neutral (0)) state of the molecules.

This measure is related to the entropy of the system as it provides the number of possible interaction states, which will be a key component of the drug molecule binding to the surface. The energy component is then included in the form of interaction energies to complete the picture. The individual interaction energies are currently not available through theory alone, and, therefore, these quantities are obtained by running computational chemistry simulations such as density functional theory calculations. Application of the Theoretical Model to Drug Delivery. We have applied the above theoretical model to understand the pH-responsive drug binding and release of sulfasalazine drug with TA-functionalized mesoporous silica nanoparticle surfaces. The sulfasalazine drug (Figure 2A) is reported to have two pKa values, the carboxylic acid group has a pKa value of 3.3 and sulfonamide nitrogen, −SO2NH, has a pKa value of 6.24 (Figure 2A). 57 Sulfonamide contains the sulfonamido (−SO2NH) group, which is weakly acidic (pKa value approximately 5−8) because of the powerful electron-withdrawing effect of −SO2− substituent and stabilization of the resulting anion by resonance.58 For practical reasons, explained in the next section, we divided sulfasalazine into two fragments, namely N-(2-pyridyl) benzene sulfonamide (PBS) (Figure 2B) and salicylic acid (SA) (Figure 2C), each of which have single pKa value. PBS is reported to have a pKa value of 6.89,59 whereas SA is reported to have a pKa value of 2.97.60,61 At a neutral pH, a certain percentage of surface hydroxyl or silanol groups (Si−O−H) of silica is deprotonated, thereby forming a negative surface charge.42 A decrease of pH will cause a predictable decrease in the surface charge as the deprotonated silanol groups become gradually protonated in an acidic solution.62

Silica has been reported to have two separate pKa values of 4.5 and 8.5.40,41,63 In the present study, we tested our theoretical model on both pKa values of 4.5 and 8.5 (Figure 2D). Density Functional Theory Calculations. To calculate the individual interaction energies, periodic density functional theory (DFT) calculations were performed in the CASTEP64 module of BIOVIA Materials Studio 2016.65 The generalized gradient approximation was performed using the exchangecorrelation functional of Perdew, Burke and Ernzerhof.66,67 Ultrasoft pseudopotentials in the Vanderbilt’s formulation68 along with a medium plane wave cutoff energy of 489.8 eV were used throughout the calculations. Generally, MSNs contain amorphous silica and their surfaces are covered with silanol groups. As a first approximation, we opted to use the most stable and well-characterized (101) surface of quartz.69 We consider the quartz (101) surface as an adequate model because it contains the requisite silanol groups, which can be protonated or deprotonated depending on the pH. For periodic DFT calculations, to mimic the mesoporous silica surface, we used the silica quartz (101) surface with 2 × 3 × 1 periodic supercells, with five atomic layers and vacuum thickness of 20 Å in the z-direction. The structure of the sulfasalazine molecule is too large to fit onto the surface within a reasonable sized cell. If we were to make the cell large enough to fit flat onto the surface and adequately sample conformational space, it would require more than 800 atoms, which would be computationally too expensive to be practical. We have, therefore, divided sulfasalazine into two fragments, PBS (Figure 2B) and SA, (Figure 2C) and placed them onto the surface. Geometry optimizations were performed using CASTEP with maximum SCF cycles set at 500 cycles. The TA functional group (FG) is attached to the silica surface on the middle Si 12282

DOI: 10.1021/acs.jpcc.8b02794 J. Phys. Chem. C 2018, 122, 12279−12290

Article

The Journal of Physical Chemistry C

Figure 4. Variations with pH in the fractions for the various combinations, S0D0, S0D−1, S−1D0, and S−1D−1, of SA and PBS with the surface at different pKa values.

species and most of the variations in the fraction occur within 1 pH unit either side of the pKa; that is, within 2 pH units, a species will convert from being mostly protonated to mostly neutral. Next, we consider how the fractions of the various combinations, S0D0, S0D−1, S−1D0, and S−1D−1, vary with the pH, where the drug molecule is either PBS or SA and the surface has either a pKa value of 4.5 or 8.5 (Figure 4). This measure provides the relative probability that a drug−surface interaction will be the interaction defined by a particular combination. It can be noticed that the FS−1D−1 values for all drug molecule/ surface pairs behave like that of FD−1 or FS−1 that has the higher value; for example, FS−1D−1 for the interaction between PBS and the surface at a pKa value of 8.5 (PBS6.89―surface8.5) is the same as FS−1 but FS−1D−1 for PBS6.89―surface4.5 is the same as FD−1. Likewise, the FS0D0 values for all drug molecule/surface pairs behave like that of FD0 or FS0 that has the lower value. These observations are significant because they provide insight into the possible interactions. At a low pH, all of the interactions are S0D0, which should be reasonably attractive interactions because of the positively charged functional group; these are largely responsible for loading. Conversely, at high pH values, S−1D−1 are the main interactions, which should be more repulsive and responsible for the drug release. At intermediate pH values, either S0D−1 or S−1D0 is dominant and these arise as peaks where the maximum is at the pH that is exactly the center point between the two pK a values; for example, for PBS6.89―surface4.5, the peak of FS−1D0 is 4.5 + ((6.89 − 4.5)/2) = 5 .695. It s hould be noted t hat f or PBS6.89―surface4.5 pairs the S−1D0 interaction is dominant whereas for all other systems the S0D−1 interaction is dominant. The reason for this difference is that for PBS6.89―surface4.5 the surface pKa value is lower than that of the drug molecule

atom in the second row replacing the hydroxyl group from the top, whereas for deprotonating the silica surface, H atom was removed from the hydroxyl group of the middle Si atom in the last row (Figure S2). The interactions of drug molecules with the silica surface, in the absence and presence of FG, were calculated using the CASTEP calculations. For the silica surface with a positively charged FG, the interaction energy (Eint) is calculated using eq 7, whereas for the system without any FG, eq 8 is used.44 E int = E(silica with FG + drug) − E(silica with FG) − E(drug)

(7)

E int = E(silica +drug) − E(silica) − E(drug)

(8)

To monitor all possible ways of the interaction of a drug with the silica surface, we have placed the drug fractions in six varied positions (Figures S3 and S4). The PBS molecule was placed on the silica surface with FGs at various positions, as shown in Figure S3A−F. Although the drug fraction SA was also treated in a similar fashion, it was placed on the silica surface with FGs in different positions, as shown in Figure S4a−f. We report 224 CASTEP calculations, where we have varied the placement of each drug fragment on the functionalized and nonfunctionalized silica surfaces in six positions each.



RESULTS AND DISCUSSION To provide an insight into the behavior of the Henderson− Hasselbalch equation (eq 1), the variations in pH, in the fractions of neutral (D0) and deprotonated states (D−1) of PBS and SA, are shown in Figure 3A,B, respectively. Similarly, the variations in the fractions of neutral (D0) and deprotonated states (D−1) of silica with pH, at pKa values of 4.5 and 8.5 (surface4.5 and surface8.5, respectively), are shown in Figure 3C,D, respectively. When the pH is the value of the pKa, there is an equal amount of the neutral (0.5) and deprotonated (0.5) 12283

DOI: 10.1021/acs.jpcc.8b02794 J. Phys. Chem. C 2018, 122, 12279−12290

Article

The Journal of Physical Chemistry C

Table 1. Lowest Interaction Energies (Eint) in kcal mol−1 Calculated for the Drug Fractions PBS and SA with a Silica Surface combination

S0D0

S0D−1

S−1D0

S−1D−1

−1

lowest lowest

PBS on S-FG (Eint (kcal mol )) −32.19 −9.03 −30.13 SA on S-FG (Eint (kcal mol−1)) −5.41 −17.00 −38.70

S0D0

S0D−1

S−1D0

S−1D−1

−4.64

86.8204

−29.39

11.41

−1

9.43 4.88

PBS on S (Eint (kcal mol )) −5.88 68.61 SA on S (Eint (kcal mol−1)) −20.73 55.49

Figure 5. Lowest energy configurations for all of the different protonation combinations of PBS and SA on the quartz surfaces. Corresponding interaction energy in kcal mol−1 is given below for each combination.

For PBS, there are few H-bonds, with the exception of S0D−1, where there are short H-bonds from the sulfonamide oxygen atoms to the silanol hydrogen. As the energy of S0D−1 is less attractive than that of S0D0, however, it can be inferred that Hbonds have a less contribution to the binding than that of the ion/dipole interactions between the positively charged TA and the PBS drug fragment. As for S−1D0, it is likely that the negative interaction energy (−30.14 kcal/mol) is due to a cooperative effect of the electronegative sulfonamide oxygen atoms that are facing toward the positively charged TA coupled with the negatively charged deprotonated silanol oxygen atom. For SA, the effect of H-bonding and proton transfer is more significant than for PBS. For S−1D0, there is a proton transfer from the phenol hydroxyl hydrogen atom of SA to the deprotonated silanol oxygen, providing the highly negative interaction energy (−38.70 kcal/mol). The carboxyl shares its proton with the newly deprotonated phenol hydroxyl oxygen. A similar proton transfer occurs for the nonfunctionalized surface, providing a similar negative interaction energy (−29.3 kcal/ mol, Figure S5). For S0D−1, there is proton sharing between the oxygen atom of the deprotonated carboxyl and the silanol hydrogen atom. As the proton sharing was also available for the nonfunctionalized surface but resulting in a large positive interaction energy (55.79 kcal/mol, Figure S5), it indicates that the positive charge of the TA group electrostatically interacting with the negatively charged deprotonated SA was significant for the negative interaction energy (−17.00 kcal/mol). The interaction energy of S0F0 is lower for SA on the nonfunctionalized surface than on the functionalized surface, and this is because of favorable H-bonds between the hydrogen atom of the phenyl hydroxyl with silanol oxygen atoms and silanol hydrogen atoms with the hydroxyl oxygen (Figure S5)

but for all other systems the surface pKa value is higher than that of the drug fragment. We will call the FS−1D0 or FS0D−1 peaks the modulator peaks as they are not always responsible for loading or release but should modulate the pH-responsive release profile by affecting how sensitive to pH the binding is in the mid pH range; this will become more apparent when we input the interaction energies into the theoretical model below. It is notable to point out that as the difference in pKa values becomes greater the modulator peak becomes both higher and wider; for instance, compare SA 2.97 ―surface 4.5 with SA2.97―surface8.5. Interaction Energies of Each Combination. We have calculated the interaction energies (Eint) for each of the combinations using periodic density functional theory calculations and employing eq 7 for drug fractions on silica with a FG and eq 8 for drug fractions on silica in the absence of a FG. The lowest interaction energies (Eint) for the four possible combinations, S0D0, S0D−1, S−1D0, and S−1D−1, at all of the drug positions are shown in Table 1. In the presence of the functional group, both the PBS and SA drug fragments have negative, attractive interaction energies for all combinations, except for the S−1D−1, which is slightly positive. In the absence of the functional group, S−1D−1 is highly positive and the difference between the absence and presence of the functional group highlights the role of the positively charged TA group in increasing the binding interaction in mid pH ranges. Figure 5 shows the lowest energy configurations for all four combinations, for both PBS and SA on the functionalized surface. Short H-bonds that are less than 2.5 Å are shown for understanding the factors that contribute to the differences in Eint among the different combinations. 12284

DOI: 10.1021/acs.jpcc.8b02794 J. Phys. Chem. C 2018, 122, 12279−12290

Article

The Journal of Physical Chemistry C

Figure 6. pH-Dependent binding energies (EpH int ) calculated for (A) PBS, (B) SA with surface (pKa 4.5), and (C) the average of both PBS and SA.

Figure 7. pH-Dependent binding energies (EpH int ) calculated for (A) PBS, (B) SA with surface (pKa 8.5), and (C) the average of both PBS and SA.

Figure 8. Contributions of the different combinations to the pH-dependent binding energies (EpH BE ) calculated for PBS and SA with the functionalized surface at both pKa values.

Prediction of the Drug Loading and Release Behavior. Now as we have the energies of each combination, we can input

the lowest energies into eq 6 and plot the pH-responsive interaction energy. The resultant EpH int versus pH plots for PBS 12285

DOI: 10.1021/acs.jpcc.8b02794 J. Phys. Chem. C 2018, 122, 12279−12290

Article

The Journal of Physical Chemistry C

In contrast to PBS, for SA, the effect of the change in pKa value of the surface is drastic because the combination that acts as the modulator remains the same but the difference between the pKa value of the drug fragment and the surface gets larger. In this case, the effect of the modulator becomes greater as the difference between the pKa value of the drug fragment and the surface increases. It is notable to point out that the energy of the modulator (S0D−1, −17.00 kcal/mol) in the case of SA is much more negative than in the binder combination (S0D0, −5.4 kcal/mol) and results in the modulator curve being very influential to the overall EpH int curve. In the case of the nonfunctionalized surface (Figure S6), the modulator tends to act in an opposite fashion than for the functionalized surface because the S0D−1 values are highly positive (55.49 and 68.61 kcal/mol for SA and PBS, respectively). The modulator is responsible for the large repulsive EpH int values in the mid pH ranges. For PBS6.89―surface4.5, the modulator has little effect because it is S−1D0, which has only a slightly negative interaction energy (−4.64 kcal/mol). Application for the Design of Drug Delivery Systems. Currently, the main basis for deciding on an appropriate functional group and, therefore, the design of the pHresponsive drug delivery system is to consider whether the drug will interact with the functional group more than it would with the silica surface. The strategy considered by experimental researchers does indeed consider the charge of the silanol or functional group, how that matches to the charge of the drug, and how these charges are dependent on the pH. At a particular pH, Lee et al.43 considered whether the functional group or silanol groups are deprotonated, protonated, or neutral and considered the drug molecule to be always deprotonated. The choice of a particular functional group is dependent on a qualitative description of attraction and repulsion that relies on the chemical intuition yet ignores many detailed factors. The current work, however, has uncovered other factors such as hydrogen bonding interactions are at play and that to understand the drug loading and release behavior we should consider other protonation states of the drug. Furthermore, the current model provides us with tools for predicting the release rates of a range of functionalized MSNs (FMSNs). The following rules apply to the theoretical model: (1) If the pKa value of the surface is greater than that of the drug fragment, the modulator will be S−1D0; if pKa value of the surface is lower than that of the drug fragment, the modulator will be S0D−1. (2) If the interaction energy of the modulator is much more negative than that of the loading combination, S0D0, the modulator peak will make the curve much more complicated with less binding at very low concentrations and a plateau of binding in mid range pH values. The same applies for large positive energies that will affect the modulator in an equal but opposite way. (3) The difference in pKa values of the surface and drug fragment will accentuate the effect of the modulator. We can use these rules when choosing the ideal drug−surface combination; for example, usually, it is desirable for a very sharp pH response and in this case one should choose a drug− surface combination that has very close pKa values or low modulator interactions that would reduce the effect of the modulator. Furthermore, researchers should decide on whether they desire the surface or the drug to have the higher pKa value that will allow a different modulator to take effect. Armed with

and SA with functionalized and nonfunctionalized silica surfaces with pKa values of 4.5 and 8.5 are shown in Figures 6A,B and 7A,B, respectively. The loading and release profiles for the overall drug molecule, sulfasalazine, were obtained by the taking the average of the pH-responsive binding energies EpH int obtained for both the fragments PBS and SA and are shown in Figures 6C and 7C. More negative values for EpH int are indicative of binding, reflecting the loading of the drug molecule on the surface, pH values show the lack of attractive whereas positive Eint interaction or repulsion, indicating the release of the drug fragment from the surface. Overall, results show that for the functionalized silica at both surface pKa values and both drug fragments, there is binding (negative EpH int values) at lower pH (1−5) and at pH values greater than 7 the molecule is released from the surface (positive EpH int values); this is in quantitative agreement with experimental drug loading results by Lee et al.43 It should be 0 0 noted that the EpH int values at pH 1 are equal to the Eint of S D pH −1 −1 and Eint values at pH 14 are equal to the Eint of S D . PBS shows a very predictable drug release curve, moving from a plateau of negative EpH int values at low pH to a plateau of slightly values at high pH. Conversely, SA is less attractive positive EpH int at very low pH values than at pH ∼ 4. The pKa value for the surface does not really affect the profiles for PBS, except for making the release steeper at pH 8.5. For SA, however, a change of the pKa value for the surface from 4.5 to 8.5 causes a larger attractive plateau and release at higher pH values. As for the nonfunctionalized surfaces, quite different behavior to the functionalized surface is observed using our theoretical model. There is a similar shape of curve for PBS on functionalized and nonfunctionalized surfaces (both pKa pH values), but nonfunctionalized Eint values are much less attractive at low pH values and much more repulsive at high pH values. For SA, the EpH int values are only attractive at very low values but become highly repulsive at pH 4 before becoming less repulsive at high pH values. It can be noticed from the average curves (Figures 6C and 7C) that for the nonfunctionalized surfaces there is binding only at very low pH values. Our theoretical model, therefore, agrees with the experimental results that there was no significant binding of sulfasalazine at pH 5.5 for the nonfunctionalized surface where our model predicts positive EpH int values. To provide more details about the factors that affect the EpH int values, we have individually plotted the contribution of each combination, S0D0, S−1D0, and S−1D−1 (Figure 8). We have also plotted S0D0 + S−1D−1, the loading and release contributions, to understand the effect of the modulator peak, S−1D0 and S0D−1, on the overall EpH int values. First, it should be recalled that the modulator for SA 2.97 ―surface 4.5 is S −1 D 0 but for SA2.97―surface8.5 the modulator is S0D−1. As the interaction energy for S0D−1 is more negative (−30.13 kcal/mol) than for S−1D0 (−9.03 kcal/mol), the S0D−1 modulator contribution to EpH int is much more pronounced, resulting in a deeper negative well. For PBS6.89―surface4.5, therefore, the overall EpH int value is quite different to S0D0 + S−1D−1 because of the large modulating effect of S 0 D −1 ; this is in contrast to PBS6.89―surface8.5, which has very little modulation. Interestingly, for PBS, there is very little difference between the overall EpH int curves for the two different pKa values of the surface because the modulation, in the case of PBS6.89―surface4.5, shifts the EpH int curve to be comparable to that of PBS6.89―surface8.5. 12286

DOI: 10.1021/acs.jpcc.8b02794 J. Phys. Chem. C 2018, 122, 12279−12290

Article

The Journal of Physical Chemistry C

combination existing at a particular pH; this can then be applied for catalytic reactions at enzyme active sites,71 understanding the pH dependence of catalysis for homogeneous catalysis or any acid-/base-catalyzed organic mechanisms. Electronic properties such as band gap or fluorescent properties that may depend on pH may also be calculated in a similar way. The present model is not constrained to QM calculations but can be adopted for any molecular modeling technique and more complicated, multipoint data. For example, four classical MD simulations can be performed at different protonation states; a computer program can be written to input each point on a graph, for any MD analysis such as a radial distribution function (RDF) of a functional group with water, into eqs 2, 3, 4, or 5. All points on the RDF of each of the four simulations can be added to provide the RDF graph at the desired pH. Free energy profiles obtained by the potential of mean force or binding free energies can be treated in a similar manner. It should be noted that if a system has more than two acid− base species or species that have several protonation states, such as protonated, neutral, and deprotonated, our theoretical model should be extended. For more complicated systems, the number of combinations would increase exponentially and this provides the opportunity for complex theoretical models that may be best derived computationally and whose theoretical consequences can be explored through linear algebra53 or Monte Carlo techniques.72 Our model can be used to analyze the effect of the pHresponsive interactions on more realistic systems using classical molecular dynamics. As our model only requires that researchers input their data into eqs 2−6, it can easily be implemented without any specialist software and, therefore, the present study makes pH-responsive molecular modeling accessible for any researcher. Furthermore, we propose that any molecular simulation that only treats single protonation state at a defined pH will benefit from running a series of simulations and using our equations to provide their pHresponsive property. As most systems of interest, such as bio and nano systems, involve acid−base behavior, we envisage that our theoretical model will have a wide applicability across the vast landscape of computational chemistry and molecular simulation studies.

our theoretical model and with computational chemistry techniques, we are in a position to design drug delivery systems with improvements in the desired properties of drug loading and release and pH response. As researchers make educated guesses about which functional groups to use but without fully understanding the molecular basis behind their choice, there are limitations in how far they can progress in their molecular design. Moreover, there are practical limitations, which usually result in researchers using commercial precursors to synthesize the FMSNs. Although it is possible to modify the functional groups using organic synthesis techniques, these techniques are expensive and time-consuming. Without meaningful molecular insight into the factors that will affect the desired properties, this technique will be trial and error and may not yield any meaningful results. With the current method, however, it is possible to model any molecule that will accentuate pH-responsive binding and release even if they are difficult to synthesize. Such a strategy will allow us to tune the functional group until we find the optimal functional group/drug match that will bind with a particular drug molecule most strongly at the loading pH yet will release the drug molecules at the desired pH. For instance, the functional group can be modified with electron-withdrawing or -accepting substituents to tune the pKa value or its interaction with the surface at different protonation states. As alluded to, current drug delivery systems by Lee et al.43 used reasonably large pH changes from acidic pH of the stomach to neutral pH of the intestine to affect the drug release. However, one of the aims of advanced drug delivery systems is to manipulate much smaller pH changes for the controlled release of an anticancer drug into a tumor where there are slightly acidic pH conditions (pH 6.5−6.8),15 bacterial infections (pH 5−6),16 and particular organelles such as lysosomes (pH ∼ 4−5).17,18 Once we have a deep understanding of the factors that allow the loading and release of a drug at pH 3 and 7, respectively, it is conceivable that one can design functional group molecules that will load at pH 7 but release at the slightly acidic pH of tumors and bacterial infections. Potential for Other Applications. Apart from drug delivery applications, our model can be applied to calculate the adsorption of a single molecule on a surface or in fact any interaction between two molecules where there is acid−base behavior. One will now only have to input their four QMderived interaction energies, with their correct pKa values, into eq 6 in any graphical software such as gnuplot or excel, as used in our study, to obtain the plot of EpH int . First, we envisage that our method can be used for the calculation of the pH-responsive interaction of drugs and ligands with biomolecules such as the active sites of enzymes/ receptors. Generally, these interactions are studied for a single drug molecule and enzyme. QMMM calculations are usually performed where the active site is calculated as QM and the rest of the molecule is treated by a force field. As amino acid residues involved within enzyme active sites often have acidic or basic side chains,13,70 the method presented here should be useful. The present method is not limited to calculating interaction energies, but the pH response of any property can be determined by scaling the property by the probability of the combination. For example, the activation barrier in a reaction mechanism can be determined for the different combinations of protonation states and can be scaled by the probability of that



CONCLUSIONS We have provided a novel theoretical model, which combined with computational chemistry calculations will aid in understanding the probability and strength of interactions of two molecular species with acid−base behavior at different pH values. The relative probability of combinations of a pair of molecules with two different protonation states is derived (S0D−1, S−1D0, S0D0, and S−1D−1), and density functional theory simulations of each combination are then used to calculate the pH-responsive binding energy. In the present study, we applied the theoretical model for the prediction of the experimentally observed pH-dependent loading and release of the sulfasalazine drug molecules to positively charged TA-functionalized mesoporous silica nanoparticles, which has consequences for pH-responsive drug delivery. Overall, we found that our model, combined with periodic functional theory calculations on model TA-functionalized quartz surfaces, showed a strong agreement with experimental results that sulfasalazine was bound to MSN surfaces at low pH values but released at high pH values. Furthermore, in agreement with the experimental results, we 12287

DOI: 10.1021/acs.jpcc.8b02794 J. Phys. Chem. C 2018, 122, 12279−12290

Article

The Journal of Physical Chemistry C

(4) Hu, C.; Yu, L.; Zheng, Z.; Wang, J.; Liu, Y.; Jiang, Y.; Tong, G.; Zhou, Y.; Wang, X. Tannin as a Gatekeeper of pH-Responsive Mesoporous Silica Nanoparticles for Drug Delivery. RSC Adv. 2015, 5, 85436−85441. (5) Yoshida, T.; Lai, T. C.; Kwon, G. S.; Sako, K. pH- and IonSensitive Polymers for Drug Delivery. Expert Opin. Drug Delivery 2013, 10, 1497−1513. (6) Yuan, L.; Tang, Q.; Yang, D.; Zhang, J. Z.; Zhang, F.; Hu, J. Preparation of pH-Responsive Mesoporous Silica Nanoparticles and Their Application in Controlled Drug Delivery. J. Phys. Chem. C 2011, 115, 9926−9932. (7) Liu, J.-Q.; Li, X.-F.; Gu, C.-Y.; da Silva, J. C. S.; Barros, A. L.; Alves-Jr, S.; Li, B.-H.; Ren, F.; Batten, S. R.; Soares, T. A. A Combined Experimental and Computational Study of Novel Nanocage-Based Metal-Organic Frameworks for Drug Delivery. Dalton Trans. 2015, 44, 19370−19382. (8) DeMuth, P.; Hurley, M.; Wu, C.; Galanie, S.; Zachariah, M. R.; DeShong, P. Mesoscale Porous Silica as Drug Delivery Vehicles: Synthesis, Characterization, and pH-Sensitive Release Profiles. Microporous Mesoporous Mater. 2011, 141, 128−134. (9) Heffernan, M. J.; Murthy, N. Polyketal Nanoparticles: A New pH-Sensitive Biodegradable Drug Delivery Vehicle. Bioconjugate Chem. 2005, 16, 1340−1342. (10) Lee, C.-S.; Na, K. Photochemically Triggered Cytosolic Drug Delivery Using pH-Responsive Hyaluronic Acid Nanoparticles for Light-Induced Cancer Therapy. Biomacromolecules 2014, 15, 4228− 4238. (11) Gao, C.; Zheng, H.; Xing, L.; Shu, M.; Che, S. Designable Coordination Bonding in Mesopores as a pH-Responsive Release System. Chem. Mater. 2010, 22, 5437−5444. (12) Zhu, D.; Wen, H.-M.; Li, W.; Cui, X.-B.; Ma, L.; Kang, A. pHResponsive Drug Release from Porous Zinc Sulfide Nanospheres Based on Coordination Bonding. RSC Adv. 2014, 4, 33391−33398. (13) Bisswanger, H. Enzyme Assays. Perspect. Sci. 2014, 1, 41−55. (14) Liu, L.; Yao, W.; Rao, Y.; Lu, X.; Gao, J. pH-Responsive Carriers for Oral Drug Delivery: Challenges and Opportunities of Current Platforms. Drug Delivery 2017, 24, 569−581. (15) He, N.; Chen, Z.; Yuan, J.; Zhao, L.; Chen, M.; Wang, T.; Li, X. Tumor pH-Responsive Release of Drug-Conjugated Micelles from Fiber Fragments for Intratumoral Chemotherapy. ACS Appl. Mater. Interfaces 2017, 9, 32534−32544. (16) Padan, E.; Bibi, E.; Ito, M.; Krulwich, T. A. Alkaline pH Homeostasis in Bacteria: New Insights. Biochim. Biophys. Acta, Biomembr. 2005, 1717, 67−88. (17) Appelqvist, H.; Wäster, P.; Kågedal, K.; Ö llinger, K. The Lysosome: From Waste Bag to Potential Therapeutic Target. J. Mol. Cell Biol. 2013, 5, 214−226. (18) Zhou, J.; Tan, S.-H.; Nicolas, V.; Bauvy, C.; Yang, N.-D.; Zhang, J.; Xue, Y.; Codogno, P.; Shen, H.-M. Activation of Lysosomal Function in the Course of Autophagy Via mTORC1 Suppression and Autophagosome-Lysosome Fusion. Cell Res. 2013, 23, 508−523. (19) Bunker, A.; Magarkar, A.; Viitala, T. Rational Design of Liposomal Drug Delivery Systems, a Review: Combined Experimental and Computational Studies of Lipid Membranes, Liposomes and Their PEGylation. Biochim. Biophys. Acta, Biomembr. 2016, 1858, 2334− 2352. (20) Beg, S.; Rahman, M.; Jain, A.; Saini, S.; Midoux, P.; Pichon, C.; Ahmad, F. J.; Akhter, S. Nanoporous Metal Organic Frameworks as Hybrid Polymer−Metal Composites for Drug Delivery and Biomedical Applications. Drug Discovery Today 2017, 22, 625−637. (21) Henry, S. M.; El-Sayed, M. E. H.; Pirie, C. M.; Hoffman, A. S.; Stayton, P. S. pH-Responsive Poly(Styrene-Alt-Maleic Anhydride) Alkylamide Copolymers for Intracellular Drug Delivery. Biomacromolecules 2006, 7, 2407−2414. (22) Car, A.; Baumann, P.; Duskey, J. T.; Chami, M.; Bruns, N.; Meier, W. pH-Responsive PDMS-b-PDMAEMA Micelles for Intracellular Anticancer Drug Delivery. Biomacromolecules 2014, 15, 3235− 3245.

predict that there should be no binding of sulfasalazine to the nonfunctionalized silica surface at pH 4.5. We found that S0D0 and S−1D−1 contributed to binding and release at low and high pH, respectively, but we considered the S0D−1 and S−1D0 combinations as modulators because they affected the mid pH ranges and, therefore, the sensitivity and pH of the release profile. In turn, the main factors that affected the modulators were the difference in pKa values of the functionalized surface and drug molecule and the interaction energy of the modulators. As it should be possible to predict and test functional group/drug pairs against the aforementioned factors, we propose that our model will be instrumental in the chemical design of the pH-responsive drug delivery systems. The theoretical model and its implementation can be trivially adapted to various computational chemistry problems, and, therefore, we consider the work as a new way of approaching molecular simulations involving species with acid−base properties. We assert that our basic equations can be an integral part of any future molecular simulation study that considers how properties change as a function of pH, as well those that want to accurately simulate those properties at the desired pH.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.8b02794. Overall equation, statistics of binding, various figures describing the position of function group (FG) on the silica surface, varied positions of drug fractions on functionalized and nonfunctionalized silica surfaces, details of H-bonding observed, contribution of different combinations and the table indicating all of the interaction energies obtained (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel: +27 31-2608520, +27 31-26084. ORCID

Adam A. Skelton: 0000-0003-0155-8287 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work was supported by a grant from College of Health Sciences (CHS) at UKZN at Durban in South Africa. We are also grateful for the helpful support of CHPC (www.chpc.ac.za) and to UKZN CHPC cluster for providing computational resources.



REFERENCES

(1) Khandogin, J.; Brooks, I. C. L. Toward the Accurate FirstPrinciples Prediction of Ionization Equilibria in Proteins. Biochemistry 2006, 45, 9363−9373. (2) Socher, E.; Sticht, H. Mimicking Titration Experiments with MD Simulations: A Protocol for the Investigation of pH-Dependent Effects on Proteins. Sci. Rep. 2016, 6, No. 22523. (3) Swails, J. M.; York, D. M.; Roitberg, A. E. Constant pH Replica Exchange Molecular Dynamics in Explicit Solvent Using Discrete Protonation States: Implementation, Testing, and Validation. J. Chem. Theory Comput. 2014, 10, 1341−1352. 12288

DOI: 10.1021/acs.jpcc.8b02794 J. Phys. Chem. C 2018, 122, 12279−12290

Article

The Journal of Physical Chemistry C (23) Mao, J.; Li, Y.; Wu, T.; Yuan, C.; Zeng, B.; Xu, Y.; Dai, L. A Simple Dual-pH Responsive Prodrug-Based Polymeric Micelles for Drug Delivery. ACS Appl. Mater. Interfaces 2016, 8, 17109−17117. (24) Caminade, A.-M.; Turrin, C.-O. Dendrimers for Drug Delivery. J. Mater. Chem. B 2014, 2, 4055−4066. (25) Ambrogio, M. W.; Thomas, C. R.; Zhao, Y.-L.; Zink, J. I.; Stoddart, J. F. Mechanized Silica Nanoparticles: A New Frontier in Theranostic Nanomedicine. Acc. Chem. Res. 2011, 44, 903−913. (26) Ayad, M. M.; Salahuddin, N. A.; El-Nasr, A. A.; Torad, N. L. Amine-Functionalized Mesoporous Silica KIT-6 as a Controlled Release Drug Delivery Carrier. Microporous Mesoporous Mater. 2016, 229, 166−177. (27) Miriyala, N.; Ouyang, D.; Perrie, Y.; Lowry, D.; Kirby, D. J. Activated Carbon as a Carrier for Amorphous Drug Delivery: Effect of Drug Characteristics and Carrier Wettability. Eur. J. Pharm. Biopharm. 2017, 115, 197−205. (28) Yang, K.; Luo, H.; Zeng, M.; Jiang, Y.; Li, J.; Fu, X. Intracellular pH-Triggered, Targeted Drug Delivery to Cancer Cells by Multifunctional Envelope-Type Mesoporous Silica Nanocontainers. ACS Appl. Mater. Interfaces 2015, 7, 17399−17407. (29) Natarajan, J. V.; Nugraha, C.; Ng, X. W.; Venkatraman, S. Sustained-Release from Nanocarriers: A Review. J. Controlled Release 2014, 193, 122−138. (30) Datt, A.; El-Maazawi, I.; Larsen, S. C. Aspirin Loading and Release from MCM-41 Functionalized with Aminopropyl Groups via Co-Condensation or Postsynthesis Modification Methods. J. Phys. Chem. C 2012, 116, 18358−18366. (31) Slowing, I. I.; Trewyn, B. G.; Giri, S.; Lin, V. S. Y. Mesoporous Silica Nanoparticles for Drug Delivery and Biosensing Applications. Adv. Funct. Mater. 2007, 17, 1225−1236. (32) Jin, S.; Li, D.; Yang, P.; Guo, J.; Lu, J. Q.; Wang, C. Redox/pH Stimuli-Responsive Biodegradable PEGylated P(MAA/BACy) Nanohydrogels for Controlled Releasing of Anticancer Drugs. Colloids Surf., A 2015, 484, 47−55. (33) Rimola, A.; Costa, D.; Sodupe, M.; Lambert, J. F.; Ugliengo, P. Silica Surface Features and Their Role in the Adsorption of Biomolecules: Computational Modeling and Experiments. Chem. Rev. 2013, 113, 4216−4313. (34) Vallet-Regi, M.; Rámila, A.; del Real, R. P.; Pérez-Pariente, J. A New Property of MCM-41: Drug Delivery System. Chem. Mater. 2001, 13, 308−311. (35) Vallet-Regí, M.; Balas, F.; Arcos, D. Mesoporous Materials for Drug Delivery. Angew. Chem., Int. Ed. 2007, 46, 7548−7558. (36) Horcajada, P.; Rámila, A.; Pérez-Pariente, J.; Vallet-Regí, M. Influence of Pore Size of MCM-41 Matrices on Drug Delivery Rate. Microporous Mesoporous Mater. 2004, 68, 105−109. (37) Muñoz, B.; Ramila, A.; Perez-Pariente, J.; Diaz, I.; Vallet-Regi, M. MCM-41 Organic Modification as Drug Delivery Rate Regulator. Chem. Mater. 2003, 15, 500−503. (38) Heikkilä, T.; et al. Evaluation of Mesoporous TCPSi, MCM-41, SBA-15, and TUD-1 Materials as API Carriers for Oral Drug Delivery. Drug Delivery 2007, 14, 337−347. (39) Doadrio, A. L.; Sousa, E. M.; Doadrio, J. C.; Perez Pariente, J.; Izquierdo-Barba, I.; Vallet-Regi, M. Mesoporous SBA-15 HPLC Evaluation for Controlled Gentamicin Drug Delivery. J. Controlled Release 2004, 97, 125−132. (40) Leung, K.; Nielsen, I. M. B.; Criscenti, L. J. Elucidating the Bimodal Acid−Base Behavior of the Water−Silica Interface from First Principles. J. Am. Chem. Soc. 2009, 131, 18358−18365. (41) Liu, X.; Cheng, J.; Lu, X.; Wang, R. Surface Acidity of Quartz: Understanding the Crystallographic Control. Phys. Chem. Chem. Phys. 2014, 16, 26909−26916. (42) Kroutil, O.; Chval, Z.; Skelton, A. A.; Předota, M. Computer Simulations of Quartz (101)−Water Interface over a Range of pH Values. J. Phys. Chem. C 2015, 119, 9274−9286. (43) Lee, C.-H.; Lo, L.-W.; Mou, C.-Y.; Yang, C.-S. Synthesis and Characterization of Positive-Charge Functionalized Mesoporous Silica Nanoparticles for Oral Drug Delivery of an Anti-Inflammatory Drug. Adv. Funct. Mater. 2008, 18, 3283−3292.

(44) Fox, S. J.; Dziedzic, J.; Fox, T.; Tautermann, C. S.; Skylaris, C. K. Density Functional Theory Calculations on Entire Proteins for Free Energies of Binding: Application to a Model Polar Binding Site. Proteins 2014, 82, 3335−3346. (45) Adeowo, F. Y.; Honarparvar, B.; Skelton, A. A. Density Functional Theory Study on the Complexation of NOTA as a Bifunctional Chelator with Radiometal Ions. J. Phys. Chem. A 2017, 121, 6054−6062. (46) Adeowo, F. Y.; Honarparvar, B.; Skelton, A. A. The Interaction of NOTA as a Bifunctional Chelator with Competitive Alkali Metal Ions: A DFT Study. RSC Adv. 2016, 6, 79485−79496. (47) Skelton, A. A.; Agrawal, N.; Fried, J. R. Quantum Mechanical Calculations of the Interactions between Diazacrowns and the Sodium Cation: An Insight into Na+ Complexation in Diazacrown-Based Synthetic Ion Channels. RSC Adv. 2015, 5, 55033−55047. (48) Skelton, A. A.; Fried, J. R. The Insertion of Gas Molecules into Polyhedral Oligomeric Silsesquioxane (POSS) Cages: Understanding the Energy of Insertion Using Quantum Chemical Calculations. Phys. Chem. Chem. Phys. 2013, 15, 4341−4354. (49) Baik, M. H.; Friesner, R. A.; Lippard, S. J. Theoretical Study on the Stability of N-Glycosyl Bonds: Why Does N7-Platination Not Promote Depurination? J. Am. Chem. Soc. 2002, 124, 4495−4503. (50) Skelton, A. A.; Walsh, T. R. Interaction of Liquid Water with the Rutile TiO2 (110) Surface. Mol. Simul. 2007, 33, 379−389. (51) Mongan, J.; Case, D. A.; McCammon, J. A. Constant pH Molecular Dynamics in Generalized Born Implicit Solvent. J. Comput. Chem. 2004, 25, 2038−2048. (52) Aronson, J. N. The Henderson-Hasselbalch Equation Revisited. Biochem. Educ. 1983, 11, 68. (53) Bashford, D.; Karplus, M. Multiple-Site Titration Curves of Proteins: An Analysis of Exact and Approximate Methods for Their Calculation. J. Phys. Chem. 1991, 95, 9556−9561. (54) Antosiewicz, J.; Briggs, J. M.; Elcock, A. H.; Gilson, M. K.; McCammon, J. A. Computing Ionization States of Proteins with a Detailed Charge Model. J. Comput. Chem. 1996, 17, 1633−1644. (55) Skelton, A. A.; Fenter, P.; Kubicki, J. D.; Wesolowski, D. J.; Cummings, P. T. Simulations of the Quartz(101̅1)/Water Interface: A Comparison of Classical Force Fields, Ab Initio Molecular Dynamics, and X-Ray Reflectivity Experiments. J. Phys. Chem. C 2011, 115, 2076−2088. (56) Skelton, A. A.; Wesolowski, D. J.; Cummings, P. T. Investigating the Quartz (101̅0)/Water Interface Using Classical and Ab Initio Molecular Dynamics. Langmuir 2011, 27, 8700−8709. (57) Graham, G. G.; Pile, K. D. Sulfasalazine and Related Drugs. In Encyclopedia of Inflammatory Diseases; Parnham, M., Ed.; Springer: Basel, 2015; pp 1−5. (58) Cairns, D. Essentials of Pharmaceutical Chemistry, 2nd ed.; Pharmaceutical Press: London, 2003. (59) Marvinsketch. https://chemaxon.com/products/marvin (accessed Nov 10, 2017). (60) McCubbin, J. A.; Voth, S.; Krokhin, O. V. Mild and Tunable Benzoic Acid Catalysts for Rearrangement Reactions of Allylic Alcohols. J. Org. Chem. 2011, 76, 8537−8542. (61) Baba, T.; Matsui, T.; Kamiya, K.; Nakano, M.; Shigeta, Y. A Density Functional Study on the pKa of Small Polyprotic Molecules. Int. J. Quantum Chem. 2014, 114, 1128−1134. (62) Bolt, G. H. Determination of the Charge Density of Silica Sols. J. Phys. Chem. 1957, 61, 1166−1169. (63) Ong, S.; Zhao, X.; Eisenthal, K. B. Polarization of Water Molecules at a Charged Interface: Second Harmonic Studies of the Silica/Water Interface. Chem. Phys. Lett. 1992, 191, 327−335. (64) Clark, S. J.; Segall, M. D.; Pickard, C. J.; Hasnip, P. J.; Probert, M. I. J.; Refson, K.; Payne, M. C. First Principles Methods Using CASTEP. Z. Kristallogr. - Cryst. Mater. 2005, 220, 567−570. (65) Accelrys, Materials Studio, 2016. http://accelrys.com/products/ collaborative-science/biovia-materials-studio/Molecular (accessed Nov 19, 2017). (66) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. 12289

DOI: 10.1021/acs.jpcc.8b02794 J. Phys. Chem. C 2018, 122, 12279−12290

Article

The Journal of Physical Chemistry C (67) Meyer, B.; Rabaa, H.; Marx, D. Water Adsorption on ZnO(1010): From Single Molecules to Partially Dissociated Monolayers. Phys. Chem. Chem. Phys. 2006, 8, 1513−1520. (68) Vanderbilt, D. Soft Self-Consistent Pseudopotentials in a Generalized Eigenvalue Formalism. Phys. Rev. B 1990, 41, 7892−7895. (69) Kubicki, J. D.; Sofo, J. O.; Skelton, A. A.; Bandura, A. V. A New Hypothesis for the Dissolution Mechanism of Silicates. J. Phys. Chem. C 2012, 116, 17479−17491. (70) Iqbal, A.; Clifton, I. J.; Chowdhury, R.; Ivison, D.; Domene Nunez, C.; Schofield, C. J. Structural and Biochemical Analyses Reveal How Ornithine Acetyl Transferase Binds Acidic and Basic Amino Acid Substrates. Org. Biomol. Chem. 2011, 9, 6219−6225. (71) Bauler, P.; Huber, G.; Leyh, T.; McCammon, J. A. Channeling by Proximity: The Catalytic Advantages of Active Site Colocalization Using Brownian Dynamics. J. Phys. Chem. Lett. 2010, 1, 1332−1335. (72) Huang, Q.; Herrmann, A. Calculating pH-Dependent Free Energy of Proteins by Using Monte Carlo Protonation Probabilities of Ionizable Residues. Protein Cell 2012, 3, 230−238.

12290

DOI: 10.1021/acs.jpcc.8b02794 J. Phys. Chem. C 2018, 122, 12279−12290