Breast cancer is the most diagnosed cancer in women and the second leading cause of cancer-related mortality in women . Heritage is the most critical risk factor, and 15 to 20% of breast cancer is familiar . One of the characteristics of breast cancer is that it can be wholly cured given an early diagnosis . The mortality rate from breast cancer has been reduced by 1.9% annually from 2002 to 2011 and 1.3% from 2011 to 2020 . Diagnostics and treatments have continuously improved through the years. However, the situation is different in each country considering the costs and technological advances in each country. In the United States, 300,590 cases of breast cancer had been estimated for the year 2023, with a total of 43,700 deaths . Latin America has over 210,000 new cases and around 60,000 deaths yearly . For the year 2020, it was estimated that about 2.3 million breast cancer cases were diagnosed in women globally, and about 685,000 died from this disease . A recurrent problem with standard treatments are the side effects. Regarding the use of chemotherapeutic drugs, such issues are nephrotoxicity of cisplatin, cardiotoxicity of doxorubicin, and pulmonary fibrosis from the use of bleomycin . Besides, in the case of radiotherapy, fibrosis, atrophy, and neuronal damage caused by irradiation can occur . Consequently, novel treatments try to reduce the secondary effects while retaining the benefits of standard approaches.

Chemotherapy is one of the most extensively applied treatments for breast cancer, with different drug targets depending on the type of cancer. Progesterone- or estrogen-receptor-positive tumors are related to cancers with low mortality . Another common target in chemotherapy is human epidermal growth factor receptor 2 (HER2). Only 15 to 20% of all tumors are HER2-positive, overexpressing Erb-B2 receptor tyrosine kinase 2 (ERBB2) in the cell membrane. HER2 tumors are usually more aggressive than other ones, but the advantage is that their treatment is very effective . Another chemotherapy target is the chemokine C-X-C motif receptor 7 (CXCR7) . This G-protein is targeted because studies show a possible positive effect on inhibiting the metastasis of cervical cancer cells . However, more clinical and preclinical studies on CXCR7 and its co-player CXCR4 are required since alterations have been detected in diseases such as cancer, central nervous system and cardiac disorders, and autoimmune diseases .

In recent years, nanomaterials have attracted the attention of different scientific communities by providing them with new solutions for drug delivery . These nanotechnological applications have made it possible to obtain treatments that release substances at specific sites of interest, reducing the required drug amount and side effects. Nanostructures to form these drug delivery systems can be divided into organic and inorganic , with the latter one being the less extensively studied. One option currently considered in pharmacy and medicine is carbon-based nanomaterials because of their physicochemical, mechanical, electrical, thermal, and optical properties , as well as their capacity to modify existing drugs. Fullerene derivatives have been proposed recently, particularly those obtained from fullerene C60. The unmodified fullerene C60 is known as a “free radical sponge” because its double bonds tend to accept free radicals . Because of its size, surface area, and capacity to extinguish or generate reactive oxygen species, C60 is very promising in medicine and clinical therapy . It is also possible to modify pristine fullerenes by adding polar functional groups (e.g., –COOH, –OH, or –NH2), to improve water solubility, antioxidant properties, and even biological activity . For instance, polyhydroxy fullerenes (PHFs) exhibit properties suitable for biomedical applications, such as water solubility, biodegradability, biocompatibility, and hypoallergic response. It has been shown that PHFs can inhibit cancer tumor growth and positively regulate the immune system . The same is valid for carboxylated fullerenes ; for instance, C60[C(COOH)2]3 is well known for its high biological activity in plants and within mitochondrial dynamics .

Since the evaluation of novel drugs is a task that requires significant human and material resources, innovative strategies have been formulated as alternatives. Quantitative structure–activity and quantitative structure–property relationships (QSAR/QSPR) are a paradigm that can be useful in choosing promising molecules, considering the information on inactive and active compounds, through in silico approaches. According to the QSAR/QSPR paradigm, a given activity/property, f, can be modeled using a set of quantitative descriptors, x1, x2, x3,..., xn, theoretically determined or measured by experiments . A relationship f(x1, x2, x3,..., xn) can be defined to predict the activity or property of molecules after the evaluation of their quantitative descriptors. However, the QSAR/QSPR paradigm does not explain how to select the descriptors or how to build the mathematical function. Consequently, the following paragraphs discuss basic concepts about selecting descriptors and regression techniques implemented in this manuscript.

Lipinski’s rule of five is a compendium of guidelines commonly used to determine if a molecule can be proposed as an orally delivered drug according to its physicochemical properties. According to this rule, a drug compound should have a molecular weight below 500 g/mol, a octanol–water partition coefficient (LogP) below 5, less than five hydrogen bond donor sites, and less than ten hydrogen bond acceptors sites. It is possible to add two other conditions, namely polar surface area (PSA) ≤ 140 Å2 and less than ten rotatable bonds . Taking advantage of the readiness of these quantities in public datasets, the current study proposes some of these quantities as potentially suitable descriptors for predictive models. Besides, Pearson’s hard–soft acid–base (HSAB) theory suggests other descriptors to describe and predict the interactions between chemical species, such as those between a drug molecule as a ligand and a protein . These quantitative values are based on the vertical ionization energy (I) and electron affinity (A). According to Koopmans’ theorem, both can be approximated by I = −EHOMO and A = −ELUMO, where EHOMO is the energy of the highest occupied molecular orbital (HOMO), and ELUMO is the energy of the lowest unoccupied molecular orbital (LUMO). It is advantageous to combine these properties to find out if an interaction between two species will occur and to obtain new quantitative relationships. Another helpful descriptor is the global electrophilicity, calculated as ω = χ2/2η . Electrophilicity is related to the energetic stabilization that a species gains by obtaining an additional electron.

First, 42 drugs related to chemotherapy treatments for breast cancer were proposed. Although the most notable fullerene derivatives for biological applications are those with several hydrophilic groups, the carboxylic acid derivative C60–COOH has been studied as well. Baglayan and coworkers carried out a conformation analysis within DFT to obtain the ground state structure for C60–COOH . In addition, they discussed its usage as a potential drug carrier for the antimetabolic and anticancer drug 5-fluoruracil . Similarly, Parlak and Alver reported a theoretical study on the interactions and stability of paracetamol complexes with C60–COOH . Consequently, this work proposes the interaction of C60–COOH fullerene with anticancer drugs. As a complement, a water-soluble fullerene predicted as stable at the normal human body temperature was proposed to study the interactions with doxorubicin and gemcitabine . The water-soluble fullerene is introduced to avoid known mutagenic reactions related to breast cancer . It was also studied as a potential carrier for bedaquiline, an agent against tuberculosis . The current study only considered molecules and complexes formed with up to 100 atoms to be affordable with our computational resources.

For all calculations within DFTB3, the 3OB parameter set was used . To carry out the global optimization procedure, Balloon 1.8.2 and DFTB+ 17.1 were used for the initial conformational study by genetic algorithms and final optimization at the DFTB3 level, respectively. London dispersion forces were considered in the DFTB3 and global optimization procedures by Lennard-Jones potentials, as implemented in UFF and MMFF94 force fields, respectively. The solvent effect was included by the Born solvation model within DFTB3. The study considered the chemotherapy drugs isolated and interacting with pristine C60 fullerene as well as its carboxylic acid derivative C60–COOH. Eight initial drug–fullerene structures were proposed to obtain their global optimization by means of DFTB3. The drugs were initially set at 1.5 Å of minimal distance from the fullerene. Once the global optimization was done, the same steps as for the isolated drugs were carried out for the molecular docking. The datasets were modified to take into account the effect of the fullerenes. Also, the validation set was reduced because of the large size of the complexes.

The atypical chemokine receptor 3, also known as CXCR7 or G-protein-coupled receptor 159 (GPR159) , was selected as the target protein for molecular docking. The iterative assembly refinement server (I-Tasser) was used to produce an initial structure for the CXCR7 protein by the homology approach. The sequence was extracted from the UniProtKB/Swiss-Prot dataset. From all homology structures produced by the I-Tasser server, the one with the highest confidence coefficients was selected to produce a reliable initial structure . The lowest-energy structure, as in the study of Muthiah and coworkers , was validated using PROCHECK to check the quality of the protein structure. The PDB produced with the previous step was subsequently optimized by an energy minimization through Amber force fields using the USCF Chimera 1.14 toolkit . The secondary structural features were stabilized by TMpred and HMMTOP during energy minimization. Last, the protein was prepared by setting atomic charges and hydrogen atoms and merging the nonpolar groups. Once the structures were optimized, molecular docking was performed with the CXCR7 protein, using Autodock Vina 1.1, to obtain the docking score, established hydrogen bonds, and the binding site (pocket). The above was done for all drugs in the dataset and an external validation set.

IBM Watson AI was used to build the models and to predict the docking score through the Extra Trees regressor algorithm . It was also used to obtain the most significant quantum descriptors used in each model. Extra Trees, an abbreviation of “extremely randomized trees”, is a mathematical method used to estimate a relationship between data and the covariates . The Extra Trees algorithm creates many decision trees , but the sampling of each one is random. Thus, a dataset for each tree contains unique samples. The optimization of the hyperparameters associated with the decision trees obtained was performed by the derivative-free global search algorithm known as RBfOpt, which fits a radial basis function mode to accelerate the discovery of the hyperparameters . All the above was used through the AutoAI tool within IBM Watson, an automatized routine to select the model with the best performance among those available in the platform. Since this method does not produce exportable mathematical models, another approach was used as detailed below .

Multiple linear regression (MLR) could be a tool to solve the problem in a complementary way to Extra Trees regression. MLR is a mathematical model that can be seen as an extension of linear regression. In terms of n input variables, x1, x2,…, xn, the outcome y can be expanded by the following linear expansion :

Figure 1: Workflow diagram of the stages during modeling.

Table 1: Quantitative descriptors proposed to model the docking score of the isolated drugs, as well as of those modified with fullerenes C60 and C60–COOH, interacting with the protein CXCR7.

Table 2: Docking score, number of established H-bonds, and protein–ligand interacting residue in three-letter symbol, up to 3 Å distance. Drugs marked with an asterisk were used as external validation set.

The docking scores ranged from −10.1 to −4.6 kcal/mol for the training set. The molecule with the most significant bond strength, according to its docking score, was olaparib, whereas the one with the lowest bond strength was fluorouracil. The number of hydrogen bonds was computed ranging from 0 to 5. It is important to note that the number of hydrogen bonds is not directly related to the docking score since there are weak and strong hydrogen bonds. This assumption was proved by the analysis performed to obtain the predictive models, as discussed below.