The interactions of folate with the enzyme furin: a computational study

Entrance of coronavirus into cells happens through the spike proteins on the virus surface, for which the spike protein should be cleaved into S1 and S2 domains. This cleavage is mediated by furin, a member of the proprotein convertases family, which can specifically cleave Arg-X-X-Arg↓ sites of the substrates. Here, folate (folic acid), a water-soluble B vitamin, is introduced for the inhibition of furin activity. Therefore, molecular insight into the prevention of furin activity in the presence of folic acid derivatives is presented. To this aim, molecular docking, molecular dynamics (MD) simulations, and binding free energy calculations were performed to clarify the inhibitory mechanism of these compounds. In this regard, molecular docking studies were conducted to probe the furin binding sites of folic acid derivatives. The MD simulation results indicated that these drugs can efficiently bind to the furin active site. While the folic acid molecule tended to be positioned slightly towards the Glu271, Tyr313, Ala532, Gln488, and Asp530 amino acids of furin at short and long ranges, the folinic acid molecule interacted with Glu271, Ser311, Arg490, Gln488, and Lys499 amino acids. Consequently, binding free energy calculations illustrated that folic acid (−27.90 kcal mol−1) has better binding in comparison with folinic acid (−12.84 kcal mol−1).


Introduction
Coronaviruses, a family of Coronaviridae, can cause signicant human pathologies such as respiratory tract infections in humans and other mammals. 1 Coronavirus infections are usually mild, but some beta coronaviruses including Middle East respiratory syndrome coronavirus (MERS-CoV) and severe acute respiratory syndrome coronavirus (SARS-CoV) may induce critical symptoms. 2,3 In December 2019, an outbreak of lower respiratory tract infections was reported in Wuhan, China. 4 The pathogen was recognized as a novel RNA beta coronavirus, later named as SARS-CoV-2. 5 The infection caused by this virus, COVID-19, has been declared by the World Health Organization (WHO) as a pandemic. 6 In view of the novelty of SARS-CoV-2, further studies are required to obtain more insights about its pathogenesis.
The coronavirus (CoV) genome encodes four structural proteins, comprising spike (S), membrane (M), envelope (E), and nucleocapsid (N). The spike (S) protein of coronaviruses mediates receptor binding and fusion of the virus with the target cells. 7 Each class of coronavirus attaches to a specic cellular receptor to facilitate virus entrance into cells. The angiotensin-converting enzyme 2 (ACE2) and CD209L are shown responsible for SARS-CoV entrance. 8,9 It is reported that SARS-CoV-2 enters the respiratory tract by interacting with the ACE2 receptor. 10 The spike protein comprises an amino (N)-terminal S1 subunit and a carboxyl (C)-terminal S2 subunit. The entrance of the virus is facilitated by cleavage of S protein to S1/S2 subunits. The S1 subunit binds to the ACE2 receptor, while the S2 site interacts with the cell membrane to mediate receptordependent endocytosis, 11 as shown in Fig. 1a. The coronavirus spike protein is cleaved into S1 (receptor binding subunit) and S2 (membrane fusion subunit) by a proteolytic activation at the furin consensus motif RRRR 537 YS (R ¼ arginine, Y: cleavage site) in virus-infected cells. Additionally, the S2 subunit of the S protein is further cleaved at the second furin site (RRRR 690 YS) in the infected cells expressing S constructs. [12][13][14][15][16] Mutations of one basic residue in the RRRR 690 YS motif and use of non-furin cleavable PRRRYS sequence demonstrated that furin may play an important role in furin-dependent entry. 17 The working protease is a cellular proprotein convertase that circulates between plasma membrane, early endosome, and trans-Golgi network (TGN), by participation in endocytic and exocytic paths. 18,19 This proprotein convertase is a major candidate for processing the surface glycoproteins of pathogenic viruses. 20,21 Furin can cleave precursor proteins with specic motifs to produce mature proteins with biological activity. The rst (P1) and fourth (P4) amino acids at the N-terminus of the substrate cleavage site must be arginine "Arg-X-X-Arg Y" (R-X-X-R, X ¼ any amino acid, Y: cleavage site). If the P2 position is basic lysine or arginine, the cleavage efficiency could be improved by about 10 times. 22 The results of a series of analyses have proposed that one of the important reasons for the high infectivity of COVID-19 is a redundant furin cut site in the virus spike protein. 23 Our aim is to suggest folic acid as a potential inexpensive, safe, and non-immunogenic drug candidate for the prevention or treatment of early stages of respiratory disease associated with COVID-19 ( Fig. 1). Folic acid is a type of B vitamin normally found in foods such as spinach, broccoli, asparagus, dried beans, lentils, peas, and oranges. Folic acid helps the body produce and maintain new cells and also prevent changes to DNA that may lead to cancer. Noticeably, folic acid deciency is associated with a variety of human malignancies, including colorectal cancer. The over-expression of folate receptors in the early stages of malignant cell formation can be due to folic acid deciency. Besides, folate malnutrition can cause a high incidence of adenomatous polyps and premalignant lesions of the colon. 24 To this aim, the molecular dynamics (MD) simulations of the interactions of furin enzyme with folic acid and one of its active metabolites, folinic acid, was performed here for the rst time to evaluate the interplay of these molecules with furin.

Molecular dockings
In this study, docking of two folate analogs, folic acid and folinic acid against human furin was performed using molegro virtual docker (MVD) version 6.0 soware. 25 An X-ray crystal structure of human furin was used for docking studies taken from the Protein Data Bank (PDB code 5MIM) accessed at the URL (http://www.rscb.org/pdb) with reasonable resolution (#1.9Å). The optimized coordinates (see Fig. 2) of folic acid and folinic acid were obtained by DFT calculations at the B3LYP/6-311+G** level of theory implemented in Gaussian soware. 26 Two drugs were docked against furin protein and 100 independent runs were performed with the guided differential evolution algorithm. In case of furin-drug complexes, the program generally identied ve different binding sites. Among these ve predicted cavities, the one with the volume of 75.776 A 3 was selected as the potential binding site for investigation. For each docking of protein-ligand, 10 docked poses were generated and listed in Table 1. The Moldock score of the best drug-furin complex of folic acid, and folinic acid were À140.40, and À136.90 kcal mol À1 , respectively.

Molecular dynamics simulations in water
The best congurations from the above-mentioned procedures were selected as the initial structures for MD simulation in aqueous media. In this case, the MD simulations were performed with the GROMACS 4.5.4 program using the GRO-MOS96 53A6 force eld employing periodic boundary conditions in three dimensions. 27 The simulations consisted of  three different systems: (I) furin in water (furin) as control system, (II) furin-folic acid in water (furin/folic acid), and (III) furin-folinic acid in water (furin/folinic acid). The partial atomic charges were calculated by using natural population analysis as implemented in Gaussian 09 program. All simulations were performed in the presence of water by using simple point charge (SPC) model. 28 The net charge of the systems was neutralized by addition of the same amount of sodium or calcium ions. Aer energy minimization, an equilibration NVT run has been carried out over 500 ps, while restraining the position of furin by force constant of 1000 kJ mol À1 nm À2 to their initial position. The systems were simulated for 1 ns using an NPT ensemble. Finally, production runs have been conducted over 500 ns with time steps of 2 Â 10 À3 ps at 310 K. The rst 100 ns of trajectories were set aside for equilibration and were discarded during analysis.

Binding free energy calculations using MM-PBSA
The binding free energies between drug molecules and furin were calculated by using Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) method to comparatively evaluate binding ability of inhibitors to furin using the g_mmpbsa tool of GROMACS. 29 In this regard, binding free energy was evaluated based on the following equations: in which G complex , G protein and G ligand indicate free energies of the docked complex, protein, and drug, respectively. Molecular mechanical (MM) energy, DE MM represents the sum of electrostatic and van der Waals interactions of inhibitors with proteins, independently (DE MM ¼ DE elec + DE vdW ). The solvation free energy, DG solv was identied based on the equation of DG solv ¼ DG pol + DG nonpol , in which DG pol and DG nonpol are the polar and nonpolar solvation free energies, respectively. In agreement with a number of previous computational investigations, 30-32 the contribution of conformational entropy of molecules was ignored in the reported relative DG binding .

Results
Evaluation of the binding modes using molecular docking The ligand-furin docking was performed to predict the major binding sites of folic acid and folinic acid molecules in the active site of furin protein. In general, the interaction energies between protein and ligands are obtained by docking studies. The results showed that folic acid and folinic acid molecules interacted well with the active site residues of furin by formation of hydrogen bonds (H-bond). Different atom sites of the two drug molecules established hydrogen bonding interactions with various amino acids of furin as shown in Fig. 3. Interestingly, the binding sites of folinic acid and furin were clearly different. The interactions of Gly307, Glu271, Tyr313, Gln488, Ala532, Arg490, and Asp530 residues with hydrogen, oxygen, and nitrogen atoms of folic acid were dominant; whereas, strong H-bonds were established through Ser311, Glu271, Arg490, Lys449, and Gln488 residues of folinic acid molecules.   Table 1, the folic acid and folinic acid drugs showed the highest H-bond energy of À14.27 kcal mol À1 and À13.85 kcal mol À1 , respectively. Furthermore, the obtained interaction energies were À161.60, and À159.20 kcal mol À1 for folic acid and folinic acid, respectively.

Elucidation of inhibitors/furin interactions by MD simulations
The MD simulations of folate derivatives interaction with furin were also conducted to gain additional insight into the specic mechanism by which folic acid and folinic acid molecules can exert their potential inhibitory actions in the furin active sites of COVID-19 patients. Atoms labeling for the two drug molecules under study are shown in Fig. 2. The snapshots of the simulated systems aer 400 ns (water molecules were removed for clarity) are depicted in Fig. 4. The position and orientation of drug molecules in the simulations are in line with docking results. The intermolecular interactions between different residues of furin and drug molecules were studied here.
To obtain the inuence of drug molecules on the overall stability of furin structure, a comparative structural evaluation of furin as control system and furin-folic acid and furin/folinic acid as complex systems has been performed. Throughout the MD simulations, both drugs affect the furin secondary structures (see Fig. 5 and Table 2). Furin in all systems is largely dominated by the coil, b-sheet and a-helix structures. In this regard, the average of number of H-bonds between folic acid and furin were obtained to be 6.25 AE 1.20 per time frame, Fig. 4 One mode of drugs binding to the furin protein taken from a snapshot of the simulation at 400 ns. Water molecules were removed for clarity. relative to the value of 2.86 AE 0.96 H-bonds per time frame for folinic acid (see Fig. 6a).
The changes in the structural stability of furin is measured for each system by root-mean-square deviations (RMSD), rootmean-square uctuations (RMSF), to observe uctuations in amino acid residues, and radius of gyration (R g ). The R g was calculated to examine the compactness of protein structure in the presence of inhibitors.
The results are shown in Fig. 6b indicate that R g value is slightly higher for complex systems as compared with control system. For instance, the R g uctuates to the value $2.26 nm for furin/folic acid system, and $2.24 nm for furin system. Thus, ligands induces furin structure to adopt less compact conformations in water in contrast to control system. RMSD of furin in all three systems is shown in Fig. 6c. We see that RMSD of control system converges around $0.26 nm. RMSD uctuates to a lower value for furin-folic acid complex ($0.21 nm) as compared with furin-folinic acid complex ($0.23 nm). In order to evaluate structural exibility of furin, RMSFs of residues are calculated and shown in Fig. 6d. For better clarity, RMSFs of C a atoms of residues in furin were computed for all systems. Generally, the RMSFs of the whole regions are decreased relative to the furin in control system, which shows that the structural exibility of residues are diminished by binding of inhibitors. Thus, the furin-drug complexes are likely to remain more stable than that of free furin.

Distribution functions
To obtain statistically reliable structural data, the radial distribution functions (RDFs) are calculated by averaging over trajectories of long production runs. The RDFs between furin's residues and drug molecules are shown in Fig. 7. For clarity, visual inspection of interaction sites are also demonstrated in Fig. 8 (see Fig. 2 for atom's labels). The RDFs between the center of mass of folic acid and different amino acids of furin indicates that Ala532, Tyr313, and Glu271 interacted with folic acid with a relatively high probability and small dynamics both at short and long distances. In this case, the atom.atom RDFs demonstrated that position of the rst peak obtained for H 19 / O(Glu271) and O 3 /H(Tyr313) were smaller than that of O 4 / H(Ala532). Interestingly, MD simulations showed that the distance of O 3 atom of folic acid with Tyr313 and Gln488 hydrogen atoms were calculated to be 1.66Å, and 1.96Å, respectively (see Fig. 2 for atom labeling and Fig. 8a). The main interaction of O 4 atom of folic acid with Ala532, was located at 1.97Å. The RDF peaks in Fig. 7 (right panel) showed that the a Helix is the sum of a-helix, p-helix and 3 10 -helix. folinic acid tended positioning slightly towards Ser311, Glu271, Arg490, Gln488, and Lys499 at short and long ranges. As shown in Fig. 8b, the main atomic interactions in this case were H 2 / O(Glu271), and N 5 /H(Ser311), which were located at 1.97Å. More details can be obtained from the spatial distribution function of the hydrogen bonding between the furin's amino acids and different atom sites of the drug molecules by the calculation of combined radial/angular distribution function (CDF) as a powerful tool for dening H-bond criteria. 33 Fig. 8c (up panel) manifests the most favorable hydrogen bonding interaction between the polar hydrogen atoms of the Tyr313 and the O 3 atom of carbonyl group of folic acid in the angle range 170 < q < 180 at the distance of around 1.66Å. Particularly, the CDF of Fig. 7c (down panel) indicates that there are interactions between the hydroxyl group of Ser311 and the folinic acid with the angle range 170 < q < 180 and 1.97Å distance. These interactions occurred when the drug molecules tilted substantially to directly interact with furin. Based on these ndings, the folic acid.furin interactions were more probable than folinic acid.

Binding free energies
To get an insight into molecular interactions between inhibitors and furin protein in aqueous media, binding free energy for each protein-ligand complex was evaluated using MM-PBSA method. As shown in Table 3, the binding free energy calculations indicate that the binding is an energetically favorable process for both drug molecules. Binding free energies clearly depicted that folic acid had higher binding energy (DG binding ¼ À27.90 kcal mol À1 ) as compared to folinic acid (DG binding ¼ À12.84 kcal mol À1 ).
From the data reported in Table 3, non-bonded van der Waals, non-bonded electrostatic interactions, and non-polar component to solvation are noted to be favorable for both complexes.
DE elec values are stronger than DE vdW values for both systems. This result highlight that electrostatic interactions between folic acid (folinic acid) and furin dominate over the van der Waals interactions. Thus, intermolecular nonpolar solvation and electrostatic interactions are the main forces involved in the binding of drugs with furin structure. The binding affinity (dissociation constant (K d )) of folate for human folate receptor via isothermal titration calorimetry measurements has been obtained to be $10 Â 10 À12 M. 34 Also, the binding affinity of folate receptor toward folate measured by biolayer interferometry was $1.14 Â 10 À9 M. 35 Folic acid binding to folate receptor determined by saturation radioligand-binding assay has been examined to be $1.90 Â 10 À10 M. 36 The binding free energy of protein-ligand complex formation, can be associated with the dissociation constant through DG exp ¼ RT ln K d C 0 equation, in which T is the absolute temperature, R is the ideal gas constant, C 0 is generally set to 1 with the moles per liter dimension and K d C 0 is accordingly dimensionless. By using this equation, binding free energies of folate for human folate receptor were determined to be À15.05, À12.23, and À13.30 kcal mol À1 , experimentally. Herein, the binding free energy of furin-folate complex predicted by MM-PBSA method is higher than the experimentally determined values of folate receptor-folate complex. It is encouraging that the rank of our predicted binding free energies is in agreement with the one determined by the experimental data.

Discussion
Furin enzyme is associated with a great number of pathologies, including bacterial and viral infections, cancer, and metastasis. Hence, this protein is extremely considered as a drug target. 37,38 A characteristic feature of furin in the protease family is its very limited reactivity toward typical covalent inhibitors due to spatial restrictions. 39 Previous studies have documented that furin could promote the activation of coronavirus by sequence-specic cleavage of the spike protein. Furin cleaves a wide variety of protein precursors in the preferred consensus motif RXR(K)R/R (R ¼ arginine, K ¼ lysine, X¼ any amino acid). 12,40 Therefore, furin protein appears to be a promising target for the infection treatment. The present study identied folic acid as a furin-binding protein by using an all atom MD simulation study. It is well known that folate receptors are mainly expressed in lungs and kidney in normal conditions. 41 Interestingly, the ACE2 receptors are also mostly expressed in lung. Recent studies have proposed that furin inhibition can have a substantial role in the prevention of COVID-19 infection progress. 42 Folic acid is small and stable over a broad range of temperatures and pH values, and it retains its ability to bind to the folate receptor aer conjugation with drugs or diagnostic markers. 43,44 The present study introduces the ability of folic acid to interact and inhibit furin proprotein.
In this study, structural parameters such as radial distribution functions and combined radial-angular distribution functions were used to analyze the intermolecular interactions between folic acid (or folinic acid) and furin protein. The RDF  presents the probability of nding a particle at a certain distance from another reference particle, therefore, contains the information of the average nearest neighbors' distance. 45,46 The rst peaks of simulated RDFs were positioned at very short distances with very high probabilities, which could be attributed to the strong intermolecular interaction between folic acid molecule and furin enzyme. The combined distribution functions conrmed the H-bond character of folic acid-furin interactions. The results indicated that the interactions between folinic acid and furin were high, but substantially lower than that of folic acid. Binding affinities of inhibitors to furin were computed and the results show that binding ability of folic acid is stronger than folinic acid. In this way, folic acid could block the access of COVID-19 spikes to furin and prevent the cell entry and consequently turn-over of the virus. In summary, our results suggest that folic acid could be used to inhibit the furin enzyme. The association of folic acid with furin would affect the structure of the protein and consequently interfere with its proteolytic capability. Thus, folic acid, as a safe drug, could be useful in the prevention or management of COVID-19-associated respiratory disease in the early stages of the disease.

Conclusions
In the present study, the effect of folate derivatives as safe active inhibitors on the structural stability of furin protein was investigated using molecular docking and equilibrated trajectories of MD simulations. The results indicate that binding of folate derivatives leads to the conformational changes of the protein and also affects its internal dynamics. The radial and combined distribution functions reveal that the main interaction between furin and drug molecules is through hydrogen bonding formation. The binding free energy analysis between inhibitors and furin structure with MM-PBSA method inferred that folic acid has better binding affinity as compared to folinic acid. This insights into the underlying inhibitory mechanism of folic acid that show potential inhibitory activity against furin will be benecial for the current and future COVID-19associated respiratory disease therapeutic studies.

Conflicts of interest
There are no conicts to declare.