Sage Journals: Discover world-class research

Abstract

Since late 2019 to early 2020, an outbreak caused by severe acute respiratory syndrome-coronavirus 2 (SARS-CoV-2) has become a worldwide health emergency due to its rapid infection and mortality of millions of people around the world. As the main protease M^pro or 3CLpro produced by the virus plays an important role in coronavirus survival and proliferation, it becomes an excellent drug target to identify COVID-19 inhibitors. Lawsonia inermis L. (henna) is a medicinal plant that has been used for a long time for the treatment of many fungal and bacterial infections. In the search for new anti-COVID agents from medicinal plants, we report the results of our study into the potential inhibition of M^pro by the compounds isolated from the extracts of L. inermis roots and leaves using molecular docking and molecular dynamics simulation. The molecular modeling results showed that all isolated compounds bonded spontaneously into the catalytic pockets of M^pro with binding energies <0. The docking and calculated pharmacokinetic results of the compounds (1-3, 6-8) were similar to and even better than those of the commercial COVID-19 inhibitor remdesivir. In particular, the triterpenoid glycoside suavissimoside R1 (8) showed the best binding to SARS-CoV M^pro, with the lowest binding energy ΔG and IC_50,calc. values of −8.19 kcal/mol and 0.98 μM, respectively. Furthermore, the calculations of ADMET (absorption, distribution, metabolism, excretion, and toxicity) showed that it had the lowest toxicity, with a predicted LD₅₀ value of 3320 mg/kg. These triterpenoids are worthy of further study to evaluate their actual bioactivity against SARS-CoV-2 in vitro and in vivo in the hope of contributing valuable scientific data for natural resources for the development novel drug formulations for either the prevention or treatment of COVID-19.

Keywords

henna protease of SARS-CoV-2 COVID-19 ADMET suavissimoside R1

Introduction

A novel deadly coronavirus known as severe acute respiratory syndrome-coronavirus 2 (SARS-CoV-2) has led to a worldwide health emergency.¹ The infection caused by this virus to the lungs and respiratory system is so severe that it has already killed more than 5.9 million people around the world till 24 February 2022 and infected more than 428 million people since its first case was reported in December 2019 in Wuhan city of China.² The virus produces many kinds of enzymes to enable it to invade the lung cells of humans. There are several druggable enzymes/proteins of SARS-CoV-2 that have attracted the attention of scientists around the world as research targets, including the nonstructural protein (nsp 1-16), RNA polymerase, helicase, the main structural proteins (spike [S], membrane [M], envelope [E], nucleocapsid phosphoprotein [N]), and papain-like protease (PLpro).³ Out of these drug targets, the main focus has been on M^pro or 3CLpro because they act as essential multitasking polyproteins for virus survival, to facilitate the vital processes like replication, transcription, and translation of viral RNA into functional non-structural proteins.^1,4 Inhibiting the activity of this M^pro protease would block viral proliferation and invasion.⁵ Hence, this protease has become an excellent drug target among all the proteins of the coronavirus to identify COVID-19 inhibitors. In addition, natural products from medicinal plants have become an important source for the discovery of therapeutic agents, especially in the fight against COVID-19.⁶ By using computer-aided techniques for the prediction of ligand–target interactions,⁷ many natural compounds of different classes, such as baicalin, chrysin (Scutellaria baicalensis) (flavonoids), sugetriol-3,9-diacetate (Cyperus rotundus) (sesquiterpenoid), phaitanthrin D, 2,2-di (3-indolyl)-3-indolone (Isatis indigotica) (alkaloids), and platycodin D (Platycodon grandiflorus) (triterpenoid) have been reported to be potential new scaffolds in the treatment of SARS-CoV-2 as they exhibited high binding affinity to M^pro protein.^6,8

Lawsonia inermis L. (syn. L. alba Lamk.), known as henna or hina, is a flowering plant belonging to the family Lythraceae that was first scientifically described by Carl von Linné (Linnaeus) in 1753. The species is named after the Scottish physician Isaac Lawson, a good friend of Linnaeus.⁹

The henna plant has been widely grown in many countries and continents around the world for more than three thousand years, with a variety of uses, including ornamental, dyeing, and medicine.¹⁰ Henna has also been used for a long time as a colorant in many countries. Henna is used as a body dye, especially in India and the Middle East, for tattoos, hair, beard, eyebrows, skin decoration, and nails. Henna powder is mixed with shampoo and other hair dyes in different proportions to produce colors ranging from bronze to reddish. It is also used for fabric dying, including silk, textiles, wool, and leather. In folk medicine, the bark is traditionally used in the treatment of jaundice, spleen enlargement, renal calculus, leprosy, and obstinate skin diseases. In Morocco, the infusion of leaves of L. inermis is used against diarrhea, renal lithiase, and gastric pains.¹¹ In Vietnamese traditional medicine, it is used to treat menstrual disorder, edema, rheumatism, bronchitis, and hemorrhoids.¹² Henna is used in the treatment of hair and scalp problems and for hair loss.¹³ The plant has properties that effectively prevent the growth of bacteria,¹⁴ prevent and treat fungal infections, have analgesic and anti-inflammatory effects on the human skin, and are effective against lice and sunburn.¹⁰^,15 Phytochemical studies of leaves, barks, and flowers of L. inermis revealed the presence of many kinds of compounds including flavonoids, alkaloids, terpenoids, sterols, and especially triterpenoids.^10,16

In the search for new inhibitory scaffolds from Vietnamese medicinal plants for the treatment of the SARS-CoV-2 triggered global health coronavirus disease 2019, this study was designed to evaluate the potential biological activity of natural compounds (1-12) isolated from different parts of L. inermis (roots and leaves) against the main protease M^pro from SARS-CoV-2 using several docking methods. The hope was that these docking results would help find new anti-CoV-2 lead candidates for further investigations relating to in vitro, in vivo, and pre-clinical antiviral experiments.

Materials and Methods

General Experimental Procedures

¹H-NMR (500 MHz) and ¹³C-NMR (125 MHz) spectra were measured on a Bruker AVANCE 500 MHz spectrometer with tetramethylsilane (TMS) as an internal standard, and chemical shifts are expressed in ppm. ESI-MS were obtained from a Varian FT-MS spectrometer and MicroQ-TOF III (Bruker Daltonics Inc.). Optical rotations were measured on DIP-2000 KUY polarimeter (JASCO, Tokyo, Japan), and UV and IR spectra on a UV-VIS spectrophotometer (Varian Cary, Australia) and a FT-IR spectrometer 1650-Perkin Elmer (USA), respectively. Column chromatography (CC) was carried out on silica gel (Si 60 F₂₅₄, 230-400 mesh, Merck). All solvents were redistilled before use. Analytical and preparative TLC were performed on pre-coated Kieselgel (Si 60 F₂₅₄) or RP18 plates (0.25 mm, Merck). Compounds were visualized under UV radiation (254, 365 nm) and by spraying plates with 10% H₂SO₄, followed by heating with a heat gun.

Plant Materials

L. inermis whole plants (roots, leaves, and stems) were collected from in Ngoc Truc Village, Dai Mo Commune, Ha Dong District, Hanoi. Specimens was identified at the National Institute of Medicinal Materials (NIMM), and a voucher specimen (HNIP/16698/08) was deposited at the Herbarium of the Department of Botany at Hanoi University of Pharmacy (HUP), Hanoi, Vietnam.

Extraction and Isolation

The roots of L. inermis were dried in the shade at room temperature and ground into powder. Dried powdered roots (3.0 kg) were extracted with MeOH (3 × 5 L) over 5 days at room temperature. The methanol extracts were collected and concentrated under reduced pressure to yield a black crude MeOH extract (20 g). This was suspended in hot water (1:1, v/v) and successively partitioned with n-hexane (H), chloroform (C), ethyl acetate (E), and n-butanol (B). The resulting fractions were concentrated under reduced pressure to give the corresponding n-hexane (2.1 g), chloroform (2.9 g), EtOAc (8.8 g), and n-butanol (3.2 g) fractions, which were kept at 4 °C in the dark until further analysis.

The chloroform fraction (2.8 g) was subjected to normal phase silica gel CC eluting with n-hexane/CH₂Cl₂/MeOH (gradient mixtures 7:3:0.5-3:7:0.5) to afford three fractions (CR1-3). Fraction CR1 was repeatedly chromatographed over a silica gel column and eluted with n-hexane/CH₂Cl₂ (8:2 v/v) to obtain compounds 2 [(9(11),12-oleanadien-3β-ol,] and 3 [11,13(18)-oleanadien-3β-ol] (17 mg)). The EtOAc fraction (8.5 g) was further chromatographed by silica gel CC with n-hexane/ EtOAc (8:2, v/v) as eluent to yield compounds 1 (rubinaphthin B, 21.0 mg) and 4 (catechin, 20.0 mg).

The same extraction procedure was applied to the dried leaves of L. inermis. Briefly, the leaf powder (5.0 kg) was extracted with EtOH 85% (3 × 5-10 L) (3 × 24h). The ethanol extracts were collected and concentrated under reduced pressure to yield a black crude extract (31 g), which was then suspended in hot water (1:1, v/v) and successively partitioned with n-hexane (H), chloroform (C), ethyl acetate (E), and n-butanol (B). The resulting fractions were concentrated under reduced pressure to give the corresponding n-hexane (1.4 g), chloroform (7.0 g), EtOAc (5.4 g), and n-butanol (9.2 g) fractions.

The chloroform fraction (7.0 g) was chromatographed over a silica gel column eluted with gradient mixtures (CHCl₃/MeOH [0%-100%-50/50%, 2L/each gradient]) to obtain 4 sub-fractions (CL1-4). Fraction CL4 was placed in a fridge at −4 °C overnight. The yellow precipitate that formed was re-crystallized in cold EtOAc and re-filtered to afford compound 5 (afzelin, 18 mg). Fraction CL1 (1.2 g) was separated by silica gel CC eluting with CHCl₃ /MeOH (10:1) to afford compound 6 (augustic acid, 10.0 mg). Fraction CL2 (200 mg) was purified by RP-18 CC eluting with MeOH/H₂O/acetone (1:1:0.5, v/v/v) to obtain compound 7 (1β, 2α, 3α, 19α-tetrahydroxy-12-ursen-28-oic acid, 11.0 mg). Fraction CL3 (3.1 g) was chromatographed by RP-18 CC and eluted with methanol/water (1/1, v/v) to obtain compound 8 (suavissimoside R1, 13 mg).

The n-hexane fraction H (2.1 g) was chromatographed over silica-gel CC, eluting with n-hexane/EtOAc (6:4-2:6) to obtain compounds 9 (lawson, 5 mg), 11 (1-tridecanol, 5 mg), 12 (1-pentadecanol, 5 mg), and 10 (β-sitosterol, 22 mg). The structures of all isolated compounds were elucidated based on infrared (IR), mass spectrometry (MS) and nuclear magnetic resonance (NMR) spectroscopic data.

Molecular Docking Analysis

Molecular docking was performed with AutoDock v4.2 and AutodockVina v1.2.3.^16,17 The computational software was downloaded from the website http://scripps.edu, operated under Microsoft Windows 8, installed on an Intel i7 PC with a 3.2 GHz processor and 16 GB RAM.

Phytochemicals Optimization

Ligand structures (1-10) were prepared using ChemBioDraw Ultra 12.0 and PUBCHEM. The energies of the prepared ligands were minimized using MM2 in Chem3D Ultra software, while the pdbqt format essential for docking simulation was generated by AutoDock Tools.

Protein Preparation

The main protease of SARS-CoV-2 in the form of three-dimensional crystal structures bound to a co-crystal ligand (M^pro + O6K, in orthorhombic form) were downloaded from the Protein Data Bank (www.rcsb.org) (PDB ID: 6y2g) with a resolution of 2.20 Å.⁵ All crystallized water molecules, heteroatoms, and the bound inhibitor present in the protease crystal structure were removed by BIOVIA Discovery Studio Visualizer 2021¹⁸ and AutoDockTools.¹⁶ The torsions were fixed for the ligand. The rigid grid box was attained using Autogrid. A grid box (21.1 Å × 21.4 Å × 24.3 Å) with coordinates x = −19.9, y = −6.2, z = − 26.8 and exhaustiveness = 400 was set around the active site, where the protein acid amines interacted with its co-crystal O6K (13b¹⁹). This was followed by AutoDock with a Lamarckian Genetic Algorithm to obtain the best docking conformation. The pose with best binding affinity was visualized using AutoDock Tools.⁷ The co-crystal O6K was also used to validate the docking parameters by redocking it to the protein target, and the validated parameters (root mean square deviation [RMSD] < 2 Å) were used for docking all our tested compounds. The formed interactions between the protein and the ligand in the best docked complexes were analyzed using BIOVIA Discovery Studio Visualizer 2021.¹⁸ The docking was performed in two trials, and the docking scores are presented as the average binding energy of two Gibbs free energy of binding (ΔG) values.^7,17 The small variance in ΔG values and binding posed inside the pocket may be attributed to the differences in the position of the functional groups in the selected compounds. The ΔG values were further converted to pKi, pred, where K_i (inhibitory constant) was calculated by the formula K_i = IC₅₀/1 + [S]/K_m = exp(ΔG/R × T), so IC₅₀ = exp (ΔG/R × T) × (1 + [S]/K_m)/K_i (Table 2), where K_m is the Michaelis constant; R is the gas constant, T = 298.15K, and [S] is the substrate concentration (1 M). Using this method, the binding affinities of the ligand-protease were determined and reported in kcal/mol unit.

Table 2.

Calculation of Binding Free Energy Using MM-PBSA and MM-GBSA and Average Values of RMSD (Backbone, Ligands, and C-Alpha), Rg (Radius of Gyration), SASA, and RMSF of M^pro Complexes with Tested Compounds (1-3 and 6-8).

Complexes	ΔG_MM-PBSA _{bind energy} (kcal/mol)	ΔG_MM-GBSA _{bind energy} (kcal/mol)	Average backbone RMSD (nm)	Average ligand RMSD (nm)	Average C-alpha RMSD (nm)	Average radius of gyration (nm)	Average SASA (nm)²	RMSF (nm)
M^pro	-	-	0.140536706	-	0.142376113	2.589052	264.3415101	0.0383-0.2225
M^pro-(1)	−23.8389	−22.3129	0.164272215	0.152515	0.1661877	2.586192	266.5263925	0.0412-0.2325
M^pro-(2)	−23.7150	−24.8406	0.189309436	0.05737689	0.153186	2.597401	264.802156	0.0396-0.7004
M^pro-(3)	−21.3930	−17.9548	0.140214217	0.03967340	0.142428	2.593204	268.820502	0.0422-0.2455
M^pro-(6)	−19.4095	−22.2180	0.140218816	0.05752583	0.142071	2.596201	265.713442	0.0413-0.242
M^pro-(7)	−23.5426	−23.4781	0.147829626	0.07985207	0.149772	2.594259	265.084083	0.2404-0.0385
M^pro-(8)	−27.3538	−29.6151	0.191380753	0.14710293	0.192947	2.605063	267.926171	0.416-0.242
M^pro-Remd	−26.3570	−28.0022	0.137393529	0.28573854	0.13981056	2.5884745	263.4853581	0.0389-0.2807

Molecular Dynamics Simulation

From the results of molecular docking simulations, we conducted molecular dynamics (MD) simulation of the selected compounds with the potential to inhibit SARS-CoV-2 M^pro, including 1-3 and 6-8 to verify their binding stability. Before starting the MD simulation, the post-docking compounds were supplemented with hydrogen using Avogadro software.^20,21 Then, the force field parameters AMBER (generalized AMBER force field [GAFF]) were generated for each molecule under consideration.²² Quantum mechanics calculations at level B3LYP/6-31G(d,p) were performed by GAMESS-US software,²³ followed by limited electrostatic potential energy calculations (restrained electrostatic potential (RESP)) to determine partial charges^24,25 via Antechamber,^22,26 available in AmberTool.²⁷ For proteins, the protonation states of the ionizable amino acid residues of SARS-CoV-2 M^pro had been determined using the H + + program.²⁶ The required GROMACS input files were created with ACPYPE.py scripts.²⁸ We used the Amber99SBILDN force field for protein and ionic parameterization in conjunction with the TIP3P water model to simulate an explicit water molecule. The SARS-CoV-2 M^pro–ligand complex was set at the center of a cubic box that was immersed in an aqueous solvent and the system charge was neutralized with counter ions. The energy of the systems had been minimized through the steepest descent algorithm,²⁹ and the maximum minimization step was 50 000 steps. The systems were then equilibrated with NVT and NPT equilibria. In the NVT ensemble (ie, with a constant number of molecules, volume, and temperature), there was about 0.1 ns. The NPT ensemble was performed at a temperature of 310K (37°C, Nosé-Hoover thermostat^30,31, and a pressure of 1.0 bar (Parrinello-Rahman barostat³²) for a total of 2 ns. The LINCS algorithm³³ restricts bond lengths involving hydrogen atoms (the LINCS algorithm constrained the bond lengths involving hydrogen atoms). To integrate the equations of motion, a leapfrog stochastic dynamics integrator was used with a time step of 2 fs. The coordinates were saved every 10 ps. The MD (MD production) simulation was performed with a time of 50 ns to determine the stability of the combination.

The free binding energy calculation in MD was performed by the gmx_MM-PBSA protocol based on the single trajectory of GROMACS with AMBER force field.³⁴ The binding free energy ΔG_bind of each complex in water was estimated via the commonly used Molecular Mechanics-Poisson − Boltzmann Surface Area (MM-PBSA) and Molecular Mechanics-Generalized Born Surface Area (MM-GBSA) methods.³⁵ A total of 5000 snapshots were the output from the trajectory every 10 ps. The free binding affinities and the RMSD were analyzed. Other parameters, like RMSD score (backbone, ligand, C-α), root-mean-square fluctuations (RMSF), radius of gyration (Rg), and solvent accessible surface area (SASA) calculations, were carried out with GROMACS tools.

ADMET Analysis

The toxicity profiles were verified for each candidate using ProTox (http://tox.charite.de/protox_II/) where data like hepatotoxicity, carcinogenicity, immunotoxicity, mutagenicity, cytotoxicity, toxicity class, and possible lethal dose (LD₅₀) value were calculated.³⁶ Furthermore, the overall drug-likeness characteristics were performed by in silico ADMET study (absorption, distribution, metabolism, excretion, and toxicity) of each compound measured using tools such as Molsoft (https://molsoft.com/mprop/), Admetlab 2.0 (https://admetmesh.scbdd.com/), SwissADME (https://www.swissadme.ch/), and Molinspiration (https://molinspiration.com/) to identify safe and effective drug candidates in the early stage of drug discovery, pharmacokinetic study, and data analysis.

Results and Discussion

Chemical Compositions of L. inermis

From the methanolic extract of L. inermis (roots, leaves, and stems), after fractionation with solvents of different polarities: n-hexane, chloroform, ethyl acetate, and n-butanol, and repeated chromatography on silica gel, dianion HP-20 and reversed phase RP-18, 12 known compounds were isolated (1-12), including five triterpenoids (2, 3, 6, 7, 8), two flavonoids (4-5), two naphthoquinone derivatives (1 and 9), β-sitosterol (10), and two long-chain alcohols (11-12). Their structures were identified by comparison of their melting points, optical rotation values, mass spectrometric, and spectroscopic (NMR, UV, IR, MS) data with those previously published. Specifically, from the chloroform fraction of L. inermis roots, rubinaphthin B (1, 4-naphthohydroquinone-1, 4 di-O-β-D-glucopyranoside) (1),³⁷ 9(11),12-oleanadien-3β-ol (2), 11,13(18)-oleanadien-3β-ol (3)^38,39 and catechin (4)⁴⁰ were isolated. From the chloroform and n-hexane fractions of L. inermis leaves, 3-O-α-L-rhamnopyranosyl kaempferol (afzelin) (5),⁴¹ angustic acid (6),⁴² 1β, 2α, 3α, 19α-tetrahydroxy-12-ursen-28-oic acid (7),^43,44 suavissimoside R1 (8),⁴⁵ 2-hydroxy-1,4-naphthoquinone, commonly named lawsone (9),³⁸ β-sitosterol (10), 1-tridecanol (11), and 1-pentadecanol (12) were isolated (Figure 1). The compounds like naphthoquinone derivative (1) and triterpenoid glycoside (8) were reported as new record, isolated for the first time from L. inermis species.

Figure 1.

The structures of isolated compounds 1-12 elucidated on the basis of spectroscopic and mass spectrometric data (HR-MS, 1D/2D-NMR).

Molecular Docking Analysis

Compounds 1-10 isolated from L. inermis were further analyzed by molecular docking (in silico experiments) with the main protease (non-structural M^pro) from SARS-CoV-2 (COVID-19) to determine their binding energies and molecular interactions with M^pro.

Identification of Active Binding Site

As a first approach to identify the best poses within the binding pocket of M^pro 6Y2G, an active site prediction by FTSite server was performed.⁴⁶ The results showed that three possible active binding sites were found in M^pro 6Y2G (displayed as green, pink, and blue lines in Figure 2A).

Figure 2.

(A) Cartoon representation of 3 proposed active binding sites (green, pink, and blue lines) of M^pro protease 6y2g by FTSite server. (B) The result of re-docking ligand O6K (co-crystal) into the active site region of M^pro to validate docking parameters. The redock ligand (blue layer) overlapped the original ligand (native) (green). (C) The superimposition of isolated compounds (1-10) from L. inermis and the co-crystal inhibitory ligand O6K (13b) into the active binding site (the blue pocket) of M^pro.

Next, in order to validate the docking parameters and to examine the active binding site, we investigated the crystallographic structures of M^pro complexed with inhibitors and redocked the co-crystal O6K with a known conformation and orientation into the target's active site of the protein target M^pro. The results using PyMol (or BDV) software showed that the re-docked ligand was superimposed on the original co-crystallized ligand of the M^pro-O6K complex and right in the blue active site (Figure 2B). The validation parameter RMSD value was obtained as 1.82055 Å, less than 2 Å, which means that the use of the molecular docking simulation program was suitable for this study.⁴⁷ Additionally, the orientation of the crystal ligand to the blue active binding site between domains I and II (Figure 2C) was similar to the published results of Zhang et al.⁵

Binding Parameters

After strengthening the binding site, the 10 compounds were docked to the M^pro protein and the docking simulation results highlighted that the isolated compounds (1-10) posed inside the active pocket with network interactions (hydrogen bond, Pi- bond) (Figure 2C). The docking parameters, including binding energy ΔG (kcal/mol), hydrophobic interaction, π-π stacking, and IC_50,calcd. of the 10 compounds, with remdesivir as positive control, are presented in Table 1 and Figure 3.

Figure 3.

Network of binding interactions of compounds 1-10 with amino acid residues of SARS-CoV-2 M^pro. (Domain I [residues 8-100]; Domain II [residues 101-183]; Domain III [residues 200-302]; long loop [joining DI and DII] [residues 184-199]).

Table 1.

Free Binding Energy (ΔG, kcal/mol) and IC_50,calcd. (μM [Micromolar]) and Network of Binding Interactions Calculated by Molecular Docking.

Compounds (ligands)	Estimated free binding energy, ΔG (kcal/mol)	IC_{50, calcd}. (μM)^a	Hydrogen interaction (R (Å))	π-π stacking
1	−7.30	4.43	Gln A 192 (2.52)	Met A 165
2	−7.73	2.13		Met A 49 Met A 165 Pro A 168 Ala A 191
3	−7.27	4.64	Ser B 1 (3.07) Glu A 166 (2.3)	Cys A 145 His A 41
4	−6.40	20.40	His A 163 (2.48) Ser A 144 (2.60) Cys A 145 (3.03) Gly A 143 (3.26)	His A 41 Met A 165
5	−7.45	3.45	Glu A 166 (2.85) Asn A 142 (2.8) Cys A 145 (3.3)	Met A 165
6	−7.08	6.44	Gln A 189 (3.27)	Met A 165 Met A 49 Cys A 145 His A 41 Leu A 27
7	−7.37	3.96	Glu A 166 (2.18) Gly A 143 (3.09)	Met A 49 Met A 165 His A 41 Cys A 145
8	−8.19	0.98	Gly A 143 (3.05) Thr A 26 (2.84) Glu A 166 (2.64) Phe A 140 (2.02) Cys A 145 (3.41) His A 163 (2.63) His A 164 (3.02)	Met A 49
9	−5.49	94.85	Cys A 145 (3.50) Leu A 141 (2.18) Ser A 144 (2.30)	GLU A 166
10	−6.38	20.78	Thr A 26 (2.0)	Met A 165 Met A 49 His A 41 Cys A 145
O6K (13b)	−6.66	13.08	Cys A 145 (2.49) His A 41 (2.81) Phe A 140 (3.01) HIS A 164 (2.96) SER A 144 (2.97) GLY 143 (3.34) Glu A 166 (2.99)	Met A 49 Met A 165 Pro A 168
Remdesivir	−6.43	19.13	Ser A 144 Cys A 145 His A 41	Met A 165 Met A 49
N3	−6.716		His 41 (2.54)Met A 49 (2.64)Cys A 145 (3.53)	His A 163Met A 165

At temperature = 298.15 K.

The result of the molecular docking analysis indicated that the free binding energies of the compounds to M^pro ranged between − −8.19 and −5.49 kcal/mol (Table 1). The negative values of ΔG indicate that the compounds spontaneously bind to M^pro and, therefore, could have potential inhibitory activities on the M^pro enzyme (Figure 4). The molecules occupied the same cavity, with only some differences in the involvement of amino acid residues. The variance of the ΔG values and the placement of the bond inside the active pocket could be caused by the functional groups presented in the compounds. In Table 1, the docking results of the isolated compounds (ligands) from L. inermis, especially of the five triterpenoids, with the main protease (M^pro) of SARS-CoV-2 were remarkable; the calculated binding energy of 8 < 2 < 7 < 3 < 6 were the lowest, with ΔG values of −8.19, −7.73, −7.37, −7.27, and −7.08 kcal/mol, respectively.

Figure 4.

Hydrophobic interaction and hydrogen bonds between M^pro and the selected triterpenoids (2, 3, 6, 7) OK6 and N3.

Correspondingly, the inhibition concentrations induced by these triterpenoids were with IC_50,calcd. values of 0.98, 2.13, 3.96, 4.64, and 6.44 μM, respectively. In particular, compound 8, suavissimoside R1, a triterpenoid glycoside, showed the most potential inhibitory activity against M^pro protein among the 10 tested compounds, with binding energy ΔG and IC_50,calc. values of −8.19 kcal/mol and 0.98 μM, respectively. These values were even lower than those of the positive compound, remdesivir, with ΔG and IC_50,calc values of −6.43 kcal/mol and 19.13 μM, respectively (Figures 4 and 5; Table 1). Compound 9, a naphthoquinone, namely lawson with the smallest size and molecular weight (of only 174.15 amu) among tested compounds, showed the largest IC₅₀ value of 94.85 μM (Table 3).

Figure 5.

2D and 3D presentations of interactions between compound 8 and acid amino residues inside the M^pro active site.

Table 3.

Pharmacokinetic Parameters and Toxicity Prediction of the Isolated Compounds (1-10) from Lawsonia inermis (Leaves, Stems).

Network Binding Interactions

Additionally, a network of binding interactions of the tested compounds (1-10) with the amino acid residues of SARS-CoV-2 M^pro is displayed in Figure 3. Genome sequencing analysis revealed that M^pro has three domains I, II, and III and a long loop of 15 amino acid joining domains II and III.^8,48 Our studied compounds (1-10) were found to interact with amino acidic residues belonging to both domain I (residues 8-100) and domain II (residues 101-183), including methionine (MetA 49, MetA 165), histidine (HisA 41, 163, 164), cysteine (CysA 145), glutamate (GluA 166), glutamine (GlnA 166, 189), glycine (GlyA 143), and serine (SerA 46), with different interaction distances ranging from 2.0 to 3.4 Å (Figure 3). No amino acid residues of domain III (residues 200-302) were found to bond to the tested compounds. However, the residues ThrA 190, AlaA 191, and GlnA 189 and 192, belonging to the long loop (comprising residues 184-199) connecting to the C-terminal of domain III, formed either H-bonds or π-π interactions with compounds 1, 2, and 6 (Figure 3). According to Anand et al⁴⁹ and Wu et al,⁸ the catalytic site of M^pro employs conserved residues cysteine (Cys145) and histidine (His41), namely CysHis catalytic dyad in it.⁴⁸ Importantly, a deep cleft lined by hydrophobic residues between domains I and II make up the substate-binding site, where the catalytic process takes place.⁴⁹ An in-depth interaction diagram of docking results of the selective triterpenoids and the active site residues, presented in 2D and 3D cartoons using BDV and PyMol tools, is shown in Figure 4. Compound 8 was found to form several hydrogen bonds to amino acid residues of domain I and II like Thr A 26 (in distance of 2.84 Å), Phe A 140 (2.02 Å), Gly A 143 (3.05 Å), Cys A 145 (3.41 Å), His A 163 (2.63 Å), His A 164 (3.02 Å), and Glu A 166 (2.64 Å) inside the catalytic site of M^pro (Figure 5). These data suggested that our test compounds were docked right at the catalytic center of M^pro where they could interact by forming H-bonds and π-bonds with acid amine residues playing key roles in the catalytic enzyme reactions.

MD Simulation

The binding energy was less than 0 kcal/mol, indicating that the ligand and receptor can bind spontaneously and the complex formed between the active compounds and the target might have greater stability. Therefore, based on the results of molecular docking analysis, the compounds 1-3 and 6-8 with the binding energy ΔG values < 7.08 kcal/mol were subjected to MD simulation study to verify and understand more about the conformational stability and dynamics of the protein M^pro–ligand complex.

After MD simulation for a period of 50 ns (Figure 6), the complexes were refined and their overall energy and RMSD score, RMSD of backbone, RMSD ligands, RMSF, RMSD C-α, Rg, SASA, and H-bond were obtained (Figure 6A to H and Table 2).

Figure 6.

Time plot of (A) RMSD backbone, (B) RMSD ligand, (C) 2D interaction of compound 8 with M^pro protein at the end of the simulation time of 50 ns, (D) RMSF, (E) RMSD C-alpha, (F) Rg, (G) SASA, and (H) H-bond of 5 M^pro–inhibitor complexes.

The overall binding free energies (ΔG bind energy) of the M^pro-inhibitor complex systems were calculated using the MM-PBSA and MM-GBSA methods from 500 snapshots and from the whole MD trajectories, in which ligand and protein can assume different conformations.³⁵ The estimated values of binding energy ΔG_MM-PBSA calculated for M^pro-(1), -(2), -(3), -(6), -(7), and -(8) complexes ranged from −19.4 to −27.4 (kcal/mol) (Table 2). The protein–ligand complex of M^pro with compound 8 was the most stable with the lowest ΔG_MM-PBSA binding energy value of −27.4 kcal/mol (Table 2). This value was lower than those of protein M^pro with hypericin (−26.6 kcal/mol), amentoflavone (−22.1 kcal/mol), and glycyrrhizic acid (−5.0 kcal/mol) complexes.⁵⁰ The ΔG_MM-GBSA of M^pro-(8) complex was −29.6 kcal/mol.

The RMSD backbone is the measure of the average distance between the atoms of the tested compounds to the M^pro protein backbone. The average RMSD backbone values of M^pro-(3), M^pro-(6), and M^pro-(7) were around 0.1402-0.1478 nm. Those of M^pro-(1), M^pro-(2), and M^pro-(8) were around 0.1642, 0.1893, and 0.1913 nm, respectively (Figure 6A; Table 2), which were higher than those of M^pro alone (0.1405 nm) and M^pro-remdesivir (0.1374 nm).

During the simulation time of 50 ns, the average RMSD ligand value of the Mpro-inhibitors systems was measured to be in the average range of 0.039-0.147 nm (Table 2). RMSDs of ligands in complexes M^pro-(2), -(3), -(6), and -(7) were stable throughout the simulation. However, the RMSD ligand values of M^pro-(8) fluctuated unstably in the time range of 0-8 ns. From 10 ns onwards, the simulation started to stabilize back to the RMSD ligand value of 0.1650 nm. With an average simulation time of 50 ns, the ligand RMSD value of the M^pro-(8) complex remained 0.1471 nm (Figure 6B).

As the value of RMSF represents the flexibility of conformational residues, the fluctuation in the RMSD backbone and ligand of the M^pro-(8) complex can be explained by the conformation of compound 8 changing during the MD simulation, but the structure of the M^pro-(8) complex was predicted to be stable from the eighth ns of the simulation process (total 50 ns).The interaction sites with amino acid residues also changed over time, but the hydrogen bonds from compound 8 to amino acid residues Glu 166 and Gly 143 remained well-preserved until the end of the 50 ns simulation (Figure 6C) (and supporting information). Its RMSF values were in range of 0.416-0.242 nm (Figure 6D; Table 3).

The RMSD trajectories of the C-α atom were described for the dynamics and convergence of each M^pro-inhibitor system (Figure 6E). It was shown that the C-α RMSD values of M^pro-(1)-(8) were in the range from 0.1420 to 0.1929 nm (Table 2).

Additionally, in order to evaluate the compactness of the M^pro-inhibitor system after MD simulation, we further calculated the radius of gyration (Rg). Figure 6F showed that the Rg of M^pro is observed to be nearly stable in the consistency of the oscillations throughout the simulation. The Rg value of the complex M^pro-(1) and -(8) fluctuated around the mean value of 2.5973 nm (Table 2). This was evidence that there was no significant change in the M^pro protein structure after binding to the tested natural compounds (Figure 6F; Table 2).

The interactions between the M^pro-inhibitor complex system and the aqueous solvent were studied by SASA during the time of MD simulation. This was aimed to investigate the degree of structural change that occurred after the interaction. Figure 6G shows the plot of the SASA value versus time for all Mpro-inhibitor complexes. The mean SASA values gradually decreased from 264.802156 nm² (M^pro-2) to 268.820502 nm² (M^pro-3). In general, the variation in these values was small, indicating that the protein structure after binding to the studied compounds interacting with the solvent is quite stable in the 50 ns simulation period.

The number of intermolecular hydrogen bonds in the ligand-protein complex was determined, which contributes to the stability of the structure. It is shown in Figure 6H that hydrogen bonds are present in all M^pro-inhibitor complexes. The highest number of hydrogen bonds was observed for the M^pro-8 complex over the 50 ns simulation period. Compound 8 also formed the most hydrogen bonds with the M^pro protein, of up to 7 bonds. The other compounds also showed strong interactions with M^pro. As shown in Figure 6H, the number of hydrogen bonds in M^pro-(2), -(3), -(6), and -(7) complexes was 2, 2, 4, and 6, respectively. An increase in residue fluctuations was observed in the M^pro-(8) system, notably among residues 40-70, 130-150, and 160-175, possibly as a result of strong molecular interactions between compound 8 and the protein residues (supporting information). This was similar to the result of molecular docking, where M^pro-(8) formed some hydrogen bonds with several important amino acid residues, such as, Phe A 140, Cys A 145, His A 163, Asn A 142, Glu A 166, and Ser A 46. Hydrogen-binding interactions have been shown to play a key role in stabilizing protein-ligand complexes, and the M^pro-(8) system has been shown to interact strongly with the lowest overall binding free energy among the tested compounds. This molecular docking simulation study suggested that triterpenoid compounds 1-3, 6, 7, and especially triterpenoid glycoside 8 exhibited effective interactions against SARS-CoV-2 M^pro (Table 2).

ADMET Analysis

Additionally, pre-clinical toxicity studies are important to establish the margin of safety and to consider efficiently the risk-benefit of a proposed drug⁷ so that pharmacophoric features influencing the behavior of a molecule in a living organism need to be calculated, including bioavailability (LogP—octanol/water partition coefficient), transport properties (total molecular polar surface area [TPSA]), Caco-2 permeability, enzyme and protease inhibitor potential, toxicity predicted LD₅₀ (mg/kg), and predictive toxicity class of the compounds (Table 3). Generally, the tested compounds were predicted to be in toxicity class from 3 (medium-toxic) to 6 (non-toxic). Compound (4) (catechin) was considered to non-toxic with LD₅₀ value of 10 000 mg/kg. Both compounds (8) and (2) are triterpenoids, but Caco-2 permeability and protease inhibitory potentials of (8) were higher than those of (2) (−6.147 vs −6.052 and 0.06 vs 0.03, respectively). It is thought that compound (8) might be absorbed more easily through intestinal cells and might have more potential in protease inhibition. However, by using a tool for toxicity prediction, triterpenoid (2) (9(11),12-oleanadien-3β-ol) was calculated to be the most toxic among the tested compounds, with a predicted LD₅₀ value of 288 mg/kg. Compound (8) was classified as low-toxic, with an LD₅₀ value of 3220 mg/kg. These results are very interesting and need further study (in vitro, in vivo) to evaluate the factual bioactivity of these natural compounds from L. inermis against SARS-CoV-2 (Table 3).

Conclusion

From the different parts of L. inermis, 12 compounds were isolated, including 5 triterpenoids (2, 3, 6, 7, and 8), 2 flavonoids (4, 5), 2 naphthoquinone derivatives (1 and 9), β-sitosterol (10), and 2 long-chain alcohols (11, 12). Their structures were elucidated by IR and NMR spectroscopic and mass spectrometric methods. Docking analysis found that the triterpenoids had the most significant binding affinity with the M^pro protein of SARS-CoV-2. The triterpenoids (2, 3, 6, 7, 8) showed remarkable binding energies to M^pro with ΔG values even lower than those of the control antiviral compound (remdesivir) (−6.43 kcal/mol) and integrated ligand O6K (−6.66 kcal/mol). By this molecular docking, suavissimoside R1 (8) was shown to be the best compound, binding to SARS-CoV M^pro with binding energy ΔG and IC_50,calc. values of −8.19 kcal/mol and 0.98 μM, respectively. In MD simulation, the protein–ligand complex of M^pro with compound 8 was shown to be the most stable with the lowest ΔG_MM-PBSA and ΔG_MM-GBSA binding energy values of −27.4 and −29.6 kcal/mol, respectively (Table 2). With a simulation time of 50 ns, the average ligand RMSD value of M^pro-(8) complex was stable from the eighth ns (0.1471 nm). The compound also showed the lowest toxicity, with a predicted LD₅₀ value = 3320 mg/kg. These results show that suavissimoside R1 (8) is a potential natural drug against M^pro COVID protease and is worthy of further study to evaluate its actual bioactivity against SARS-CoV-2 in vitro and in vivo. These scientific data contribute to the study of natural resources of Vietnam and can aid in developing a novel drug formulation for either the prevention or treatment of COVID-19.

Supplemental Material

sj-docx-1-npx-10.1177_1934578X221125161 - Supplemental material for Compounds Isolated from Lawsonia inermis L. Collected in Vietnam and Evaluation of Their Potential Activity Against the Main Protease of SARS-CoV-2 Using In silico Molecular Docking and Molecular Dynamic Simulation

Supplemental material, sj-docx-1-npx-10.1177_1934578X221125161 for Compounds Isolated from Lawsonia inermis L. Collected in Vietnam and Evaluation of Their Potential Activity Against the Main Protease of SARS-CoV-2 Using In silico Molecular Docking and Molecular Dynamic Simulation by Nguyen Anh Tuan, Pham Ngoc Khanh, Nguyen Xuan Ha, Ta Chi Binh, Nguyen Duy Khanh and Tran Thị Oanh in Natural Product Communications

Footnotes

Acknowledgements

This work was financial supported by the Ministry of Science and Technology for the Project “Study on production of antifungal drug for treatment of trichophyton, epidermophyton, dermatophytes from Lawsonia inermis L. in Vietnam.” The authors would like to give special thanks to Prof. Dr Nguyen Manh Cuong, INPC, who helped us to finish this paper.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Ministry of Science and Technology (grant number: MOST/19-22).

ORCID iDs

Pham Ngoc Khanh

Nguyen Xuan Ha

Supplemental Material

Supplemental material for this article is available online.

References

Sharmaa

Vijayana

Pantb

, et al. Identification of potential drug candidates to combat COVID-19: a structural study using the main protease (mpro) of SARS-CoV-2. J Biomol Struct Dyn. 2021;37(17):6649-6659.

Worldometer. https://www.worldometers.info/coronavirus/ (2021).

Chen

Liu

Guo

. Emerging coronaviruses: genome structure, replication, and pathogenesis. J Med Virol. 2020;92(4):418-423.

Islam

Parves

Paul

, et al. A molecular modeling approach to identify effective antiviral phytochemicals against the main protease of SARS-CoV-2. J Biomol Struct Dyn. 2021;39(9):3213-3224.

Zhang

Lin

Sun

, et al. Crystal structure of SARS-CoV-2 main protease provides a basis for design of improved α-ketoamide inhibitors. Science 2020;368(6489):409-412.

Verma

Twilley

Esmear

, et al. Anti-SARS-CoV natural products with the potential to inhibit SARS-CoV-2 (COVID-19). Front Pharmacol. 2020;11:561334.

Trott

Olson

. Autodock vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J Comput Chem. 2010;31(2):455-461.

Liu

Yang

, et al. Analysis of therapeutic targets for SARS-CoV-2 and discovery of potential drugs by computational methods. Acta Pharm Sin B. 2020;10(5):766-788.

Gull

Sohail

Aslam

Amin Athar

. Phytochemical, toxicological and antimicrobial evaluation of Lawsonia inermis extracts against clinical isolates of pathogenic bacteria. Ann Clin Microbiol Antimicrob 2013;12:36-36.

10.

Semwal

Combrinck

Cartwright-Jones

Viljoen

. Lawsonia inermis L. (henna): ethnobotanical, phytochemical and pharmacological aspects. J Ethnopharmacol. 2014;155(1):80-103.

11.

Boubaya

Marzougui

Yahia

Ferchichi

. Chemical diversity analysis of Tunisian Lawsonia inermis L. populations. Afr J Biotechnol. 2011;10(25):4980-4987.

12.

Cuong

Nhiem

Thao

, et al. Inhibitors of osteoclastogenesis from Lawsonia inermis leave. Bioorg Med Chem Lett. 2010;20(16):4782-4784.

13.

Ezzat

El Bishbishy

El Kersh

Zayed

Salem

Salama

. Chapter 7 - herbal cosmeticology. In: Egbuna

Mishra

Goyal

, eds. Preparation of Phytopharmaceuticals for the Management of Disorders. Academic Press; 2021: 129-168.

14.

Pasandi

Farahbakhsh

. Lawsonia inermis L. leaves aqueous extract as a natural antioxidant and antibacterial product. Nat Prod Res. 2020;34(23):3399-3403.

15.

Suleiman

Mohamed

. In Vitro activity of Lawsonia inermis (Henna) on some pathogenic fungi. Hindawi J. Mycol. 2014;2014:1-5.

16.

Morris

Huey

Lindstrom

Sanner

Belew

Goodsell

. Autodock 4 and AutoDockTools 4: automated docking with selective receptor flexibility. J Comput Chem. 2009;30:2785-2791.

17.

Eberhardt

Santos-Martins

Tillack

Forli

. Autodock vina 1.2.0: new docking methods, expanded force field, and python bindings. J Chem Inform Modell. 2021;61(8):3891-3898.

18.

BIOVIA Discovery Studio Visualizer. v21.1.0.20298. Dassault Systèmes, 2021.

19.

RCSB. (2021) https://www.rcsb.org/ligand/O6K.

20.

Avogadro. Version 1.2.0 Avogadro: an open-source molecular builder and visualization tool. 2022. http://avogadro.cc/

21.

Hanwell

Curtis

Lonie

Vandermeersch

Zurek

Hutchison

. Avogadro: an advanced semantic chemical editor, visualization, and analysis platform. J Cheminform. 2012;4:1-17.

22.

Wang

Wolf

Caldwell

Kollman

Case

. Development and testing of a general amber force field. J Comput Chem. 2004;25(9):1157-1174.

23.

Barca

Bertoni

Carrington

, et al. Recent developments in the general atomic and molecular electronic structure system. J Chem Phys. 2020;152(15):154102.

24.

Bayly

Cieplak

Cornell

Kollman

. A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model. J Chem Phys. 1993;97(40):10269-10280.

25.

Cornell

Cieplak

Bayly

Kollman

. Application of RESP charges to calculate conformational energies, hydrogen bond energies, and free energies of solvation. JACS 1993;115(21):9620-9631.

26.

Wang

Kollman

Case

. Automatic atom type and bond type perception in molecular mechanical calculations. J Mol Graph Model. 2006;25(2):247-260.

27.

Case

Ben-Shalom

Brozell

, et al. Amber 2021. University of California; 2021.

28.

Sousa da Silva

Vranken

. ACPYPE – AnteChamber PYthon parser interfacE. BMC Res Notes. 2012;5:367.

29.

Meza

. Steepest descent. Wiley Interdiscip Rev: Comput Stat. 2010;2(6):719-722.

30.

Nosé

. A molecular dynamics method for simulations in the canonical ensemble. Mol Phys. 1984;52(2):255-268.

31.

Hoover

. Canonical dynamics: equilibrium phase-space distributions. Phys Rev A Gen Phys. 1985;31(3):1695-1697.

32.

Parrinello

Rahman

. Polymorphic transitions in single crystals: a new molecular dynamics method. J Apppl Phys. 1981;52:7182-7190.

33.

Hess

Bekker

Berendsen

HJC

Fraaije

JGEM

. LINCS: a linear constraint solver for molecular simulations. J Comput Chem. 1997;18:1463-1472.

34.

Valdés-Tresanco

Valiente

Moreno

. Gmx_MMPBSA: a new tool to perform end-state free energy calculations with GROMACS. J Chem Theory Comput. 2021;17(10):6281-6291.

35.

Poli

Granchi

Rizzolio

Tuccinardi

. Application of MM-PBSA methods in virtual screening. Molecules 2020;25(8):1971-1990.

36.

Banerjee

Eckert

Schrey

Preissner

. ProTox-II: a webserver for the prediction of toxicity of chemicals. Nucleic Acids Res. 2018;46(W1):W257-W263.

37.

Liou

. Constituents of the roots of Rubia yunnanensis. Chem Pharm Bull (Tokyo) 2002;50(2):276-279.

38.

Tanaka

Matsunaga

. Triterpene dienols and other constituents from the bark of Phyllanthus flexuosus. Phytochemistry 1988;27(7):2273-2277.

39.

Tanaka

Tabuse

Matsunaga

. Triterpenes from the stem bark of Phyllanthus flexuosus. Phytochemistry 27(11):3563-3567.

40.

Kim

Byun

Bandi

AKR

Hyun

Lee

. Compounds with elastase inhibition and free radical scavenging activities from Callistemon lanceolatus. J Med Plants Res. 2009;3(11):914-920.

41.

Agrawal

. ¹³C-NMR of Flavonoids. Elservier Science Publishers; 1989.

42.

Alam

Chopra

Ali

Niwa

. Oleanen and stigmasterol derivatives from Ambroma augusta. Phytochemistry 1996;41:1197-1200.

43.

Khan

Ali

Singh

Dhawan

. Occurrence and identification of a new antiviral saponin from Lawsonia alba lam. Nat Acad Sci Lett (India) 1996;19(7&8):145-148.

44.

Hewawasam

Jayatilaka

KAPW

Pathirana

Mudduwa

LKB

. Hepatoprotective effect of Epaltes divaricata extract on carbon tetrachloride induced hepatotoxicity in mice. Indian J Med Res. 2004;120(1):30-34.

45.

Abe

Yamauchi

. Trachelosperosides, glycosides of 19-alpha-hydroxyursane-type triterpenoids from Trachelospermum asiaticum (Trachelospermum. III). Chem Pharm Bull. 1987;35(5):1748-1754.

46.

Ngan

C-H

Hall

Zerbe

Grove

Kozakov

Vajda

. FTSite: high accuracy detection of ligand binding sites on unbound protein structures. Bioinformatics 2011;28(2):286-287.

47.

Cole

Murray

Nissink

Taylor

. Comparing protein-ligand docking programs is difficult. Proteins 2005;60(3):325-332.

48.

Isa

Mustapha

Qazi

, et al. In silico molecular docking and molecular dynamic simulation of potential inhibitors of 3C-like main proteinase (3CLpro) from severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) using selected African medicinal plants. Adv Trad Med. 2020;22:1-17.

49.

Anand

Palm

Mesters

Siddell

Ziebuhr

Hilgenfeld

. Structure of coronavirus main proteinase reveals combination of a chymotrypsin fold with an extra alpha-helical domain. EMBO J. 2002;21(13):3213-3224.

50.

Saravanan

Zhang

Senthil

, et al. Structural basis for the inhibition of SARS-CoV2 main protease by Indian medicinal plant-derived antiviral compounds. J Biomol Struct Dyn. 2022;40(5):1970-1978.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.03 MB