Molecular dynamic simulations analysis of ritonavir and lopinavir as SARS-CoV 3CL(pro) inhibitors.

Since the emergence of the severe acute respiratory syndrome (SARS) to date, neither an effective antiviral drug nor a vaccine against SARS is available. However, it was found that a mixture of two HIV-1 proteinase inhibitors, lopinavir and ritonavir, exhibited some signs of effectiveness against the SARS virus. To understand the fine details of the molecular interactions between these proteinase inhibitors and the SARS virus via complexation, molecular dynamics simulations were carried out for the SARS-CoV 3CL(pro) free enzyme (free SARS) and its complexes with lopinavir (SARS-LPV) and ritonavir (SARS-RTV). The results show that flap closing was clearly observed when the inhibitors bind to the active site of SARS-CoV 3CL(pro). The binding affinities of LPV and RTV to SARS-CoV 3CL(pro) do not show any significant difference. In addition, six hydrogen bonds were detected in the SARS-LPV system, while seven hydrogen bonds were found in SARS-RTV complex.

Introduction

Reported first in the Guangdon province of China, in late 2002, as an atypical pneumonia, severe acute respiratory syndrome (SARS) is of concern, since it subsequently rapidly spread to over 25 countries and displayed an estimated mortality rate of 3–6% in infected humans, leading to more than 700 human deaths (Donnelly et al., 2003). In March 2003, four months after the first case was reported, the causative agent of SARS was identified as a new coronavirus (CoV) (Kuiken et al., 2003; Peiris et al., 2003). The complete genomic sequence (29,751 nucleotides) of the +ssRNA virus was reported a month later (Marra et al., 2003), providing the opportunity for more in-depth studies and the identification of potential targets for drug and vaccine development. The genome of SARS-CoV contains 11–14 major annotated open reading frames including those predicted to encode for known viral proteins, including the replicase polyproteins, S (spike protein), polymerase, M (membrane protein), N (nucleocapsid protein) and E (small envelope protein). The viral proteinase and replicase are the preferred targets for searching and designing of antiviral compounds. Indeed, the SARS-CoV main proteinase, 3CLpro (also called Mpro), exhibits a key role in proteolytic processing of the replicase polyproteins and, as an essential protein for viral replication and function, is viewed as a key target.

Whilst phylogenetic analysis of the SARS genome reveals little resemblance to any of the three known groups of cornoaviruses, several of the specific viral genes, including the main protease gene (SARS-CoV 3CLpro), do show significant homology. Thus, sequence alignment of the SARS-CoV 3CLpro with orthologues from other coronaviruses indicated that the enzyme is highly conserved, with 40% and 44% sequence identity, respectively, to human CoV (HCoV) 229E proteinase and procine transmissible gastroenteritis virus (TGEV) proteinase (Anand et al., 2003). The crystal structure of the 3CLpro protein has been determined (Yang et al., 2003; Hsu et al., 2005; Lee et al., 2005), revealing that it is similar to those of other CoV proteinases. SARS-CoV 3CLpro forms a homodimer, comprised of three domains. Domains I (residues 8–101) and II (residues 102–184) are β-barrels and together resemble the structure of chymotrypsin, respectively, whereas domain III (residues 201–306) consists mainly of α-helices. Domains II and III are connected by a long loop (residues 185–200). The active site of the Mpro is comprised of a catalytic dyad that consists of the conserved residues H41 and C145, which are located at the cleft between domains I and II (Huang et al., 2004). The availability of protein structures and the known biological characteristics of 3CLpro both encourage and help facilitate it to be a key target for searching for drugs directed against SARS.

Although neither an effective antiviral drug nor a vaccine against SARS is currently available, several reports have indicated that HIV-1 protease inhibitors have potential for designing SARS-CoV proteinase inhibitors (Zhang and Yap, 2004; Yamamoto et al., 2004; Jenwitheesuk and Samudrala, 2003). In particular, the proteinase inhibitor Kaletra, which is a mixture of the proteinase inhibitors lopinavir and ritonavir (Fig. 1 ), exhibits encouraging signs of being partially effective against the SARS virus (Vastag, 2003). Therefore, a deeper understanding of known proteinase inhibitors as anti-SARS drugs can be used as a starting point to design and discover better or more specific inhibitors for the treatment of SARS, as well as any mutations by predicting resistance mechanisms and/or new drugs for such mechanisms in advance.

Fig. 1

Chemical formula of HIV-1 PR inhibitor, ritonavir and lopinavir, which exhibits signs of effectives against SARS.

In the present research, the evidence for the potential benefits of HIV-1 protease inhibitors against the SARS coronavirus (SARS-CoV) is provided. Molecular dynamic simulations were performed in order to investigate the structure and dynamics behaviours of the free enzyme and its complexes with the lopinavir and ritonavir inhibitors.

Methods

Starting structure and protein preparation

In order to understand ligand–enzyme interactions as well as the effects of inhibitors on the structure and dynamic properties of the free enzyme, three molecular dynamics (MD) simulations were performed: SARS-CoV 3CLpro free enzyme (free SARS) and its two complexes with the lopinavir (SARS–LPV) and ritonavir (SARS–RTV) inhibitors. The X-ray structure of the SARS-CoV 3CLpro (PDB code 1UK3; 2.4 Å resolution) was used as the initial model to construct the SARS–RTV and SARS–LPV complexes (Fig. 2 ). Prior to the simulations, residues SerA1 and GlyA2, which had not been visible in the electron density maps, were modeled using the LEaP module in the AMBER 7 software package (Case et al., 2002). The modeled conformation at the beginning of each simulation was similar to that of residues SerB1 and GlyB2 of the other polypeptide chain in the dimer, which exhibits a well-defined electron density.

Fig. 2

Structures of (a) SARS-free enzyme, (b) SARS–lopinavir complex and (c) SARS–ritronavir complex.

Flexible docking of SARS and inhibitors

The flexible molecular docking (Muegge and Martin, 1999) with Genetic Algorithm was performed using BioMedCache 2.0 Software to find the most favorable binding interaction. Two drugs, ritonavir and lopinavir, were separately docked into the SARS-CoV 3CLpro binding site. The residues associating in the binding pocket were defined by the conserved amino acid residues, namely H41 and C145, which are also known to be part of the catalytic site of the SARS-CoV proteinase (Anand et al., 2002, Anand et al., 2003; Yang et al., 2003). Additionally, their neighbor residues within a radius of 3 Å of these conserved amino acids were selected and considered as the members of the binding pocket.

The structure of SARS-CoV 3CLpro was obtained from X-ray crystallographic structure (PDB ID:1UK3) (Anand et al., 2003). All missing hydrogen atoms were added into the enzyme using tleap in the AMBER 7 simulation package. The geometries of the inhibitors were generated from X-ray structured and all hydrogen atoms were added using Gaussian03, with the use of antechamber to assign atom types. The geometry of inhibitors was optimized by the semiempirical AM1 method. The PMF scores of the drugs were evaluated by the genetic algorithm with a population size of 50, crossover rate of 0.80, elitism of 5, mutation rate of 0.2, with the maximum generation cycle set to be 40,000. The size of the grid box was set at 25 Å×25 Å×25 Å. Since, as expected, no significant difference was found in terms of their orientations, among the complex structures within the low-energy clusters, the conformations with the lowest PMF energy score were used as the initial model for MD simulations. Finally, the complex structures with the lowest PMF energy score were constructed and used as an initial model for MD simulations.

Molecular dynamics simulations

Three MD simulations (free SARS, SARS–RTV and SARS–LPV) were carried out using the AMBER 7 simulation package. Energy minimization and MD simulations were achieved using the SANDER module of AMBER 7 with the Cornell force field (Cornell et al., 1995). The parm99 force field, which describes the structure and conformation of each compound, was applied. All MD runs reported here were achieved under an isobaric–isothermal ensemble (NPT) using a constant pressure of 1 atm and a constant temperature of 298 K. The SHAKE algorithm (Ryckaert et al., 1977) was used to constrain all bonds involving hydrogen atoms. Energy minimizations were carried out to relax the system prior to MD runs. To incorporate the solvent and counter ions into the models, each system was neutralized with the counter ions (Na+) and solvated with the TIP3P water model (Jorgensen et al., 1983). The total number of atoms for the three simulated systems, free SARS, SARS–LPV and SARS–RTV, were 77,866, 77,966 and 86,788, respectively. The crystallographic waters were also kept in the simulation box. The simulation time step was set at 2 femtosecond (fs). The simulation steps consist of thermalization, equilibration, and production phases. Initially, the temperature of the system was gradually raised to 298 K for the first 60 picosecond (ps) and then kept constant according to the Berendsen algorithm with a coupling time of 0.2 ps. The systems were maintained at 298 K, until 600 ps. Finally, the production phase started from 600 ps to 2 ns. The 1400 ps trajectories of the production phase were used to calculate the average structure. All MD simulations were carried out for 2 ns. The quality of the geometry and the stereochemistry of the protein structures were validated using PROCHECK. MD trajectories were evaluated in terms of the root-mean-square displacement (RMSD), changes of the distances, torsion angles, H–bond and binding free energies using the CARNAL, Ptraj and MM-PBSA modules of the AMBER7 simulation package.

Results and discussion

Overall enzyme and inhibitor structure

The 2 ns MD trajectories of the three simulated systems were generated. The RMSDs of the overall structure and two inhibitors with respect to the initial configuration were analyzed and plotted (Fig. 3 ). The overall RMSDs for the three systems appeared to have reached equilibrium after 600 ps. Whilst no significant difference was found for the overall (black) and for the back bone (dark grey) RMSDs of all plots, that of the RTV in the SARS–RTV complex (light grey) was found to change substantially. This is due to the rotation of the RTV in the SARS-CoV 3CLpro pocket, i.e., initial (docking) conformation of the RTV is totally different from that after equilibration.

Fig. 3

RMSDs of SARS-CoV 3CLpro free enzyme ((a) SARS) and its complexes with lopinavir ((b) SARS–LPV) and ritonavir ((c) SARS–RTV).

Flexibility of the inhibitor-binding loop regions

Protein flexibility is essential to understand the ligand-binding affinity (Carlson, 2002; Carlson and McCammon, 2000; Lin et al., 2002; Teague, 2003). Conformational changes of the catalytic binding domain play a critical role in ligand binding. Here, the ligand-binding site of the SARS-CoV 3CLpro was evaluated focusing on the amino acid sequence around the active sites (H41 and C145, see Fig. 2) of the CoV proteinase (Anand et al., 2002, Anand et al., 2003; Yang et al., 2003). The binding of the inhibitor is likely to affect the relative motion between the apex of the loop with respect to the D48 and Q189 loops upon the insertion of inhibitor. The probability distributions of the associated interatomic distances from the center of mass of the following residues for the three systems, d1:H41-C145 and d2:D48-Q189, were evaluated and plotted (Fig. 4 ). Clearly, SARS-CoV 3CLpro shows different dynamic flexibilities in the binding region when the SARS-selective inhibitors are bound to their active sites. In the SARS-CoV 3CLpro free enzyme, the most probable d1 distance between H41 and C145 takes place at about 10.0 Å, but this is contracted to ∼7.1–7.3 Å in the enzyme–inhibitor complexes. For the distance between D48 and Q189, d2, in the loop region, no significant differences were observed between free enzyme and two enzyme–inhibitor complexes. Rather the optimal distances, represented by the maximum of the plots, for both SARS-CoV 3CLpro and its two inhibitor complexes are between 8.64 and 10.02 Å (Fig. 4b).

Fig. 4

The probability distributions of the associated interatomic distances between the distances (a) d1 and (b) d2 from the center of mass of the following residues, H41 to C145 and D48 to Q189, respectively.

This could be due to the fact that the SARS-CoV 3CLpro free form was initially generated from the complexed SARS-CoV 3CLpro structure, which was identified as having a closed loop conformation. Similarly, this kind of flap (loop) closing was also observed when inhibitors bind in the active site of HIV-1 protease (Hornak et al., 2006a, Hornak et al., 2006b). In addition, the SARS–RTV distribution plots are narrower and higher than those of the SARS–LPV, indicating a lower flexibility of the catalytic pocket and the flap loop of the SARS–RTV complex compared to that of the SARS–LPV complex. This fact relates directly to the greater flexibility of the longer side chain of the RTV than that of the LPV.

Conformations of the catalytic residues

Another structural feature common to enzyme–inhibitor complexes in general is the conformation of the enzyme-binding site. To investigate the enzyme-binding site conformation of SARS-CoV 3CLpro, we analyzed the changes in the four dihedral angles of the key interactions of the enzyme, involving the amino acids at the catalytic residues and the tip of the loop upon inhibitor insertion. In this study, we present the conformational analysis of the two torsion angles of catalytic residues (Tor1: CA-CB-CG-CD of H41 and Tor2: N-CA-CB-S of C145) and the two amino acids at the tip of the loop (Tor3: CA-CB-CG-OD of D48 and Tor4: CB-CG-CD-CA of Q189), during 600–2000 ps of the simulations (summarized in Fig. 5 and Table 1 ).

Fig. 5

(a) Definitions and (b) probability distributions of the torsion angles of the two catalytic residues (Tor1: CA-CB-CG-CD of H41 and TorR2: N-CA-CB-S of C145) and two amino acids at the tip of the loop (Tor3: CA-CB-CG-OD of D48 and Tor4: CB-CG-CD-CA of Q189) during 600–2000 ps of the simulations.

Table 1

The four dihedral angles (°), Tor1: CA-CB-CG-CD of H41, Tor2: N-CA-CB-S of C145, TorR3: CA-CB-CG-OD of D48, and TorR4: CB-CG-CD-CA of Q189 taken from the maxima of the distribution plots shown in Fig. 5b

	Tor1	Tor2	Tor3	Tor4
SARS	−77	−52, 52	−57, 52	−137, 132 (broad)
SARS-LPV	−77, 187	−57 (sharp)	−147, 37	132 (broad)
SARS-RTV	77 (sharp)	−62 (sharp)	127	82

In free form, Tor1 generates a sharp peak indicating rigidity of the H41 side chain. However, this is not the case for the Q189 side chain, where Tor4 shows a broad peak ranging from −200° to 0°. In contrast, C145 and D48 residues generate two peaks (Tor2 and Tor3), representing two preferential conformations. Changes of their conformations due to complex formations were clearly observed for all four residues. For H41, C145 and Q189, where their orientations were represented by Tor1, Tor2 and Tor4, respectively, their side chains were rotated via complex formations, to the same direction for both inhibitors. This is in contrast to D48, where its side chain was induced by RTV and LPV to turn to the opposite side. Tor3 increases for the RTV complex (dotted line in Fig. 5b), whilst only the first peak of −57° decreases to −147° for the LPV complex (dashed line in Fig. 5b).

Inhibitor–enzyme interaction

The SARS–LPV and SARS–RTV interactions were analyzed using hydrogen bond interactions and the molecular mechanical energy obtained from the electrostatic and the van der Waals interactions within the systems. The results are summarized in Table 2 .

Table 2

Percentage occupation of hydrogen bonding between inhibitor and specific residues of SAR-CoV enzyme

Residue	Type	%Occupation
		SARS-LPV	SARS-RTV
T26	OG1_HG1---O3_LPV	18.7	–
T26	LPV_N4_H47---OG1	22.8	–
T26	LPV_N4_H47---O	36.9	–
H41	RTV_N2_H5---NE2	–	1.7
Y118	OH_HH---O5_LPV	2.8	–
N119	LPV_N4_H47---OD1	3.3	–
C145	RTV_N2_H5---SG	–	56.3
C145	SG_HG---O3_RTV	–	4.4
E166	N_H---N1_RTV	–	29.9
Q189	NE2_HE21---O1_LPV	9.9	–
Q189	NE2_HE21---N1_RTV	–	11.6
Q189	NE2_HE21---S1_RTV	–	7.6
Q189	NE2_HE21---O2_RTV	–	1.3

To analyze the hydrogen bond interactions, the percentage occupancy and number of hydrogen bonds between the inhibitor and binding pocket residues were determined based on the following criteria: (i) a proton donor–acceptor distance of ⩽3.5 Å, and (ii) a donor–H–acceptor bond angle of ⩾120°. It was found that the numbers of hydrogen bonds for the two systems are almost comparable (Table 2). Six hydrogen bonds were detected in the SARS–LPV system, whilst seven were found in the SARS–RTV complex. However, all the hydrogen bonds show a rather low % occupation. Interestingly, excluding Q189, the two inhibitors were found to form hydrogen bonds to different residues, RTV binds to H41, C145 and E166, while LPV binds to T26, Y118 and N119.

In terms of the gas-phase MM interaction energy (ΔE), summarized in Table 3 , the molecular mechanical interaction energy for the SARS–LPV complex (−52.9 kcal/mol) is lower than that of the SARS–RTV complex (−48.4 kcal/mol). Taking into account the hydrophobic and electrostatic interactions for the MM interaction energy, a higher contribution from van der Waals energy (−48.1 kcal/mol) and a lower contribution from electrostatic energy (−4.8 kcal/mol) were observed for the SARS–LPV complex than that for the SARS–RTV complex.

Table 3

Averaged energy contributions (kcal/mol) of the binding free energy for the two simulated systems, SARS-LPV and SARS-RTV

	SARS-LPV	SARS-RTV
ΔEelea	−4.8±4.2	−8.8±3.3
ΔEvdwb	−48.1±3.3	−39.6±2.9
ΔEMMc	−52.9±5.4	−48.4±4.6
ΔGsolnonpolar d	−3.8±0.3	−2.9±0.2
ΔGsolele e	9.6±1.7	6.0±0.8
ΔGsolf	5.9±1.4	3.1±0.8
ΔGsolele+ΔE_ele	4.9±3.4	−2.8±3.6
ΔGbindingg	−47.2±5.3	−45.3±4.3

Electrostatic energy.

van der Waals energy.

Electrostatic + van der Waals energy.

Nonpolar contribution to salvation.

Electrostatic contribution to salvation.

Total solvation energy.

Binding free energy in the absence of entropic contribution.

Inhibitor–enzyme binding free energy

The average binding free energies (ΔG) for the two complexes are shown in Table 3. In this study, the contribution of the entropic term (TΔS) was excluded due to the reasons as in our previous works that the two system models are very similar and, therefore, the contribution of TΔS to ΔG in terms of relative scale is negligible.

From the results, the total electrostatic contribution term, , of SARS–RTV (−2.8 kcal/mol) was found to be more favorable than that of SARS–LPV (4.9 kcal/mol). This can be understood by the higher polarity of the RTV side chain than that of the LPV. However, the van der Waals interactions and nonpolar solvation terms of the SARS–RTV complex are unfavorable by 4.3 kcal/mol relative to those of the SARS–LPV complex. Summing the above contributions, the total binding free energy of the SARS–LPV complex is only slightly more favorable than that of the SARS–RTV complex at −47.2 and −45.3 kcal/mol, respectively. This close similarity, with only −1.9 kcal/mol difference, is due to structural similarity of both complexes. Indeed, the inherent suggestion that the binding affinities of LPV and RTV are not significantly different for the SARS-CoV 3CLpro enzyme is in accordance with the fact that the LPV inhibitor was developed from RTV, so the structures of both inhibitions are quite similar.

Solvent-accessible region

Ligand solvation, especially in the catalytic pocket, is known to play an essential role in the catalytic mechanism of the enzyme. To ascertain basic information on ligand solvation, the radial distribution function (RDF)—g(r): the probability of finding a particle of type i in a spherical radii, r, around the particle of type j—centered on the selected atoms of the inhibitors was calculated. The evaluation was done separately for atom sets of both inhibitors, which lay either (i) in equivalent (O2, N2, O3, N3 and N4) or (ii) in different (O1, N4, O5, N5 and N6) environments (see Fig. 1 for atomic label). The results are summarized in Fig. 6 .

Fig. 6

Radial distribution functions, g(r), centered on the inhibitor atoms (see Fig. 1 for definitions) to oxygen atoms of modeled water of the two complexes, SARS–LPV and SARS–RTV. The plots are classified into two categories where the central atom types for the two inhibitors are topologically (a) equivalent and (b) different.

For both of the two complexes, at the central part of the enzyme–inhibitor complex cavity where the atoms of the two inhibitors are topologically equivalent (Fig. 6a), the plots for O2, O3 and O4 show a pronounced first peak at about 3 Å. This indicates that these three atoms were tightly solvated by water molecules. In addition, the three plots for the SARS–RTV complex (dashed line) all show lower minima, indicating that the residence time for the first shell solvation of the SARS–RTV complex is longer than that of the SARS–LPV one. In contrast, the N2 and N3 atoms of both inhibitors are solvated to a significantly lower extent, implying a high hydrophilicity around these regions.

For the outer part (O1, N4, O5, N5 and N6) of the inhibitors (see also Fig. 1), where the RDFs are given in Fig. 6b, O1 of the SARS–RTV complex (dashed line) and O5 of the SARS–LPV complex (solid line) were firmly solvated. Excluding these two regions, the two ends of the two inhibitors are much less solvated, the first peaks are broad and locate far from the referent atoms, indicating a highly hydrophobic characteristic of these two drug molecules.

Conclusion

MD simulations and analysis of free SARS, and the SARS–RTV and SARS–LPV complexes give an insight into the dynamic characteristics in terms of the flexibility, conformation, and the inhibitor-enzyme interactions. The analysis of the energetic binding affinity between the HIV-1 proteinase inhibitors has been extensively performed based on free energy calculations. We have found that the interaction of LPV and RTV does not show a significant difference for the SARS-CoV 3CLpro enzyme.