Molecular Dynamic Simulation of Centella asiatica Compound as an Inhibitor of Advanced Glycation End Products

Lili Legiawati1, Fadilah Fadilah2, Kusmarinah Bramono1, Siti Setiati3*, Em Yunir3 1Department of Dermatology and Venereology, Faculty of Medicine, Universitas Indonesia, Cipto Mangunkusumo Hospital, Jakarta, Indonesia. 2Department of Pharmacy, Faculty of Medicine, Universitas Indonesia, Depok Beji, Indonesia. 3Department of Internal Medicine, Faculty of Medicine, Universitas Indonesia, Cipto Mangunkusumo Hospital, Jakarta, Indonesia.


INTRODUCTION
The glycation process happens naturally in the human body as the process of aging. However, in diabetes mellitus (DM) patient, glycation process becomes uncontrollable, resulting in the accumulation of advanced glycation end products (AGEs) as firstly identified in 1986 (Fournet et al., 2018;Singh et al., 2014;Rahbar, 1968). The accumulation of AGEs and interaction of AGEs/receptor AGEs (RAGE) are the basis for the development of complications in DM, including cardiovascular disease, neurodegenerative disorder, chronic kidney disease, and ophthalmology disease (Sarwar and Gao, 2010). As in 2017, it was reported that 429 million patients suffered from DM with a financial burden of USD 727 billion, which estimated to grow larger in 2045, reaching 628.6 million patients (IDF, 2017).
The formation of AGEs occurs via three different parts, which are nonenzymatic Maillard reaction, polyol pathway, and lipid peroxidation (Nowotny et al., 2015). Then, AGEs will accumulate in various organs leading to the disruption of intracellular and extracellular protein structure and function (Singh et al., 2014). AGEs also initiate a signaling pathway mediated by several membrane receptors. The most studied AGEs receptor is the type 1 cell surface RAGE. The binding of AGEs and RAGE is responsible for reactive oxygen species generation and inflammatory responses through the expression of mitogen-activated protein kinase pathway, nuclear factor-kB, interleukin-6, tumor necrosis factor-α, intercellular adhesion molecule-1, and vascular cell adhesion molecule-2 (Grimm et al., 2018;Ott et al., 2014). Several studies of drugs that inhibit AGEs directly or indirectly by disrupting AGEs and RAGE have shown some promising results as the prevention of DM complication. Elparestat lowered the accumulation of AGEs effectively in peripheral neuropathy, whereas hesperidin lowered AGEs in the retina (Amore et al., 1997;Schmidt et al., 1994). The therapy using a synthetic inhibitor of AGEs, although promising, has the downside risk, which is a more prominent side effect. Therefore, the studies had leaned into searching for alternative traditional inhibitor for AGEs. Several studies on alternative traditional inhibitors, including grape seed proanthocyanidins extract (Cui et al., 2008), curcumin (Stefanska, 2012), Litsea japonica (Sohn et al., 2013), and Korea ginseng extract (Quan et al., 2013), had shown that these agents were comparable as AGEs inhibitors.
Centella asiatica (CA), a widely used traditional medicinal plant in Indonesia, was reported to have some clinical uses for preventing DM complications. Active compounds of CA consist of pentacyclic triterpenes, mainly asiatic acid, madasiatic acid, madecassic acid, asiaticoside, and madecassoside. Moreover, there were several other triterpenes, including betulinic acid, isothankunic acid, terminolic acid, and centellasapogenol A (Orhan, 2012). These active compounds have been proven to possess anti-inflammation and anticancer properties, which are beneficial for attenuating glycative process. Pentacyclic triterpenes also had anti-AGE and anti-RAGE agents based on several in vivo studies (Yin, 2015). Although the previous studies not explicitly focused on a specific active compound of CA, an in vitro study of CA proved its ability of lowering glycative activity. The CA purified extract had been shown to inhibit the N ε -carboxymethyllysine formation totally, which is one of the significant constituents of AGEs, induced by methylglyoxal (Maramaldi et al., 2014). This inhibition potential could be explored via in silico method to understand the mechanism of the action of CA better. The in silico method provides current development of the method in testing pharmacology hypothesis via measuring the affinity of a new molecule to a target and assessing the physiochemical properties (absorption, distribution, metabolism, and excretion) of new drugs (Ekins et al., 2007).
This study aims to explore the molecular docking of CA active compound as an inhibitor of AGEs and RAGE. Furthermore, molecular dynamic (MD) simulations were used in the study to assess the stability of binding complex. We also analyzed the probable molecular mechanisms of the action of CA active compounds on the AGEs and RAGE.

Ligand and protein preparation
The protein targets of AGEs and RAGE were generated from three-dimensional structure modeling in.pdb format from Protein Bank Data (PDB) (www.pdb.org) code with PDB ID: 3DKK for AGE and PDB ID: 4YBH for RAGE. The AGE protein structure (PDB ID 3DKK) named Aged Form of Human Butyrylcholinesterase inhibited by Tabun as the native ligand. The RAGE protein structure belongs to the crystal structure of the human RAGE ectodomain (VC1C2 fragment) in complex with human S100A6 as the native ligand. The three-dimensional structures were obtained from www.pdb.org. AGEs and RAGE underwent an optimization process using docking preparation, which consisted of hydrogen addition or current addition through Amber ® 14 method. The residue preparation was facilitated according to default in Python Molecular Viewer ® (PMV). The co-crystal ligand chain of AGEs and RAGE in the target was separated from the protein chain. The co-crystal ligand was selected as crystallography ligand according to the interaction model from PDB. The co-crystal ligands were prepared by hydrogen addition or charge addition. Nine of CA compounds were selected from PubChem ® database search engine. The threedimensional structure was then optimized by PMV ® .

Docking validation
The prepared protein target and crystallography ligand were later put into PyRx-Autodock-Vina. Gridbox was set in protein receptor location with coordinate X:1,291 Y:7,504 Z:37,569 and dimension X:30,031 Y:28,790 Z:27,798.
Crystallography ligand proteins from PDB were optimized for the docking process. The optimization of proteins and ligands was conducted four times. After the docking process, ligands were saved and then measured against crystallography ligand. The measurement resulted in the root-mean-square deviation (RMSD). Ligands with RMSD < 2 were eligible for the next process. Proteins from docking process and derivatives from ligand were put into PyRx ® program. The results were saved in.csv format and.sdf format. The analytical process was based on binding energy, binding affinity, and interaction residue from molecular docking. This study used molecular docking to investigate the possible mechanisms of interaction of CA compounds. The compound with a significant in silico score was deemed possible as an inhibitor.
Chimera program was used for optimization with Dunbrack Rotamer Library, which served for repairing incomplete protein residue. An improvement of protein residue charge was based on AMBER and also Chimera program. The purpose was to adjust the protein charge for biomolecular simulation.

Dynamic simulation
This study performed MD simulation using the Amber ® 14 software package (Case et al., 2014). MD simulation was the prevailing sampling algorithms to investigate the interaction between a ligand and a protein target. The docked structure of AGEs (PDB ID: 3DKK) was reported as the fourth most active compound and novel chemicals with high activity among in silico studies. Therefore, we used these docked structures as the initial structures for MD simulations. The missing hydrogen atoms of the AGEs were added by the LEaP module of the Amber ® 14 package to maintain the electroneutrality. Each complex was immersed into a cubic periodic box of TIP3P water model with at least 10˚A distance around the complex (Wang et al., 2004). For the ligands, the GAFF parameter assignments were made by using Antechamber program, and then, the partial charges were assigned by using the AM1-BCC method (Hou et al., 2011;Klebanov et al., 1998). Afterward, four 10-ns results from MD simulations were carried out with the PMEMD program without limit in the isothermal, isobaric ensemble (NPT, P = 1 atm, and T = 310 K MD. The time step was set at 2 fs. The cutoff of 10˚A was applied to treat nonbonding interactions.

RESULTS AND DISCUSSION
The compounds were screened based on the energy variation and the formation of the ligand-receptor structure, given by the binding constant and the Gibbs free energy (∆G) values. We analyzed ligand-protein binding by evaluating their binding affinity, binding structure, binding energy, and residue component.
Binding energy values were calculated as the sum of intermolecular energies (kcal/mol), comprising hydrogen bond energy, Van Der Waals energy, desolvation energy, and electrostatic energy. Binding energy represents the energy released due to the interaction between ligand and receptor. The ligand with the least amount of binding energy had the best possibility to interact with protein target. The value of pKi represents the level of affinity of the compound to the molecular target.
This study used AGE target model from Protein Data Bank named Aged Form of Human Butyrylcholinesterase inhibited by Tabun with the PDB ID: 3DKK, which contained with 529 amino acid. This three-dimensional structure has resolution 2.31 Å, R-value work 0.197, and R-value observed 0.200. The active side or receptor site for this enzyme was formed from the incorporation of both molecule chains with the primary residue. This protein target's binding site was located in Asn106, Asn188, Lys190, Se191, Arg219, Asn241, Asn245, Phe278, Pro281, Gln316, Asp324, Ser338, Arg341, Gly413, Tyr420, Arg465, Asn485, and Tyr500.
In visual observation through the docking simulation (Figure 1), the three observed ligands (asiaticoside, madasiatic acid, and madecassic acid) made similar conformations. Asiaticoside, madasiatic acid, and madecassic acid had similar conformations in the core chain and side chain, and the only differences are asiaticoside glucose chain and kiral compound. Brahmol showed the most distinctive conformation from other ligands due to Brahmol's low molecular weight and simpler side chain compared to the other compounds.
To verify whether the studied systems reached equilibrium, the RMSD of all residues of the active site (residues within 5 Å around ligand) and the heavy atoms of the ligand from the initial structure were monitored to examine the dynamic stability of the systems and plotted against time, as shown in Figure 2. The three RMSDs had small fluctuations after 12 ns, implying that the studied systems had reached the stability. We used the last 20 ns to analyze the energy and binding modes for the four complexes.
Asn188, 278, 413, and 465 as a binding site which occupied Gly413 in madasiatic acid exhibited the inhibitory effect toward the protein activity. Nonetheless, MD simulation data annexed the inhibition level based on the number of occupancy binding to the active site (Asp324). Madasiatic acid, asiaticoside, and native ligand had occupied the active site Asp324 with a binding distance of 6.80, 20.90, and 28.90, respectively. Compared to native ligand, asiaticoside possessed a more stable inhibition level than madasiatic acid.
To obtain a detailed interaction between the two ligands and AGEs, we investigated the key residues during the binding process. The key residues were identified with the binding free energy decomposition calculated through the Molecular Mechanics/ Generalized Born Surface Area (MM/GBSA) method. Energy decomposition was calculated from Van der Waals, electrostatic, solvation-free energy, and total energy contribution terms. The residues with energy contributions more than 1.5 kcal/mol were selected. The results are shown in Table 2. Asn106, Asp324, Asp376, Tyr420, and Tyr500 residues of AGEs made a significant contribution to the asiaticoside AGE complex. Meanwhile, Asn118 and Tyr500 residues made a significant contribution to the madasiatic acid AGE complex. The vast majority of the crucial residues of AGE was nonpolar, so it was reasonable to speculate that these residues were capable of forming more significant Van der Waals interactions with hydrophobic ligand. Therefore, the nonpolar attributes also exhibited a more favorable contribution of nonpolar interaction to the binding free energy.
Based on the result of docking with RAGE (Table 3), there were three ligands with the least binding energy for RAGE which are asiaticoside, isothankunik acid, and asiatic acid with a binding energy of −10.6125, −9.4469, and −9.1015 kcal/ mol, respectively. The most common residues found from these interactions were Thr55, Arg57, and Glu103.
To explore the possible binding mode of the selected compounds, thus molecular docking was engaged to study the interaction in Figure 4 (especially, asiaticoside and isothankunik acid), together with the reported most active compound native ligand as a comparison, receptor of advanced glycation end products (RAGE) (PDB ID: 4YBH). The three-dimensional crystal structure of the human RAGE ectodomain (VC1C2 fragment) in complex with human S100A6 (PDB ID:4YBH) has a resolution of 2.40 Å with R-value work 0.195 and R-value observed 0.197. The Protein Data Bank has shown that RAGE was Aged Form of Human covered with 304 amino acid sequence length and no mutation. This target has binding site of Thr28, Thr55, arg57, Trp72, Arg78, Leu79, Asp93, Glu132, Lys169, His180, Asp201, His218, Arg228, and His270.  Based on molecular dynamic evaluation (Table 2), we have discovered that asiaticoside has very close similarity in hydrogen bond occupancy with native ligand (L-peptide linking Tabun) in amino acid residues: Asp324, Val377, and Tyr500. Madasiatic acid has a strong similarity with the native ligand (L-peptide linking Tabun) in amino acid residues: Asn188 and Tyr420.   The interaction identified between CA and AGEs showed the CA's role as AGE inhibitor. The individual substances tend to have different mechanisms as an AGEs inhibitor. Asiaticoside, madasiatic acid, and madecassic acid bound to the binding site of AGEs; thus, they are capable of preventing AGEs to bind to RAGE directly. The mechanism of modifying AGEs binding site to prevent an interaction with RAGE is similar to the mechanism of known AGEs inhibitor, carnosine, a dipeptide consisting of b-alanine and histidine. As noted in the carnosine study, the modification of the binding site also prevented cross-linking between AGE or AGE formation (Klebanov et al., 1998). In the previous study of CA compound, specifically studying asiaticoside, the result showed that asiaticoside lowered the level of AGEs (Schmidt et al., 1994).
In this study, we also discovered that asiaticoside, asiatic acid, and isothankunic acid bound to the binding site of RAGE. The mechanism of direct inhibition to the binding site of RAGE was previously studied in azeliragon. The capability of a ligand to bind to RAGE binding domain blocked the capability of another ligand to bind to RAGE (Burstein et al., 2014). In the study by Wang, the asiatic acid had lowered the level of RAGE in human skin keratinocyte (HaCat) (Wang, 2014).
Currently, the CA studies in animal model have provided evidence of CA role on oxidative stress. CA prevents cognitive deficit by the attenuation of oxidative stress and anti-inflammation activity, as well as restoring cholinergic function (Chiroma et al., 2019). Specifically, CA ability to enhance antioxidant activity, e.g., superoxide dismutase (Chiroma et al., 2019) catalase, and glutathione peroxidase (Ahmad Rather et al., 2019) while decreasing oxidative stress, e.g., malondialdehyde (MDA) (Chamaco-Alonso et al., 2019) and reactive nitrogen species (RNS) (Visnawathan et al., 2019) balances the antioxidant-prooxidant condition which is usually impaired in cognitive diseases, such as Alzheimer Disease. Centella asiatica was found to enhance antioxidant pathway through promoting nuclear transcription factor (erythroid-derived 2)-like 2 (NRF2) (Matthews et al., 2019) and activate Akt/GSK3β pathway (Ahmad Rather et al., 2019).

CONCLUSION
In this study, we applied the docking studies to CA active compounds with AGEs and RAGE. AGEs had been found to react with significant binding energy to asiaticoside, madasiatic acid, and madecassic acid with a binding energy of −11.8253, −10.6724, and −10.1462, respectively. Through MD simulation, asiaticoside and madasiatic acid bound into a similar site of native ligand with significant analyzed binding energy decomposition, thus confirming their role as an inhibitor of AGEs. CA active compound had also been proven to interact with AGEs and RAGE. The docking studies showed that the ligands RAGE bound to asiaticoside, isothankunik acid, and asiatic acid with a binding energy of −10.6125, −9.4469, and −9.1015 kcal/mol, respectively. Although this study proved the inhibition AGEs and RAGE by CA via in silico study, the findings were preliminary. These ligands can be further assessed for their effectiveness in clinical trials.

CONFLICT OF INTEREST
Authors declared that they do not have any conflicts of interest.

FINANCIAL SUPPORT
None.