Fullerene Derivatives for Drug Delivery against COVID-19: A Molecular Dynamics Investigation of Dendro[60]fullerene as Nanocarrier of Molnupiravir

In this paper, a theoretical investigation is made regarding the possibility of using a water-soluble derivative of C60 as a drug delivery agent for treating Coronavirus disease 2019 (COVID-19). Molnupiravir is chosen as the transporting pharmaceutical compound since it has already proved to be very helpful in saving lives in case of hospitalization. According to the proposed formulation, a carboxyfullerene known as dendro[60]fullerene is externally connected with two molnupiravir molecules. Two properly formed nitrogen single bonds (N−N) are used as linkers between the dendro[60]fullerene and the two molnupiravir molecules to create the final form of the C60 derivate/molnupiravir conjugate. The energetics of the developed molecular system and its interaction with water and n-octanol are extensively studied via classical molecular dynamics (MD) using the COMPASS II force field. To study the interactions with water and n-octanol, an appropriate periodic amorphous unit cell is created that contains a single C60 derivative/molnupiravir system surrounded by numerous solvent molecules and simulated via MD in room conditions. In addition, the corresponding solvation-free energies of the investigated drug delivery system are computed and set in contrast with the corresponding properties of the water-soluble dendro[60]fullerene, to test its solubility capabilities.


Introduction
The COVID-19 pandemic has driven the majority of recent research efforts in the rapid development of new vaccines [1], antibody drugs [2], and antiviral medicines [3] against the SARS-CoV-2 virus. Nanotechnology and its discoveries in the last few decades could not be absent in this challenge [4].
In the field of drugs, perhaps one of the most promising medicines is the so-called molnupiravir, which seems that it may contribute to the healing from COVID-19 and its symptoms when given orally in solid-dose forms such as tablets and capsules [5], which provide a good degree of stability and deliver precise dosage. The molnupiravir reduces the ability of SARS-CoV-2, the virus that causes COVID-19, to multiply in the body [6]. Very recently, Sharov et al. [7] provided valuable density functional theory (DFT) calculations that may contribute to the understanding and development of new therapies against SARS-CoV-2 based on molnupiravir. In their in-depth computational study, they performed molecular docking simulations of three plausible tautomeric forms of molnupiravir with a series of the SARS-CoV-2 proteins to predict the corresponding structural, electronic, and optical responses. Given the action mechanism of molnupiravir, it becomes evident that the more controllable and accurate the delivery of the specific drug within the human body, the more successful the provided treatment against COVID- 19. In this effort, the drug delivery approaches and technologies via the use of nanoparticles [8] may help to achieve the desired therapeutic results against COVID-19 since they may provide safe storage and transport of the required pharmaceutical compound to its target location with a high degree of precision.
It is now well understood that the most critically affected human organ by SARS-CoV-2 is the lung [9]. Therefore, there is already an open discussion about the optional strategies to prevent or fight the action of the virus in the deep lung and to target the most important host cells via more focused drug delivery strategies [10]. In any case, there are certain areas within the human body where the virus-cell populations are denser or spots where the potential damage and risk are much higher [11] and, thus, it seems that drug delivery methods that are based on nanostructured systems [12] may offer very promising treatment solutions against the recent Coronavirus disease. In an extensive review, Abd Elkodous et al. [13] gave insights into several nanomaterial-based methodologies for increasing the effectiveness of the drugs against COVID-19 by their attachment inside or onto the molecular surface of five promising drug delivery nanostructured agents, i.e., nano-emulsions, liposomes, polymeric nanoparticles and micelles, dendrimers, and nano-suspensions. Similarly, Witika et al. [14] provided a state-of-the-art regarding the possible use of biomimetic nanocarriers for the efficient treatment of COVID-19. They have indicated that there is a great potential for adopting a variety of nanomaterials such as nanospheres, nanocrystals, solid lipid nanoparticles, nanosponges, etc. to enable effective targeted delivery of chemical substances to specific cells. In a more focused attempt, Thakur et al. [15] discussed the concept of respiratory delivery of favipiravir-tocilizumab combination through mucoadhesive protein-lipidic nanovesicles as a therapeutic tactic against COVID-19.
The carbon-based nanomaterials due to their exceptional structural and physical properties have attracted much of the research interest for use in biomedicine and, especially, in drug delivery technologies. Above the various allotropes of carbon, fullerenes [16] have been proposed as promising candidates for drug carriers, mainly because of their unique spherical molecular shape. There are numerous experimental [17][18][19][20][21][22][23] as well as theoretical reports [24][25][26][27][28][29][30][31][32][33], mainly grounded on DFT, associated with the analysis of drug delivery systems made of fullerene derivatives. Nevertheless, the relative computational approximations and predictions regarding the feasibility of using carbon nanomaterials as vehicles for anti-SARS-CoV-2 drugs are much fewer due to the freshness of the corresponding pandemic. In the field, Shahabi and Raissi [34], using Molecular Dynamics (MD) as a computational tool, conducted simulations of the transfer and delivery of the anti-SARS-CoV-2 drug Carmofur with the assistance of graphene oxide quantum dots. It was found that graphene oxide quantum dots offer efficient drug delivery around the catalytic region of the COVID-19 protein while promoting good drug penetration at the target location. Based on DFT instead, Bibi et al. [35] investigated the fullerene C 60 doped with metals Cr, Fe, and Ni as the carrier of the favipiravir anti-COVID-19 drug. In an analogous attempt, Yao et al. [36] explored the potential of using fullerene-like metal oxide nanocages as drug-carrying systems for favipiravir by estimating their drug adsorption efficiency. On the other hand, Parlak et al. [37] investigated the possibility of using a silicon-doped C 60 as a drug delivery carrier of molnupiravir by using DFT simulations. Specifically, they developed and optimized a SiC 59 molecular structure and, then, assumed that its silicon atom is externally bonded with an atom of a single molnupiravir molecule. They examined the bonding of the silicon atom of the doped fullerene with six different atoms of the drug molecule in order to find the ideal intermolecular interaction and optimum linking.
In the present study, a molecular system based on an appropriate cyclopropane derivative of C 60 fullerene was proposed and tested as a carrier agent of the molnupiravir drug. The idea is to use a functionalized molecular structure of C 60 fullerene, which is well soluble in water. In this way, not only the excellent and already proven drug-carrying capabilities of C 60 will be taken advantage of, but also the obtained system is expected to present good solubility, which is a key property and a required characteristic for drug delivery nanosystems. The MD simulations and corresponding energetic calculations are carried out to characterize the structural behavior and stability of the investigated water-soluble C 60 derivative/molnupiravir system. A carboxyfullerene known as dendro[60]fullerene is chosen as the basic structure of the delivery system, a molecule that was synthesized several years ago [38] and has already proved its good solubility in water [38,39]. In order to develop the proposed system, two molnupiravir molecules are attached to the two major building blocks of the dendro[60]fullerene via N-N single bonding. To test if the new drug delivery system has the ability to form a solution with water, the Gibbs free energy of solvation [40] is calculated using suitable MD modeling and simulations. In order to reach clear conclusions, the dendro[60]fullerene alone is also tested in water and n-octanol environments to be compared with the proposed dendro[60]fullerene/molnupiravir molecules drug delivery system. Despite the fact that there are several reports [28,35,37,[41][42][43][44][45] on the exploitation of the use of C 60 -based molecular systems against COVID-19, to the author's best knowledge, this is the first time that dendro[60]fullerene has been proposed and tested as a possible nanovehicle of molnupiravir.

Solubility and Partition Coefficient
A crucial design factor for many biochemical procedures and operations is the relative solubility of a substance in an aqueous and an organic solvent, which is usually expressed by the partition coefficient P. The partition coefficient P describes the tendency of an uncharged molecule to better dissolve in water or lipids and consequently to have a hydrophilic/lipophobic or lipophilic/hydrophobic nature, respectively. More simply, it defines the amount of a solute that dissolves in the water portion versus an organic portion. The most commonly followed tactic is to make estimations of the partitioning between water and n-octanol. In addition, it is more convenient to express the resulting partitioning by providing the constant logP, which is equal to the decimal logarithm of the partition coefficient P as follows: where [octanol] and [water] denote the concentration of solute in the n-octanol and water, respectively. It should be noted that obtaining a negative logP constant implies that the investigated compound is more hydrophilic, more polar, and has poorer lipid bilayer permeability. In contrast, a positive value for logP leads to the conclusion that the compound is more lipophilic, more nonpolar, and has poor aqueous solubility. Finally, when logP is close to zero, the molecule under investigation is similarly partitioned between the lipid and aqueous phases. The above short analysis makes clear that logP may be used extensively as a prognosticator of the interactions of drug components and products in the human body, which is separated into aqueous as well as lipid regions. Therefore, the logP constant may define the delivery and distribution of the drug molecules between the aqueous phase outside the human cell's membrane and the lipid phase inside the cell's membrane.

Solubility and Solvation Free Energy
The solubility of a compound may be expressed with respect to free energy by the difference of the Gibbs free energies [40] of the neutral form of a compound species in two different solvents, α, and β, as follows: or ∆G = RT ln(10) log(P α→β ) where R = 1.987 × 10 −3 kcal/mol/K is the gas constant, T is the temperature in K, and P α→β is the partition coefficient defined as the ratio of the number concentrations of the solute in the phases α and β, respectively. A very common method for approximating the free energy between subsequent states is to computationally calculate the derivative of the total energy, including both the potential and kinetic energy of the studied system. This theoretical approach is known as "thermodynamic integration" and is commonly used for free energy calculations based on computer simulations [46]. According to this numerical technique, the influence of the interaction between the solute and the solvent is progressively increased from zero to full scale in a finite number of steps. A complete molecular simulation must be performed separately for each step in order to compute the corresponding derivative of the solvationfree energy. When all simulations are finalized, the free energy may be defined by the numerical integration of all calculated derivatives. A very brief mathematical representation of the method is given below.
Let us assume a pair of end states α and β and that their corresponding Hamiltonians, i.e., total energies including both potential and kinetic energies, are denoted as H α (q, p; λ) and H β (q, p; λ), where q and p represent all of the positions and momenta of the system while λ is a coefficient taking values from 0 to 1. A total Hamiltonian H may be defined as follows: where f and g are functions of λ that are utilized to mix the two Hamiltonians, which are generally chosen such that H = H α for λ = 0 and H = H β for λ = 1. The Gibbs free energy difference between α and β states may now be calculated by the following relationship: The above integral may be approximated by performing a series of simulations regarding a discrete set of λ values and then adopting a numerical quadrature approach.
The calculation of the solvation-free energy is performed in three simulation stages. First, the ideal contribution is computed. During this stage, the solute molecule is assumed to be in a vacuum while the charges are steadily reduced to zero, whereas the rest of the interactions are kept the same. Second, the van der Waals contribution is calculated. In this kind of simulation process, the chargeless molecule is coupled to the solvent by activating the van der Waals interaction in several steps. The final stage is related to the approximation of the electrostatic contribution. Here, the charges of the solute molecular system are regenerated while the solute compound is progressively immersed in the solvent. Evidently, the total solvation free energy is equal to the sum of the ideal, van der Waals, and electrostatic aforementioned contributions.

Simulation Details
The present theoretical investigations are performed by utilizing force-field-based MD. The specific numerical technique has been proved to be reliable, and efficient, followed by a reasonable computational cost when used for studying pharmaceutical and drug applications at the molecular level [47]. The proposed MD simulations are grounded on the COMPASS II (Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies II) [48,49] based force-field, which consists of terms for bond stretching (b), angle bending (θ), dihedral angle rotating (ϕ), out-of-plane angle distortion (χ) as well as crossterms that account for interactions between the above four types of internal coordinates, and two non-bonded functions, a Coulombic function for electrostatic interactions due to the charges q and a 9-6 Lennard-Jones potential for van der Waals (vdW) interactions. The general form of the utilized potential is given by the following: or in a more analytical form as follows: In the above equation, r ij = r i − r j is the distance between atoms I and j having the atomic positions r i and r j , respectively. The subscripts 0 in the above parameters represent the corresponding reference values. The stiffness-like coefficients k, k 1 , k 2 , k 3 and k 4 are the COMPASS-II force field constants. The parameters q i and q j are the charges of atoms i and j, respectively. Finally, the functional terms r ij0 and ε ij are given, respectively, by the following: where ε i and ε j are the van der Waals well depths of the atoms i and j, correspondingly.
For the purpose of the present investigation, the Andersen thermostat and Berendsen barostat are adopted to control, respectively, the temperature and pressure of the system during the simulations [50]. In addition, the group-based summation method with cut-off radii of 12.5 Å and long-range corrections are assumed for the estimations of the vdW and the electrostatic interactions [50]. In all numerical tests, it is considered that the environment of the unit cells is characterized by a room temperature of T = 298 K and, if required, a pressure of p = 1 atm. The optimization of the unit cells is achieved by using the steepest descent algorithm under a convergence tolerance of 2 × 10 −5 kcal/mol, 10 −3 kcal/mol/Å, and 10 −5 Å regarding the energy, force, and distance variation, respectively.

Proposed Drug Delivery System
Because of their great physical properties, fullerenes have proved to be ideal nanomaterials in many biological applications [51]. Especially the Buckminsterfullerene C 60 , due to its symmetric hollow spherical shape, shares unique attributes and, thus, has gained much attention in the field of biosciences. The compact size of the fullerene C 60 makes it an ideal nanovehicle for drug delivery applications. Moreover, this nanomaterial, because of its physicochemical properties and size, is optimally suited to enter deep into the lung through the airways.
Regarding the possibility of using C 60 for carrying drug molecules in the human body, one important disadvantage is the fact that it presents significant hydrophobic behavior. In view of this limitation, here, a water-soluble C 60 derivative is proposed as a nano-carrier. Specifically, in order to enhance solubility and avoid aggregation during delivery, the dendro[60]fullerene is adopted as the basic structure of the drug delivery system, which is shown in Figure 1. The specific dendritic monoadduct of C 60 has 18 carboxylic groups in the periphery and may be synthesized by nucleophilic cyclopropanation of C 60 with a bis(polyamide)-malonate dendrimer and subsequent deprotection of the terminal t-butyl groups [38]. Generally, the higher the number of carboxylic acid groups, the higher the water solubility.
Regarding the possibility of using C60 for carrying drug molecules in the human body, one important disadvantage is the fact that it presents significant hydrophobic behavior. In view of this limitation, here, a water-soluble C60 derivative is proposed as a nano-carrier. Specifically, in order to enhance solubility and avoid aggregation during delivery, the dendro[60]fullerene is adopted as the basic structure of the drug delivery system, which is shown in Figure 1. The specific dendritic monoadduct of C60 has 18 carboxylic groups in the periphery and may be synthesized by nucleophilic cyclopropanation of C60 with a bis(polyamide)-malonate dendrimer and subsequent deprotection of the terminal t-butyl groups [38]. Generally, the higher the number of carboxylic acid groups, the higher the water solubility. The next step is the selection of an efficient drug against COVID-19, capable of being carried and delivered by dendro[60]fullerene. The drug molecule that has been preferred is molnupiravir, an antiviral drug officially indicated for the treatment of COVID-19 [5,6]. Special attention is required regarding the link between the proposed nanovehicle and the drug molecule/es. The use of smart linkers between the two compounds that may respond and react due to a stimulus is very crucial for the effective on-site and on-purpose release of the drug [52]. The ideal linker selection requires substantial experimental effort using a variety of environmental conditions, representing the different human body regions. In this study, the main idea is to adopt input positions for the drug that offer simple linkage as well as unique and sensitive response properties in the connection region. The most convenient way to achieve these two features is to use a simple bond as a linker, which is characterized by a distinctively small dissociation energy. The introduction of a linkage of small strength in comparison with the rest of the bonds in the system will allow the easy and controllable release of the drug molecules, avoiding structural modifications of any of the two coupled basic compounds at the drug release instant. By observing in Fig  The next step is the selection of an efficient drug against COVID-19, capable of being carried and delivered by dendro[60]fullerene. The drug molecule that has been preferred is molnupiravir, an antiviral drug officially indicated for the treatment of COVID-19 [5,6]. Special attention is required regarding the link between the proposed nanovehicle and the drug molecule/es. The use of smart linkers between the two compounds that may respond and react due to a stimulus is very crucial for the effective on-site and on-purpose release of the drug [52]. The ideal linker selection requires substantial experimental effort using a variety of environmental conditions, representing the different human body regions. In this study, the main idea is to adopt input positions for the drug that offer simple linkage as well as unique and sensitive response properties in the connection region. The most convenient way to achieve these two features is to use a simple bond as a linker, which is characterized by a distinctively small dissociation energy. The introduction of a linkage of small strength in comparison with the rest of the bonds in the system will allow the easy and controllable release of the drug molecules, avoiding structural modifications of any of the two coupled basic compounds at the drug release instant. By observing in Figure 2a the dendro[60]fullerene in contrast with the molecular structure of the molnupiravir, one may notice that both compounds contain similar −N−H (−Nitrogen−Hydrogen) subunits. Thus, it becomes evident that the single N−N bond could be proposed and tested as a possible simple linker at these locations, as depicted in Figure 2b. To retain the molecular symmetry of the carrier, two corresponding −N−N− links are created with two drug molecules in opposing positions. Although the dendro[60]fullerene contains 6 more −N−H subunits in its periphery, drug molecule inputs in these sites are avoided, in order not to disturb the dendritic formations at the two edges of the system, which would probably harm its overall solubility.
Although the N-N bond has a relatively small dissociation energy of 240 kJ/mol, this value is only 20% smaller than that corresponding to the C-N bond, which is the next weakest link existing in the proposed system [53]. In addition, it is worth mentioning at this point that over the last few years, plenty of reactions capable of forming the nitrogen−nitrogen bond in natural product biosynthesis have been revealed, overcoming the technical difficulties due to the inherent nucleophilic nature of nitrogen elements [54]. Moreover, some nitrogen-nitrogen-bond-containing bioactive products have already been used in clinical practice, evidencing their promising long-life capabilities [54]. Although the N-N bond has a relatively small dissociation energy of 240 kJ/mol value is only 20% smaller than that corresponding to the C-N bond, which is the weakest link existing in the proposed system [53]. In addition, it is worth mentioni this point that over the last few years, plenty of reactions capable of forming the n gen−nitrogen bond in natural product biosynthesis have been revealed, overcomin technical difficulties due to the inherent nucleophilic nature of nitrogen elements Moreover, some nitrogen-nitrogen-bond-containing bioactive products have alr been used in clinical practice, evidencing their promising long-life capabilities [54].

Results and Discussion
The main scope of the computational study is the characterization of the solubil the proposed dendro[60]fullerene/molnupiravir drug delivery system through solva free energy MD computations. For comparison, the pure water-soluble dendro[60]fu ene is computationally investigated in parallel. Both compounds are tested in water n-octanol (1-octanol) to enable calculations via Equation (3).
Firstly, all the molecular units are optimized at room temperature using the ste descent algorithm and COMPASS II potential to obtain the corresponding lowest en stable structures. Figure 3 illustrates the optimized molecular structure of the dendro[60]fullerene. The obtained optimized potential energy of this C60 derivati

Results and Discussion
The main scope of the computational study is the characterization of the solubility of the proposed dendro[60]fullerene/molnupiravir drug delivery system through solvationfree energy MD computations. For comparison, the pure water-soluble dendro[60]fullerene is computationally investigated in parallel. Both compounds are tested in water and n-octanol (1-octanol) to enable calculations via Equation (3).
Firstly, all the molecular units are optimized at room temperature using the steepest descent algorithm and COMPASS II potential to obtain the corresponding lowest energy stable structures. Figure 3 illustrates the optimized molecular structure of the pure dendro[60]fullerene. The obtained optimized potential energy of this C 60 derivative is E C60der = −1882.09 kcal/mol. In addition, Figure 4 shows the equilibrated three-dimensional atomistic structure of a single molnupiravir molecule, which leads to a potential energy value of E Molnup = −39.28 kcal/mol. Two optimized molnupiravir molecules are appropriately attached to the two major building blocks (see Figure 2) of the adopted optimized nanovehicle. A new optimization process is carried out for the entire system this time, which leads to the molecular structure of Figure 5. The arisen optimized drug delivery molecular system that is created via the connection of the dendro[60]fullerene and two molnupiravir molecules via nitrogen-nitrogen bonding has a total energy of E drugsys = −2008.80 kcal/mol. According to the above energy values, the binding energy between the nanocarrier and the two drug molecules is equal to E bind = E drugsys − (E C60der + 2E Molnup ) = −48.15 kcal/mol, the negative value of which indicates that the interaction between the nanocarrier and the drug is reasonably stable.
value of EMolnup = −39.28 kcal/mol. Two optimized molnupiravir molecules are appropriately attached to the two major building blocks (see Figure 2) of the adopted optimized nanovehicle. A new optimization process is carried out for the entire system this time, which leads to the molecular structure of Figure 5. The arisen optimized drug delivery molecular system that is created via the connection of the dendro[60]fullerene and two molnupiravir molecules via nitrogen-nitrogen bonding has a total energy of Edrugsys = −2008.80 kcal/mol. According to the above energy values, the binding energy between the nanocarrier and the two drug molecules is equal to Ebind = Edrugsys − (EC60der + 2EMolnup) = −48.15 kcal/mol, the negative value of which indicates that the interaction between the nanocarrier and the drug is reasonably stable.   value of EMolnup = −39.28 kcal/mol. Two optimized molnupiravir molecules are appropriately attached to the two major building blocks (see Figure 2) of the adopted optimized nanovehicle. A new optimization process is carried out for the entire system this time, which leads to the molecular structure of Figure 5. The arisen optimized drug delivery molecular system that is created via the connection of the dendro[60]fullerene and two molnupiravir molecules via nitrogen-nitrogen bonding has a total energy of Edrugsys = −2008.80 kcal/mol. According to the above energy values, the binding energy between the nanocarrier and the two drug molecules is equal to Ebind = Edrugsys − (EC60der + 2EMolnup) = −48.15 kcal/mol, the negative value of which indicates that the interaction between the nanocarrier and the drug is reasonably stable.   Next, representative unit cell models are constructed for the evaluation of the interactions of the dendro[60]fullerene ( Figure 3) and the dendro[60]fullerene/molnupiravir system ( Figure 5) in water and n-octanol solvents. As expected, for the solute molecules presented in Figures 3 and 5, the optimized state of the solvent molecules should be ob- Next, representative unit cell models are constructed for the evaluation of the interactions of the dendro[60]fullerene ( Figure 3) and the dendro[60]fullerene/molnupiravir system ( Figure 5) in water and n-octanol solvents. As expected, for the solute molecules presented in Figures 3 and 5, the optimized state of the solvent molecules should be obtained as well. The steepest descent optimization leads to the equilibrium molecular structures of water and n-octanol, which are illustrated in Figure 6a,b, respectively. Next, representative unit cell models are constructed for the evaluation of the interactions of the dendro[60]fullerene ( Figure 3) and the dendro[60]fullerene/molnupiravir system ( Figure 5) in water and n-octanol solvents. As expected, for the solute molecules presented in Figures 3 and 5, the optimized state of the solvent molecules should be obtained as well. The steepest descent optimization leads to the equilibrium molecular structures of water and n-octanol, which are illustrated in Figure 6a,b, respectively. The representative amorphous simulation box for each of the investigated solvent/solute systems is generally constructed as described below. Initially, a periodic and symmetric cubic unit cell with a centrally positioned single solvent molecule is developed by assuming a large enough lattice length of 50 nm. The cubic unit cells containing the dendro[60]fullerene and the drug delivery molecular system are illustrated in Figures 7a  and 8a, respectively. Then, a Connolly surface is created to define the boundary between the solvent molecular structure and its environment, which is going to be filled with solute molecules. The developed Connolly surfaces for the dendro[60]fullerene and the dendro[60]fullerene/molnupiravir compound may be seen in Figure 7b and Figure 8b. The representative amorphous simulation box for each of the investigated solvent/ solute systems is generally constructed as described below. Initially, a periodic and symmetric cubic unit cell with a centrally positioned single solvent molecule is developed by assuming a large enough lattice length of 50 nm. The cubic unit cells containing the dendro[60]fullerene and the drug delivery molecular system are illustrated in Figures 7a and 8a, respectively. Then, a Connolly surface is created to define the boundary between the solvent molecular structure and its environment, which is going to be filled with solute molecules. The developed Connolly surfaces for the dendro[60]fullerene and the dendro[60]fullerene/molnupiravir compound may be seen in Figures 7b and 8b. The remaining empty volume in the simulation box is packed with a number of solvent molecules. The population of the solvent molecules is progressively increased in the empty unit cell volume until the theoretical density of the solvent under consideration is reached. Note that a density value of 1 g/cm 3 and 0.824 g/cm 3 is assumed for the water and n-octanol phase, respectively. In addition, a geometry optimization procedure is carried out followed by energy minimization until a stable structure of the representative simulation box is achieved. To further relax the simulation boxes and to obtain the true density of the solute/solvents systems, a dynamics analysis is performed using the NPT ensemble and a time step of 1 fs. An external pressure of 1 atm and a temperature of 298 K are applied for 100 ps to equilibrate the simulation boxes under normal environmental conditions. Figure 9 illustrates the density variation of the unit cells over 100 ps for the four tested solute/solvent combinations. The results show that both dendro[60]fullerene and the drug delivery molecular system when immersed in water and n-octanol converge to a density value of about 1.15 g/cm 3   The remaining empty volume in the simulation box is packed with a number of solvent molecules. The population of the solvent molecules is progressively increased in the empty unit cell volume until the theoretical density of the solvent under consideration is reached. Note that a density value of 1 g/cm 3 and 0.824 g/cm 3 is assumed for the water and noctanol phase, respectively. In addition, a geometry optimization procedure is carried out followed by energy minimization until a stable structure of the representative simulation box is achieved. To further relax the simulation boxes and to obtain the true density of the solute/solvents systems, a dynamics analysis is performed using the NPT ensemble and a time step of 1 fs. An external pressure of 1 atm and a temperature of 298 K are applied for 100 ps to equilibrate the simulation boxes under normal environmental conditions. Figure 9 illustrates the density variation of the unit cells over 100 ps for the four tested solute/solvent combinations. The results show that both dendro[60]fullerene and the drug delivery molecular system when immersed in water and n-octanol converge to a density value of about 1.15 g/cm 3  To compute the Gibbs free energies at room temperature through Equation (5), a series of simulations for the four unit cells of Figure 7c,d and 8c,d is performed by increasing the coupling coefficient λ from 0 to 1 using a stable increment of Δλ = 0.05. Note that for each new increment of the coupling coefficient, a full independent MD analysis is performed for 100 ps under the NVT ensemble, to achieve the required equilibration of the system under consideration. A time step of 1 fs is adopted for all simulations. At the end of each NVT analysis, an optimization process is carried out to further minimize, if possible, the potential energy of the system. The simulations for each unit cell are associated with the separate computation of the ideal energy, van der Waals, and electrostatic contribution. As already mentioned, the sum of these contributions expresses the corresponding Gibbs free energy. Given the above simulation details, it becomes evident that the free energy calculations via the MD method require a substantial computational effort. Figures 10 and 11 illustrate the free energy contributions (total, ideal, van der Waals, electrostatic) of the two considered solutes, i.e., the drug delivery system as well as the water-soluble dendrimer fullerene, with respect to the coupling parameter increase in the water and the n-octanol, respectively. Evidently, the final free energy values correspond to λ = 1. To compute the Gibbs free energies at room temperature through Equation (5), a series of simulations for the four unit cells of Figure 7c,d and Figure 8c,d is performed by increasing the coupling coefficient λ from 0 to 1 using a stable increment of ∆λ = 0.05. Note that for each new increment of the coupling coefficient, a full independent MD analysis is performed for 100 ps under the NVT ensemble, to achieve the required equilibration of the system under consideration. A time step of 1 fs is adopted for all simulations. At the end of each NVT analysis, an optimization process is carried out to further minimize, if possible, the potential energy of the system. The simulations for each unit cell are associated with the separate computation of the ideal energy, van der Waals, and electrostatic contribution. As already mentioned, the sum of these contributions expresses the corresponding Gibbs free energy. Given the above simulation details, it becomes evident that the free energy calculations via the MD method require a substantial computational effort. Figures 10 and 11 illustrate the free energy contributions (total, ideal, van der Waals, electrostatic) of the two considered solutes, i.e., the drug delivery system as well as the water-soluble dendrimer fullerene, with respect to the coupling parameter increase in the water and the n-octanol, respectively. Evidently, the final free energy values correspond to λ = 1.
As implied by the data of both diagrams in Figures 10 and 11, for all the investigated solute/solvent system cases, the ideal contribution of the solvation-free energy is very close to the absolute value of the electrostatic contribution. This occurs because the ideal contribution is actually the free energy cost of removing charges from the solute molecular system. On the other hand, van der Waals's contribution gets lower values while it presents a mixed ascending/descending behavior as the coupling coefficient is increased. The Gibbs free energy of the drug delivery system is about −51.15 kcal/mol and −53.14 kcal/mol in water and n-octanol, respectively, meaning that it is somewhat more soluble in n-octanol than in water. On the other hand, the dendro[60]fullerene presents solvation-free energy of −59.52 kcal/mol and −62.33 kcal/mol in water and n-octanol respectively. The hydration free energy of dendro[60]fullerene (−59.52 kcal/mol) is more negative than the −51.15 kcal/mol value, a fact that proves that it has a slightly higher water-solubility compared to the same compound with two additional molnupiravir molecules on its structure. By substituting the above computed Gibbs free energy values into Equation (3), one may find that the constant logP for the proposed drug delivery system and the dendro[60]fullerene is equal to 1.46 and 2.06, respectively. These values prove that the proposed drug-delivery system not only retains the good water-solubility of the nanocarrier alone but seems to be even less lipophilic, suggesting that it is a good candidate for delivering the antiviral product molnupiravir in certain regions of the human body. Note that the unmodified C 60 fullerene molecule, which is highly lipophilic, according to experimental evidence presents an octanol-water partition coefficient of about 6.67 [55]. Unfortunately, since dendro[60]fullerene joined with molnupiravir is studied here for the first time, more focused comparisons with other relevant available studies [28,35,37,[41][42][43][44][45]    As implied by the data of both diagrams in Figures 10 and 11, for all the investigated solute/solvent system cases, the ideal contribution of the solvation-free energy is very close to the absolute value of the electrostatic contribution. This occurs because the ideal contribution is actually the free energy cost of removing charges from the solute molecular system. On the other hand, van der Waals's contribution gets lower values while it presents a mixed ascending/descending behavior as the coupling coefficient is increased. The  As implied by the data of both diagrams in Figures 10 and 11, for all the investigated solute/solvent system cases, the ideal contribution of the solvation-free energy is very close to the absolute value of the electrostatic contribution. This occurs because the ideal contribution is actually the free energy cost of removing charges from the solute molecular system. On the other hand, van der Waals's contribution gets lower values while it presents a mixed ascending/descending behavior as the coupling coefficient is increased. The

Conclusions
A more lipophilic drug delivery system against COVID-19 may be sequestered by fatty tissue and therefore difficult to excrete and reach its target for drug release. For this reason, a new idea for a drug delivery system has been proposed here. According to the proposed concept, a highly water-soluble derivative of the C 60 , i.e., the dendro[60]fullerene has been used as a nanocarrier and two molnupiravir molecules as the drug agents. The hydrophilic nanocarrier has been linked with two drug molecules via simple N−N bonds, which may easily break under a stimulus due to their low dissociation energy. This would permit the effective synchronized release of the drug to the required target organs. The MD-based computations of the solvation free energy of the proposed drug delivery system, in contrast to the dendro[60]fullerene alone, proved the effective solubility of the proposed nanovehicle. More theoretical and experimental research is required in order to validate and adjust the proposed C 60 derivative-based drug delivery scheme.