The nanomolar affinity of C-phycocyanin from virtual screening of microalgal bioactive as potential ACE2 inhibitor for COVID-19 therapy

2.1 Software and data collection

Molecular docking simulations were carried out using AutoDock VINA software, downloaded from the website of The Scripps Research Institute (https://vina.scripps.edu/). Visualization of the docking results was performed with BIOVIA Discovery Studio Visualizer 2020. The protein’s crystal structure was downloaded from the Protein Data Bank website (https://www.rcsb.org). ChemDraw Ultra 8.0 (PerkinElmer Inc.) was used to draw 2D and 3D structures and optimize ligand geometry. Data collection strategies and structure of compounds from potential microalgae and cyanobacteria as antivirals were carried out through MarinLit database searches (https://pubs.rsc.org/marinlit/) and publications related to active compounds from marine organisms from 1970 to 2020.

2.2 Ligand preparation and target receptors for molecular docking

The 2D structure of the compound was prepared with ChemDraw Ultra 8.0 and then converted to a 3D structure using the Chem 3D program. All ligands were saved in PDB format. The PDB format of the crystal structure of ACE2 receptor in complex with their inhibitors was downloaded from the RCSB Protein Data Bank (https://www.rcsb.org). The complex for ACE2 and its inhibitor, MLN-4760, was available under PDB ID 1R4L. The MLN-4760 is a cell-permeable, bioavailable, highly potent inhibitor of ACE2 (IC₅₀ = 440 pM against soluble human ACE2) that exhibits good selectivity over bovine carboxypeptidase A and porcine ACE. Afterward, the Biovia Discovery Studio 2020 (BDS) was used to separate the protein from the inhibitor. Concomitantly, undesired molecules were also eliminated. Both protein and inhibitor molecules were then saved into two different PDB files (target receptor and control ligand). Finally, the inhibitor molecule file was put aside, whereas the protein file was loaded into Autodock Tools 1.5.6 software for further preparation steps. These steps cover the addition of polar hydrogen and Kollman charges, removing non-polar hydrogen atoms, and converting into PDBQT format following (Afriza et al., 2018; Khayrani et al., 2021). The file conversion into PDBQT format allowed the file to be loaded in Autodock Vina for molecular docking simulation (Huey et al., 2012).

2.3 Preparation of grid and molecular docking parameters

The grid parameter files were prepared with the help of AutoDockTools 1.5.6. As for its input and output, Vina uses a similar PDBQT molecular structure file format used by AutoDock. The pdbqt files were generated interactively and viewed using MGLTools. The lowest binding energy and Root Mean Square Deviation (RMSD) were further analyzed.

2.4 Molecular docking simulations

Positive control docking was performed using the Autodock Vina program to find the coordinates and parameters that are most appropriate on the active sites. The observed docking parameters were amino acid residues and hydrogen bonds, and the RMSD value should not be more than 2 Å. As for target receptors that are not complex with ligands, blind docking was performed to determine the binding site of the ligand with the receptor.

2.5 Virtual screening of potential compounds and molecular interaction analysis

Virtual screening was performed using PaDEL-ADV software (https://www.yapcwsoft.com/dd/padeladv/) following the protocol from (Hardianto et al., 2018). This program reads each ligand file in pdbqt format, previously prepared using AutodockTools (see sections 2.2 and 2.3). Autodock Vina was then used to dock the ligand with the receptor. Individual binding modes were extracted from the output PDBQT file using vina_split. Results for each binding mode were extracted from the log file and saved into the csv file. The log file and all the related PDB and PDBQT files were then compressed into a zip file. The docking results were sorted by their free energy binding of compounds and by the most histogram groups. Afterward, the compound with the lowest binding energy or maximum multitarget receptors was chosen for further analysis of the interaction between the ligand compound and its receptor.

Autodock tools integrated with the Autodock Vina generated grid maps for each ligand atom. The grid boxes were made to include one site at a time and perform docking. Each coordinate used in the grid box is described in Table 1, with 1 Å spacing on the grid box. Autodock Vina produces energy binding for each molecule and ten models with different conformations and energy. The lowest binding energy was considered as the ligand with the maximum binding affinity. Additionally, Visualization from docking analysis was performed using Biovia Discovery Studio by examining the number and type of interactions between the ligand and the target protein and its native ligand that was used as the control.

2.6 Molecular dynamic simulations

The molecular dynamics simulation was carried out on the ACE2 protein system with potential inhibitors from screening results on the ACE2 binding site. This simulation used the GPU with the program pmemd.cuda on AMBER20 (Case et al., 2021). Initial minimization was performed using 1000 steps of steepest descent. Afterward, 2000 steps of conjugate gradient minimization were carried out by applying 5 kcal mol⁻¹ harmonic force resistance. The force resistance was gradually released from all heavy atoms (including side chains) to only atoms in the mainframe with the same minimization protocol. Additionally, 5000 steps of unimpeded conjugate gradient minimization were performed to eliminate spatial collisions.

Afterward, the system was set to 310 K in increments (0–100 K; 100–200 K; 200–310 K for 20 ps each) for 60 ps. Moreover, the system density and pressure were equilibrated for 1000 ps. The resistance force was released gradually during the equilibrium phase. Then the production stage was carried out for 200 ns. The time-step in the production phase was set to 2 fs. The MD simulation trajectories were then analyzed using the cpptraj module on AmberTools 20.

2.7 CPC isolation from Arthrospira platensis

Dry powder of A. platensis was obtained from the Jepara Aquaculture Research Facility (BBPBAP Jepara), Indonesia. A total of 20 mg of A. platensis dry biomass with 5 mL of sodium phosphate buffer (0.1 M, pH 7) used and subjected to repeated maceration extraction (ME) and then vortexed (Vortex V1 plus, BOECO, Germany). The cell debris were removed by centrifugation (Tomy MX-307, Tomy Seiko Co.ltd., Japan) at 3140 × g for 5 min. The supernatant was then pooled (500 mL) and marked as crude extract and the extraction solvent was colorless. The phycocyanin crude extract was then purified following protocol from (Liao et al., 2011) and the purity of resulting CPC was calculated. The CPC was then freeze dried (BIOBASE BK-FD10PT, China) for in vitro experiments.

2.8 Surface plasmon Resonance-Based competition assay

3.1 Molecular docking analysis for ACE2 receptor

In this study, docking analysis between the control ligand and protein ACE2 showed five hydrogen bonds with several amino acids, such as Arg²⁷³, Thr³⁷¹, Tyr⁵¹⁵, Pro³⁴⁶, and Thr³⁷¹ (Table 2). Apart from the main hydrogen bonds, other interactions were also found, two types of electrostatic bonds with Arg⁵¹⁴ and Arg⁵¹⁸, as well as hydrophobic interaction between ligand and amino acids Tyr⁵¹⁰ and His³⁴⁵. From the docking analysis on several bioactive compounds from microalgae and cyanobacteria, and the various interactions that occur between these compounds and the ACE2 receptor, a total of three potential candidate compounds were selected based on the number of interactions that resemble the native ligand one, as well as the number of overall interactions involved. Phycocyanobilin, phycocyanin, and 7-dehydroporiferasterol peroxide have 5, 4, and 3 hydrogen bonds, respectively, with several hydrophobic interactions. Although not many interactions were formed compared to the native ligand, some of these compounds have approximately similar potential as the native ligand, especially for the phycocyanobilin compound with five hydrogen bonds and three hydrophobic interactions. Based on theoretical and experimental correlations between H-bond pairings and their effects on ligand binding affinity, H-bonds could enhance receptor-ligand interactions when both the donor and acceptor have either significantly stronger or significantly weaker H-bonding capabilities than the hydrogen and oxygen atoms in water (Chen et al., 2016). Thus, our virtual screening revealed that the binding affinity of the screened compounds was similar to or higher than that of the control ligand of ACE2 (MLN-4760).

The docking prediction using Autodock VINA software revealed that among the different molecules targeted, phycocyanobilin could have various interactions with the ACE2 receptor. In this receptor, phycocyanin has a docking score of −9.7 kcal mol⁻¹, which is lower than the ACE2′ native ligand with four hydrogen bonds formed in their amino acids’ residues (Asn¹³¹, Arg²⁵⁵, Glu³⁸⁸, and Asp²⁵¹; Fig. 1a). Additionally, a total of seven hydrophobic bonds were also observed in the amino acid residues (Arg²⁵⁵, Phe²⁵⁶, Ala¹³⁵, Phe⁴⁸⁶, His⁴⁸⁷, Tyr⁴⁹², and Tyr⁴⁹⁷; Fig. 1a). On the other hand, phycocyanobilin has a slightly higher docking score (−9.5 kcal mol⁻¹) compared to that of phycocyanin by forming six hydrogen bonds with the residues of amino acids Lys³⁶³, Thr³⁷¹, Arg⁵¹⁸, His³⁴⁵, Glu⁴⁰⁶, and Glu¹⁴⁵ (Fig. 1b). Apart from the hydrogen bonds, phycocyanobilin also forms electrostatic interactions with the ACE2 receptor at the residue Arg²⁷³ as well as hydrophobic interactions phi – sigma (Tyr⁵¹⁰ and Phe²⁷⁴) and phi – alkyl (His³⁷⁸ and Pro³⁴⁶). From docking prediction, it can thus be assumed that the binding affinity of phycocyanin and phycocyanobilin is lower than those of the native ligand. This result suggests that these compounds can potentially be used as ACE2 inhibitors.

Close

Fig. 1

A detailed view of low-energy binding conformations of phycocyanin (a) and phycocyanobilin (b) bound to the active site of ACE2 generated by molecular docking (Autodock VINA software). Image was created using BIOVIA Discovery Studio 2020. The green dashed line indicates the hydrogen bonds with amino acid residues around the ligand, while the pink dashed line indicates hydrophobic interaction between the ligand and the amino acid residues. The orange dashed line indicates electrostatic interaction.

Export to PPT

Phycocyanobilin (PCB, C₃₃H₃₈N₄O₆) and phycocyanin (CPC, C₁₆₅H₁₈₅N₂₀O₃₀) could thus be good candidates for developing anti-COVID drugs, and more tests were conducted on these compounds. PCB is a blue phycobilin commonly found in cyanobacteria and the chloroplasts of red algae, glaucophytes, and some cryptomonads (Hayes, 2012b; Santoyo et al., 2012). This molecule has a molecular weight of 586.68 g mol⁻¹. It is present only in the phycobiliproteins, allophycocyanin (APC), and phycocyanin (CPC), of which it is the terminal acceptor of energy (de Morais et al., 2015). A previous study by (Radibratovic et al., 2016) revealed that PCB could bind to human serum albumin (HSA), a protein found mainly in the blood of humans, increasing its ability to prevent the proteolytic activity of other proteins.

3.2 Molecular dynamics simulation of PCB and CPC on ACE2 receptor

The molecular dynamics simulations of the ACE2 receptor that binds to PCB, CPC, and native ligands showed a significant change in distance (Å) from the initial conformation when compared to the structure of ACE2 Apo (Fig. 2). The dynamics of the receptor that binds to the ligand was more stable than the Apo structure. This result is indicated by the root mean square deviation (RMSD) value, which increased to 3.5 Å after 150 ns of simulation (Fig. 2). In the receptor bound to the control, the RMSD value ranged from 2 to 2.5 Å during the simulation. On the other hand, the receptor bound to the phycocyanin and phycocyanobilin compounds, the RMSD profile showed a very stable value of 1.5 Å. This stable conformation indicates favorable interaction between the ligand and the receptor, thereby changing the flexible nature of the ACE2 structure to become more rigid.

Close

Fig. 2

The RMSD analysis of 200 ns-length MD simulation results for ACE2 receptor and each of its ligands (i.e control, phycocyanin, and phycocyanobilin).

Export to PPT

The RMSD value of the control ligand experienced a relatively high progression compared to the phycocyanin and phycocyanobilin compounds, where the shift occurred up to 15 Å at 150 ns of the simulation (Fig. 3a-c). After this period, the native ligand again dropped its RMSD value to the similar level as the phycocyanin ligand (Fig. 3). This result indicates lower bond stability in the native ligand, which causes the RMSD value to increase to 15 Å.

Close

Fig. 3

The conformational changes of each ligand which is control (a) phycocyanin (b), and phycocyanobilin (c) at 100 ns interval throughout the simulations.

Export to PPT

The ACE2 protein has 597 amino acids. From the root mean square fluctuation (RMSF) profile, the ACE 2 receptor that is not bound to the ligand has a higher fluctuation than the system bound to the ligand. In contrast, as for phycocyanin and phycocyanobilin ligands, RMSF values produced a more stable profile than the native system, suggesting that these two compounds have a potential interest as inhibitors for ACE2 protein (Fig. 4).

Close

Fig. 4

The RMSF analysis of the 597 amino acid residues of ACE2 receptor during the 200 ns simulation.

Export to PPT

The molecular mechanics/generalized Born surface area (MM/GBSA) method evidenced that during the MD simulation, the phycocyanin had the lowest total binding free energy (–32.8158 kJ mol⁻¹) compared to other ligands, such as phycocyanobilin (-17.162 kJ mol⁻¹) and the native ligands (-5.0989 kJ mol⁻¹, Table 3). Further hydrogen bond analysis results are also consistent with these findings, where the average number of hydrogen bonds of phycocyanin was the highest among all, followed by phycocyanobilin and the control, respectively (Fig. 5). It is also indicating that the hydrogen bond contributes to most of the interactions between ligand and ACE2 receptor. Therefore, this finding suggests that the interaction of phycocyanin and phycocyanobilin compounds with the ACE2 receptor is relatively strong, and they have stable binding free energy. To our knowledge, little is known about the potential of phycocyanin and phycocyanobilin as potential SARS-CoV-2 inhibitors by targeting the ACE2 receptor. Recent studies from Al-Khafaji et al. (2020) demonstrated that the molecular docking calculations and molecular dynamic analysis support the fact that the phycocyanin could be a potential candidate inhibitor for the ACE2 receptor towards SARS-CoV-2. On the other hand, Pendyala et al. (2021) demonstrated that phycocyanobilin potentially inhibits the binding of other SARS-CoV-2 receptors, such as the main protease (M_pro) and papain-like protease (PL_pro). The in silico study revealed that the binding affinity score of phycocyanobilin with M_pro and PL_pro were −8.6 and −9.8 kcal mol⁻¹, respectively.

Close

Fig. 5

The average number of hydrogen bonds between each ligand throughout 200 ns simulations. The angle of 135° and 3.0 Å cutoff were applied to the calculations.

Export to PPT

3.3 Surface plasmon Resonance-Based competition assay

The SPR method is a powerful tool for discovering protein–protein interaction (PPI) inhibitors. Moreover, this method can quantitively explain highly potent PPIs and weak fragment or low molecular weight (LMW) ligand–protein interactions. In this study, the function of phycocyanin (CPC) as an inhibitor of SARS-CoV-2 was evaluated to further validate the molecular docking and molecular dynamics studies using the SPR technique. With this method, first, we confirmed that after the dissociation of ACE2 and BSA 1 %, the dynamic response did not return to the baseline (Fig. 6). This result indicates that ACE2 and BSA 1 % well bonded on the surface of the SPR chip so that the measurement of CPC binding is ready to be carried out. This SPR technique can provide information on the affinity and kinetics of CPC binding by ACE2. The ACE2 − CPC binding SPR response data at different concentrations (5–125 μM) are shown in Fig. 7a-b.

Close

Fig. 6

The SPR dynamic responses of ACE2 and BSA1% immobilization processes.

Export to PPT

Close

Fig. 7

(a) The SPR dynamic response resulting from the interaction of ACE2 – CPC at a concentration of 5 – 125 μM, (b) kinetic adsorption analysis on each CPC association response by ACE2 with different concentrations to obtain the k_obs values, and (c) the linear regression of C vs k_obs to obtain the k_on value, and (d) Non-linear regression curve of the BS isotherm adsorption model on the concentration effect to the dynamic response data.

Export to PPT

To determine the efficacy of the developed inhibitor, it is necessary to evaluate the residence time (τ), where the longer the time, the more potent the drug developed (Zhu et al., 2021). The τ was obtained through the value of the dissociation rate (k_off) (Copeland et al., 2006). Fig. 7b shows the SPR curve of the association of CPC with ACE2, which was then analyzed using an adsorption kinetic model. It was demonstrated that the association curve of SPR concentration of 5–125 μM follows the pattern of different adsorption kinetic models. Concentrations of 5 and 10 μM followed the Pseudo second-order model. In contrast, 25 – 125 μM concentrations follow Avrami's kinetic model with relative coefficient values (R²) shown in Table S1. The obtained linear regression parameters are summarized in Table S1.

The linear regression between C and k_obs revealed the linearity relationship $k_{\mathit{obs}}=0.00035145\left[C\right]-0.000701142$ with R² = 0.99214 (Fig. 7c-d). From the linear equation, the computed values k_on and k_off were 0.00035145 s⁻¹ μM⁻¹ and 0.000701142 s⁻¹, respectively. On the other hand, the affinity constant (K_A) and dissociation constant (K_D) for ACE2 – CPC interactions were obtained using K_D =  $\frac{k_{\mathit{off}}}{k_{\mathit{on}}}$ and K_A = $1/K_D$ , thus the value of K_D and K_A were 1.995 μM and 0.501 μM⁻¹, respectively (Swinney et al., 2014). Additionally, the value of τ = 1/k_off obtained was 0.396 h, while the half-life value (t_1/2 = 0.693/k_off) was 0.275 h (Copeland et al., 2006). It appeared that the range of residence time values obtained in this study was relatively narrow. This result is likely caused by the lack of strength of ACE2 in binding CPC. Nevertheless, this argument needs to be validated. Therefore, it is necessary to calculate the binding energy ( $\mathrm{\Delta }$ G⁰ = R.T. Ln K_D) involved during the adsorption process (Karoui et al., 2020). As a hypothesis, if the gas constant (R) is 1.987 cal K⁻¹ mol⁻¹ and the temperature (T) is 298 K, the value of $\mathrm{\Delta }$ G⁰ obtained is − 7.77 kcal mol⁻¹.

To predict the binding mechanism of CPC by ACE2, the data were evaluated using a non-linear regression isotherm adsorption model. The curve pattern follows the Brouers-Sotolongo (BS) model with a relative coefficient value (R²) of 0.9996. The parameters obtained from the C vs ${(\mathrm{\Delta }RU/Ro)}_m$ non − linear regression are ${(\mathrm{\Delta }RU/Ro)}_m$  = 76.087, K_BS of 0.00337, and α of 1.548. In the present study, the value of α was greater than 1, suggesting that the CPC binding process by ACE2 takes place through slow kinetic biosorption, and each ACE2 binds to CPC molecules with unequal energy (Karoui et al., 2020).

3.4 Potential microalgae candidates for developing SARS-CoV-2 inhibitor

In the present study, our docking and molecular dynamics analysis to search for potential SARS-CoV-2 inhibitor from natural products demonstrated that phycocyanobilin, with CPC in particular, could be used as a potential inhibitor, as this compound has an excellent binding affinity and significant number type of interaction, both characteristics that resemble those of its natural ligands. Moreover, the results of these in silico experiments were also confirmed in vitro with the SPR experiment. CPC can be obtained from Arthrospira platensis (previously known as Spirulina platensis), a cyanobacteria species that is extensively cultured (Olaizola, 2003; Singh et al., 2011), mainly as a source of food, feed supplement, and natural blue pigment (Belay et al., 1996; McClane et al., 2006; Singh et al., 2011). It is worth noting that CPC from several Arthrospira species has also therapeutic potential but without any adverse effect to the living organisms (Chamorro et al., 1997, 1996), in contrast to other natural blue pigments and bioactive compounds such as marennine-like pigments from the marine diatoms Haslea ostrearia and H. nusantara (Prasetiya et al., 2017, 2019; Prasetiya et al., 2020a; Prasetiya et al., 2020b).

In this study, CPC demonstrated potential activity as a candidate for the SARS-CoV-2 inhibitor. A previous study from (McCarty, 2007) demonstrated that CPC constitutes up to 14 % of the total dry weight of A. platensis, in which phycocyanobilin represents 4.7 % of the mass of CPC. Regarding the production of phycocyanin from A. platensis biomass, several studies revealed that extraction and purification is relatively simple as it can be completed with four major steps that cover; crude extract preparation, ammonium sulfate precipitation, dialyses, and anion exchange chromatography (Kumar et al., 2014) or through a hexane extraction process combined with high pressure (Seo et al., 2013).

To complete this search of compounds that could be potential candidates of SARS-CoV-2 inhibitors, our approach highlighted that the genus Chlorella also presents a significant interest. According to our docking analysis, several bioactive compounds from the database such as oxocholesterol, ergosterol, dehydroporiferasterol, and triterpenoid are produced by Chlorella vulgaris and they display inhibitory activity against SARS-CoV-2. To our knowledge, most Chlorella species mainly inhabit freshwater and are commonly found in very nutrient-rich waters. However, it is possible to find these microalgae in marine environment since a few marine species are also known (Levine and Fleurence, 2018). Chlorella strains or species are the most cultivated microalgae since they are widely used as nutraceutical and feed supplements, as well as in the pharmaceutical and cosmetics industry (Borowitzka, 1997; Li et al., 2014; Silva et al., 2018; Yasuhara-Bell and Lu, 2010). Additionally, previous studies demonstrated that several species of Chlorella were found to contain ergosterol as their major sterol (de Morais et al., 2015; Rahal et al., 2014; Yasuhara-Bell and Lu, 2010). This molecule also has been reported to have the inhibitory activity to zoonotic influenza viruses, such as H5N1, H7N9, and H9N2 (Silva et al., 2018).