Fluorinated derivatives of tetrahydroaltersolanol molecule on COVID-19, HIV, and HTLV protease by DFT and molecular docking approaches

Structural, optoelectronic, and biological properties of tetrahydroaltersolanol (C16H20O7) and fluorinated derivatives are calculated using density functional theory (DFT) and molecular docking approaches. It is shown that the pure C16H20O7 molecule has a direct HOMO–LUMO energy gap about 3.1 eV. The substitution of F atom at A category decreases the electronic energy gap, while it is constant at B category. In A category, the behavior of the pure molecule changed from insulator to semiconductor with various substitution of F atom. The electronic properties were depended on the F sites in the pure molecule. The molecular electrical transport properties and charge-transfer possibilities increase with decreasing energy gap. The pure C16H20O7 molecule with high energy gap has low chemical reactivity and substitution of F atom at considered molecule increases chemical reactivity. Obtained results show that F-O bonds in trifurcation bonds of C16H19O7(F14), C16H19O7(F16), and C16H19O7(F17) molecules play a key role in confronting with COVID-19, HIV, and HTLV proteases, respectively. Optical spectra, such as the dielectric functions, electron energy-loss spectroscopy, refractive index, extinction coefficient, and reflection spectra show that fluorinated derivatives of C16H20O7 at B category can be used in the new drugs. Graphical abstract


Introduction
Fungi are found in all environments with a universal spread and grow in different habitats [1][2][3]. Hydroanthraquinones, alkaloids, polyketides, and peptides are bioactive materials which produce by species of fungi and some of them are toxic to plants, insects and animals. Fungal metabolites have various applications in medicine as antiviral, antibacterial, and anticancer compounds [4][5][6]. The reduced species of anthraquinones that are formed in fungi such as Alternaria solani [7,8], A. porri [9], Dactylaria lutea [10], Phomopsis juniperivora [11], and Stemphylium botryosum [12] are called altersolanol, which contains many different types and exhibits interesting properties. For instance, altersolanols exhibited significant antimicrobial, antiprotozoal, and cytotoxic activities [13]. Marine fungi have recently emerged as a newfound source of fungal natural products and maybe have a potentially crucial point in drug discovery.
Genetic material (DNA and RNA) of the virus encodes viral protease as essential enzymes to complete the viral infections cycle. These enzymes cleave selectively the scissile peptide bonds with different catalytic mechanisms in viral or cellular protein. Serine, cysteine, or aspartic acid residues in the active site of viral proteases recognize the specific peptide bond that is often located within conserved sequence motifs. The 2019-nCOV as novel coronavirus has main cysteine protease named (M pro , 3C) [14], HIV and HTLV protease [15]. The aspartic acid protease family have potential drug targets because of their essential activities. One of the main methods commonly used to study the physical and chemical properties of materials is density functional theory (DFT) approach [16] as implemented in the Wien2k code [17]. Fluorine with high electronegativity in carbon-fluorine bond plays a key role to create a Maliheh Azadparvar and M. Kheirabadi contributed equally to this work. strongly polarized bond with the highest energies. A hydrogen atom or a hydroxyl group can displace with fluorine [18,19]. Fluorine is rarely found in natural compounds, and monofluorinated compounds are produced by some plants and bacteria [20,21]. They often play a key role in defense strategy due to their toxic properties. Nowadays, fluorinated compounds have various applications in agriculture and medicine [22]. On the other hand, hydroanthraquinone derivatives as bioactive metabolites were obtained from a soft coral-derived Alternaria sp. fungus [23]. One of them is a tetrahydroaltersolanol with C 16 [23].
In this study, we have investigated the structural, optoelectronic, and biological properties of tetrahydroaltersolanol (C 16 H 20 O 7 ) [23] and fluorinated derivatives by using density functional theory (DFT) and molecular docking approaches. The 19 fluorinated derivatives of tetrahydroaltersolanol are constructed via hydrogen substitution by fluorine. The optoelectronic and ADME properties were predicted by the Wien2k code and Molsoft web server, respectively. Finally, binding mode of tetrahydroaltersolanol and designed compounds to HIV, HTLV protease, and COVID-19 protease were obtained by molecular docking and compared with the obtained DFT results.

A. The FP-LAPW method as implemented in WIEN2k
In the first principal calculations, the full-potential linearized augmented plane wave method (FP-LAPW) is used for simulation of the structural and optoelectronic properties of pure C 16 [16,17,24]. The exchange-correlation potentials are calculated using the generalized gradient approximation (GGA). The tetrahydroaltersolanol (C 16 H 20 O 7 ) compound is crystalized at orthorhombic structure with the "P 21 21 21" space group, as shown in Fig. 1. The lattice constants are a = 5.6702 Å, b = 10.8579 Å, and c = 23.9550 Å. The number of k-points is selected about 1000 for self-consistent calculations. In the FP-LAPW calculations, the cut-off energy, which defines the separation of core and valence states is also selected − 6 Ry with the convergence parameter of 0.0001 Ry.

B. Molecular docking and molecular dynamic
Molecular docking is a computation process by which ligands are moved in a 3D space so that, finally, the best space structures of the target and ligand positions are found by getting the highest scores in the functions. Before docking process, to predict drug-likeness of compounds, ADME properties of the compounds were predicted using Molsoft web services (https:// molso ft. com/ mprop/) which estimates the absorption, distribution, metabolism, and excretion of candidate compounds. For docking process, proteins and ligands optimized. At the first step, primary ligand optimization was performed by Hyperchem software [25], a salvation box was created for them, and molecular dynamic simulation was done by amber force field [25] under the temperature of 300 K in 10 picoseconds in order to achieve the most sustainable state and the lowest energy of the structure of inhibitors for being used in the future studies. Then, these structures were saved as PDB files to be used in the docking process [26]. At the second step, the crystallography structure of HIV, HTLV, and COVID-19 protease enzymes were selected from the protein data bank with the respective PDB codes of 1YT9, 3LIN, and 6M03. Then, all the nonamino acid molecules (water, ion, ligand, inhibitor, etc.) were removed. Then, the structure of the HTLV and HIV protease's optimization was done by Gromacs 5.1.2 software [27] in a quasi-natural condition in terms of temperature, atmosphere pressure, pH, and the existence of water as a solvent to get their balanced structure. First, the topology of protein was built by using the GROMOS96 53a5 force field. Then, a cubic box was considered between the protein and the edge of the box with the distance of 1.0 nm, and the water molecules were added to the box as solvent by using simple point charge (SPC) model. Then, four chlorine ions were added for neutralizing the system load. In the following, energy minimization was done by the steepest descent algorithm. Then, NVT equilibration and NPT equilibration were done each for 100 picoseconds, and finally, molecular dynamic simulation was done for 10,000 picoseconds (10 ns). The docking process was done by Auto Dock 4.2 software. In the docking process, LGA (Lamarckian genetic algorithm) was used and all the non-rotatable amide bonds were selected. The cubic box with the size of 126 × 126 × 126 Å was set to cover all the protein-ligand complexes. The number of GA runs for each 250 docks was selected by 25,000,000 energy measurement [28]. The content of hydrogen bonds and hydrophobic interactions between ligand and protein in the complex was analyzed by PLIP server. Data were drawn with PyMol software [29].

Results and discussion
Today, in medicine, plants and fungi are used as rich sources of various metabolites to treat various diseases. The resulting metabolites can have inhibitory effects on the cancer cells or viral proteins. On the other hand, rare natural f luorinated compounds are produced by plants and bacteria that play a major role in defense mechanisms. These compounds have attracted special attention as effective drug compounds. Fluorinated derivatives of chemical compounds are usually toxic; still it should also be borne in mind that the enzymes that break down fluorinated compounds are found only in small numbers in some organisms. Therefore, the necessary precautions should be taken in their therapeutic uses. One of the hydroanthraquinone derivatives extracted from marine fungi with antiviral activity against porcine reproductive and respiratory syndrome virus (PRRSV) with IC 50 values 65 micromolar (μM) [23] was selected to design its f luorinated derivatives. The hydrogen atoms in hydroanthraquinones (Cif code: 1,503,469) were substituted with fluorine (Fig. 2). As the first step to design new compounds in lab, they should be energetically stable. To estimate stability of designed fluorinated substance, electronic and optical properties of them was calculated. After that, drug likeness and ADME properties of fluorinated substances was predicted. Finally, binding mode of compounds to HTLV, HIV, and COVID-19 protease was obtained.

A. Optoelectronic results
The physical and chemical properties of materials are studied by investigating the behavior of electrons near the fermi level. The electrons in fermi level play a key role in the reactive properties of drugs. In this part, an attempt has been made to obtained correlation between DFT calculations and pharmaceutical properties. It is shown that several computational parameters such as highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) band gap and absorption spectra were predicted in order to describe the pharmaceutical properties of tetrahydroaltersolanol (C 16 H 20 O 7 ) and fluorinated derivatives. The calculated  Table 1. For better presentation of the obtained results, all compounds are classified into one of two groups: category A or category B. Category A involves the fluorine atoms that joined to the oxygen atom, while in category B, the fluorine atom is attached to carbon that are located outside benzene rings with three hydrogen bonds or on benzene rings. The results show that the C 16 H 20 O 7 compound has a direct bandgap about 3.1 eV (Fig. 3) Fig. 4. It is seen that the contribution of s-states is small and p-states have a remarkable contribution. The s-states of hydrogen atom and p-states of carbon atom have a small contribution in the total density  (Fig. 5). Figure 6 shows that the p-states of the fluorine atom have the high contribution near the Fermi level at HOMO and  (F14 to F19), the p-states of the oxygen and fluorine atoms play a key role in the pharmaceutical properties of compounds.
To analyze the optical properties of a material, we must first obtain its dielectric function, which consists of both real and imaginary parts (ε 1 (ω), ε 2 (ω)) [30]. It is the linear response of a system to electromagnetic waves. This function is calculated by considering the transition between the occupied and unoccupied states in the energy band or the transition between HOMO and LUMO states. The imaginary and real parts of the dielectric function are defined as [31]: where P represents the Cauchy integral part. � � c k ⟩ and � � v k ⟩ are the electronic states at the conduction band and the valence band, respectively. In Fig. 7a and b, the real and imaginary parts of the dielectric function for some selected compounds such as C 16 (F18) are plotted at x, y, and z directions. It is also seen from the ε 1 (ω) spectra that obtained values are positive for all photon energy ranges (Fig. 7a). Hence, the dielectric nature of these compounds is identified. As seen, calculated ε 1x (ω) of all compounds have small variations between 2 and 4, except C 16 H 19 O 7 (F4), while ε 1y (ω) and ε 1z (ω) vary between 0 and 6. So, there is an anisotropy between three directions of the ε 1 (ω) spectrum corresponding to different polarizations. This polarization is due to the response of material to the electrical field of electromagnetic wave. The maximum value of polarized ε 1x (ω), ε 1y (ω), ε 1z (ω) for the C 16 H 20 O 7 and selected B category of compounds occurred at 3.1 eV. The first peaks are related to maximum values of variations ε 2 (ω) spectrum (Fig. 7b) (F4) in A category of compounds, a maximum polarization of ε 1 (ω) spectra is observed at the x-direction. The first peak of ε 2x (ω), ε 2y (ω), and ε 2z (ω) spectra is occurred around 1 eV due to the electronic transitions between the F-p state in HOMO and the O-p state in LUMO. For C 16 H 19 O 7 (F13) compound, the first peaks of ε 1y (ω) and the ε 2y (ω) spectra at energy rang of 0-3 eV appears due to the small bandgap of compound and electronic transition from valence band to the unoccupied states in the conduction band. The calculated static dielectric constants at zero energy ε 1 (0) are presented in Table 1.
For B category of compounds, the maximum static dielectric constant ε 1 (0) is observed at the y-direction, while for C 16 H 19 O 7 (F4) compound, the maximum value of the ε 1 (0) is obtained at the x-direction and it is about 6.35. It is found that the fluorination at B category of compounds does not effect on ε 1 (0) of compounds. As a result, a more detailed analysis showed that these Fig. 7 a The real part of dielectric spectrum in terms of energy at three directions: x, y, z. b The imaginary parts of dielectric spectrum at x, y, and z directions compounds exhibit the same mechanism of reaction between the drug molecules and the diseases proteases.
Electron energy-loss spectroscopy (EELS) is an analytical technique that evaluates the losing energy of an electronic beam during interaction with the substance. EELS corresponded to the HOMO and LUMO transitions of electrons and their elastic scattering [30]. It is an optical spectrum [32] used to get Plasmon energy. EELS is applied to identify [33] the different (organic) compounds used in biological systems or pharmaceuticals and to detect any impurities present [34] even during the synthesis or formulation process [35]. Also, EELS is used as a technique for quantification of penetration to cell walls [36]. EELS defined as: Calculated EELS spectra for C 16 Fig. 8 at the x-, y-, and the z-directions. According to this figure, the prominent peaks in ELOSS spectra are located at the energy of 24 eV for pure and fluorinated compounds of B category. The maximum values of peaks are about 2.3, 2.4, and 2.2 at the x-, y-, and the z-direction, respectively. The ELOSS spectrum for B category of compounds is the same as pure (Fig. 8b). The value of main peak of C 16 H 19 O 7 (F4) is 1.9 at energy of 26 eV.
(3) ELOSS = −Im 1 Absorption spectroscopy [37] is defined as the measurement of the interaction between electromagnetic radiation, and a substance is another important analytical chemistry technique to ensure the safety, quality, and efficacy of compounds used in the pharmaceutical industry. Especially, ultraviolet/visible (UV/Vis) light spectroscopy can be used in the manufacturing process to identify contaminants within a substance or measure the kinetics of a reaction. In this study, we have investigated the absorbance spectra of pure and fluorinated C 16 H 20 O 7 compounds, as in the following: Figure 9 displays the absorbance spectra of the pure and selected fluorinated compounds at the x-, y-, and the z-directions. The maximum absorption is in the wavelength range of 50 to 125 nm that is at the extreme ultraviolet radiation (EUV) region. As shown, all compounds have the considerably same behavior in these regions. In ultraviolet (UV) regions in the wavelength range of 200 nm to 400 nm, the absorption coefficients are less than the EUV region (inset panels of Fig. 4). Also, there are some slight differences among the absorbance spectra of compounds. High similarity in absorbance spectrum of fluorinated compounds with pure is crucial from the viewpoints of pharmaceutical usages.    Fig. 10 a and b The calculated refraction index and extinction coefficient spectra for pure C 16 H 20 O 7 at x, y, z directions. The inset plots are related to the comparison of refraction index and extinction coefficient spectra among selected six compounds at the energy range 0 to 7 Ev One of the most known physicochemical properties of substances is the refractive index [38]. The refractive index, that means the ratio of the speed of light passing in the vacuum to its passing through a substance, is complex quantity n( ) = n( ) + k( ), where n(ω) is the refractive index for the real part and k(ω) is the extinction coefficient for the imaginary part. We have calculated as follows: Refractive index is mainly applied to identify a particular substance, confirm its purity, or measure its concentration. Generally, it is used to measure the concentration of a solute in an aqueous solution. So, it has the most applications to the determination of drug concentration in the pharmaceutical industry.
We have presented the refractive index and extinction coefficient spectra in terms of photon energy for C 16 H 20 O 7 compound, as shown in Fig. 10a and b. The calculated n( )andk( ) spectra exhibit some anisotropic behavior at the energy range 0 to 7.5 eV. For high-energy photons (higher than 16 eV), n( ) is lesser than unity which means that the nature of the compound changes from linear to non-linear, and the superluminal effect occurs. The inset graphs are related to comparing of refraction index and extinction coefficient spectra among selected compounds at energy range 0 to 7 eV. Generally, where n and k are refractive index and extinction coefficient, respectively. As seen in Fig. 11, the main peak about 0.33 at 0.5 eV, is related to the C 16 H 19 O 7 (F4) compound from category A of compounds at the x-direction. Other compounds have not considerable reflectivity at the x-direction. For the other compounds, the maximum value of reflectivity of 0.3 at 3.25 eV is observed at the z-direction. For pure C 16 H 20 O 7 compound, large optical reflectivity anisotropy is observed at 6.5 eV in the x-, y-, z-directions. On the other hand, the behaviors of reflectivity spectra for the pure and the other B category of compounds, except C 16 H 19 O 7 (F13), are quite similar at the y-direction.
In summary, in comparison with the experimental results, there is not the available experimental Fig. 11 The calculated reflectivity spectra at x, y, and z directions   [39] providing stronger evidence that it makes these materials suitable for pharmaceutical applications. In Fig. 12, we have compared obtained absorption spectrum for Altersolanol L with the experimental results for the same family of compounds for Altersolanol A and Altersolanol J compounds [39]. There is a good agreement between obtained absorption spectra at 224 and 275 nm with the other experimental results.

B. Prediction of drug-likeness
To predict drug-likeness of designed compounds, Molsoft software was used to estimate the pharmacological similarity of fluorinated derivatives with the five famous Lipinski's rules [40]. Lipinski's rule of five: As per Lipinski's rule of five [40], an orally active drug has (i) not more than five hydrogen bond donors (OH and NH groups), (ii) not more than ten hydrogen bond acceptors, (iii) a molecular weight under 500, and (iv) a partition coefficient log P under 5. The capacity of the hydrogen bond is an essential factor in the permeability of the drug to the body [41]. Drug permeation or absorption in the body decreases with increasing the number of hydrogen bond donors to more than five and the number of hydrogen bond acceptors to more than 10. The designed compounds have 7 H-bond donors and 5 to 7 H-bond acceptors (Table-2). Lipophilicity and partition coefficient (log P) values are an indicator of the tendency of the compound to dissolve in the membrane. The values describe the uptake, distribution, biotransformation, and excretion of organic chemicals in biological systems [42]. High log P values represent high solubility and good penetration of lipid membranes but low solubility in water and cytoplasm. Log P values for C 16 H 20 O 7 and its fluorinated derivatives are approximately − 4 ( Table 2) confirming higher affinity to the aqueous phase.
Blood-brain barrier (BBB) is a physicochemical barrier that maintains homeostasis of the central nervous system (CNS) to regulate the content of the interstitial fluid in the brain [42]. The blood-brain barrier prediction score is greater than four indicates the chemical pass the BBB. This score was calculated in Molsoft server by Gupta Fig. 12 Comparison of calculated absorption spectrum Altersolanol L at the x-direction (in unit of 10 -4 /Cm) with the experimental absorption spectra of Altersolanol A and Altersolanol J (without unit) [39] Fig. 13 RMSD and total energy for HTLV protease during 10 ns molecular dynamics simulations. a Root mean square deviations (RMSDs) of HTLV-1 protease (purple) and HIV protease (light purple). b Total energy of HTLV-1 protease (light purple) and HIV protease (purple). Both parameters confirm structural stability of protease during simulation et al. [43]. In this study, the BBB scores for C 16 H 20 O 7 and its fluorinated derivatives are 2.3 to 2.4 (Table 2) confirming them as unpenetrated chemicals to the brain.
The aqueous solubility of compounds is significant to oral absorption and distribution characteristics. All designed compounds have aqueous solubilities of Mol-LogS from − 0.7 to − 0.5. Based on ADME prediction, the C 16  In them, fluorine is substituted with hydrogen in the methyl group (C16). To get more optimal and reliable results from the docking process, first, the 3D structure of HTLV-1 protease (PDB codes: 3LIN) and HIV protease (PDB code: 1YT9) were optimized by Gromacs in 10 ns. The RMSD results (root mean square deviation) and total energy showed that the protein structure during running time is stable (Fig. 13). HTLV-1 and HIV have a homodimer activity structurally, and each monomer of these enzymes has 125 and 99 amino acids, respectively. The active site of these proteases is located in a valley between two monomers ASP25, and ASP32 in HIV and HTLV-1 protease, respectively, in each monomer act as the main catalytic amino acids [44]. The sequence alignment evidence shows 28% identity between the primary sequence of HTLV-1 protease and HIV protease enzyme, while they are more conserved in the active site (45%) [45].
Ruker et al. [46] Tables 3 and 4 (2), and ALA99B in HTLV protease. Nowadays, the world has been confronted the novel coronavirus (2019-nCOV)-dependent acute respiratory illness. It seems that imbalanced immune cells' activity can lead to the self-destruction of the host, resulting in high severity of COVID-19. Lim et al. [47] showed that a load of the virus decreases significantly in Korean patients diagnosed with COVID-19 when treated with lopinavir and ritonavir as two potent inhibitors of HIV protease. Dipti Mothay et al. [48] analyzed some potent inhibitors of COVID-19 protease that they are able to treat the HIV patients. Based on in silico study, Arg141of protease model of COVID-19 had molecular interactions with lopinavir, ketoamide, and ritonavir; Gln175 showed molecular interactions with nelfinavir and Thr75 and His176 from COVID-19 protease bind to remdesivir via noncovalent interactions. Günther et al. [49] introduced an active site and two allosteric binding sites for potent inhibitors of COVID-19 protease: The hydrophobic pocket is created by Ile213, Leu253, Gln256, Val297, and Cys300 within the C-terminal dimerization domain; the deep groove is located between the catalytic domains (Cys145 and His 141) and the dimerization domain (Arg4, Ser10, Gly11, Glu14, Asn28, Ser139, Phe140, Ser147, Glu290, Arg298) [50]. Günther, et al. showed that AT7519 is the only compound in their screen that bound to loop 107-110, Val202, and Thr292, and the chlorinated benzene ring in AT7519known as a potent compound for human cancer treatment [49] is engaged in various van der Waals interactions to loop 107-110, Val202, and Thr292 [50,51]. Muramatsu et al. introduced substratebinding residues His41, Met49, Gly143, Ser144, His163, His164, Met165, Glu166, Leu167, Asp187, Arg188, Gln189, Thr190, Ala191, and Gln192 for COVID-19 protease [52]. Based on our results, C 16 H 19 O 7 (F14) has a similar binding mode to COVID-19 protease with the AT7519 compound mentioned above. Two compounds are halogenated, which may effect on binding affinity to the special site of the COVID-19 protease. This site may be a target to inhibitor binding for COVID-19 protease. The C 16 H 19 O 6 (F14) has hydrogen bonds and hydrophobic interactions with Arg298 in dimerization domain and loop 107-110. It may introduce as potent inhibitor to COVID-19 protease (Table 5).

Conclusions
First principle calculations based on DFT and molecular docking approaches are used for investigation of optoelectronic and biological properties of pure C 16 (F17), have a considerable contribution on their stability and chemical reactivity. It is shown that fluorination at B category of compounds does not effect on intrinsic static dielectric constant of compounds. So, these compounds have the same nature in the reaction between the drug molecules and the diseases proteases. The anisotropic optical properties