Interaction of Spike protein and lipid membrane of SARS-CoV-2 with Ursodeoxycholic acid, an in-silico analysis

Numerous repositioned drugs have been sought to decrease the severity of SARS-CoV-2 infection. It is known that among its physicochemical properties, Ursodeoxycholic Acid (UDCA) has a reduction in surface tension and cholesterol solubilization, it has also been used to treat cholesterol gallstones and viral hepatitis. In this study, molecular docking was performed with the SARS-CoV-2 Spike protein and UDCA. In order to confirm this interaction, we used Molecular Dynamics (MD) in “SARS-CoV-2 Spike protein-UDCA”. Using another system, we also simulated MD with six UDCA residues around the Spike protein at random, naming this “SARS-CoV-2 Spike protein-6UDCA”. Finally, we evaluated the possible interaction between UDCA and different types of membranes, considering the possible membrane conformation of SARS-CoV-2, this was named “SARS-CoV-2 membrane-UDCA”. In the “SARS-CoV-2 Spike protein-UDCA”, we found that UDCA exhibits affinity towards the central region of the Spike protein structure of − 386.35 kcal/mol, in a region with 3 alpha helices, which comprises residues from K986 to C1032 of each monomer. MD confirmed that UDCA remains attached and occasionally forms hydrogen bonds with residues R995 and T998. In the presence of UDCA, we observed that the distances between residues atoms OG1 and CG2 of T998 in the monomers A, B, and C in the prefusion state do not change and remain at 5.93 ± 0.62 and 7.78 ± 0.51 Å, respectively, compared to the post-fusion state. Next, in “SARS-CoV-2 Spike protein-6UDCA”, the three UDCA showed affinity towards different regions of the Spike protein, but only one of them remained bound to the region between the region's heptad repeat 1 and heptad repeat 2 (HR1 and HR2) for 375 ps of the trajectory. The RMSD of monomer C was the smallest of the three monomers with a value of 2.89 ± 0.32, likewise, the smallest RMSF was also of the monomer C (2.25 ± 056). In addition, in the simulation of “SARS-CoV-2 membrane-UDCA”, UDCA had a higher affinity toward the virion-like membrane; where three of the four residues remained attached once they were close (5 Å, to the centre of mass) to the membrane by 30 ns. However, only one of them remained attached to the plasma-like membrane and this was in a cluster of cholesterol molecules. We have shown that UDCA interacts in two distinct regions of Spike protein sequences. In addition, UDCA tends to stay bound to the membrane, which could potentially reduce the internalization of SARS-CoV-2 in the host cell.

www.nature.com/scientificreports/ Then, density balancing was done over 750,000 steps with a time step of 0.002, 8.0 Å nonbonded cut-off, 1.0 pressure relaxation time, Langevin thermostat for temperature control, a collision frequency of 2.0/ps (gammaln = 2.0) and positional constraints. Finally, a molecular dynamic of 4,000,000 steps (total 10 ns) was executed with a time step of 0.002 ps, a cut-off radius of 8.0 Å, Langevin thermostat for temperature control (ntt = 3) at 300 K, a collision frequency of 2.0 (gamma_ln = 2.0) and positional constraints. For the "SARS-CoV-2 Spike protein-6UDCA", we surround the trimeric structure of the Spike protein with six UDCA residues through MD. The same parameters were used for minimization, state solvation, heating, density balance and production of molecular dynamics. All sugar residues attached to Spike protein were removed. The trajectory analysis was visualized with VMD 1.9.3 37 and the distance and hydrogen bond calculations were carried out with CPPTRAJ 39 .
Minimizations of 50,000 steps were performed at constant volume, no pressure scale, no restriction on hydrogen atoms and a nonbonded cut-off of 10.0 Å. After this, a preheating from 0 to 100 °K, 250,000 steps, with a 0.002-time step, hydrogen bond length restrictions, no evaluation for hydrogen bonds, relative geometric tolerance for coordinates of 0.0000001, Langevin thermostat and a collision rate of 1.0/ps. Then, we simulated a warm-up from 100 to 300 °K with the same parameters as in the previous preheat but for 250,000 steps, periodic boundary conditions at constant pressure, anisotropic pressure coupling, and a pressure relaxation time of 2.0 ps, all according to the AMBER14 instructions 43 . For the density balance, periodic limit conditions were used at constant pressure, Langevin's thermostat, a collision frequency of 1.0/ps, relative geometric tolerance for coordinates of 0.0000001, a cut-off radius of 10.0 Å during 1, 750,000 steps with a time step of 0.002 (for a total of 3.5 ps). The molecular dynamics was done during 60,000,000 steps with a time step of 0.002 (total 120 ps), with restrictions for hydrogen bonds, relative geometric tolerance for coordinates of 0.0000001, Langevin thermostat, with a collision frequency of 1.0/ps at temperature 300 °K, 10 Å cut-off, anisotropic pressure coupling, and periodic boundary conditions at constant pressure.
We used VMD 1.9.3 36 for the visualization of the trajectories and CPPTRAJ 39 for the calculation of distances, RMSD, and the calculation of hydrogen bonds.

Results
Interaction between the trimeric structure of the Spike protein and UDCA by docking. In the "SARS-CoV-2 Spike protein-UDCA" analysis, we employed the trimeric structure of the Spike protein and UDCA using the HEX 6.3 program. UDCA showed an affinity energy of − 386.35 kcal/mol with the residues conforming the three alpha central helices region. This set of three alpha helices is formed by the contribution of one alpha helix of every monomer (residues K986 to C1032). We think that this set of three alpha helices hold together the three monomers of the Spike protein before a post-fusion conformation.
When applying MD to this structure (spike protein-UDCA), we observed that the UDCA remained in the same position during the whole trajectory, and there were even moments when hydrogen bonds were formed with residues close to it ( Table 2). UDCA had affinity mainly with residues R995 and T998 of monomers A, B, and C. Residues R995 and T998 were found in the central helices from subunit S2 (S2) of the protein (Fig. 2).
The oxygen atom OG1 of T998 of the monomer A interacts with the hydrogen HOA1 of UDCA. The hydrogen atom HGB of T998 of monomer B interacts with the oxygen atom O4 of the UDCA, and the oxygen atom OG1 of the residue T998 of monomer C interacts with the hydrogen H222 of UDCA with average distances of 2.67 ± 0.41 Å. On the other hand, the distances between monomers A, B, and C in the pre-fusion state were 5.6 and 8.6 Å for OG1 and CG2 of every T998, respectively, and in the post-fusion state, it is 9.58 and 9.23 Å for the same atoms (Fig. 3A,B). In presence of UDCA, these distances are not modified and remained at 5.93 ± 0.62 and 7.78 ± 0.51 Å, respectively ( Fig. 2A). This leads us to consider that UDCA could act as a connecting bridge that would not allow an adequate conformational shape for the post-fusion state. The RMSD, RMSF as well as Rg of monomers A, B, and C were maintained stable as shown in Fig. 2B-E, indicating an equilibrated system. The monomer C showed the smallest RMSD compared with monomer A and B (RMSD (A), 4.03 ± 0.6 Å; RMSD (B), 3.14 ± 0.42 Å and RMSD (C), 2.89 ± 0.32 Å). RMSF average of monomer C was 2.25 ± 0.56 and the average Rg of monomer C was 44.31 ± 0.28.
Likewise, the distances between UDCA and residues R995, T998 and D994, R995 and T998 of the second and third monomers were preserved. In Table 2, we observed that ARG 995, THR 998 and GLY 999 atoms showed a Table 2. Interactions between UDCA and T998 of central helices of S2 of Spike protein of SARS-CoV-2. These distances and angles correspond to H bonds formed during the trajectory of molecular dynamics. As observed in Fig. 2B, once these atoms are at a close distance, they remain together along the trajectory. In molecular dynamics, most of the H-bonds are at the beginning of the trajectory. H-bonds are formed between UDCA (denoted as IU5_2906 according to AMBER files) and T998, R995, D994 and G999 from three monomers (H-bonds obtained with CPPTRAJ).  36 . We observed that in the closed conformation (7DF3.pdb) (Fig. 3A), the OG1 and CG2 atoms of every T998 are separated for 5.6 and 8.6 Å from each other,  www.nature.com/scientificreports/ respectively. However, for the open conformation (7BNO.pdb) (Fig. 3B), those same atoms had a distance of 9.58 and 9.23 Å, respectively. This suggests that in post-fusion conformation, these alfa-helices are slightly separated, and this separation is important for correcting fusion between viral and host membranes; although, in the presence of UDCA, we observed that the conformation remains closed, as seen in Fig. 2A. In accordance with the observed distances, we think that if UDCA binds to this region and serves as a bridge between these residues, it could prevent the Spike protein from folding correctly to its open conformation, and thus not carry out membrane fusion, Fig. 2A.
UDCA showed affinity towards a region between HR1 and HR2. In the "SARS-CoV-2 Spike protein-6UDCA", we placed six UDCAs around the trimeric structure of the Spike protein in its closed conformation (6VSB.pdb). We visually analyse the trajectory using VMD 1.9.3. We found that three of the six residues were attracted to three different regions of the Spike protein in different areas to the first analysis. One of them was attracted to an ASN536 residue of the C monomer for a very short period of 50 ps, while the rest of the trajectory was kept repelled. Another UDCA was attracted to a residue of THR of chain C with a very short period of 25 ps. The third one was attracted by CYS1126 of chain C with a longer period. The last interaction UDCA-CYS1126 is conserved in different periods of the trajectory. Also, these residues interact for long periods of time, for up to 375 ps, which suggests a strong attraction. The distances between UDCA and N536 was 10.94 ± 4.69 Å, between UDCA and T602 was 10.94 ± 4.69 Å, and between UDCA and C1126 was 6.99 ± 3.15 Å. In monomer C, we found 1.39 of RMSF with CYS1126, less than the average (1.81 ± 0.55) of the whole monomer C; monomer A with 4.03 ± 0.6, monomer B with 3.14 ± 0.42 and monomer C with 2.89 ± 0.32. The RMFS's and Rg of each monomer were similar (Fig. 4A-D). A graphical representation is shown in Fig. 4E, where we see the UDCA interacting with the region between HR1 and HR2.
Lipid membranes and UDCA. In "SARS-CoV-2 membrane-UDCA", we used dioleoyl-or palmitoyl-as controls and observed that the electron density remained stable on both sides of the membrane, in a range of ± 20 Å along all trajectories (Figs. 5A and 6), indicating balanced simulations. Three UDCAs were attracted to the virion-like membrane, which was retained for a long time (30 ns), as was observed in the distances of the centres of mass of the UDCAs with respect to the PC and PE. UDCA bound to these residues from the beginning of the trajectory and remained together along the trajectory (Fig. 5B). In the Supplementary files, we show the simulation of a membrane virion-like with UDCA.
In the ERGIC-like membrane, only one leaf contained a cholesterol residue, the other side contained 62.5% PC, 28.2% PE and 8% PS. Possibly because of cholesterol, the UDCAs that were attracted to the membrane only remained attached for a short time and then were detached, as shown in Fig. 5C,D.
In a eukaryotic PLASM-like membrane, there was 37,5% CHL, 33% PC, 17% PE, and 12,5% PS. We found that two UDCA residues maintained a close distance to the phospholipids, Fig. 5E,F. However, they did not stay together along the membrane. The behaviour of PC and PE is the same, they tend to attract UDCA residues. It is notorious that in the PLASM-like membrane, a cluster of three cholesterols is formed (charmm-gui.org server aleatory constructed that way), and only this way, cholesterol attracted UDCA, since if cholesterol is found alone, it does not attract UDCA.
It is noteworthy that UDCA does not show attraction to membranes containing palmitoyl (phospholipids with 20 carbons), whatever its constitution.

Discussion
UDCA interacts in the central helix region of S2 and HR1 and HR2. We observed that UDCA interacts in two regions of Spike protein. One of the regions which shown attraction is in the area with three alpha-helices of three monomers of the Spike protein (residues R995, T998, V991, D994, T998, T998, and R995), and the other, with a CYS1126 residue in a zone between the HR1 and HR2 regions.
Also, we compared the pre-fusion and post-fusion conformation of the alpha-helices of the Spike protein 34,36,43 and observed that in the pre-fusion conformation (7DF3.pdb, Fig. 3A) and the post-fusion conformation (7BNO. pdb, Fig. 3B), those same atoms have a distance of 9.58 and 9.23 Å, respectively. As was measured before 44 in post-fusion conformation, these alfa-helices are slightly separated; this separation is important for a correct fusion between viral and host membranes, but in the presence of UDCA, we observed that the conformation remains closed.
This leads us to consider that UDCA could join in two regions of the Spike protein, in the S2 region, at residue T998 found in the central helices formed by residues K986-C1032. Precisely in this region, there are two mutations of interest, the K986 and V987, these are already being induced with proline substitutions in several vaccines that are in development 45 . In the joint of the central helices of S2, UDCA could prevent conformational changes in the Spike protein once it binds to ACE2. UDCA could act as a link between the three monomers, since as we observed in Fig. 2A,B, UDCA forms hydrogen bonds with the T998 residues of each of the monomers.

UDCA interacts in the central helix region of S2 and HR1 and HR2.
A key step in the coronavirus infection process is membrane fusion 46 . In this stage, the virus protein fusion core is formed by the union of three copies of HR1 and three of HR2. This stage had been verified in a certain way by the inhibition of fusion through the introduction of peptides, analogous to HR1 and HR2. According to the current model, HR1 and HR2 are joined to form the coiled-coil bundle necessary for fusion 46,47 .
The post-fusion stage takes place after the receptor-binding domain (RBD) has contact with its target protein (ACE2). After this contact, the Spike protein has a series of conformational changes, among these, the union of the three HR1 and three HR2 48 . Another region where the UDCA has a high probability of binding is an   www.nature.com/scientificreports/ intermediate region located between the monomers HR1 and HR2, since it is more exposed, as seen in Fig. 4C. In this region, UDCA interacts with the residue C1126. This residue is part of the C1082-C1126 disulphide bridge, which could modulate the flexibility of the protein 49,50 .
The residue C1126 is in the region between this HR1 and HR2 and in our results of "SARS-CoV-2 protein-6UDCA", the oxygen atoms O4 and O4A of UDCA stays bound to the hydrogen atoms HB2 and HB3 of C1126, respectively, up to 375 ps. UDCA has a β-hydroxylation that provides it with greater hydrophilicity and thus could interact with hydrogens 51,52 . This causes that UDCA and cysteine C1126 to suddenly separate, however, nanoseconds later, both molecules tend to bound, as shown in Fig. 4A (red and magenta line).

UDCAs interact in a stable time interval with the virion-like membrane and PLASM-like membrane.
In membranes, UDCA interacts with PE and PC of virus-like and ERGIC-like membranes. Although with the virion-like membranes the interaction show more affinity, due to three of four residues are attracted to it, whereas in ERGIC-like membrane, only show that one UDCA is attracted.
The finding that three cholesterol molecules are present and UDCA is easily attracted, is related to published studies stating that eukaryotic membranes are rich in cholesterol, which is used by SARS-CoV-2 to enter cells, as it depends on a cholesterol-rich lipid raft 53 . Previously, another study found that cholesterol increases interactions with viral membranes due to its electrostatic and solvation properties 54 , moreover, cholesterol is also known to promote oligomerization of the SARS-CoV-2 FP 55 . Therefore, we think that if UDCA binds to cholesterol lipid rafts in eukaryotic cells, UDCA could hinder the interaction between cholesterol and FP and avoid oligomerization. According to Fig. 5A,B, UDCA may form a cluster with cholesterol-rich lipid rafts in eukaryotic membranes, and thus prevent that the fusion peptide of the Spike protein of SARS-CoV-2 from oligomerizes.
In Fig. 5A,B, we observed that in virion-like membranes, the distance of the centres of mass between UDCA and PE, at the beginning of the simulation are far apart, although they tend to attract each other and once they are joined, they no longer separate. Whereas ERGIC-like membrane, the distance from centres of mass of UDCA towards PE, shows that they are not attracted, except when three cholesterols are together, the UDCA tends to be attracted (PLASM-like membrane, Fig. 5A-F).
Likewise, bile acids are endogenous inhibitors of the NLRP3 inflammasome via NLRP3 ubiquitination, through the TGR5-cAMP-PKA axis. Bile acids could act through this pathway and regulate the activation of the NLRP3 inflammasome 56 and attenuate COVID-19 57 . Furthermore, it has been reported that UDCA presents anti-apoptotic activity by modulating mitochondrial membrane permeability transition 58 .
However, within the limitations of this work, we were unable to build a more realistic membrane due to the lack of force fields of various lipids. In addition, the sugars that accompany the Spike protein could serve as protection against various compounds including UDCA. Therefore, removing them for simulation reduces reliability.
In conclusion, the interaction between UDCA and various residues of the Spike protein and its membrane indicates that UDCA has an affinity and remains attached to the central helix of S2 of the Spike protein monomers. UDCA interacts in a zone between HR1 and HR2 regions, it also tends to remain membrane-bound in the presence of cholesterol molecules. The data showed that UDCA could interact with the SARS-CoV-2 Spike protein and a membrane model like SARS-CoV-2, destabilizing the interaction of the virus its target cells.
Received: 4 July 2021; Accepted: 1 November 2021 Figure 6. Area per lipid for PLASM-IU5-like membrane. As shown in the graph, it tends to be stable along the trajectory. This parameter was similar for all of the rest of the membranes simulated. Obtained with CCPTRAJ.