e-Graphene: A Computational Platform for the Prediction of Graphene-Based Drug Delivery System by Quantum Genetic Algorithm and Cascade Protocol
Suqing Zheng1,2* 
Fu Lin- 1School of Pharmaceutical Sciences, Wenzhou Medical University, Wenzhou, China
- 2Chemical Biology Research Center, Wenzhou Medical University, Wenzhou, China
- 3Institute of Frontier and Interdisciplinary Science, Shandong University, Qingdao, China
- 4Center of Chemical Biology, Guangzhou Institute of Biomedicine and Health, Chinese Academy of Sciences, Guangzhou, China
Graphene, as a novel category of carbon nanomaterials, has attracted a great attention in the field of drug delivery. Due to its large dual surface area, graphene can efficiently load drug molecules with high capacity via non-covalent interaction without chemical modification of the drugs. Hence, it ignites prevalent interests in developing a new graphene/graphene oxide (GO)-based drug delivery system (GDDS). However, current design of GDDS primarily depends on the prior experimental experience with the trial-and-error method. Thus, it is more appealing to theoretically predict possible GDDS candidates before experiments. Toward this end, we propose to fuse quantum genetic algorithm (QGA) and quantum mechanics (QM)/semi-empirical quantum mechanics (SQM)/force field (FF) to globally search the optimal binding interaction between the graphene/GO and drug in a given GDDS and develop a free computational platform “e-Graphene” to automatically predict/screen potential GDDS candidates. To make this platform more pragmatic for the rapid yet relatively accurate prediction, we further propose a cascade protocol via firstly conducting a fast QGA/FF calculation with fine QGA parameters and automatically passing the best chromosomes from QGA/FF to initialize a higher level QGA/SQM or QGA/QM calculation with coarse QGA parameters (e.g., small populations and short evolution generations). By harnessing this platform and protocol, systematic tests on a typical GDDS containing an anticancer drug SN38 illustrate that high fabrication rates of hydroxyl, epoxy, and carboxyl groups on a pristine graphene model will compromise the stability of GDDS, implying that an appropriate functionalization rate is crucial for the delicate balance between the stability and solubility/biocompatibility of GDDS. Moreover, automatic GDDS screen in the DrugBank database is performed and elicits four potential GDDS candidates with enhanced stability than the commonly tested GDDS containing SN38 from the computational point of view. We hope that this work can provide a useful program and protocol for experimental scientists to rationally design/screen promising GDDS candidates prior to experimental tests.
Introduction
Pristine graphene, discovered by Geim and Novoselov, is a two-dimensional (2D) nanomaterial consisting of a single layer of carbon atoms (Novoselov et al., 2004). It has gained tremendous interests in various research areas (e.g., energy storage and sensors) due to its unique and impressive electrical, thermal, and mechanical properties (Novoselov et al., 2012; Lightcap and Kamat, 2013; Wassei and Kaner, 2013; El-Kady et al., 2016; Georgakilas et al., 2016; Aditya et al., 2017; Yu et al., 2017). Meanwhile, graphene has large dual surface area and high loading capacity, excellent chemical and mechanical stability, and good solubility and biocompatibility after simple fabrication on the surface, consequently graphene and its common derivatives, such as graphene oxides (GO) have exhibited considerable potential in the drug delivery (Chung et al., 2013; Mao et al., 2013; Yang et al., 2013; Goenka et al., 2014; Xianfeng and Feng, 2014; Reina et al., 2017; Ghosal and Sarkar, 2018; Mohammad Omaish et al., 2019).
In 2008, Dai group reported the first case of using GO as a nanocarrier to transport an anticancer drug SN38 (Liu et al., 2008). In their work, they demonstrated that SN38 indeed binds non-covalently with GO to form a stable GO/SN38 complex that affords remarkable activity with IC50 of about 6 nM for HCT116 cell line; hence, this pioneering work opens up a novel application of graphene/GO in the field of drug delivery system. Since then, more and more research works on the non-covalent graphene/GO-based drug delivery system (GDDS) are surging and summarized in the comprehensive reviews (Sun et al., 2008; Yang et al., 2008, 2016; Pan et al., 2012; Liu et al., 2013; McCallion et al., 2016; Shim et al., 2016; Zhang et al., 2017; Yi et al., 2019). However, current design of GDDS still highly relies on the traditional trial-and-error approach. Hence, it will be more efficient if GDDS can be evaluated or predicted with the computational methods before experiments.
Regarding the computational predictions of GDDS, molecular dynamics (MD) simulation and quantum mechanics (QM) programs are powerful to accurately investigate the non-bonded interaction between the graphene/GO and ligand for a given GDDS (Gráfová et al., 2010; Ramraj and Hillier, 2010; Calero et al., 2013; Cho et al., 2013; Guo et al., 2013; Vovusha et al., 2013, 2018; Mudedla et al., 2014; Vincent and Hillier, 2014; Wang et al., 2015; Mahdavi et al., 2016, 2020; Krepel and Hod, 2017; Safdari et al., 2017; Ajala et al., 2019; Alkathiri et al., 2019; Azhagiya Singam et al., 2019; Mason et al., 2019). However, both MD and QM methods are usually time-consuming, need to know a reasonable initial binding mode between the graphene/GO and ligand before calculations and tend to entrap the graphene/GO and ligand in a local minimum. To increase the chance of escaping from the local minimum, meta-dynamics (MTD) simulation method (Laio and Parrinello, 2002) is one of the most popular enhanced sampling methods (Bernardi et al., 2015) and extensively employed to promote the crossing of high energetic barriers via gradual additions of the Gaussian potential over time. According to the MTD method, the Grimme group developed a very useful conformation search program CREST (Bannwarth et al., 2019), where the root-mean-square deviation (RMSD) of a small molecule is adopted as a reaction coordinate (RC) in the MTD simulation, but this RC cannot capture the relative position between the small molecule and graphene/GO. Alternatively, semi-classical MD simulation method was proposed by the Ceotto group to explore various minima on a potential energy surface, which also is not a computationally cheap method and mainly developed for the prediction of molecular spectroscopy (Conte and Ceotto, 2020; Gandolfi et al., 2020). As a result, these methods may be not very suitable for the rapid prediction/screen of GDDS.
Alternatively, the fast molecular docking programs, such as AutoDock (Morris et al., 2009), which were originally designed to globally search the optimal conformations of a small molecule within the pocket of a biological target, can be borrowed for the quick docking of a small molecule to a graphene/GO model if appropriate parameters for the graphene/GO can be provided. Nevertheless, native scoring functions in these docking programs, initially derived from protein/ligand binding affinities, are also not dedicated to the prediction of drug/graphene (or GO) binding interaction. In short, these molecular docking programs seem not specialized for the accurate prediction/screen of GDDS. Therefore, a single method, such as pure MD, QM, or molecular docking, probably is not the best choice for the fast yet relatively accurate prediction/screen of GDDS, and combining the advantages of some of these methods would be more promising.
So far there are no dedicated solutions for the prediction/screen of GDDS. Toward this aim, in this work, we propose to fuse quantum genetic algorithm (QGA) (Narayanan and Moore, 1996) and QM/semi-empirical quantum mechanics (SQM)/force field (FF) to predict the optimal drug/graphene (or GO) interaction and further develop a specialized computational platform “e-Graphene” (Figure 1) to predict/screen the potential GDDS in an automatic fashion, which was implemented on the basis of our previous visual FF derivation toolkit “VFFDT” (Zheng et al., 2016). To well balance the speed and accuracy of GDDS prediction/screen, we further suggest an efficient cascade protocol: (1) firstly perform a cheap QGA/FF calculation with fine QGA parameters (e.g., larger populations and longer evolution generations); and (2) automatically pass the best chromosomes from QGA/FF to evoke a higher level QGA/SQM or QGA/QM calculation with coarse QGA parameters (e.g., smaller populations and shorter evolution generations) and also implement this protocol in e-Graphene. To the best of our knowledge, the QGA method and cascade protocol were never used to predict/screen GDDS.
Figure 1. The main interface of e-Graphene program for the automatic prediction/screen of graphene-based drug delivery system (GDDS) by quantum genetic algorithm (QGA) and cascade protocol.
Based on this proposed computational platform and protocol, we organize “Methods” and “Discussion” sections as follows. First, the implementation of main functions in e-Graphene program is briefly introduced in “Methods” section. Second, extensive tests on an exemplary GDDS containing an anticancer drug SN38, which was reported in Liu et al. (2008), are conducted to investigate the impact of different graphene-functionalization rates and sites on the stability of GDDS in “Discussion” section. Finally, drug molecules from DrugBank database (https://go.drugbank.com/) (Wishart et al., 2017) are used to computationally screen potential GDDS candidates for a pristine graphene model, and their results are fully summarized and analyzed in “Discussion” section.
Methods
In e-Graphene program, six main functions have been implemented and listed as follows: (1) generate graphene/GO models with different sizes and shapes, which can be intuitively visualized in three-dimensional (3D) viewer of e-Graphene; (2) inspect whether graphene/GO models satisfy the Hückel's rule (4n + 2) and full connectivity of Pi-bonds; (3) remove the duplicated GO models; (4) minimize graphene/GO models or ligands by the use of the natively implemented Tripos FF or external QM/SQM programs; (5) develop GDDS prediction/screen functions based on QGA and FF/SQM/QM; (6) implement an automatic and efficient cascade protocol for the GDDS prediction/screen; and (7) visualize GDDS prediction/screen results (Clark et al., 1989). For the sake of conciseness, each function is succinctly introduced one by one.
Generation and Inspection of Circular/Rectangular Graphene/Graphene Oxide Models
For different computational methods, users have to generate distinct graphene models with different sizes and shapes. e.g., in the QM studies, users often adopt a small circular graphene model satisfying the Hückel's rule (4n + 2). For the convenience of users, four options for the automatic generation of a pristine graphene model have been implemented in e-Graphene.
To improve the biocompatibility and solubility, pristine graphene is usually subjected to the chemical modification to become the GO material. However, the precise control of exact functionalization sites/rates on the pristine graphene is still quite challenging in the experiment. Therefore, an ensemble of GO models is generated to approximately stand for the fabricated GO material in the computational study. In our current implementation, common chemical modification groups (hydroxyl, epoxy, and carboxyl groups) can be automatically added to a pristine graphene model. Furthermore, the number of decoration groups (hydroxyl, epoxy, and carboxyl groups), functionalization direction (upper/lower/dual-surface), modification manner (random/systematic manner), and maximum number of GO models (Figure 2) can be specified by users. It is worth mentioning that for each GO model, e-Graphene can automatically inspect both the Hückel's rule (4n + 2) and full connectivity of Pi-bonds for GO models to ensure that the produced GO models still preserve the aromatic property.
Figure 2. Automatic generation of an ensemble of graphene oxide (GO) models for a given pristine graphene model.
Elimination of Duplicated Graphene Oxide Models
An ensemble of GO models can be automatically enumerated in a random or systematic manner by e-Graphene; however, some produced GO models are actually the same structure after certain translations and rotations. Thus, it is necessary to remove those duplicated GO models from the ensemble before subsequent calculations. In this work, we firstly calculate a large-diameter circular fingerprint “2048bit-ECFP10,” which is one type of commonly used extended connectivity fingerprint (ECFP) (Rogers and Hahn, 2010) to represent each GO model and then eliminate the duplication if some GO models share the same ECFP fingerprint, which has been implemented on the basis of our previous work about machine learning-guided bitterant prediction (Zheng et al., 2018).
Minimization of Graphene/Graphene Oxide Models
Before the GDDS prediction/screen, all the generated graphene/GO models are usually subjected to the minimization. For the rapid optimization, we have implemented Tripos FF in e-Graphene to automatically minimize the models by Powell minimization method that is natively coded in e-Graphene. At the same time, we also have developed an interface to call the external QM/SQM programs, such as Gaussian03/09/16 (Frisch et al., 2016), ORCA4.2 (Neese, 2012), MOPAC2009/2012/2016 (Stewart, 2016), and XTB6.3 (Bannwarth et al., 2019) to further minimize the graphene/GO model in a higher level of accuracy. In order to minimize all the graphene/GO models in batch mode, automatic calculation of the total formal charge for each model is developed in e-Graphene, whereas the multiplicity for each model is set to 1 by default due to our current focus on the ground state of the graphene/GO model. It is of note that this module in e-Graphene can also be utilized to minimize general small molecules.
Implementation of QGA-Guided Prediction/Screen of GDDS
Quantum genetic algorithm was proposed by Narayanan and Moore for the first time in 1996 (Narayanan and Moore, 1996) and introduced the concepts from quantum computation to the genetic algorithm. However, QGA was never used to predict the GDDS before. Thus, we combine QGA and QM/SQM/FF for this purpose. The whole implementation workflow of QGA is displayed in Figure 3, and the key steps are detailed as follows.
Figure 3. Implementation workflow of the quantum genetic algorithm (QGA) guided prediction of graphene-based drug delivery system (GDDS).
First, all the freely rotatable bonds (except the single bonds in the terminal groups, such as methyl group) for a small molecule are automatically detected and recorded, and all the torsion angles of the freely rotatable bonds, besides the overall translation (three variables) and rotation (quaternions with four variables) of the small molecule, are firstly encoded in a chromosome with the form of real number. Then this real-coded chromosome, reflecting one specific conformation of the small molecule, is further converted to a binary-coded chromosome for the subsequent QGA.
In QGA, each bit in each chromosome is referred to as a qubit, which is the basic unit of information in a quantum computer, and the corresponding state of each qubit is formulated by Dirac notation |ϕ > = α|0> +β|1> with the constraint of |α|2+|β|2=1, where α and β refer to the probability amplitudes of a qubit (Figure 4). For the simple description, each qubit is denoted by (α, β) or (cosθ, sinθ). Obviously, each qubit can be in the ground state (|0> or |1>), or the quantum superposition state of both ground states. Once a qubit is measured, it will be collapsed to the ground state |0> (or |1>) with a probability of |α|2 (or |β|2). If |α|2 (or |cosθ|2) is larger than |β|2 (or |sinθ|2), the qubit will be degraded to the ground state |0>, otherwise, it will be collapsed to the ground state |1>. Therefore, the concept “qubit” introduced in QGA provides the vast searching space.
Figure 4. Schematic illustration of a qubit. For the convenience of description, a qubit |ϕ > is simply denoted by (α, β) or (cosθ, sinθ).
Before the evolution of each chromosome, each qubit (cosθ, sinθ) in each chromosome is initialized to [cos(π/4), sin(π/4)], indicating that all the qubits initially have the same probability of being collapsed to either ground state |0> or |1>. If a randomly generated floating number between 0.0 and 1.0 is less than | cos(π/4) |2, the qubit will be degraded to the bit “1” (or |1>), otherwise, it will be collapsed to the bit “0” (or |0>). After the measurement of all the qubits in a chromosome, each chromosome is transformed to a conventional binary string containing only the bit “1/0,” which can be further decoded for the fitness calculation to update each chromosome for the next round of evolution.
Quantum rotation gate operator (QRGO) is utilized to update each chromosome with the rotation angle (Δθ) by comparing with the best chromosome in the whole population (Li et al., 2018), which has been defined by Equation (1). All the qubits in each chromosome are updated according to Equation (2), and the detailed update scheme of quantum rotation angle (Δθ) is tabulated in Supplementary Table 1 (Li et al., 2018). Similarly, quantum NOT gate operator (QNGO) is exerted to transform the qubits (Equations 3, 4) for the mutation of a chromosome. In addition, crossover operator is also applied to two different chromosomes to exchange parts of qubits in one chromosome with the corresponding qubits in the other chromosome. Moreover, catastrophe operator is used to reinitialize all the qubits in all the chromosomes if the whole population is trapped in a local minimum without any change of average fitness score in certain successive steps. After applying different operators with certain probabilities, all the newly updated chromosomes are decoded for the fitness calculation and the best chromosome are updated accordingly.
where QRGO refers to quantum rotation gate operator and Δθ is the variation of quantum rotation angle.
where each qubit is updated according to the QRGO; α and β are the probability amplitudes of the original qubit; and α′ and β′ are the probability amplitudes of the updated qubit.
where QNGO refers to quantum NOT gate operator.
where each qubit is updated according to the QNGO; α and β refer to the probability amplitudes of the original qubit; and α′ and β′ are the probability amplitude of the updated qubit.
Fitness of each chromosome is calculated from the potential energy of GDDS (Equation 5) or binding energy between the graphene (or GO) and ligand by Equations (6, 7). The potential energy or binding energy is computed by using the natively implemented Tripos FF in e-Graphene for the fast prediction or evoking the external QM/SQM programs for more accurate prediction.
where EGDDS, EG, and Elig are the potential energies of a given GDDS conformation, graphene/GO, and ligand, respectively. EDelta refers to the binding energy between the graphene/GO and ligand in a given GDDS conformation. The potential energy can be calculated by the internally implemented Tripos FF or externally evoked SQM/QM programs, such as Gaussian03/09/16, ORCA4.2, MOPAC2009/2012/2016, and XTB6.3. It is worth noting that Equation (5) is used to calculate the fitness by default, although both options (Equations 5, 6) are supported in e-Graphene.
QGA run can be terminated if the maximum number of evolution generation is reached or the average fitness difference between certain successive steps are consistently lower than 0.001, which is the default convergence criterion and can be customized by users. Once completed, all the best chromosomes are decoded to return the final conformations with customized Tripos Mol2 format files containing the energy information, which can be visualized by another module in e-Graphene as shown below.
Implementation of an Efficient Cascade Protocol for the GDDS Prediction/Screen
To make e-Graphene more efficient, we further develop a cascade protocol for the GDDS prediction/screen. Even though the default two-layer cascade protocol (QGA/FF-QGA/SQM or QGA/FF-QGA/QM) is strongly recommended and extensively evaluated in this work, more computationally expensive three-layer (QGA/FF–QGA/SQM-QGA/QM) cascade protocol is also supported in e-Graphene program (Figure 5). Users can select either protocol, define the relevant parameters for each layer, and set the number of the best chromosomes that are automatically migrating from the previous layer to the next layer for the population initialization of QGA in the next layer. It should be noted that the coarse QGA parameters, such as smaller populations and shorter evolution generations, should be given for the layer with QGA/SQM or QGA/QM in the practical prediction/screen; moreover, solvation and dispersion effects should be considered in the layer with QGA/SQM or QGA/QM.
Figure 5. A two-/three-layer cascade protocol implemented in e-Graphene.
Visualization of the GDDS Prediction/Screen Result
After GDDS prediction/screen, we are usually very interested in the optimal binding mode/energy between the graphene/GO and ligand for a given GDDS. For this purpose, we have implemented a convenient function to intuitively and synchronously visualize the graphene (or GO)-ligand binding mode/energy, and box information in 3D viewer of e-Graphene (Figure 6) so that users can easily examine whether QGA-guided prediction/screen gives reasonable results before further experimental tests.
Figure 6. Visualization of graphene-based drug delivery system (GDDS) prediction/screen result.
Results and Discussion
In this work, we propose to combine QGA and QM/SQM/FF for the intelligent prediction/screen of GDDS and further put forward an efficient cascade protocol by seamlessly integrating consecutive runs of QGA/FF–QGA/SQM (or QGA/FF–QGA/SQM-QGA/QM). For the convenience of users, we have developed a relatively comprehensive program “e-Graphene” with a graphic user interface to prepare and perform automatic prediction/screen of GDDS with/without the cascade protocol. The whole package including program, manual, and tutorial files is freely available on the Dropbox and Baidu Cloud public shared-folders with the respective download links (https://www.dropbox.com/sh/xj8wnc08kmw8b7a/AAAyXFbObUsk28QwFVfLjXCka?dl=0) and (https://pan.baidu.com/s/15_N7FNPrfjjoaA8f2s5gwQ with the four-letter extraction code: urdf). It is worth mentioning that e-Graphene has been fully tested and directly run on the local machine with 64-bit version of Windows10 operating system (OS) and can also prepare GDDS prediction/screen jobs including the input files, scripts, and a binary package for a remote computer with 64-bit version of Redhat5/6/7, CentOS5/6/7, or Ubuntu16.04/18.04/20.04-LTS OS. To demonstrate the applicability of this program and protocol, two typical examples are provided as below.
Different Chemical Functionalization Rates and Sites Impacting on GDDS Stability
In the work of Liu et al., GO material was demonstrated to effectively carry an anticancer drug SN38 by non-covalent interactions (Liu et al., 2008). However, the GO material is generally assumed to be a mixture from the structural point of view due to the fact that the oxidation process of pristine graphene is quite challenging to be precisely manipulated in terms of the exact functionalization sites and rates; hence, distinct functionalization sites and rates may have certain influence on the stability of GDDS, which was never probed in the original work of Liu et al. (2008). For the clarity, the stability mentioned in this work specifically refers to the binding energy between the graphene/GO and ligand for a given GDDS.
To tentatively explore this problem from the computational perspective, appropriate models should be produced to stand for the pristine and fabricated graphene materials. More concretely, a circular graphene model with the diameter of 21.760 Å and the total atom number of 168 is generated via clicking on the menu (Figure 1), and then this model is randomly functionalized by N hydroxyl group(s), N epoxy group(s), and N carboxyl group(s) (N = 1, 2, 3, 4, and 5) to produce five sets of 20 GO models, which represent five different rates of chemical modification on the pristine graphene model and are simply denoted by GO-1, GO-2, GO-3, GO-4, and GO-5, respectively, for the convenience of following discussion. Subsequently, one pristine graphene model and 100 GO models are combined to form a graphene/GO model database, which is systematically minimized in batch by evoking an external and efficient SQM program XTB6.3 with the keyword “–opt –gfn 2 –gbsa water reference” including the solvation and dispersion effect. Similarly, the 3D structure of SN38 with the neutral form is retrieved from the PubChem database (https://pubchem.ncbi.nlm.nih.gov) and also optimized by XTB6.3 with the same keyword.
In order to investigate the effect of different graphene-functionalization rates and sites on the stability of GDDS including SN38, the two-layer cascade protocol (QGA/FF–QGA/SQM) will be applied to all 101 graphene/GO-SN38 pairs due to our limited computational resource. In this protocol, the full parameter settings for current tests are shown in Supplementary Figures 1–3, while the key parameters will be highlighted as follow: the evolution generation number and population size in the fast QGA/FF calculation are set to 5,000 and 100, respectively, whereas the counterparts in the slow QGA/SQM calculation are set to 500 and 20, respectively. Ten best chromosomes from the cheap QGA/FF calculation are automatically passed to QGA/SQM. SQM is set to XTB6.3 with the keyword “–gfn 2 –gbsa water reference” including the solvation and dispersion effect as well. The respective probabilities of QRGO, QNGO, and crossover operation are set to 0.5, 0.1, and 0.7, respectively. Catastrophe operator is applied if the average fitness difference between two consecutive generations is always lower than 0.001 within 50 successive steps during the QGA evolution. The fitness is calculated by Equation (5) for the current test. According to these parameters, this cascade protocol was repeated three times for each GDDS because of the intrinsic randomness in QGA. Therefore, 303 runs with this cascade protocol were conducted, and all the corresponding results for the best chromosomes are listed in Supplementary Table 2.
Before evaluating the impact of different functionalization rates and sites on the stability of GDDS, the convergence analysis was first conducted to inspect whether QGA/FF in our cascade protocol could assist the convergence of QGA/SQM. For the convenience of discussion, one QGA run for the GDDS containing the pristine graphene model and SN38 was taken as an example. The average fitness of the whole chromosomes along the evolution of QGA/SQM is monitored because our final optimal solutions are actually derived from QGA/SQM rather than QGA/FF. In Figure 7, it is clearly shown that catastrophe operator is automatically applied in the 68th and 325th steps, respectively, where the average fitness declines dramatically due to the fact that all the chromosomes are completely reinitialized in order to escape from the previous local minimum. However, the average fitness in the 1st step is much larger than that in the 68th or 325th step because the best chromosomes from QGA/FF are used to seed the initial population for QGA/SQM. Apparently, it is much easier to reach the plateau when QGA/SQM run starts from the 1st step, whereas it takes longer evolution to achieve the convergence if QGA/SQM run is completely randomly reinitialized from the 68th or 325th step. From this perspective, QGA/FF indeed can accelerate the convergence of QGA/SQM, and in the practical predictions, users can turn off the option of removing the catastrophe operator and turn on the option of activating the automatic inspection of convergence in QGA/SQM for the faster calculation, while still preserving the catastrophe operator in QGA/FF to produce diverse optimal chromosomes for the subsequent QGA/SQM in the cascade protocol.
Figure 7. Average fitness for the whole population along the quantum genetic algorithm (QGA)/semi-empirical quantum mechanics (SQM) evolution: (A) original curve and (B) local enlarged curve.
To scrutinize the influence of different functionalization rates on the GDDS stability, the binding energy between the graphene/GO model and SN38 is derived from the best chromosome after running the two-layer (QGA/FF–QGA/SQM) cascade protocol. All the results are listed in Supplementary Table 2 and Table 1. From Table 1, it is shown that the 95% confidence intervals of binding energy for graphene/SN38, GO-1/SN38, GO-2/SN38, GO-3/SN38, GO-4/SN38, and GO-5/SN38 are −25.8118 ± 1.2778, −25.7589 ± 0.5120, −24.1041 ± 1.8056, −23.3682 ± 0.8053, −21.9780 ± 1.2295, and −20.3435 ± 1.0982 (kcal/mol), respectively, which suggest that GO-4/SN38 and GO-5/SN38 are obviously less stable than GO-1/SN38 and GO-2/SN38. Therefore, high functionalization rates compromise the stability of GDDS, indicating that an appropriate rate of fabrication on the pristine graphene is very important to achieve a stable GDDS with good solubility and biocompatibility.
Table 1. The 95% confidence intervals of binding energies between the graphene/GO model and SN38 derived from three runs with the two-layer cascade protocol.
To further evaluate the impact of different functionalization sites on the stability of GDDS with the same fabrication rate, the histogram of binding energies between the GO models and SN38 is as plotted in Figure 8. It is manifested that the binding energies of GOs/SN38 with the same functionalization rate, but different fabrication sites, have large variations. Here, we take the GO-1 models as examples, the lowest and highest binding energy between GO-1 and SN38 are −27.7450 and −22.3940 kcal/mol, respectively. Thus, different functionalization sites have a remarkable impact on the GDDS stability, but unfortunately the exact fine-tuning of functionalization sites on the pristine graphene still remains challenging from the experimental perspective, which indicates that it is more reasonable to harness the GO models with different functionalization sites to stand for the GO material with a certain average functionalization rate.
Figure 8. The histogram of average binding energies (kcal/mol) between GO models and SN38. For each GO-X (X = 1, 2, 3, 4, and 5), all the GO-X models have the same functionalization rate with different functionalization sites.
GDDS Screen for a Given Graphene/GO Model and a Given Chemical Database
Our cascade protocol can also be used to computationally screen the potential GDDS candidates for a given graphene (or GO) model and a given chemical database. For the demonstration purpose, we adopt the same circular pristine graphene model (diameter: 21.760 Å) as mentioned above and choose the DrugBank database (Wishart et al., 2017) for this screen. More specifically, 3D structures of 8,820 drug molecules with the neutral form are retrieved from the web server of DrugBank and then are subjected to the filtering of ChemAxon LogP larger than 5.0, which affords 200 highly hydrophobic drug molecules. Furthermore, if the longest distance within the 3D structure of a drug molecule is larger than the diameter of our current circular graphene model, this drug molecule will be excluded in our screen because our current graphene model in this demonstration is not large enough and its edge effect will be severe for this molecules. According to this criterion, 121 drug molecules are obtained and further minimized by XTB6.3 program with the keyword “–opt –gfn 2 –gbsa water reference,” which are finally utilized for the current screen by the two-layer cascade protocol (QGA/FF–QGA/SQM). SQM is also set to XTB6.3 with the keyword “–gfn 2 –gbsa water reference” that takes into account the solvation and dispersion effect. All the parameters for this protocol used in the screen are exactly the same as those in the previous section. This screen campaign was also performed for three times due to the inherent stochastic nature of QGA.
Once three runs of GDDS screen in e-Graphene were completed, all the conformations for the best chromosomes were harvested and prioritized according to the average binding energy over three repeats (Supplementary Table 3 and Figure 9). Only 12 GDDS candidates (or 12 graphene/drugs) had lower binding energies (−36.3613 ~ −27.6150 kcal/mol) than the 95% CI of binding energy for the graphene model/SN38 (−27.0896, −24.5340) kcal/mol as shown in our previous test. Further manual scrutinization of 12 compounds in the DrugBank database suggested that four of them were anticancer drugs (midostaurin, nilotinib, tucatinib, and arzoxifene) that usually need to be delivered to the specific cancer cells, and their average binding energies with the graphene model were −31.8856, −30.8095, −29.6934, and −27.8303 kcal/mol, respectively (Table 2), whereas most of the remaining drugs for the external usage were not necessary for the delivery by the graphene/GO. Thus, only these four anticancer drugs are proposed to be used for the potential GDDS candidates from the computational perspective.
Figure 9. The histogram of average binding energies (kcal/mol) between the graphene model and 121 drugs after screen with the two-layer (QGA/FF–QGA/SQM) cascade protocol.
Table 2. Four predicted anticancer drugs for the potential GDDS candidates.
Limitation and Outlook of Our Computational Platform and Protocol
Although e-Graphene could provide a pragmatic computational platform and protocol for the automatic prediction/screen of GDDS based on QGA and QM/SQM/FF, it indeed had several drawbacks: (1) QGA belongs to the stochastic method, which produces different results for different runs. Thus, multiple runs are strongly recommended to achieve more convergent results, which undoubtedly brings some extra computational burden. (2) The speed and accuracy of two-layer (QGA/FF–QGA/SQM) or three-layer (QGA/FF–QGA/SQM–QGA/QM) cascade protocol was largely dependent on the computationally expensive SQM or QM with different keywords and parameters. Careful benchmarks of different choices are highly advised to obtain the comparatively accurate yet fast combination for the practical prediction/screen of GDDS. Thus, it is suggested for users to clearly understand the weakness of e-Graphene before applying it to their own projects.
As recent booming of quantum computation in various areas, QGA possesses intrinsic quantum operators, thus we tentatively envision that if QGA can be ported to the quantum computer simulator and subsequent real quantum computer, QGA might release its tremendous native capacity. However, it still has a very long way to go before QGA-based protocol can be fully deployed in the quantum computer due to the following two main obstacles: (1) the current quantum computers support a small number of qubits, e.g., 53 bits from Google quantum computer (Arute et al., 2019); (2) the accessibility of quantum computers is quite limited, which undoubtedly imposes a steep learning curve, even though quantum computer simulators could alleviate this circumstance to some extent. Despite of those difficulties, we believe that this emerging area still draws more and more researchers to foster its sustainable development.
Conclusion
In this work, we present a practical computational platform “e-Graphene” for the automatic prediction/screen of GDDS by the seamless integration between QGA and QM/SQM/FF and further propose a cascade protocol for the efficient prediction/screen of GDDS. In this platform, a graphene/GO model can be automatically generated according to the shape and size specified by users, subsequently minimized by the internally implemented Tripos FF or externally evoked QM/SQM programs, and finally used for the preparation of GDDS prediction/screen job with/without the cascade protocol, which can be ran on the local machine with Windows10 and on the remote machine with Red Hat/CentOS/Ubuntu. By using this platform, two representative test cases are adopted to demonstrate its pragmatic functions. One is systematic tests on a extensively used GDDS containing an anticancer drug SN38, which manifests that high rates of functionalization on the pristine graphene reduces the stability of GDDS and also indicates that an appropriate fabrication rate is crucial to attain an optimal GDDS with good stability, solubility, and biocompatibility. The other is the GDDS screen in the DrugBank database for a pristine graphene model, which harvests four potential GDDS candidates with better stability than the common GDDS containing SN38. At last, the limitations of this computational platform are also addressed, although it is envisioned that this free platform could enable experimental scientists to automatically predict possible GDDS candidates before experimental tests.
Data Availability Statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.
Author Contributions
SZ and JX performed the calculations and analysis. LW and DZ tested the program. YX and FL revised the manuscript. SZ and FL conceived the workflow, wrote the manuscript, and developed the program. All authors contributed to the article and approved the submitted version.
Funding
The authors acknowledge the support of Zhejiang Provincial Natural Science Foundation of China (Grant no. LY17B030007), the basic research funding (Grant no. KYYW201928), and startup funding from Wenzhou Medical University.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2021.664355/full#supplementary-material
References
Aditya, D. G., Abhay, R. S., and Vilasrao, J. K. (2017). Graphene: a comprehensive review. Curr. Drug Targets 18, 724–733. doi: 10.2174/1389450117666160709023425
Ajala, A. O., Voora, V., Mardirossian, N., Furche, F., and Paesani, F. (2019). Assessment of density functional theory in predicting interaction energies between water and polycyclic aromatic hydrocarbons: from water on benzene to water on graphene. J. Chem. Theory Comput. 15, 2359–2374. doi: 10.1021/acs.jctc.9b00110
Alkathiri, F. A., McCallion, C., Golovanov, A. P., Burthem, J., Pluen, A., and Bryce, R. A. (2019). Solvation of pristine graphene using amino acids: a molecular simulation and experimental analysis. J. Phys. Chem. C 123, 30234–30244. doi: 10.1021/acs.jpcc.9b06648
Arute, F., Arya, K., Babbush, R., Bacon, D., Bardin, J. C., Barends, R., et al. (2019). Quantum supremacy using a programmable superconducting processor. Nature 574, 505–510. doi: 10.1038/s41586-019-1666-5
Azhagiya Singam, E. R., Zhang, Y., Magnin, G., Miranda-Carvajal, I., Coates, L., Thakkar, R., et al. (2019). Thermodynamics of adsorption on graphenic surfaces from aqueous solution. J. Chem. Theory Comput. 15, 1302–1316. doi: 10.1021/acs.jctc.8b00830
Bannwarth, C., Ehlert, S., and Grimme, S. (2019). GFN2-xTB—an accurate and broadly parametrized self-consistent tight-binding quantum chemical method with multipole electrostatics and density-dependent dispersion contributions. J. Chem. Theory Comput. 15, 1652–1671. doi: 10.1021/acs.jctc.8b01176
Bernardi, R. C., Melo, M. C. R., and Schulten, K. (2015). Enhanced sampling techniques in molecular dynamics simulations of biological systems. Biochim. Biophys. Acta Gen. 1850, 872–877. doi: 10.1016/j.bbagen.2014.10.019
Calero, C., Martí, J., Guàrdia, E., and Masia, M. (2013). Characterization of the methane–graphene hydrophobic interaction in aqueous solution from Ab initio simulations. J. Chem. Theory Comput. 9, 5070–5075. doi: 10.1021/ct400554q
Cho, Y., Min, S. K., Yun, J., Kim, W. Y., Tkatchenko, A., and Kim, K. S. (2013). Noncovalent interactions of DNA bases with naphthalene and graphene. J. Chem. Theory Comput. 9, 2090–2096. doi: 10.1021/ct301097u
Chung, C., Kim, Y.-K., Shin, D., Ryoo, S.-R., Hong, B. H., and Min, D.-H. (2013). Biomedical applications of graphene and graphene oxide. Acc. Chem. Res. 46, 2211–2224. doi: 10.1021/ar300159f
Clark, M., Cramer III, R. D., and Van Opdenbosch, N. (1989). Validation of the general purpose tripos 5.2 force field. J. Comput. Chem. 10, 982–1012. doi: 10.1002/jcc.540100804
Conte, R., and Ceotto, M. (2020). “Semiclassical molecular dynamics for spectroscopic calculations,” in Quantum Chemistry and Dynamics of Excited States, 595–628. doi: 10.1002/9781119417774.ch19
El-Kady, M. F., Shao, Y., and Kaner, R. B. (2016). Graphene for batteries, supercapacitors and beyond. Nat. Rev. Mater. 1:16033. doi: 10.1038/natrevmats.2016.33
Frisch, M. J., Trucks, G. W., Schlegel, H. B., Scuseria, G. E., Robb, M. A., Cheeseman, J. R., et al. (2016). Gaussian 16 Rev. C.01. Wallingford, CT. Available online at: https://gaussian.com/citation/
Gandolfi, M., Rognoni, A., Aieta, C., Conte, R., and Ceotto, M. (2020). Machine learning for vibrational spectroscopy via divide-and-conquer semiclassical initial value representation molecular dynamics with application to N-methylacetamide. J. Chem. Phys. 153:204104. doi: 10.1063/5.0031892
Georgakilas, V., Tiwari, J. N., Kemp, K. C., Perman, J. A., Bourlinos, A. B., Kim, K. S., et al. (2016). Noncovalent functionalization of graphene and graphene oxide for energy materials, biosensing, catalytic, and biomedical applications. Chem. Rev. 116, 5464–5519. doi: 10.1021/acs.chemrev.5b00620
Ghosal, K., and Sarkar, K. (2018). Biomedical applications of graphene nanomaterials and beyond. ACS Biomater. Sci. Eng. 4, 2653–2703. doi: 10.1021/acsbiomaterials.8b00376
Goenka, S., Sant, V., and Sant, S. (2014). Graphene-based nanomaterials for drug delivery and tissue engineering. J. Control. Release 173, 75–88. doi: 10.1016/j.jconrel.2013.10.017
Gráfová, L., Pitonák, M., Rezáč, J., and Hobza, P. (2010). Comparative study of selected wave function and density functional methods for noncovalent interaction energy calculations using the extended S22 data set. J. Chem. Theory Comput. 6, 2365–2376. doi: 10.1021/ct1002253
Guo, Y.-n., Lu, X., Weng, J., and Leng, Y. (2013). Density functional theory study of the interaction of arginine-glycine-aspartic acid with graphene, defective graphene, and graphene oxide. J. Phys. Chem. C 117, 5708–5717. doi: 10.1021/jp310088e
Krepel, D., and Hod, O. (2017). Selectivity of a graphene nanoribbon-based trinitrotoluene detector: a computational assessment. J. Phys. Chem. C 121, 21546–21552. doi: 10.1021/acs.jpcc.7b07540
Laio, A., and Parrinello, M. (2002). Escaping free-energy minima. Proc. Natl. Acad. Sci. U.S.A. 99, 12562–12566. doi: 10.1073/pnas.202427399
Li, Y., Gong, G., and Li, N. (2018). A parallel adaptive quantum genetic algorithm for the controllability of arbitrary networks. PLoS ONE 13:e0193827. doi: 10.1371/journal.pone.0193827
Lightcap, I. V., and Kamat, P. V. (2013). Graphitic design: prospects of graphene-based nanocomposites for solar energy conversion, storage, and sensing. Acc. Chem. Res. 46, 2235–2243. doi: 10.1021/ar300248f
Liu, J., Cui, L., and Losic, D. (2013). Graphene and graphene oxide as new nanocarriers for drug delivery applications. Acta Biomater. 9, 9243–9257. doi: 10.1016/j.actbio.2013.08.016
Liu, Z., Robinson, J. T., Sun, X., and Dai, H. (2008). PEGylated nanographene oxide for delivery of water-insoluble cancer drugs. J. Am. Chem. Soc. 130, 10876–10877. doi: 10.1021/ja803688x
Mahdavi, M., Fattahi, A., Tajkhorshid, E., and Nouranian, S. (2020). Molecular insights into the loading and dynamics of doxorubicin on PEGylated graphene oxide nanocarriers. ACS Appl. Bio Mater. 3, 1354–1363. doi: 10.1021/acsabm.9b00956
Mahdavi, M., Rahmani, F., and Nouranian, S. (2016). Molecular simulation of pH-dependent diffusion, loading, and release of doxorubicin in graphene and graphene oxide drug delivery systems. J. Mater. Chem. B 4, 7441–7451. doi: 10.1039/C6TB00746E
Mao, H. Y., Laurent, S., Chen, W., Akhavan, O., Imani, M., Ashkarran, A. A., et al. (2013). Graphene: promises, facts, opportunities, and challenges in nanomedicine. Chem. Rev. 113, 3407–3424. doi: 10.1021/cr300335p
Mason, S. E., Beton, P. H., and Besley, N. A. (2019). AIRBED: a simplified density functional theory model for physisorption on surfaces. J. Chem. Theory Comput. 15, 5628–5634. doi: 10.1021/acs.jctc.9b00576
McCallion, C., Burthem, J., Rees-Unwin, K., Golovanov, A., and Pluen, A. (2016). Graphene in therapeutics delivery: problems, solutions and future opportunities. Eur. J. Pharm. Biopharm. 104, 235–250. doi: 10.1016/j.ejpb.2016.04.015
Mohammad Omaish, A., Kalamegam, G., Abdurahman, E., Sidi, A. B., and Adnan, M. (2019). Graphene and graphene-based materials in biomedical applications. Curr. Med. Chem. 26, 6834–6850. doi: 10.2174/0929867326666190705155854
Morris, G. M., Huey, R., Lindstrom, W., Sanner, M. F., Belew, R. K., Goodsell, D. S., et al. (2009). AutoDock4 and AutoDockTools4: automated docking with selective receptor flexibility. J. Comput. Chem. 30, 2785–2791. doi: 10.1002/jcc.21256
Mudedla, S. K., Balamurugan, K., and Subramanian, V. (2014). Computational study on the interaction of modified nucleobases with graphene and doped graphenes. J. Phys. Chem. C 118, 16165–16174. doi: 10.1021/jp503126q
Narayanan, A., and Moore, M. (1996). “Quantum-inspired genetic algorithms,” in Proceedings of IEEE International Conference on Evolutionary Computation (Nagoya). doi: 10.1109/ICEC.1996.542334
Novoselov, K. S., Fal′ko, V. I., Colombo, L., Gellert, P. R., Schwab, M. G., and Kim, K. (2012). A roadmap for graphene. Nature 490, 192–200. doi: 10.1038/nature11458
Novoselov, K. S., Geim, A. K., Morozov, S. V., Jiang, D., Zhang, Y., Dubonos, S. V., et al. (2004). Electric field effect in atomically thin carbon films. Science 306, 666–669. doi: 10.1126/science.1102896
Pan, Y., Sahoo, N. G., and Li, L. (2012). The application of graphene oxide in drug delivery. Expert Opin. Drug Deliv. 9, 1365–1376. doi: 10.1517/17425247.2012.729575
Ramraj, A., and Hillier, I. H. (2010). Binding of pollutant aromatics on carbon nanotubes and graphite. J. Chem. Inf. Model. 50, 585–588. doi: 10.1021/ci1000604
Reina, G., González-Domínguez, J. M., Criado, A., Vázquez, E., Bianco, A., and Prato, M. (2017). Promises, facts and challenges for graphene in biomedical applications. Chem. Soc. Rev. 46, 4400–4416. doi: 10.1039/C7CS00363C
Rogers, D., and Hahn, M. (2010). Extended-connectivity fingerprints. J. Chem. Inf. Model. 50, 742–754. doi: 10.1021/ci100050t
Safdari, F., Raissi, H., Shahabi, M., and Zaboli, M. (2017). DFT calculations and molecular dynamics simulation study on the adsorption of 5-fluorouracil anticancer drug on graphene oxide nanosheet as a drug delivery vehicle. J. Inorg. Organomet. Polym. 27, 805–817. doi: 10.1007/s10904-017-0525-9
Shim, G., Kim, M.-G., Park, J. Y., and Oh, Y.-K. (2016). Graphene-based nanosheets for delivery of chemotherapeutics and biological drugs. Adv. Drug Deliv. Rev. 105, 205–227. doi: 10.1016/j.addr.2016.04.004
Sun, X., Liu, Z., Welsher, K., Robinson, J. T., Goodwin, A., Zaric, S., et al. (2008). Nano-graphene oxide for cellular imaging and drug delivery. Nano Res. 1, 203–212. doi: 10.1007/s12274-008-8021-8
Vincent, M. A., and Hillier, I. H. (2014). Accurate prediction of adsorption energies on graphene, using a dispersion-corrected semiempirical method including solvation. J. Chem. Inf. Model. 54, 2255–2260. doi: 10.1021/ci5003729
Vovusha, H., Banerjee, D., Yadav, M. K., Perrozzi, F., Ottaviano, L., Sanyal, S., et al. (2018). Binding characteristics of anticancer drug doxorubicin with two-dimensional graphene and graphene oxide: insights from density functional theory calculations and fluorescence spectroscopy. J. Phys. Chem. C 122, 21031–21038. doi: 10.1021/acs.jpcc.8b04496
Vovusha, H., Sanyal, S., and Sanyal, B. (2013). Interaction of nucleobases and aromatic amino acids with graphene oxide and graphene flakes. J. Phys. Chem. Lett. 4, 3710–3718. doi: 10.1021/jz401929h
Wang, X., Liu, Y., Xu, J., Li, S., Zhang, F., Ye, Q., et al. (2015). Molecular dynamics study of stability and diffusion of graphene-based drug delivery systems. J. Nanomater. 2015:872079. doi: 10.1155/2015/872079
Wassei, J. K., and Kaner, R. B. (2013). Oh, the places you'll go with graphene. Acc. Chem. Res. 46, 2244–2253. doi: 10.1021/ar300184v
Wishart, D. S., Feunang, Y. D., Guo, A. C., Lo, E. J., Marcu, A., Grant, J. R., et al. (2017). DrugBank 5.0: a major update to the DrugBank database for 2018. Nucleic Acids Res. 46, D1074–D1082. doi: 10.1093/nar/gkx1037
Xianfeng, Z., and Feng, L. (2014). Application of graphene/graphene oxide in biomedicine and biotechnology. Curr. Med. Chem. 21, 855–869. doi: 10.2174/0929867320666131119124325
Yang, K., Feng, L., and Liu, Z. (2016). Stimuli responsive drug delivery systems based on nano-graphene for cancer therapy. Adv. Drug Deliv. Rev. 105, 228–241. doi: 10.1016/j.addr.2016.05.015
Yang, X., Zhang, X., Liu, Z., Ma, Y., Huang, Y., and Chen, Y. (2008). High-efficiency loading and controlled release of doxorubicin hydrochloride on graphene oxide. J. Phys. Chem. C 112, 17554–17558. doi: 10.1021/jp806751k
Yang, Y., Asiri, A. M., Tang, Z., Du, D., and Lin, Y. (2013). Graphene based materials for biomedical applications. Mater. Today 16, 365–373. doi: 10.1016/j.mattod.2013.09.004
Yi, L., Zhang, Y., Shi, X., Du, X., Wang, X., Yu, A., et al. (2019). Recent progress of functionalised graphene oxide in cancer therapy. Drug J. Target. 27, 125–144. doi: 10.1080/1061186X.2018.1474359
Yu, X., Cheng, H., Zhang, M., Zhao, Y., Qu, L., and Shi, G. (2017). Graphene-based smart materials. Nat. Rev. Mater. 2:17046. doi: 10.1038/natrevmats.2017.46
Zhang, Q., Wu, Z., Li, N., Pu, Y., Wang, B., Zhang, T., et al. (2017). Advanced review of graphene-based nanomaterials in drug delivery systems: synthesis, modification, toxicity and application. Mater. Sci. Eng. C 77, 1363–1375. doi: 10.1016/j.msec.2017.03.196
Zheng, S., Jiang, M., Zhao, C., Zhu, R., Hu, Z., Xu, Y., et al. (2018). e-bitter: bitterant prediction by the consensus voting from the machine-learning methods. Front. Chem. 6:82. doi: 10.3389/fchem.2018.00082
Keywords: graphene, graphene-based drug delivery system, GDDS, prediction, quantum genetic algorithm, cascade protocol, QGA
Citation: Zheng S, Xiong J, Wang L, Zhai D, Xu Y and Lin F (2021) e-Graphene: A Computational Platform for the Prediction of Graphene-Based Drug Delivery System by Quantum Genetic Algorithm and Cascade Protocol. Front. Chem. 9:664355. doi: 10.3389/fchem.2021.664355
Received: 05 February 2021; Accepted: 06 April 2021;
Published: 07 May 2021.
Edited by:
Nino Russo, University of Calabria, ItalyReviewed by:
Stephan Irle, Oak Ridge National Laboratory (DOE), United StatesRiccardo Conte, University of Milan, Italy
Copyright © 2021 Zheng, Xiong, Wang, Zhai, Xu and Lin. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Suqing Zheng, zsq_2016@126.com; Fu Lin, lin1449@126.com