Taxifolin Adsorption on Nitrogenated Graphenes: Theoretical Insights

: Solid-state drug delivery systems for the drug substances transport are of great importance nowadays. In the present work, the non-covalent interactions between taxifolin ( Tax ) and graphene as well as nitrogenated (N-doped) graphenes were systematically studied by using a wide set of theoretical techniques. Symmetry-adapted perturbation theory (SAPT0) calculations confirmed more favorable adsorption of Tax on N-doped graphenes compared to pristine graphene. It was established that dispersion interactions play the main role in the attractive interactions (>60%), whereas electrostatic and induction forces contribute only moderately to the attraction (~25% and 7–8%, respectively). Independent gradient model (IGM) analysis visually demonstrated the existence of dispersion interactions and hydrogen bonding in the studied Tax complexes. Ab initio molecular dynamics calculations indicated stability of these complexes at different temperatures. Our results show that N-doped graphenes with the enhanced interaction energy (E int ) toward Tax are promising candidates for the technical realization of the targeted drug delivery systems.


Introduction
In the last decade, scientific efforts have focused on the development of novel drug delivery systems (DDSs) for the transport of drug substances with marginal negative effects on the human organism.The need for highly efficient, flexible and controlled loading/release systems is driven by significant advances in patient compliance, clinical efficacy, and extended drug shelf life through controlled drug release.For these reasons, new DDSs have the potential to become one of the fastest growing segments of the pharmaceutical industry.It is believed that such systems can help one to overcome the limitation of bioavailability of different drugs and increase their efficiency for the introduction of novel therapies and theranostics.Current DDSs are based on nanostructural materials, and have been shown to offer substantial advantages, such as solubility, bioavailability, pharmacological activity, targeted delivery, and degradation [1][2][3].
Taxifolin (Tax, dihydroquercetin) is a natural compound that belongs to the class of flavonoids.Much attention has been paid to Tax owing to its important pharmacological properties [14].It is a powerful antioxidant with a pronounced effect in the prevention of several forms of cancer.Recently, it has shown promising suppressive activity against inflammation, malignancies, microbial infections, and cardiovascular and liver diseases.In present time, the protective role of Tax as a COVID-19 recovery is being evaluated (see, e.g., the recent review and the references therein [14]).Different works reported a successful Tax release performance by zinc oxide nanoparticles [15], collagen-acrylic hydrogels [16], zein-caseinate nanoparticles [17], and cyclodextrine [18].However, the papers dealing with interactions between Tax molecules and graphene derivatives are rather scarce.The bright examples are the work of Zhao et al. on quercetin, luteolin, and kaempferol adsorption on graphene oxide (GO) [19] and the paper of Tiwari et al. on quercetin adsorption on GO [20], and the report of Garcia et al. on a family of selected flavonoids and pristine graphene [21].
This work investigates adsorption properties of nanostructured platforms consisting of nitrogenated (N-doped) graphenes as nanocarriers loaded with the Tax molecule.Ndoped carbon materials are among the most efficient metal-free catalysts [22].It was found earlier that N-doping of graphene can slightly increase the activity of the oxygen reduction reaction, but N-doping at the graphene interface does not contribute significantly to the activity enhancement [23].Also, N-doped graphenes have been studied for electrochemical sensing applications owing to their ability to enhance conductivity [24].In addition to various catalytic applications, other applications can be found in medicine as N-doped graphene materials exhibit reduced cytotoxicity and antioxidant effects in a nitrogen content-dependent manner [25].Moreover, N-doped graphene materials show less toxicity compared to GO [26].
More specifically, our aim is to perform a thorough comparative analysis of adsorption properties of graphene (G) and N-graphenes (pyrrolic (Npyrr-G), pyridinic (Npyrid-G) and quaternary graphene (2N-G)) toward Tax by using density functional theory (DFT), symmetry adapted perturbation theory (SAPT0) [27,28], the quantum theory of atoms-inmolecules (QTAIM) [29], and the independent gradient model (IGM) method [30,31].Ab initio molecular dynamics (AIMD) simulations reveal the behavior of all complexes at room and liquid nitrogen temperatures.We suggest that the present theoretical study reveals the features of Tax/adsorbent interactions, indicating a promising future for nanoscale DDSs.

Computational Details
Geometries of the studied Tax complexes were optimized at the BLYP-D3/def2-SVP [32][33][34] level using the Orca 5.0.3 package [35].Grimme's dispersion correction (D3) was employed to properly treat van der Waals interactions [36].The BLYP-D3 functional has been shown to provide reliable results on non-covalent complexes [37,38].For each of these complexes, the total interaction energy (E int ) was calculated using the 0th-order symmetry-adapted perturbation theory (SAPT0) method [27] with the jun-cc-pVDZ basis set [39,40], which allows decomposing E int to the exchange (E ex ), electrostatic (E el ), dispersion (E disp ), and induction (E ind ) components.It was established that such a combination of methods accurately predicts adsorption energies.SAPT0 calculations were carried out using the Psi4 code (v.1.4) [41].It should be noticed that the more negative E int value for a given complex denotes the stronger adsorption.The simplest SAPT0 method is defined in the following equation [28]: In this notation, E (vw) defines the order in V and in W A + W B in the Hamiltonian (H): where the H is written as a sum of the usual monomer Fock operators (F), the fluctuation potential of each monomer (W), and the interaction potential (V); the subscript (resp) indicates that orbital relaxation effects are included.
The δ HF term takes into account the higher-order induction effects, and it is included in the definition of SAPT terms.
Normally, the adsorption energy terms are calculated as following: where We should briefly describe the components of the E int energy term [42].The electrostatic effect (E el ) describes the classical interactions between the static charge distributions of the interacting species.It may be either attractive (−) or repulsive (+).Induction effects (E ind ) arise from the distortion of a molecule in response to the electric field of all neighbor molecules, and they are always attractive (−).Dispersion interactions (E disp ) arise owing to instantaneous fluctuations of electron distributions of the molecules, and they are also attractive (−).The Pauli or exchange (E ex ) repulsion can be described from a classical point of view.The atomic approach results in a decrease in electron density in the internuclear region.As a result, the decreased screening of nuclear charges occurs, which leads to the increased nuclear-nuclear repulsion (+).
The studied model of graphene (G) consists of 48 carbon atoms.Three major cases of N-doped graphenes, namely pyrrolic (Npyrr-G), pyridinic (Npyrid-G) and quaternary graphene (2N-G) were also used as adsorbents.These models were chosen according to the experimental work on N-graphenes [43].Their stability has been reported elsewhere [44].The peripheral rims of all adsorbent models were saturated with hydrogen atoms to avoid dangling.The nanosheets of similar sizes have been recently employed as models for different adsorption studies [44,45].
The procedure for energy calculations was as follows.Firstly, to obtain the minimum energy structure, the Tax molecule was placed above the adsorbent sheet at a height corresponding the physisorbed state (~3.0 Å) and allowed to relax freely.The mutual initial configurations were randomized (see Figures S1-S4).The minimum energy configuration was used as the input configuration for SAPT0 calculations.
The initial structures for ab initio molecular dynamics (AIMD) simulations were also relaxed using the BLYP/def2-SVP method.All AIMD simulations used a time step of 1 fs, an NVT ensemble using a Nose-Hoover thermostat, and a velocity Verlet algorithm for integrating the equation of motion.The time period of 1000 fs has been used for the full AIMD computation cycle.These calculations were also done using the Orca 5.0.3 program suite.
The Multiwfn 3.8.0 program [46] has been used to obtain the independent gradient model (IGM) [30,31] isosurfaces and as well as obtaining the QTAIM parameters.VMD (v.1.9.3) [47] and Chemcraft [48] programs have been used for visualization aims.The detailed description of the IGM analysis has been given in the Supplementary Materials.

Results and Discussions
The models used for Tax, graphene and N-graphenes are demonstrated in Figures 1 and 2. Table 1 shows the calculated interaction energy (E int ) and its constituents for Tax adsorption obtained at the SAPT0/jun-cc-pVDZ level of theory.To figure out the possible influence of the size of the adsorbent on the E int values, we consider larger graphene sheets compared to the studied model.Generally, the widely used graphene model (circumcoronene) consists of 54 C atoms [49].Also, we involve into consideration the graphene model consisting of 80 C atoms (Figure S5).The total E int values are calculated to be −31.05kcal/mol for the former and −33.72 kcal/mol for the latter.These values are in a good accordance with that for Tax adsorption on G presented herein (−30.04 kcal/mol).We, hence, proceed with the initial model.
kcal/mol for the former and −33.72 kcal/mol for the latter.These values are in a goo cordance with that for Tax adsorption on G presented herein (−30.04 kcal/mol).We, h proceed with the initial model.

Tax
EPM of Tax  kcal/mol for the former and −33.72 kcal/mol for the latter.These values are in a good accordance with that for Tax adsorption on G presented herein (−30.04 kcal/mol).We, hence, proceed with the initial model.First, we analyze the E int energy term and its components for Tax adsorption.The magnitude of E int for pristine graphene is the smallest one.At the same time, the gradual increase in E int from Npyrr-G (−32.63 kcal/mol) to Npyrid-G (−37.16 kcal/mol) and further to 2N-G (−43.42 kcal/mol) is observed.It is of note, that the last case represents a much larger E int value compared to E int for graphene (~40%).To study these differences in detail, we now consider the electrostatic, inductive and dispersion contributions to the total attraction energy for all adsorbents.First, the marked increase in E el can be observed for N-graphenes, except for the case of Npyrr-G (Table 1).The most significant E el value of −28.34 kcal/mol is in the case of Tax/2N-G interactions.Besides this, only the 2N-G structure shows an increase in the relative contribution of E el to E int (26.9%).

EPM of Tax
We use electrostatic potential maps (EPMs) to further study electrostatic contributions in non-covalent interactions between the four adsorbents and Tax (Figures 1 and 2).Blue areas characterize the abundance of electrons (negative charge), whereas red areas exhibit electron depletion (positive charge).For Tax, the negatively charged regions can be observed at O atoms, whereas positively charged regions are located at H atoms of OH groups (Figure 1). Figure 2 shows that, in contrast to the pristine G structure, all N-graphenes show non-uniform EPMs.2N-G is characterized by the electron-depleted central area and electron-rich peripheral regions.At the same time, pyrrolic and pyridinic N atoms introduce marked changes at the edge of the Npyrr-G and Npyrid-G models, respectively.At a qualitative level, we may assume that the substituted benzene rings of Tax would interact with graphenes via π-π stacking.The OH-groups of Tax may preferably interact with the π-system of G and N-graphenes through H-bonding.As a whole, non-uniform EPMs of N-graphenes and especially of 2N-G will enhance intermolecular electrostatic interactions.These facts support our results obtained from SAPT0 calculations, which exhibit notable contributions of such interactions to the total attraction.The Mulliken charge transfer (CT) is calculated to be 0.031, 0.050, 0.064, and 0. 025 e-towards G, 2N-G,  Npyrid-G, and Npyrr-G, respectively.The larger CT values are attributed to the cases of larger interactions between Tax and N-graphenes ( 2N-G and Npyrid-G).
The independent gradient model (IGM) method allows the types of non-covalent interactions between interacting species to be obtained through a simple analysis of colored 3D isosurfaces.The isosurfaces obtained in the framework of IGM indicate the existence of two types of interactions between Tax and four adsorbents studied herein.
According to the IGM method, we can easily conclude that the stronger interactions (blue-colored areas of isosurfaces) are hydrogen bonding (H-bonding), and the weaker binding (green-colored areas of isosurfaces) is due to van der Waals interactions (Figure 3).The vast areas of green-colored isosurfaces are consistent with the results of SAPT0 calculations.This approves the importance of dispersion interactions for the favorable accommodation of Tax molecules on the surfaces of pristine graphene and N-graphenes.For the further examination of intermolecular interactions, we involve the Bader's QTAIM analysis [29].According to the Koch and Popelier's QTAIM criteria for the existence of H-bonding, the following conditions should be fulfilled: first, the bond critical point (BCP) (3,−1) for the proton/acceptor contact exists; second, the electron density ρ(r) value for BCP (3,−1) should be in the range of 0.002 to 0.034 a.u.; and third, the corresponding Laplacian,  2 ρ(r), should be in the range from 0.024 to 0.139 a.u.[50].Additionally, there are factors indicating the existence of non-covalent interactions, and they are the following: (i) −G(r)/V(r) > 1 and (ii) the positive value of H(r).We analyzed the voids located in between the interacting Tax, graphene, and N-graphenes, and selected the corresponding BCPs (3,−1).The QTAIM data on BCPs for possible H-bonds between Tax and G are summarized in Table 2.The QTAIM parameters for Tax/2N-G, Tax/Npyrr-G, and Tax/Npyrid-G are given in Tables 3-5, respectively.We verify three, two, four and one intermolecular H-bonds in the Tax/G, Tax/2N-G, Tax/Npyrr-G, and Tax/Npyrid-G complexes, respectively.Additionally, the existence of non-covalent bonding has been proved.To sum up, according to the QTAIM analysis, all the complexes exhibit both non-covalent interactions between species and H-bonding.No indicators of covalent bonding have been found.Green isosurfaces denote weak van der Waals interactions, whereas blue denote strong attractive interactions.Atomic color code is the same as in Figure 2.
According to the Koch and Popelier's QTAIM criteria for the existence of H-bonding, the following conditions should be fulfilled: first, the bond critical point (BCP) (3,−1) for the proton/acceptor contact exists; second, the electron density ρ(r) value for BCP (3,−1) should be in the range of 0.002 to 0.034 a.u.; and third, the corresponding Laplacian, ∇ 2 ρ(r), should be in the range from 0.024 to 0.139 a.u.[50].Additionally, there are factors indicating the existence of non-covalent interactions, and they are the following: (i) −G(r)/V(r) > 1 and (ii) the positive value of H(r).We analyzed the voids located in between the interacting Tax, graphene, and N-graphenes, and selected the corresponding BCPs (3,−1).The QTAIM data on BCPs for possible H-bonds between Tax and G are summarized in Table 2.The QTAIM parameters for Tax/2N-G, Tax/Npyrr-G, and Tax/Npyrid-G are given in Tables 3-5, respectively.We verify three, two, four and one intermolecular H-bonds in the Tax/G, Tax/2N-G, Tax/Npyrr-G, and Tax/Npyrid-G complexes, respectively.Additionally, the existence of non-covalent bonding has been proved.To sum up, according to the QTAIM analysis, all the complexes exhibit both non-covalent interactions between species and H-bonding.No indicators of covalent bonding have been found.The next component of interest is the induction interactions.The E ind term contributes only marginally to the total attraction (~7-8%).Although absolute values of E ind increase from G to N-graphenes (−4.80 (G), −7.46 (Npyrid-G), −6.21 (Npyrr-G), and −9.18 (2N-G) kcal/mol), its relative contribution remains nearly the same (we can see only small increase of 1-2%) (Table 1).It witnesses that permanent dipole moments of N-graphenes (3.18, 1.95, and 3. 28 D for 2N-G, Npyrr-G, and Npyrid-G, respectively) affect the strength of intermolecular interactions, and it leads to more favorable Tax adsorption.These effects, however, are rather small, and the E ind term contributions are smaller than those of E el .
We now turn to the discussions on London dispersion interactions, which contribute significantly to attraction in all cases (see Table 1 and Figure 3).For G and N-graphenes, the E disp relative contributions are of roughly the same magnitude, but there is a striking difference in their absolute values.Indeed, E disp increases by absolute value from −46.32 (G) to −68.02 (2N-G) kcal/mol.Moreover, its relative contribution to the attractive interactions is very large (>60%), and the order of the energy constituents is as follows: E disp > E el > E ind .From the analysis of the relative inputs of energy terms to attraction, it becomes clear that the dispersion term is roughly several per cent smaller in the 2N-G case and that this is compensated for by a simultaneous increase in the induction and electrostatic contributions.
Our previous work studied Tax interactions with a fragment of the arabinogalactan molecule [51].In contrast to the present report, the E el and E disp terms were of almost equal importance for the stabilization.In addition, we detect numerous H-bonds between two interacting species.Another theoretical work studied adsorption of different flavonoids (quercetin, luteolin, and kaempferol) on graphene oxide [19].It was established that the main driving force for favorable adsorption are π-π interactions and H-bonding.As a whole, our conclusions are similar to those obtained previously [19], but the SAPT0 method is not able to separate H-bonding from E int .Instead of that, the substantial portion of H-bonding is included in the E ind term.However, in our case, the E ind terms are rather small, and, supposedly, the E disp dominates.
We further compare behavior of all studied complexes at different temperatures (T = 77 K (liquid nitrogen temperature) and T = 300 K (room temperature)) by using results of ab initio molecular dynamics (AIMD) simulations.Although both classical MD [52] with empirically derived forces and AIMD with forces derived from DFT calculations can be suitable for this purpose, we select the AIMD method.It becomes popular owing to its advantages, e.g., transferability between systems and ability to accurately describe polarization and charge transfer [53].For the Tax/G complex and T = 77 K, we observe only minor changes in the mutual configuration of the complex (Figure 4).The arrangement of Tax and G species is preserved during the whole period of simulation.Only small changes in the framework of G can be observed.For T = 300 K, the larger changes of the framework of G can be noticed, and the sheet of G acquires the corrugated shape.At the same time, the distances between Tax and G are nearly equal during the simulation time.
Solids 2024, 5, FOR PEER REVIEW 8 main driving force for favorable adsorption are π-π interactions and H-bonding.As a whole, our conclusions are similar to those obtained previously [19], but the SAPT0 method is not able to separate H-bonding from Eint.Instead of that, the substantial portion of H-bonding is included in the Eind term.However, in our case, the Eind terms are rather small, and, supposedly, the Edisp dominates.We further compare behavior of all studied complexes at different temperatures (T = 77 K (liquid nitrogen temperature) and T = 300 K (room temperature)) by using results of ab initio molecular dynamics (AIMD) simulations.Although both classical MD [52] with empirically derived forces and AIMD with forces derived from DFT calculations can be suitable for this purpose, we select the AIMD method.It becomes popular owing to its advantages, e.g., transferability between systems and ability to accurately describe polarization and charge transfer [53].For the Tax/G complex and T = 77 K, we observe only minor changes in the mutual configuration of the complex (Figure 4).The arrangement of Tax and G species is preserved during the whole period of simulation.Only small changes in the framework of G can be observed.For T = 300 K, the larger changes of the framework of G can be noticed, and the sheet of G acquires the corrugated shape.At the same time, the distances between Tax and G are nearly equal during the simulation time.A similar situation exists for the Tax/2N-G complex and T = 77 K (Figure 5).We may observe only marginal changes in both the molecular geometry and the distances between Tax and 2N-G.For T = 300 K, the 2N-G sheet acquires a slightly curved shape, but generally smaller changes compared to the previous case can be noticed.In all considered cases, the Tax molecule is in the physisorbed state throughout the simulation.
For the Tax/Npyrr-G complex and T = 77 K, we detect very small mutual reorientations of Tax and Npyrr-G species (Figure 6).For the higher T = 300 K, the displacement of Tax is slightly larger, but it generally retains its initial location.The framework of Npyrr-G is almost intact in these two cases.For the Tax/Npyrid-G complex (Figure 7), we can detect very small displacements for the lower T = 77 K, and slightly larger displacements/distortions for the case of T = 300 K.No sufficient shifts of Tax from the initial positions can be found for Tax/Npyrid-G simulations.The distances between the center-ofmass (COM) of Tax and studied N-graphenes are summarized in Table S1.They are determined as the perpendicular line dropped from the COM to the underlying surface.While the simulation is proceeding, changes of the distances are taking place in studied systems.The comparison of the distances shows that the travelling of Tax is much A similar situation exists for the Tax/2N-G complex and T = 77 K (Figure 5).We may observe only marginal changes in both the molecular geometry and the distances between Tax and 2N-G.For T = 300 K, the 2N-G sheet acquires a slightly curved shape, but generally smaller changes compared to the previous case can be noticed.In all considered cases, the Tax molecule is in the physisorbed state throughout the simulation.
For the Tax/Npyrr-G complex and T = 77 K, we detect very small mutual reorientations of Tax and Npyrr-G species (Figure 6).For the higher T = 300 K, the displacement of Tax is slightly larger, but it generally retains its initial location.The framework of Npyrr-G is almost intact in these two cases.For the Tax/Npyrid-G complex (Figure 7), we can detect very small displacements for the lower T = 77 K, and slightly larger displacements/distortions for the case of T = 300 K.No sufficient shifts of Tax from the initial positions can be found for Tax/Npyrid-G simulations.The distances between the center-of-mass (COM) of Tax and studied N-graphenes are summarized in Table S1.They are determined as the perpendicular line dropped from the COM to the underlying surface.While the simulation is proceeding, changes of the distances are taking place in studied systems.The comparison of the distances shows that the travelling of Tax is much smaller in the case of T = 77 K, and the larger distance can be observed in the case of Npyrr-G (>4 Å).As a whole, the AIMD analysis confirms the stable behavior of Tax adsorbed on both pristine graphene and N-doped graphenes.employed in this work as the second solvent as this mixture allows Tax to be dissolved.To mimic this solvent, we add the 3:3 ratio of water and ethanol molecules.For the aim of comparison, gas phase calculations were also made.Now we calculate the E int,s values based on the supermolecular approach (Table 6).
where E comp is the energy of the studied complex, E Tax is the energy of Tax, E N-graphene is the energy of corresponding N-graphene, and E solv is the energy of solvent molecules used.When being involved, the solvents lead to a substantial stabilization of all complexes (Table 6).The E int,s change varies from a three (2N-G) to an almost five times (Npyrr-G) increase.We link this behavior with the formation of numerous H-bonds between the studied species and solvent molecules.This indicates the applicability of studied adsorbents as Tax DDSs.It should be mentioned, though, that the large magnitudes of interaction energies may complicate the drug release.To avoid this possible problem, one should modify the graphene surface or select the solvent with care.Recent papers describe the endogenous and exogenous stimuli-responsive drug delivery [56,57].The former includes pH-sensitive, redox-responsive, and ROS-responsive (reactive oxygen species) drug delivery and release.The latter deals with NIR-responsive, thermo-responsive, and electro-responsive drug delivery.Yang et al. report the strong pH dependence of loading and release of doxorubicin hydrochloride on the graphene oxide surface.They claimed that such a dependence may be caused by the stronger H-bonding interactions under basic conditions than those under acid conditions [58].We assume that the similar technique can be suitable for Tax.
To summarize, based on the present results and those obtained elsewhere, we highlight the possible ways of use for N-doped graphene.Theoretical investigations of the non-

Figure 3 .
Figure 3.The IGM isosurfaces (isovalue = 0.005) for Tax adsorption on all studied adsorbents.Green isosurfaces denote weak van der Waals interactions, whereas blue denote strong attractive interactions.Atomic color code is the same as in Figure 2.

Figure 3 .
Figure 3.The IGM isosurfaces (isovalue = 0.005) for Tax adsorption on all studied adsorbents.Green isosurfaces denote weak van der Waals interactions, whereas blue denote strong attractive interactions.Atomic color code is the same as in Figure2.

Figure 4 .
Figure 4. AIMD simulations of the Tax/G complex.Atomic color code: hydrogen-light blue, oxygen-red, and carbon-yellow.

Figure 4 .
Figure 4. AIMD simulations of the Tax/G complex.Atomic color code: hydrogen-light blue, oxygen-red, and carbon-yellow.

Table 1 .
E int and its constituents (kcal/mol) for interactions between Tax and the studied adsorbents.The values in parentheses denote the percentage contribution to attractive energy.

Table 6 .
E int,s (kcal/mol) for interactions between Tax and the studied adsorbents obtained via the supramolecular approach.