Molecular modelling studies unveil potential binding sites on human serum albumin for selected experimental and in silico COVID-19 drug candidate molecules.

Human serum albumin (HSA) is the most prevalent protein in the blood plasma which binds an array of exogenous compounds. Drug binding to HSA is an important consideration when developing new therapeutic molecules, and it also aids in understanding the underlying mechanisms that govern their pharmacological effects. This study aims to investigate the molecular binding of coronavirus disease 2019 (COVID-19) therapeutic candidate molecules to HSA and to identify their putative binding sites. Binding energies and interacting residues were used to evaluate the molecular interaction. Four drug candidate molecules (β-D-N4-hydroxycytidine, Chloroquine, Disulfiram, and Carmofur) demonstrate weak binding to HSA, with binding energies ranging from -5 to -6.7 kcal/mol. Ivermectin, Hydroxychloroquine, Remdesivir, Arbidol, and other twenty drug molecules with binding energies ranging from -6.9 to -9.5 kcal/mol demonstrated moderate binding to HSA. The strong HSA binding drug candidates consist of fourteen molecules (Saquinavir, Ritonavir, Dihydroergotamine, Daclatasvir, Paritaprevir etc.) with binding energies ranging from -9.7 to -12.1 kcal/mol. All these molecules bind to different HSA subdomains (IA, IB, IIA, IIB, IIIA, and IIIB) through molecular forces such as hydrogen bonds and hydrophobic interactions. Various pharmacokinetic properties (gastrointestinal absorption, blood-brain barrier permeation, P-glycoprotein substrate, and cytochrome P450 inhibitor) of each molecule were determined using SwissADME program. Further, the stability of the HSA-ligand complexes was analyzed through 100 ns molecular dynamics simulations considering various geometric properties. The binding free energy between free HSA and compounds were calculated using Molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) and molecular mechanics generalized Born surface area (MM/GBSA) approach. The findings of this study might be useful in understanding the mechanism of COVID-19 drug candidates binding to serum albumin protein, as well as their pharmacodynamics and pharmacokinetics.

Introduction

In December 2019, Wuhan, Hubei Province, China, reported numerous cases of a new respiratory disease. By January 2020, it had been proven that these illnesses were caused by a new coronavirus known as severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and the disease came to be known as coronavirus disease 2019 (COVID-19) (Ashour et al., 2020). SARS-CoV-2 is a new betacoronavirus that shares 79% of its genome sequence with severe acute respiratory syndrome coronavirus (SARS-CoV) and 50% similarity with Middle East respiratory syndrome coronavirus (MERS-CoV) (Lu et al., 2020). It has a genomic structure that is similar to that of other betacoronaviruses. The six functional open reading frames (ORFs) are replicase (ORF1a/ORF1b), spike (S), envelope (E), membrane (M), and nucleocapsid (N) in order from 5′ to 3′. Between the structural genes, there are additional seven potential ORFs encoding accessory proteins (Chan et al., 2020). The demand for new antiviral drugs against SARS-CoV-2 continues as the COVID-19 pandemic spreads quickly (Villamagna et al., 2020). There is currently insufficient evidence that any antiviral drugs available today may effectively treat COVID-19 pneumonia. However, numerous clinical trials on possible antiviral treatments are now underway (Zhaori et al., 2020). The antiviral treatments can be classified into two groups- the one that targets the coronavirus, either by inhibiting a key viral enzyme involved in genome replication or by preventing its entry into human cells and the second that is intended to influence the human immune system, either by enhancing the innate response against viruses or by suppressing the inflammatory processes that lead to lung damage (Tu et al., 2020). The majority of these drugs were developed for different infections before being repurposed for COVID-19. The method of repurposing existing antiviral agents authorized or under development for other viral diseases has been embraced in light of the necessity to discover an effective treatment for symptomatic individuals (Singh et al., 2020). The SOLIDARITY study, an international clinical trial, was recently initiated by the World Health Organization (WHO) to address this issue. Lopinavir and ritonavir, with interferon β and chloroquine, and remdesivir are among the drugs tested in this study (Uddin et al., 2020). These drugs like others need to bind to the drug carrier protein known as human serum albumin which influences their activity and half-life.

Human serum albumin (HSA) is a single-chain, non-glycosylated polypeptide with a molecular weight of 66,500 Da and 585 amino acids (Yamasaki et al., 2013). The structure of HSA indicates the existence of three domains, namely domains I (residues 1–195), II (196–383), and III (384–585), which are not only topologically identical but also have comparable three-dimensional structures, as anticipated by amino acid sequence comparison (He and Carter, 1992, Yang et al., 2014). The polypeptide chain has an approximate dimension of 80 × 80 × 30 Å and forms a heart-shaped structure (Sugio et al., 1999). Except for Cys34 (in domain I), HSA has 35 cysteine residues, all of which are involved in disulfide bond formation, which helps to stabilize the protein. Interdomain and inter-subdomain interactions have a crucial role in the stability of the HSA molecule, according to crystallographic evidence (Yamasaki et al., 2013). It is the most abundant protein in blood plasma that serves a variety of important physiological functions. HSA controls colloidal osmotic pressure and transports a variety of endogenous substances such as fatty acids (FA), hormones, bile acids, amino acids, metals, and toxic metabolites, among others (Yang et al., 2014). Furthermore, through binding with HSA, a wide range of drugs are transported to their target organs/tissues (Yamasaki et al., 2013). As a result, HSA not only sequesters bound pharmaceuticals from oxidation and impacts drug distribution in vivo, but it also changes drug pharmacokinetic and pharmacodynamic characteristics (Tayyab and Feroz, 2020).

HSA's remarkable ability to bind a wide range of drugs is one of its most distinguishing features (Yang et al., 2014). Given the large quantity of HSA in plasma, drug binding affinity to HSA is an essential consideration when developing novel therapies. Furthermore, the interaction of drugs that bind to HSA at the same time might alter HSA binding behaviour and potentially affect the drugs' therapeutic efficacy (Tesseromatis and Alevizou, 2008). The structures of HSA–ligand complexes have shown not only where distinct drug binding sites on HSA are located, but also reveal how several drugs interact with HSA (Yang et al., 2013, Yang et al., 2012). When studying the processes determining the pharmacological effects of these molecules, understanding drug binding characteristics to HSA is critical (Yang et al., 2014). In this study, the molecular mechanism underlying the interaction between drug candidate molecules for COVID-19 and human serum albumin was investigated. The interaction between the drug candidate molecules and HSA was evaluated in terms of binding energies and binding site identification through molecular docking studies. The pharmacokinetic properties of the COVID-19 drug candidate molecules were also computationally investigated. The dynamic behaviour of the free HSA and HSA-ligand complexes were further explored through molecular dynamics simulations in an aqueous environment.

Materials and methods

Retrieval of drug candidate molecules

A total of 55 chemical structures were downloaded from PubChem database (Kim et al., 2016) which comprises 16 experimental drug candidate molecules for COVID-19 (Set A) (Caly et al., 2020, Choy et al., 2020, Dai et al., 2020, Devaux et al., 2020, Jin et al., 2020, K.Y. Wang et al., 2020, M. Wang et al., 2020, X. Wang et al., 2020), 26 in silico drug candidate molecules for COVID-19 i.e. 15 molecules in Set-B from others’ publications (Beck et al., 2020, Hall and Ji, 2020, Ho, 2020, Ke et al., 2020, Pant et al., 2020), 11 molecules in Set C from our previous publications (Gurung et al., 2020a, Gurung et al., 2020b, Gurung et al., 2020c), and 13 control drug molecules that bind serum albumin proteins (Set D) (Varshney et al., 2010). The compounds that lacked a 3D structure were transformed using OpenBabel version 2.4.1 (O’Boyle et al., 2011) and optimized with the MMFF94 force field (Halgren, 1996).

Retrieval of protein structures

The three-dimensional X-ray crystal structure of human serum albumin was obtained from Protein Data Bank (http://www.rcsb.org/) using PDB ID: 6HSC. The structure consists of human serum albumin in complex with Aristolochic acid at a resolution of 1.9 Å.

Molecular docking studies

The binding of COVID-19 drug candidate molecules to HSA was studied using AutoDock Vina program which is based on a sophisticated gradient optimization method (Trott and Olson, 2010). The serum protein was prepared by removing the heteroatoms such as ions, water molecules and cocrystal ligands and the addition of polar hydrogen atoms and Kollman charges. The ligands were prepared by adding hydrogen atoms, Gasteiger charges and optimally defining the torsions. Blind docking was performed using a grid box centred at XYZ coordinates of −71.2957, 1.2689 and 14.0226 with dimensions of X: 98.2245 Å, Y: 51.0570 Å and Z: 90.0794 Å and exhaustiveness were set to 8. The binding conformations were clustered and ranked based on their binding affinities. The molecular interactions such as hydrogen bonds and hydrophobic interactions between HSA and compounds were evaluated using LigPlot + program version 1.4.5 (Laskowski and Swindells, 2011).

Determination of pharmacokinetic properties

Various pharmacokinetic properties (gastrointestinal absorption, blood-brain barrier permeation, P-glycoprotein substrate, and cytochrome P450 inhibitor) of each molecule were determined using SwissADME tool (Daina et al., 2017).

Molecular dynamics simulation

The trajectories of the free HSA and HSA-ligand complexes were studied through 100-ns of MD simulations using GROningen MAchine for Chemical Simulations (GROMACS) 2019.2 software (Hess et al., 2008) with GROMOS96 43a1 force field. The free HSA and the HSA-ligand complexes were prepared for MD simulation by considering a cubic box of 1 Å spacing and solvated with simple point charge (SPC216) waters. A leap-frog time integration algorithm was used for integrating Newton's equations of motion. The systems were neutralized by adding an appropriate number of counterions and subsequently energy minimized. The systems were heated to 300 K and equilibrated for 100 ps in NVT (Number of particles, Volume and Temperature) ensemble and another 100 ps in NPT ensemble (Number of particles, Pressure and Temperature). After heating and equilibration, the systems were subjected to production MD run for 100 ns in NPT ensemble. PRODRG web server (Schüttelkopf and Van Aalten, 2004) was used to generate topology of the ligand. MD analysis was performed by choosing parameters such as RMSD (root mean square deviation), RMSF (root mean square fluctuation), Rg (radius of gyration), total SASA (solvent accessible surface area) and the number of hydrogen bonds (NHBs). Graphs were plotted using Xmgrace tool.

Binding free energy analysis

The binding free energy (ΔGbind) of the compounds was calculated using LARMD program (Yang et al., 2020) which uses the following equation (Equation 1).where ΔEbind is the binding energy, TΔSsol is the solvation entropy and TΔSconf is the conformational entropy. While entropy was computed using an empirical method (Hao et al., 2009, Pan et al., 2008), the enthalpy was derived using the Molecular mechanics Poisson–Boltzmann surface area (MM/PBSA) or molecular mechanics generalized Born surface area (MM/GBSA) method (Hou et al., 2011).

Results

A total of 55 molecules comprising of experimental (N = 16, Set A1-A16) and in silico COVID-19 drug candidate molecules (N = 26) which include set B17-B31 (N = 15) and set C32-C42 (N = 11) along with the control drugs (known binders of HSA) (N = 13, set D43-D55) were docked into the human serum albumin protein using the blind docking method (Fig. 1). The molecular docking results of the compounds are represented in Table 1. On comparing the binding energy scores with the control data set (−7 to −9.3 kcal/mol), the selected COVID-19 drug candidate molecules were classified into three different categories- weak, moderate and strong binding molecules. Four drug candidate molecules viz., Chloroquine (A9), β-D-N4-hydroxycytidine (A10), Disulfiram (A11), and Carmofur (A13) demonstrate weak binding to HSA, with binding energies ranging from −5 to −6.7 kcal/mol, according to molecular docking studies. Ivermectin (A1), Hydroxychloroquine (A2), Arbidol (A3), Ebselen (A4), Remdesivir (A5), Lopinavir (A6), Emetine (A7), Homoharringtonine (A8), Shikonin (A14), Compound 11a (A15), Zanamivir (B17), Atazanavir (B20), Bedaquiline (B21), Brequinar (B22), Clofazimine (B24), Gemcitabine (B26), Tolcapone (B27), Vismodegib (B28), Cobicistat (B29), Darunavir (B31), Bonducellpin D (C34), Vincapsusine (C38), ZINC000254565785 (C41) and ZINC000726422572 (C42) with binding energies ranging from −6.9 to −9.5 kcal/mol demonstrated moderate binding to HSA. The strong HSA binding drug candidates consist of fourteen molecules Tideglusib (A12), Compound 11b (A16), Indinavir (B18), Saquinavir (B19), Celecoxib (B23), Conivaptan (B25), Ritonavir (B30), Ergotamine (C32), Dihydroergotamine (C33), Glecaprevir (C35), Daclatasvir (C36), Paritaprevir (C37), Alloyohimbine (C39) and Gummadiol (C40) with binding energies ranging from −9.7 to −12.1 kcal/mol. All these molecules bind to different HSA subdomains (IA, IB, IIA, IIB, IIIA, and IIIB) (Table 1) through molecular forces such as hydrogen bonds and hydrophobic interactions. The Carmofur (A13) binds to HSA with a binding energy of −6.7 kcal/mol and the HSA_A13 complex is stabilized by a hydrogen bond with Ser480 and hydrophobic interactions through Gly328, Lys212, Ala213, Leu327, Ala210, Val482, Leu347, Lys351, Leu331, Ala350, Val216 (Fig. 2A). The Bonducellpin D (C34) binds to HSA with a binding energy of −7.9 kcal/mol and the HSA_C34 complex interaction occurs through four hydrogen bonds with Glu153, Ser192, Gln196 and His288 and hydrophobic interactions through Lys199, Lys195, Ala291, Phe157, Glu292, Glu188, Asp451 (Fig. 2B). The compound 11b (A16) binds to HSA with a binding energy of −10.4 kcal/mol and the HSA_A16 complex is stabilized through two hydrogen bonds with Trp214 and Asp451 and hydrophobic interactions through Lys195, Lys199, Val482, Ala210, Ser480, Ser202, Leu481, Ser454, Leu457, Leu453, Glu450, Arg485, Val344, Leu198 (Fig. 2C). The control Indomethacin (D43) binds to HSA with a binding energy of −8.6 kcal/mol and the HSA_D43 interaction is facilitated through one hydrogen bond with Ser287 and hydrophobic interactions through Leu260, Ile290, Ala291, Lys199, Arg222, Val241, Trp214, Arg257, His242, Leu238 and Ala261 (Fig. 2D). The pharmacokinetic properties such as gastrointestinal absorption, blood-brain barrier permeation, P-glycoprotein substrate, and cytochrome P450 inhibitor which play a significant role in absorption, distribution, metabolism, excretion, and toxicity (ADMET) were calculated for each molecule (Table 2).

Fig. 1

The chemical structures of selected experimental and in silico COVID-19 drug candidate molecules used in the study.

Table 1

Binding energy and binding sites determined using molecular docking studies.

Molecules	Name	PubChem ID	Binding Energy (kcal/mol)	Binding strength	Binding region
A1	Ivermectin	6321424	−9.4	Moderate	IB-IIA-IIIB
A2	Hydroxychloroquine	3652	−7.1	Moderate	IIIA
A3	Arbidol	131411	−7.3	Moderate	IB-IIA
A4	Ebselen	3194	−8.5	Moderate	IIA-IIB-IIIA
A5	Remdesivir	121304016	−8.9	Moderate	IIA-IIIA
A6	Lopinavir	92727	−9.1	Moderate	IB-IIA-IIIA
A7	Emetine	10219	−8.1	Moderate	IB-IIA-IIIA
A8	Homoharringtonine	285033	−6.9	Moderate	IIA-IIB
A9	Chloroquine	2719	−6.5	Weak	IIA-IIB
A10	β-D-N4-hydroxycytidine (NHC)	197020	−6.5	Weak	IA-IIA
A11	Disulfiram	3117	−5	Weak	IIIA
A12	Tideglusib	11313622	−10	Strong	IB-IIA-IIIB-IIIA
A13	Carmofur	2577	−6.7	Weak	IIA-IIIB-IIIA
A14	Shikonin	479503	−8.9	Moderate	IIA-IIB
A15	Compound 11a	–	−9.3	Moderate	IB-IIA
A16	Compound 11b	–	−10.4	Strong	IIA-IIB-IIIA
B17	Zanamivir	60855	−7.8	Moderate	IIIA
B18	Indinavir	5362440	−9.7	Strong	IB-IIA-IIIA
B19	Saquinavir	441243	−10.2	Strong	IB-IIA
B20	Atazanavir	148192	−8.7	Moderate	IB-IIA-IIIA
B21	Bedaquiline	5388906	−8.7	Moderate	IB-IIA
B22	Brequinar	57030	−9	Moderate	IB-IIIA
B23	Celecoxib	2662	−9.8	Strong	IIA-IIIA
B24	Clofazimine	2794	−9.5	Moderate	IB-IIIA
B25	Conivaptan	151171	−11.2	Strong	IIA-IIB-IIIA
B26	Gemcitabine	60750	−7.3	Moderate	IIIA
B27	Tolcapone	4659569	−9.3	Moderate	IB
B28	Vismodegib	24776445	−8.5	Moderate	IB-IIIA
B29	Cobicistat	25151504	−9.3	Moderate	IIA-IIIA
B30	Ritonavir	392622	−10.1	Strong	IIA-IIB-IIIA
B31	Darunavir	213039	−7.9	Moderate	IB-IIIA
C32	Ergotamine	8223	−11.5	Strong	IIA-IIB-IIIA
C33	Dihydroergotamine	10531	−11.9	Strong	IB-IIA
C34	Bonducellpin D	10835061	−7.9	Moderate	IB-IIA
C35	Glecaprevir	66828839	−9.7	Strong	IB-IIA-IIIB
C36	Daclatasvir	25154714	−12.1	Strong	IB-IIIA
C37	Paritaprevir	45110509	−10.5	Strong	IB-IIA-IIIA
C38	Vincapsusine	11646359	−7.6	Moderate	IB-IIA
C39	Alloyohimbine	120716	−9.7	Strong	IIA-IIIA
C40	Gummadiol	21722930	−9.9	Strong	IA-IIA
C41	ZINC000254565785		−9.2	Moderate	IIIA
C42	ZINC000726422572		−9	Moderate	IA-IIA
D43	Indomethacin	3715	−8.6	–	IIA
D44	Furosemide	3440	−8	–	IIIA
D45	Warfarin	54678486	−8.7	–	IIA-IIB
D46	Bezylpenecillin	5904	−8.4	–	IIA-IIB
D47	Chlorpropamide	2727	−7.5	–	IIA-IIIA
D48	Phenytoin	1775	−9.3	–	IA
D49	Diazepam	3016	−7.7	–	IIA
D50	Ibuprofen	3672	−7.4	–	IIA-IIB
D51	Naproxen	156391	−8.3	–	IIA-IIB
D52	Clofibrate	2796	−7	–	IIA-IIIA
D53	Chlorpromazine	2726	−7.1	–	IIA
D54	Imipramine	3696	−7.9	–	IA
D55	Quinidine	441074	−7.8	–	IB-IIIA

Fig. 2

Binding poses and molecular interactions of the drug molecules with HSA- (A) HSA_A13 (B) HSA_C34 (C) HSA_A16 and (D) HSA_D43. The domain I of HSA consists of subdomains IA coloured green (residues 5–107) and IB coloured orange (residues 108–196); domain II comprises subdomains IIA coloured blue (residues 197–297) and IIB coloured yellow (residues 298–383); domain III has subdomains IIIA coloured magenta (residues 384–497) and IIIB coloured cyan (residues 498–582). Hydrophobic interactions are shown as semi-arcs with red eyelashes, whereas hydrogen bonds are shown as green dashed lines.

Table 2

Pharmacokinetic properties of the selected molecules.

Molecule	Lipinski #violations	Veber #violations	ESOL Class	GI absorption	BBB permeant	Pgp substrate	CYP1A2 inhibitor	CYP2C19 inhibitor	CYP2C9 inhibitor	CYP2D6 inhibitor	CYP3A4 inhibitor
A1	2	1	Poorly soluble	Low	No	Yes	No	No	No	No	No
A2	0	0	Soluble	High	Yes	No	Yes	No	No	Yes	No
A3	0	0	Moderately soluble	High	No	No	No	Yes	Yes	Yes	Yes
A4	0	0	Soluble	High	Yes	No	No	No	No	No	No
A5	2	2	Moderately soluble	Low	No	Yes	No	No	No	No	Yes
A6	1	1	Poorly soluble	High	No	Yes	No	Yes	No	No	Yes
A7	0	0	Moderately soluble	High	Yes	Yes	No	No	No	No	No
A8	1	1	Soluble	High	No	No	No	No	No	Yes	Yes
A9	0	0	Moderately soluble	High	Yes	No	Yes	No	No	Yes	Yes
A10	0	0	Very soluble	Low	No	No	No	No	No	No	No
A11	0	0	Soluble	High	No	No	Yes	Yes	Yes	No	Yes
A12	0	0	Moderately soluble	High	Yes	No	Yes	Yes	Yes	No	No
A13	0	0	Soluble	High	No	No	Yes	No	No	No	No
A14	0	0	Soluble	High	No	No	Yes	No	Yes	No	No
A15	0	1	Moderately soluble	High	No	Yes	No	No	Yes	No	Yes
A16	0	1	Soluble	High	No	Yes	No	No	Yes	Yes	Yes
B17	2	1	Highly soluble	Low	No	Yes	No	No	No	No	No
B18	1	1	Moderately soluble	High	No	Yes	No	No	No	No	No
B19	2	2	Moderately soluble	Low	No	Yes	No	No	No	No	Yes
B20	2	2	Poorly soluble	Low	No	Yes	No	No	No	No	Yes
B21	2	0	Poorly soluble	Low	No	Yes	Yes	No	No	Yes	Yes
B22	0	0	Poorly soluble	High	No	Yes	Yes	Yes	No	No	No
B23	0	0	Moderately soluble	High	No	No	Yes	No	Yes	No	No
B24	1	0	Poorly soluble	Low	No	No	No	No	No	No	No
B25	0	0	Poorly soluble	High	No	No	No	Yes	No	No	No
B26	0	0	Very soluble	High	No	No	No	No	No	No	No
B27	0	0	Soluble	High	No	No	No	No	Yes	No	Yes
B28	0	0	Moderately soluble	High	No	No	No	Yes	Yes	No	Yes
B29	2	2	Poorly soluble	Low	No	Yes	No	No	No	No	Yes
B30	2	2	Poorly soluble	Low	No	Yes	No	No	No	No	Yes
B31	1	2	Moderately soluble	Low	No	Yes	No	No	No	No	Yes
C32	1	0	Moderately soluble	High	No	Yes	No	No	Yes	Yes	Yes
C33	1	0	Moderately soluble	High	No	Yes	No	No	No	Yes	Yes
C34	0	0	Soluble	High	No	Yes	No	No	No	Yes	No
C35	2	1	Poorly soluble	Low	No	Yes	No	No	No	No	No
C36	2	2	Poorly soluble	Low	No	Yes	No	No	Yes	No	Yes
C37	2	1	Poorly soluble	Low	No	Yes	No	No	No	No	Yes
C38	0	0	Soluble	High	Yes	No	No	No	No	Yes	No
C39	0	0	Moderately soluble	High	Yes	Yes	No	No	No	Yes	No
C40	0	0	Soluble	High	No	No	No	No	No	Yes	No
C41	0	0	Moderately soluble	High	Yes	No	Yes	Yes	Yes	No	Yes
C42	0	0	Moderately soluble	High	No	No	Yes	Yes	Yes	Yes	Yes
D43	0	0	Moderately soluble	High	Yes	No	Yes	Yes	Yes	No	No
D44	0	0	Soluble	High	No	No	No	No	No	No	No
D45	0	0	Soluble	High	Yes	No	No	Yes	Yes	No	No
D46	0	0	Soluble	High	No	No	No	No	No	No	No
D47	0	0	Soluble	High	No	No	No	No	No	No	No
D48	0	0	Soluble	High	Yes	No	No	No	No	No	No
D49	0	0	Soluble	High	Yes	No	Yes	Yes	Yes	Yes	Yes
D50	0	0	Soluble	High	Yes	No	No	No	No	No	No
D51	0	0	Soluble	High	Yes	No	No	No	No	No	No
D52	0	0	Soluble	High	Yes	No	Yes	Yes	No	No	No
D53	1	0	Moderately soluble	High	Yes	No	Yes	Yes	Yes	Yes	No
D54	0	0	Moderately soluble	High	Yes	No	Yes	No	No	Yes	No
D55	0	0	Soluble	High	Yes	No	No	No	No	Yes	No

The free HSA and HSA-ligand-bound complexes (HSA_A13, HSA_C34, HSA_A16 and HSA_D43) were subjected to 100 ns of MD simulation to determine the stability of their trajectories. The geometric properties analysed in the study include RMSD, RMSF, Rg, SASA and NHB (Table 3). Root mean square deviation (RMSD) computes the average distance between the backbone atoms of starting structure (reference structure) with simulated structures (frame by frame) when superimposed. The RMSD plot of backbone atoms of the free HSA and HSA-ligand-bound complexes show stable fluctuations that sustained till the end of the simulation period (Fig. 3). The average RMSD of the backbone atoms of free HSA, HSA_A13, HSA_C34, HSA_A16 and HSA_D43 complexes were 0.448642 ± 0.048833 nm, 0.462448 ± 0.064559 nm, 0.490956 ± 0.058781 nm, 0.589058 ± 0.078123 nm and 0.478107 ± 0.0692 nm respectively (Table 4). A higher RMSD was seen in HSA-ligand complexes which may be due to structural modulation of the protein upon binding of the ligand molecules. The average RMSD of the heavy atoms of A13, C34, A16 and D43 were 1.457102 ± 0.504222 nm, 0.609009 ± 0.092277 nm, 0.322053 ± 0.045001 nm and 0.612155 ± 0.058157 nm respectively (Table 4). The RMSF plot was generated (Fig. 4) which provide further insights into the structural regions of the protein contributing to greater fluctuations. Root mean square fluctuation (RMSF) computes fluctuations (standard deviation) of atomic positions of each amino acid residue in the trajectory. The binding of A13 to HSA causes a significant increase in the amplitude of fluctuations around Ser58-Glu60, Arg114-Arg117 and Asp301 (Fig. 4A) and the region encompassing residues Arg114-Arg117, Lys564 and Glu565 witnessed enhanced amplitude of fluctuations upon binding of D43 (Fig. 4D). It is interesting to note that the residues Asp301 and Glu321 in the free HSA show increased amplitude of fluctuations compared to the HSA-ligand complexes (Fig. 4). Similarly, binding of A16 to HSA causes an increase in the amplitude of fluctuations particularly in the region encompassing Pro113 (Fig. 4C). The radius of gyration (Rg) computes the structural compactness of a molecule and the radii of gyration about the x-, y- and z-axes, as a function of time. From the Rg plot (Fig. 5) it is evident that free HSA, HSA_A13, HSA_C34, HSA_A16 and HSA_D43 complexes show a stable Rg graph from 10 ns onwards which averaged at 2.618873 ± 0.022697 nm, 2.655904 ± 0.026257 nm, 2.633126 ± 0.024887 nm, 2.587976 ± 0.030269 nm and 2.618333 ± 0.03179 nm respectively (Table 3). This indicates that the overall compactness of HSA is altered upon binding of ligand molecules. Fig. 6 represents the plot of total SASA (hydrophobic + hydrophilic) of the free HSA and HSA-ligand complexes with time. Solvent-accessible surface area (SASA) can be explained as an approximate area of a biomolecular surface that is accessible to solvent molecules as a function of time. The average total SASA value for free HSA, HSA_A13, HSA_C34, HSA_A16 and HSA_D43 was calculated to be 264.0154 ± 6.698827 nm2, 271.306 ± 7.533425 nm2, 267.8342 ± 5.879287 nm2, 265.9578 ± 8.209247 nm2 and 269.5354 ± 8.744237 nm2 respectively (Table 3). The increase in total SASA for HSA-ligand complexes indicates that the folding of the protein is altered upon binding of the ligand molecules. Fig. 7 A and B represent the plot of the number of hydrogen bonds within the protein (intramolecular hydrogen bonds) and HSA-ligand complexes respectively. The number of intramolecular protein hydrogen bonds for free HSA, HSA_A13, HSA_C34, HSA_A16 and HSA_D43 were 482.4416 ± 11.10697, 470.2488 ± 11.8971, 481.6254 ± 11.67786, 475.4775 ± 11.42339 and 475.3816 ± 11.46369 respectively (Table 3). A decrease in the number of intramolecular hydrogen bond formations may be attributed to the less folded state of the protein upon ligand binding. The number of hydrogen bonds between HSA and ligand in HSA_A13, HSA_C34, HSA_A16 and HSA_D43 were 1.082917 ± 0.867247, 0.588412 ± 0.717235, 2.310689 ± 1.061308 and 1.641359 ± 1.008091 respectively (Table 3).

Table 3

Average geometric properties of the systems calculated during 100 ns MD simulation studies.

Systems	RMSD (nm)		Rg (nm)	Total SASA (nm2)	Number of Hydrogen bonds
	Protein	Ligand			Protein	Protein-Ligand
HAS	0.448642 ± 0.048833		2.618873 ± 0.022697	264.0154 ± 6.698827	482.4416 ± 11.10697
HSA_A13	0.462448 ± 0.064559	1.457102 ± 0.504222	2.655904 ± 0.026257	271.306 ± 7.533425	470.2488 ± 11.8971	1.082917 ± 0.867247
HSA_C34	0.490956 ± 0.058781	0.609009 ± 0.092277	2.633126 ± 0.024887	267.8342 ± 5.879287	481.6254 ± 11.67786	0.588412 ± 0.717235
HSA_A16	0.589058 ± 0.078123	0.322053 ± 0.045001	2.587976 ± 0.030269	265.9578 ± 8.209247	475.4775 ± 11.42339	2.310689 ± 1.061308
HSA_D43	0.478107 ± 0.0692	0.612155 ± 0.058157	2.618333 ± 0.03179	269.5354 ± 8.744237	475.3816 ± 11.46369	1.641359 ± 1.008091

Fig. 3

The RMSD plot of backbone atoms of the free HSA and HSA-ligand complexes.

Table 4

Statistics of binding free energy calculation of HSA-ligand complexes (kcal/mol).

Protein-ligand complexes	ELE1	VDW2	GAS3	PBSOL4	PBTOT5	GBSOL4	GBTOT5	-TS6	△GPB7	△GGB7
HSA_A13	−2.87 ± 1.23	−40.49 ± 1.95	−43.36 ± 2.46	12.10 ± 2.24	−31.26 ± 2.43	9.61 ± 1.45	−33.75 ± 1.99	14.46 ± 1.65	−16.80	−19.29
HSA_C34	−6.76-±3.16	−50.79 ± 2.92	−57.54 ± 3.52	24.72 ± 3.08	−32.83 ± 2.69	18.29 ± 2.43	−39.25 ± 2.61	22.73 ± 2.74	−10.10	−16.52
HSA_A16	−6.45 ± 2.00	−65.66 ± 2.67	−72.12 ± 3.37	25.51 ± 2.48	−46.60 ± 3.44	16.39 ± 1.86	−55.72 ± 2.94	22.38 ± 2.19	−24.22	−33.34
HSA_D43	15.79 ± 5.50	−37.93 ± 2.87	−22.15 ± 6.28	−9.15 ± 5.10	−31.30 ± 3.07	−5.57 ± 5.30	−27.72 ± 2.74	15.76 ± 2.06	−15.54	−11.96

Electrostatic energy as calculated by the MM force field; 2Van der Waals contribution from MM; 3Total gas-phase energy; 4Non-polar and polar contributions to solvation based on PB/GB model; 5Final estimated binding free energy calculated from GAS and PBSOL/GBSOL; 6Entropy; 7Binding free energy with entropy

Fig. 4

RMSF plot of backbone atoms of HSA-ligand complexes-(A) HSA_A13 (B) HSA_C34 (C) HSA_A16 and (D) HSA_D43 compared to the free HSA. The labelled regions indicate high fluctuating regions.

Fig. 5

Plot of Radius of gyration versus time for free HSA and HSA-ligand complexes.

Fig. 6

Plot of total SASA for free HSA and HSA-ligand complexes versus time.

Fig. 7

Variations in the number of intramolecular hydrogen bonds with time free HSA and HSA-ligand complexes (A) intermolecular hydrogen bonds (B) protein-ligand hydrogen bonds.

Molecular mechanics Poisson–Boltzmann surface area (MM/PBSA) and molecular mechanics generalized Born surface area (MM/GBSA) binding free energies of A13 (ΔPB = −16.80 kcal/mol, ΔGB = −19.29 kcal/mol), C34 (ΔPB = −10.10 kcal/mol, ΔGB = −16.52 kcal/mol), A16 (ΔPB = −24.22 kcal/mol, ΔGB = −33.34 kcal/mol) and the control (D43) (ΔPB = −15.54 kcal/mol, ΔGB = −11.96 kcal/mol) were tabulated in Table 4. In each of the HSA-ligand complexes, van der Waals energy is the dominant force of molecular interactions. The top ten residues which contribute significantly towards the binding interaction between A13 and HSA include Arg209, Leu347, Ala213, Leu327, Ala210, Lys351, Gly328, Leu331, Ala350 and Val482 (Fig. 8A). In case of HSA_C34 complex, residues such as Gln196, Lys195, Trp214, Ser192, Arg257, Tyr150, Leu238, Glu292, Lys199 and Phe149 has higher contributions (Fig. 8B). The major residues contributing towards the binding interaction between A16 and HSA include Glu450, Trp214, Ser454, Leu481, Lys199, Val344, Leu198, Lys195, Arg485 and Leu457 (Fig. 8C). In case of HSA_D43 complex, residues such as Arg257, Leu238, Leu260, Ser287, Ile290, Leu219, Ala291, Arg222, Ala261 and Ile264 has higher contribution to the total binding energy (Fig. 8D).

Fig. 8

Heatmap showing the top ten residues contributing significantly to the total binding free energy of HSA-ligand complexes-(A) HSA_A13 (B) HSA_C34 (C) HSA_A16 and (D) HSA_D43.

Discussion

There is a high paucity of new antiviral drugs against the COVID-19 pandemic which is caused by a novel betacoronavirus known as SARS-CoV-2. Numerous clinical trials on possible antiviral treatments are still underway and the majority of these drug candidate molecules are the repurposed drugs. The drug molecules in the body are bound to the drug carrier protein know as human serum albumin which influences their activity and half-life. In the process of developing novel therapies, the binding affinity of the drugs towards HSA is an important parameter of study given the large quantity of HSA in plasma. In addition, the therapeutic efficacy of drugs is altered upon binding to HSA (Tesseromatis and Alevizou, 2008). Human serum albumin (HSA) is a non-glycosylated polypeptide with a molecular weight of 66,500 Da (Yamasaki et al., 2013). The polypeptide chain forms a heart-shaped structure with a dimension of approximately 80 × 80 × 30 Å and (Sugio et al., 1999). The presence of three domains, namely domains I (residues 1–195), II (196–383), and III (384–585) are structural characteristics of HSA (He and Carter, 1992, Yang et al., 2014). In our studies, we have characterized the binding energy and putative binding sites of selected anti-COVID-19 drug candidate molecules to HSA. Comparing the binding energies of the control drugs to HSA, we have classified the strength of their binding into weak, moderate, and strong. The weak binding molecules include four molecules such as Chloroquine (A9), β-D-N4-hydroxycytidine (A10), Disulfiram (A11), and Carmofur (A13). Molecules such as Ivermectin (A1), Hydroxychloroquine (A2), Arbidol (A3), Ebselen (A4), Remdesivir (A5), Lopinavir (A6), Emetine (A7), Homoharringtonine (A8), Shikonin (A14), Compound 11a (A15), Zanamivir (B17), Atazanavir (B20), Bedaquiline (B21), Brequinar (B22), Clofazimine (B24), Gemcitabine (B26), Tolcapone (B27), Vismodegib (B28), Cobicistat (B29), Darunavir (B31), Bonducellpin D (C34), Vincapsusine (C38), ZINC000254565785 (C41) and ZINC000726422572 (C42) exhibited moderate binding to HSA. The molecules strongly binding to HSA comprises Tideglusib (A12), Compound 11b (A16), Indinavir (B18), Saquinavir (B19), Celecoxib (B23), Conivaptan (B25), Ritonavir (B30), Ergotamine (C32), Dihydroergotamine (C33), Glecaprevir (C35), Daclatasvir (C36), Paritaprevir (C37), Alloyohimbine (C39) and Gummadiol (C40). These molecules bind to different domains of HSA through both hydrogen bonds as well as hydrophobic interactions. We have selected one drug molecule from each category-A13 (weak binding), C34 (moderate binding), A16 (strong binding), and a control (D43) for molecular dynamics simulation studies. The dynamic structural changes in the HSA upon binding of these molecules resulted in variations in the geometric properties such as RMSD, RMSF, Rg, total SASA, and the number of hydrogen bonds. Molecular mechanics Poisson–Boltzmann surface area (MM/PBSA) and molecular mechanics generalized Born surface area (MM/GBSA) binding energies of the molecules exhibited negative ΔG values indicating the spontaneous binding interaction of these molecules with HSA. Due to the computational limits in our present studies, the molecular dynamics simulations of only a few selected molecules were performed. It is worthwhile to investigate the number of binding sites, thermodynamic parameters, and dissociation constants using experimental techniques such as fluorescence quenching, isothermal titration calorimetry (ITCs) and surface plasmon resonance (SPR) which will provide further insights into the binding of COVID-19 drug candidate molecules to HSA.

Conclusion

Human serum albumin (HSA) is an important drug carrier protein that modulates the activity and half-life of the majority of drug molecules. Using molecular docking and dynamics simulation techniques, we were able to characterize the binding sites of the selected experimental and in silico COVID-19 drug candidates. Our study unravels that these drug molecules bind to different domains of HSA using intermolecular forces such as hydrogen bonds and hydrophobic interactions. The drug molecules also show different degrees of binding strength to HSA and the HSA-drug complexes were also stable throughout the molecular dynamics simulation time. Various experimental techniques such as fluorescence quenching, isothermal titration calorimetry (ITC) and surface plasmon resonance (SPR) will shed further light on the thermodynamic binding properties of COVID-19 drugs to HSA.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper