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, 3CL pro (also called M pro ), 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 3CL pro ), do show significant homology. Thus, sequence alignment of the SARS-CoV 3CL pro 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 . The crystal structure of the 3CL pro protein has been determined Hsu et al., 2005; Lee et al., 2005) , revealing that it is similar to those of other CoV proteinases. SARS-CoV 3CL pro forms a homodimer, comprised of three domains. Domains I (residues 8-101) and II (residues 102-184) are b-barrels and together resemble the structure of chymotrypsin, respectively, whereas domain III (residues 201-306) consists mainly of a-helices.
Domains II and III are connected by a long loop (residues 185-200). The active site of the M pro 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 3CL pro 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.
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.
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 3CL pro 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 3CL pro (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 welldefined electron density.
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 3CL pro 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 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 3CL pro was obtained from X-ray crystallographic structure (PDB ID:1UK3) . 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.
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.
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 3CL pro pocket, i.e., initial (docking) conformation of the RTV is totally different from that after equilibration.
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 3CL pro 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 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 3CL pro shows different dynamic flexibilities in the binding region when the SARS-selective inhibitors are bound to their active sites. In the SARS-CoV 3CL pro 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 3CL pro and its two inhibitor complexes are between 8.64 and 10.02 Å (Fig. 4b) .
This could be due to the fact that the SARS-CoV 3CL pro free form was initially generated from the complexed SARS-CoV 3CL pro 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 b) . 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.
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 3CL pro , 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 ).
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 À2001 to 01. 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 À571 decreases to À1471 for the LPV complex (dashed line in Fig. 5b) .
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 .
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 p3.5 Å , and (ii) a donor-H-acceptor bond angle of X1201. 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 (DE MM ), 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.
The average binding free energies (DG binding ) for the two complexes are shown in Table 3 . In this study, the contribution of the entropic term (TDS) was excluded due to the reasons as in our previous works that the two system models are very similar and, therefore, the contribution of TDS to DG binding in terms of relative scale is negligible.
From the results, the total electrostatic contribution term, DG ele sol þ DE ele , 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 3CL pro enzyme is in accordance with the fact that the LPV inhibitor was developed from RTV, so the structures of both inhibitions are quite similar.
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 ij (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 .
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.
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 3CL pro enzyme. 