The folic acid metabolite 7,8-dihydrofolate can decrease the interactions of ACE2 with the Wild type, but not with the Beta and Delta SARS-CoV-2 variants: A silico study.

There is evidence that high levels of folic acid in human plasma may prevent SARS-CoV-2 infections or reduce its symptoms. However, the mechanisms that enable this inhibition are still unknown. This article proves, through molecular dynamics simulation, that folic acid metabolite 7,8-dihydrofolate (DHF) has high aﬃnity to bind as ligands to the angiotensin-converting enzyme 2 (ACE2), and when this complex has been formed, the Wild type SARS-CoV-2 virus cannot bind to the ACE2 receptor, possibly inhibiting infection to the host cell. In contrast, the Beta and Delta variants of this same virus can join the ACE2 with high aﬃnity even with the presence of DHF. This results lead to the conclusion that DHF may inhibit infection from the Wild type SARS-CoV-2 virus, but not its Beta and Delta variants. These results could explain the almost double increase in severe cases of COVID-19 due to the delta variant in pregnant women compared to the Wild type SARS-CoV-2.


I. INTRODUCTION
At the moment of writing this report, the pandemic caused by the SARS-CoV-2 (COVID- 19) virus, which was originally reported in the city of Wuhan, China at the end of 2019, has infected approximately 310 million and killed 5.5 million people around the world. By the end of 2020, COVID-19 vaccines have been produced and distributed in an attempt to stop this pandemic. However, the reality is that demand far exceeds production and distribution capacity. Furthermore, current vaccines may not be as effective at generating immunity in the future, as new strains emerge. Therefore, it is imperative to research and test other preventive and therapeutic methods that could complement the vaccination efforts.

A. Folic acid and COVID-19
Recent studies support the hypothesis that high folic acid (FA) levels could be a positive factor in diminishing the severity of COVID-19. In [1], the authors found that pregnant women, which usually take FA, seemed to be less affected by the SARS-CoV-2 virus than for the A-H1N1 influenza variant. For this study, the researchers contrasted COVID-19 hospitalization data from the UK's National Institute of Health [2] against A-H1N1 data from the Department of Public Health of California [3]. The study in [4] classified 162 COVID-19 patients by their intensity of symptoms and different clinical variables that may explain this difference, including folic acid concentration. Table 1 shows their results and suggests that patients with higher concentrations of FA experienced milder symptoms. However, the molecular mechanisms by which FA could reduce the severity or prevent infection of the SARS-CoV-2 virus are still under debate. Moreover, there is a lack of research on how these mechanisms may interact with newer variants of this virus.  The SARS-CoV-2 virus initiates cell infection after binding to the ACE2 receptor [5] through the interaction with its Receptor Binding Domain (RBD). RBD is a structure that emerges from the virus' Spike as a result of the action of furin [6] and it has been the target of various attempts to combat this pandemic [7,8]. The hypothesis tested in this research is that some folate metabolites, acting through an unknown mechanism, reduce the severity of the SARS-CoV-2 infection. In order to gain insight into this mechanism, the 7,8-dihydrofolate (DHF) metabolite was chosen because this is the first metabolite to which FA is reduced after entering the bloodstream. This reduction is catalyzed by the DHF-Reductase enzyme, which is continuously regenerated by thymidylate synthase in the dihydrofolatetetrahydrofolate cycle [9].

C. Mechanisms of action of folic acid on the SARS-CoV-2 virus
Furin is an endoprotease that is found inside the cell, on its surface and in the extracellular space. This may explain its ability to process a large amount of different substrates, such as growth factors and cytokines, hormones, adhesion molecules, collagen, metalloproteinases, clotting factors, receptors, membrane channels, and albumin [10]. SARS-CoV-2, like many other viruses, requires furin to penetrate the host cells [11]. Zahra Sheybani and colleagues [12] demonstrated, using molecular simulation, that folic acid and folinic acid form a complex with furin. However, the binding place of these metabolites is not in the catalytic region of furin [13], therefore they cannot exert a competitive inhibition. This makes sense since if this binding happened in pregnant women taking FA, it could severely affect the furin physiological processes. Talia Serseg and colleagues discovered that FA can inactivate the 3CL pro protease [14], which is essential for the replication of the SARS-CoV-2 virus, and in general of all coronaviruses. The authors reached this conclusion after employing autodock vina software [15] to demonstrate that FA binds with high affinity to this protease. Gonzalez-Perez M. and colleagues analyzed the chemical interactions between FA and SARS-CoV-2 and found that FA can destabilize the spike of this virus [16]. In a work similar to the one reported here, Vipul Kumar and colleagues [17] used autodock vina [15] to determine that both, FA and its metabolites tetrahydrofolic acid and 5-methyltetraydrofolic acid, have high affinity to bind as ligand to ACE2 and therefore, they assumed that these metabolites can prevent the RBD-ACE2 binding. As in the other proposals previously cited in this subsection [12,14], they draw their conclusions based solely on the binding energy of the ligand with the protein and the stability (Root-Mean-Square Deviation) of this binding. Even if the binding energy and stability ligand-protein are high, this does not mean that the ligand modifies the functions of the protein, either by increasing or decreasing its biological activity in a competitive or non-competitive way. For example, glucose has a high affinity to bind as a ligand to hemoglobin; however, although glucose can make some conformational alterations in the structure of hemoglobin [18], it does not play a significant role in the oxygen affinity pattern of diabetics at the concentrations normally found in vivo [19]. For this reason, this research employed autodock vina [15] not only to find the highest affinity poses of FA and some of its metabolites at the ACE2 receptor, but we also used these poses to build an ACE2-metabolite complex to simulate the interactions between RBD and ACE2 in different scenarios using MD. These interactions are further described in the following sections.

II. MATERIALS AND METHODS
The first objective of this research was to employ molecular dynamics simulations (MD) to reproduce and study the bindings that may naturally occur between different SARS-CoV-2 RBD variants and ACE2 in the presence and absence of folic acid metabolites. The simulation consist of building a box, adding the RBD, ACE2 and, depending on the scenario, one of the folic acid metabolites. Then, the box is filled with water and the system is neutralized with electrolytes. The final steps are running the simulation and studying the interactions between the RBD and the ACE2 throughout the simulation time (trajectory). The final objective is to find if the folic acid metabolite prevents the binding of RBD with ACE2.
For MD the GROMACS [20] software version 2018.7 was used on a supercomputer in which each node had an Intel(R) Xeon(R) CPU E5-2680 v4 @ 2.40GHz 14 cores (28 threads) and four GPU NVIDIA(R) Tesla(R) P100-PCIE-12GB(R) cards. Four different experimental scenarios were designed, one acted as control and three more, as interventions.

A. Control scenario
In the control scenario, which allowed to validate the methodology and materials, the RBD should bind to the ACE2 correctly after a certain period of time. That is, the RBD-ACE2 complex would be obtained with binding energies, hydrogen bonds and distances between the two molecules, similar to the results obtained in vivo. However, even for very powerful computers, these simulations become excessively demanding when the molecules are left at random spatial positions before waiting for the correct binding to happen. Therefore, a decision had to be made to place the RBD and ACE2 molecules face-to-face in the precise position where they could join, but separated by a distance that would allow them to form the RBD-ACE2 complex.
In order to obtain the RBD of Wild type SARS-CoV-2 bound to the human ACE2 receptor, the crystal structure 6LZG was downloaded from the Protein Data Bank (PDB) site [21] (This complex was originally obtained by Q. Wang and colleagues by x-ray diffraction with a resolution of 2.5Ångströms). Then, using the Chimera software [22], all other than amino acids, such as water and electrolytes, were eliminated. In practice, the separation that allows the molecules to join was found by trial and error. Experiments by MD were carried out starting at a distance of 14Å and all the way down to 6Å. For the separation of RBD from ACE2, the Pymol(R) software was used. At distances greater than 6Å, the attraction between both molecules is so low that, instead of approaching, they move away from each other. At 6Å (RBD-ACE2-6A), there is a Lennard-Jones Short Range (LJ-SR) potential energy of around -30 KJ/mol, which causes the molecules to always approach each other, until they form a complex with the same properties observed in living systems. This simulation was based on a protocol published by Justin A. Lemkul [23]. This protocol verifies that the crystal proteic structure is correct, and then builds its topology using in this case the CHARM36-julio2020 force field and the TIP3P water model. Notice that during this process GROMACS establishes a neutral PH by default. Then, a single dodecahedron shaped cell is defined, the RBD-ACE2-6A complex is placed inside the cell and then is filled with water, ions are added as neutralizing agents (sodium or chlorine). In order to avoid steric clashes or inappropriate geometry, minimum energy is obtained by MD. The following step is to balance the solvent and the ions around the protein. To perform this, a two phase equilibrium process by MD is employed. The first phase is conducted to a constant number of particles, volume and temperature (NVT). The second phase stabilizes pressure and and is performed on a constant number of particles, pressure and temperature (NPT). Following this procedures, the system is stabilized at 300 • K and 0 Bar. The last step is to run the simulation. Using this protocol, the system was simulated for 100 ns.

B. Intervention scenarios
After validating the control scenario, the experimentation scenarios were performed. These experiments follow a similar procedure than the control scenario. The experiments try different FA metabolites and SARS-CoV-2 variants. However other metabolites such as (6s)-5-methyl-tetrahydrofolate, (6S)-5,6,7,8-Tetrahydrofolate, (6R)-5,10methylenetetrahydrofolate (2-) and FA itself were also tested. Although these metabolites bind with high affinity to the ACE2 receptor, MD simulation found that they do not prevent the binding of RBD with ACE2, and therefore are not reported in this work.

Intervention scenario A
In this scenario, it will be studied whether RBD of Wild type SARS-CoV-2 can bind to ACE2 when the DHF metabolite is bound to ACE2 as a ligand. This scenario uses the crystal structure built for the control scenario, RBD-ACE2-6A. The chimera software [22] was used to remove the RBD from this complex, leaving only the ACE2 with residues from SER 19 to ALA 614. The 3D structure for the DHF was obtained from the PubChem database [24], using ID 135398604. The autodock tools [25] were used to define the docking conditions for the DHF in protein ACE2: For the ligands, its corresponding torsions were established and then saved in pdbqt format. For the protein, only polar hydrogen molecules were added. Then Kolman charges were added and merged with non-polar hydrogens. The resulting protein is stored in the pdbqt format. The search space for docking the ligand was defined as a cube of 47.25 Angströms per side, which almost completely cover the ACE2 (blind mode). Using the autodock vina software [15] and a script bash, a thousand searches for the pose with the smallest union energy (higher affinity) were executed. The smallest energy pose that was able to produce inhibition was -8.3 kcal/mol (Figure 1). With this pose the ACE2-DHF complex was built.
Using the same method as the control scenario, the Wild type RBD was placed in front of the ACE2-DHF complex, but separated by 6Ångströms, and then the simulation ran for 100 nsec. The objective was to study the interaction of RBD with ACE2 when DHF acts as a ligand in ACE2. The objective of this scenario is to bound answer the question: Can RBD approach ACE2 to form a stable RBD-ACE2 complex when ACE2 has DHF as a ligand? and if affirmative, study the interactions between RBD and ACE2.

Intervention scenario B
The only difference between this scenario and scenario A is that the RBD for the Beta variant of the SARS-CoV-2 virus was used, PDB id 7LYN (B.1.351).

Intervention scenario C
This scenario follows the same methods as scenario A, but it employs the DELTA variant of the SARS-CoV-2 virus. Therefore, the PDB id 7ORB was downloaded. Then the RBD was obtained from the complex and it was prepared following the same procedure as in the control scenario. The RBD of the DELTA variant has two mutations, L452R and T478K, because the RBD isolated from PBD 70RB only has the L452R mutation and lacks the T478K mutation, therefore this mutation was added using the Pymol (R) software. The rest of the protocol is the same as in scenario A.
FIG . 1: The DHF in the ACE2. This pose was found using autodock vina. The labels show the aminoacid residues interacting with the metabolite.

III. RESULTS
For the control scenario and scenarios B and C described in the previous section, initially the RBD and the ACE2 move closer to each other until they reach a mean distance of approximately 0.17 nm (1.7Ångströms) and in that position they are maintained throughout the entire simulation trajectory, see Figure 2. To estimate the binding energy between the RBD and the ACE2, the potential LJ-SR energy was calculated using a GROMACS module. If this energy is lower, the union of the RBD with the ACE2 is stronger and more stable. Figure 3 shows that the LJ-SR energy between the RBD and the ACE2 starts high (≈ -30 kJ/mol) and, as time progresses, it decreases until it reaches a mean value of approximately -250 kJ/mol. Also, at the beginning of the simulation, there is only one hydrogen bond between RBD and ACE2 and, as time progresses, the amount of hydrogen bonds rises to an average of 7 at 100 ns, see Figure 4. This leads to the conclusion that for the control scenario and scenarios B and C, the RBD binds to ACE2 with high affinity, and form a strong and very stable complex.
In contrast, in the scenario A, the RBD and ACE2 move further away at 5 ns, and then they move closer. At 75 ns, they definitely separate, the complex is broken and the LJ-SR energy becomes zero. The hydrogen bonds also disappear, which leads to the conclusion that the RBD cannot form a stable complex with the ACE2. Therefore, the Wild type SARS-CoV-2 virus could not stay attached to the ACE2 receptor and consequently could not initiate host cell infection. It can be seen that in the control (Wild type) and B (Beta + DHF) and C (Delta + DHF) cases, the LJ-SR energy is the lowest, which can be interpreted as a binding with higher affinity. In contrast, in the A intervention (Wild + DHF) the binding energy is higher (lower affinity). In the case of control (Wild type) and intervention B (Beta + DHF) and C (Delta + DHF), it can be observed that during the trajectory the greatest number of hydrogen bonds is maintained. In contrast, it can be observed that in the A scenery (Wild + DHF) at 75 ns there are zero hydrogen bonds and this is maintained for the rest of the trajectory.

IV. DISCUSSION
The results from the molecular simulations seem to demonstrate that DHF can inhibit, non-competitively and in a reversible way, the binding of ACE2 with the Wild type SARS-CoV-2 virus RBD. However, since there is no evidence that DHF can establish covalent bonds with ACE2, it could be assumed that the DHF can spend a short time in ACE2 and later take another pose, in such case, it cannot perform an inhibition or even abandon it. Therefore, it can be said that DHF's inhibitory capacity could depend on the relation of it plasma concentration and the amount of ACE2 receptors present, that is dose-response effect.  Figure 5 shows that RBD mutations for the Beta and Delta variants have not increased the affinity of RBD for the ACE2 receptor, when compared to the Wild type SARS-CoV-2. In other words, the mutations present in the RBD of these variants do not seem to improve the binding energy of their RBD with the ACE2 receptor, so the increase in their contagious capacity may not be explained by this mechanism. On the other hand, the mutations present in the rest of the virion of the variants of this virus may contribute to this capacity. However, if we consider that the virus mutations in its RBD allow it to evade the protection offered by high levels of DHF or other folic acid's metabolites in human plasma. This might explain the higher infection rate observed with these new variants in the general population. If folate levels were an important factor in the relative resistance to infections by this virus in healthy people, this might be different in pregnant women, since their immune system is altered to avoid rejection of the fetus [26], and therefore, the prescribed supplements could no longer provide significant protection for the Delta and other variants, as it did for the Wild type SARS-CoV-2 [1]. Nicola Vousden and colleagues [27] have found that in the UK, out of 3,371 pregnant women who have been hospitalized with symptomatic COVID-19, 24% of admitted women infected with the Wild type SARS-CoV-2 had moderate or severe disease, compared to 36% with the Alpha variant and 45% with the Delta variant. Therefore, it could be said that pregnant women are up to twice as likely to develop severe symptoms when acquiring the infection of mutated SARS-CoV-2 which are able to neutralize the protection generated by FA supplements.
However an important fact that needs to be considered is that simulations are simplified abstractions of the real world phenomena. Therefore, the results presented in this paper need to be validated experimentally.

V. CONCLUSIONS
The results reported here seem to show that the DHF folic acid metabolite has a high affinity to bind to the ACE2 receptor, and that when DHF acts as a ligand, the RBD of the Wild type SARS-CoV-2 might not form a stable, low-energy complex with the ACE2 receptor and therefore cannot initiate infection of a host cell. However, the RBD for the Beta and Delta variants can form a high affinity complex with the ACE2 receptor, even when the DHF acts as a ligand. Therefore, this folic acid metabolite cannot prevent infection from these SARS-CoV-2 variants, which could explain the almost double increase in severe cases of COVID-19 due to the delta variant in pregnant women.

VI. ACKNOWLEDGMENTS
This work was partially supported by the Universidad Autónoma de San Luis Potosí. The high-performance computing resources required for this work were provided by the IPICYT's National Supercomputing Center.