Effect of pH on stability of dimer structure of the main protease of coronavirus-2

The viral main protease (Mpro) from a novel severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) is a key enzyme essential for viral replication and has become an attractive target for antiviral drug development. The Mpro forms a functional dimer and exhibits a pH-dependent enzyme activity and dimerization. Here, we report a molecular dynamics (MD) investigation to gain insights into the structural stability of the enzyme dimer at neutral and acidic pH. Our data shows larger changes in structure of the protein with the acidic pH than that with the neutral pH. Structural analysis of MD trajectories reveals a substantial increase in intersubunit separation, the loss of domain contacts, binding free energy and interaction energy of the dimer which implies the protein instability and tendency of dimer dissociation at acidic pH. The loss in the interaction energy is mainly driven by electrostatic interactions. We have identified the intersubunit hydrogen-bonding residues involved in the decreased dimer stability. These findings may be helpful for rational drug design and target evaluation against COVID-19.


Introduction
Coronaviruses (CoVs) infecting humans were first identified in 1947 [1]. CoVs are enveloped viruses consisting of single positive-strand RNA. They infect animals and human such as bats, pets, livestock, poultry, and humans. In humans, CoVs can cause not only severe respiratory disorders but also gastrointestinal, and neurological damages. CoVs belong to subfamily Coronavirinae of the family Coronaviridae. The Coronavirinae is subdivided into four genera alpha-, beta-, gamma-and deltacoronavirus (α-, β-, γ-, and δ-CoV) [2]; each genus is further divided into four lineage subgroups.
The first pandemic of severe acute respiratory syndrome from coronavirus (SARS-CoV) appeared in 2002 and ended in early 2004. The novel severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) named as COVID-19, emerged in late 2019 and rapidly caused a global pandemic. The outbreak of COVID-19 has now become a global threat and detrimentally affected public health care emergency around the world. The SARS-CoV-2 has a ~82% genome similarity to that of the SARS-CoV [3]. The SARS-CoV-2 belongs to the β-CoV group. The SARS-CoV-2 is the 7th human coronavirus (HCoV) discovered until now [4]. The other six HCoVs include HCoV-229, HCoV-OC43, HCoV-NL63, HCoV-HKU1, SARS-CoV and Middle East respiratory syndrome coronavirus (MERS-CoV) [5,6]. The SARS-CoV-2, SARS-CoV, HCoV-OC43, HCoV-HKU1 and MERS-CoV belong to the genus β-CoV while the HCoVs, HCoV-229E and HCoV-NL63, belong to the genus α-CoV [6]. The β-CoV is well known because it is highly lethal to human. The human CoVs is genetically closely related to bat CoVs such as tylonycteris bat coronavirus HKU4 (Bat CoV-HKU4) [7]. The origin of the virus is still a subject of debate. Analysis of the whole-genome sequence similarity indicates bats, a high probability of the original career, which are speculated to be primordial hosts for most CoV lineages [7][8][9][10].
There is an urgent need to explore antiviral targets and effective drugs for the COVID-19. Identification of potential protein targets and biological mechanism of the viral pathogenic process is important for drug development [11][12][13][14]. Among various protein targets of SARS-CoV-2, the main protease (M pro ) or 3-chymotrypsin (C)-like cysteine protease (3CL pro ) is an attractive therapeutic target due to its key enzyme in the replication cycle of the virus [15][16][17]. The M pro is functional as dimer [3,17,18]. The enzyme is highly conserved across various CoVs but mutation in M pro appears lethal to the virus [3] which makes M pro as a promising target for the development of a broadspectrum of antivirals against coronaviruses. Many potent inhibitors for the SARS-CoV-2 M pro have been proposed by in silico approaches [19][20][21][22]. A significant number of crystal structures of coronavirus M pro determined by X-ray diffraction provided mechanistic insights into protein function and facilitated drug development. From 3D structure of the SARS-CoV-2 M pro , the two subunits are arranged almost perpendicular to one another as shown in Fig. 1A. Each subunit comprises the N-terminus (NT) and three domains (DM I, II and III) [3,17,23,24]. The NT or "N-finger" consists of residues ~1-9. The DM I and II encompass residues ~10-100 and ~101-182, respectively. The catalytic residues are in a cleft between these two domains. The DM III possessing residues ~200-300 is reported to be involved in the dimerization of M pro . The DM II and III are linked together by a polypeptide loop. The overall structure of SARS-CoV-2 M pro is very similar to SARS-CoV M pro [25]. The SARS-CoV M pro has 96.1% of sequence identity compared to the SARS-COV-2 M pro (Table 1), suggesting a close evolutionary relationship with the SARS coronavirus from the first pandemic. The primary structure of M pro from Bat CoV-HKU4 has a relatively lower sequence identity with 49.2% but 3D structure share highly conserved structural features and conserved active site (Fig. 1B).
A number of studies proposed that the SARS-CoV-2 M pro exhibits a pH-dependent proteolytic activity that is important in dimerization and therefore essential for enzyme catalysis [26][27][28][29]. The crystal structures of the SARS-CoV M pro dimer have been solved under different pH conditions and the optimum pH for the enzyme activity was observed to be pH 7.4 [26,28]. Sharma et al. [30] demonstrated that the SARS-CoV-2 M pro is most stable at neutral pH whereas the protein destabilizes at acidic pH. Biochemical and biophysical characterization showed that the enzyme resists to thermal denaturation in a wide range of pH values (from 5 to 10) with preference to more basic pH [27]. The M pro enzyme has a high content of acidic residues (Fig. 1C). From an isoelectric point calculation (Table 1), the SARS-CoV-2 M pro molecule is negatively charged at a neutral pH. At a pH below the isoelectric point (pI) of the protein, this could influence the protein charge state which may subsequently alter the interactions and the stability of the enzyme. Since the dimerization is necessary for the enzymatic activity, changing the protonation state of the charged residues of the enzyme can significantly affect the stability of the protein dimer. However, interaction resulting from the affected residues in a reduced stability of the dimer is not fully understood. In this study, we employed a computational approach to gain insights into the effect of pH on the structure of protein dimers and to identify residues involved in dimer stability, especially those forming salt-bridge or hydrogen bond (Fig. 1C). For comparison, we carried out molecular dynamics (MD) simulations of the M pro from SARS-CoV-2 [24], SARS-CoV [31] and Bat CoV-HKU4 [32]. The details of simulated systems are summarized in Table 1.

Structure models and molecular dynamics simulations
Protein structures used in this study were taken from Protein Data Bank with PDB code as shown in Table 1. All simulations were carried out using dimer structure of CoV M pro . The simulated systems were built using visual molecular dynamics (VMD) [33]. The protein models were solvated using the water model TIP3P and neutralized using 0.1 M NaCl solutions. Based on the pI values (Table 1), the M pro protein carries a net negative charge at a neutral pH. Protonation states of the charged residues in the studied model systems were defined based on the pK a calculation of two pH conditions, in neutral, and in acidic solution. At neutral pH, the protonation states of ionizable residues in the proteins were assigned by comparing the pH values set to the system and the pK a values of amino acids predicted using PROPKA [34]. For the proteins at acidic pH, all negatively charged residues (Asp and Glu) are in protonated form since there is a difficulty in obtaining accurate structure-based pK a predictions for residues buried in protein interiors [35]. All simulated systems were run with a time step of 2 fs. Periodic boundary condition is applied with a box size of ~131 × 131 × 131 Å 3 . A distance cutoff of 12 Å was used for calculating nonbonded interactions. Electrostatic interactions were calculated with particle mesh Ewald summation via fast Fourier transform, including van der Waals interactions, with a switching distance of 10 Å. Langevin dynamics at a constant temperature of 298 K was used with a damping coefficient of 1 ps − 1 . The pressure was kept constant at 1 atm using the Nose − Hoover Langevin piston method, with a piston period of 200 fs and a damping time of 50 fs. Energy minimization was performed to remove bad contacts between atoms. Atomistic MD simulations of each M pro system were performed for 100 ns with three independent runs, each with different starting (initial) velocities; average with 90% confidence intervals are reported. The CHARMM36 force field parameters were applied for proteins [36]. The MD simulations were performed with the program NAMD version 2.12 [37]. A summary of simulation systems is shown in Table 1.

MD trajectory analysis
The structure and dynamic properties of the studied model systems were calculated as follows: root-mean-square deviation (RMSD), rootmean-square-fluctuation (RMSF) per residue, radius of gyration (R g ), the solvent-accessible surface area (SASA), secondary structure content, intersubunit distance between residue pairs, inter-domain distance between subunits, intersubunit binding free energy and interaction energy and the hydrogen bond occupancy. Analysis of the trajectory was mainly carried out using the simulation periods of 80-100 ns and averaged over 3 runs for each simulation system. The RMSDs relative to the initial structure were calculated using backbone atoms. The RMSF per residue relative to the average structure were calculated based on Cα atoms. R g and SASA of the dimer were computed using "measure" commands in VMD. The secondary structure content was obtained according to the secondary structure definition assigned by VMD. The hydrogen bond occupancy is the percentage of time that the hydrogen bond is observed during the trajectory. For analysis of hydrogen bond occupancy, the geometric criteria for distance and angle between the donor and acceptor heteroatoms were ≤ 3.5 Å and ≥ 120 • , respectively. Analysis of the structural properties were achieved using TCL scripts in VMD [33]. The average distance between heteroatoms of interacting pairs was plotted using python scripts and gnuplot.

MM-GBSA binding free energy calculation
The molecular mechanic method combined with the generalized Born and surface area continuum solvation (MM-GBSA) were employed to estimate the binding free energy of M pro dimer. This method has been widely used to assess the relative binding affinities between proteins and ligand as well as between protein molecules. In the MM-GBSA calculations, snapshots of the systems were extracted from the last 20 ns of MD trajectories, and the explicit water molecules and ions were removed. Binding free energy, ΔG bind , that is responsible for the stabilization of dimer is defined by, where solvation energy of dimer, ΔG solv Dimer , is subtracted with solvation energies of each monomer S and S ′ , ΔG solv S and ΔG solv S ′ , respectively. The difference in free energies between the dimer and monomer gives binding free energy. The free energy of each component is approximated by decomposing the energies into where ΔE(gas) contains the internal energy, ΔE(int), electrostatic energy ΔE(elec) and van der Waals energy, ΔE(vdW), calculated in the gas phase. ΔG(solvation) is composed of the polar contribution, ΔG(GB) and the nonpolar contribution, ΔG(nonpolar). The ΔG(GB) was calculated using the generalized Born implicit solvent methods with the interior and exterior dielectric constants of 1 and 78.5, respectively. The ΔG (nonpolar) was calculated from the solvent-accessible surface-area (SASA). All parameters are use according to the previous work [38]. The -TΔS term was not included in this study due to the low accuracy for estimating the entropic contribution. All the MM-GBSA calculations were carried out by NAMD version 2.12 [37,39] using CHARMM36 force field parameters [36].

Nonbond interaction energy calculation
The interaction analysis provided detailed understanding regarding the interactions within the M pro dimer at different pH conditions. To quantify the contribution of individual amino acid residues in the protein domains, nonbond interaction energy (E Inter Nonb ), including electrostatic (E Inter Elec ), van der Waals (E Inter VdW ) interaction energies and per-residue interaction energy were performed using the NamdEnergy plugin available in NAMD [37]. A distance cutoff of 12 Å was used [40]. The last 20 ns of MD trajectory for each simulation was used to compute the nonbonded interaction energy.

Principal component analysis
A principal component analysis (PCA) of the Cα atom fluctuations was performed to extract the concerted domain motions that are relevant for biological function. Prior to performing PCA, the set of conformations sampled from the last 40 ns (60-100 ns) of MD trajectories were aligned to its average structure to eliminate all translational and rotational motions. The principal components of the motions correspond to the eigenvectors of the Cα coordinate covariance matrices and the eigenvalues are the amplitude of the eigenvectors along multidimensional space. The displacements of Cα atoms of each eigenvector illustrate the concerted motions of the protein along each direction. PCA was performed using the WORDOM software [41]. The essential dynamics along the principal components were generated using VMD [33].

Root-mean-square deviation
RMSD with respect to the starting structure was computed to illustrate the influence of pH on the dimer structure in each model system during the MD simulations ( Fig. 2A). It should be noted that the RMSD has reached the plateau within the first 50 ns of the simulations for all studied systems, except for the acidic model system of Bat CoV-HKU4 Note a pI is the theoretical isoelectric point of a 3C-like proteinase monomer predicted from Expasy web server. b %identity of amino acid sequence with respect to that of the SARS-CoV-2 M pro were computed by Supermatcher (www.bioinformatics.nl). whose RMSD continuously increases throughout the simulation. As judge by convergent RMSD values, we analyzed the trajectory of the last 20 ns to calculate structural information to identify the structural changes induced by the pH change.
The backbone RMSD relative to the starting structure ( Fig. 2A) revealed that all neutral pH systems attained an equilibration with the fluctuation between 1.0 and 2.0 Å. On the other hand, the M pro dimer under low pH environment exhibited a wider fluctuation with the RMSD between 2.0 and 5.0 Å. The fluctuation in structural deviation was larger in acidic pH than that at neutral pH. The proteins tend to undergo substantial structural changes after modifying the charge state of sidechain of ionizable residues from neutral to acidic pH. All the RMSD profiles suggest the M pro dimer is less stable upon acidification. Based on our MD data, the enzyme maintained its dimeric structure better in the neutral solution with a significant loss of its dimer stability in the acidic solution. Fig. 2B shows the RMSF plots i.e. the amplitude of backbone fluctuation of individual residue with respect to its average reference position. We can see that the RMSF of the three proteins at neutral pH are similar in corresponding regions while their flexibility at acidic pH increases substantially throughout the entire protein region (Fig. 2B). Two notable regions, excluding disordered N-and C-termini residues, show a high level of flexibility compared with other regions. These two regions comprise of residues ~47-51 and ~180-190. The former is in a loop between 2 short helices of DM I while the latter is in the long linker loop between DM II and DM III. The most notable change in RMSF from neutral to acidic pH was observed for DM III in SARS-CoV-2 and Bat CoV-HKU4 and DM I in SARS-CoV. For the SARS-CoV-2 system at acidic pH, the highly flexible residues (~210-270) in DM III are located on the four helices, three of which are exposed to aqueous environment.
To identify residues and domains and quantify their relative fluctuation, RMSD per-residue was evaluated as presented in Supplementary  Fig. S1. We found that the magnitude of RMSD of almost all residues were higher on average with the acidic pH than that with the neutral pH. Consequently, the collective fluctuations of the initial positions of residues for each protein studied may result in higher fluctuations in all three domains when the pH changes from neutral to acidic pH. For example, large fluctuations (RMSD >4.5 Å) at low pH spanned across the protein particularly in the regions, including residues 47-52 in DM I, 164-172 in DM II and 186-193 in the loop (as shown by red contour in Supplementary Fig. S1). These regions contain ionizable residues (Supplementary Table S1) that are involved in formation/disruption of salt-bridges depending on the charged state of sidechain. We also notice from RMSD per-residue plots of SARS-CoV-2, SARS-CoV and Bat CoV-HKU4 that DM I and DM II have large effects at an early stage whereas DM III showed large conformational changes after ~20 ns of the simulation. While the acidity-induced structural changes appeared in DM I and DM II during entire timescale of present simulation, it showed at the later stage in DM III.
From a trajectory analysis, the total number of hydrogen bonds within each protein subunit remarkably decreased as pH change from neutral to low pH values (Supplementary Table S2). Such a change was associated with protonation of the carboxylate sidechain (glutamic acid and aspartic acid) that were engaged in inter-residue salt-bridge/ hydrogen bond interactions within the protein. Although some charged residues are accessible to the solvent, oppositely charged residues are close to enough to form ionic pairs or salt bridges. The protonated sidechain on acid residues that is electrically neutral, disrupted these interactions due to the increased electrostatic repulsion with the positively charged groups. Note that the observed conformational changes occurred in the same regions in both subunits (S and S ′ ), implying that the acid pH affects the monomer structure in similar fashion. Changes in the dimer structure and interaction energy in response to acidification will be discussed in further details later.

Structural properties
The influence of pH on the structural properties of M pro dimer were further investigated by analyzing such quantities as radius of gyration (R g ), solvent accessible surface area (SASA) and the percentage of secondary structure segments using VMD TCL command scripts. The radius of gyration (R g ) of M pro dimer was found to increase by changing to acidic pH from neutral, all variants of covid examined (Fig. 3A), i.e. conformation of the protein expanded by acidity. The solvent accessible surface area (SASA) of the M pro also increased with acidity, more in SARS-CoV and Bat CoV-HKU4 than that in SARS-CoV-2 (Fig. 3B). This trend implies that the dimer structures in acidic pH are more accessible to the solvent which is consistent with the loss of compactness of the dimer in a low pH environment. Fig. 3C and D show the partial loss of secondary structure elements (β-strand or α-helix) in the M pro dimer along the MD trajectories in acidic pH compared to that of the neutral system. We found that the partial disruption of the secondary structure was attributed mainly to β-strand. We noted that the DM I and II have the chymotrypsin-like double β-barrel fold structure whereas the DM III consists of five antiparallel α-helices. Thus, the loss of the β-strand secondary structure content occurred only in DM I and DM II but not in DM III (Fig. 1A). This implies that DM I and DM II are more susceptible to disruption of protein dimer upon acidification than DM III. Our MD results clearly show that a structure disruption of the M pro dimer is induced by pH which supports the hypothesis of "pH-dependent conformational changes" by Tan et al. [28]. Our results are also in good agreement with experimental analysis on decreased enzyme stability under acidic pH [27].

Intersubunit contact map
Intersubunit distance between residue pairs refer to contact between residues of one monomer (M pro ) with residues of another monomer of the dimer. Intersubunit contact maps showed the proximity between residues of the two monomers. They can be used to quantitatively evaluate the tendency of dimer dissociation. Contact maps presented in the 2-dimensional matrix (Fig. 4) of the protein dimer were based on the center-of-mass distance between intersubunit residues with 9 Å cut-off. At neutral pH, the dimer contacts of domains II-II, II-III and III-III are clearly visible, and their interactions are more numerous and stronger than those of domain -I-I, I-II and I-III. The pattern of the contact maps for SARS-CoV-2 and SARS-CoV is nearly identical, suggesting a very similar arrangement and spatial proximity of protein dimer. In addition, we observed strong contact for residues in the N-terminus and the DM II. The interactions between DM II and the NT of another subunit have been reported which suggests possible hydrogen bonding from the charged residues [3,24] that help control the correct orientation of binding site [42][43][44]. On the other hand, the residues of the Bat CoV-HKU4 M pro possess much smaller number of intersubunit domain contacts implying that the enzyme has a weak compact structure compared to SARS-CoV-2 and SAR-CoV enzyme. In the acidic pH systems, we observed the same contact maps with a lighter intensity for the SARS-CoV M pro . Especially, the separation distances for domains II-II, II-III and III-III became larger.
The average separation between domain subunits of the dimer was calculated using the last 20 ns MD trajectories; quantitative estimates are presented in Table 2. Overall, the separation between the same domains of the dimer, i.e. -I-I, II-II and III-III increased significantly as pH changes from neutral to acidic conditions, indicating that the two monomers moved apart further. To better understand the conformationally relevant motions of the proteins, we carried out a principal component analysis (PCA) of backbone fluctuations. MD snapshots of the last 40 ns were used to extract the dominant domain motions of the proteins. Fig. 5 illustrates the relevant collective motions of Cα atoms of SARS-Cov-2 M pro projected along the first eigenvector. PCA analysis reveals a rigid body domain motion in both subunits. At acidic pH, the three domains of each monomer rotate relative to the central body of the dimer by swinging outward away from the C2 symmetric axis of the protein. In the concerted domain rotation, DM I and DM II move outward in the same direction whereas the movement of DM III is perpendicular to that of the first two domains. This motion creates a scissor-like motion away from the C2 axis (supplementary movies S1). Based on the outward scissoring movement, the displacement of protein domains results in a larger separation between monomers.

Binding free energy of dimer
To quantify the binding energy of the dimer in the two pH conditions, we performed binding free energy calculations using MM-GBSA in NAMD. The binding free energy (ΔG bind ) that determines the stability of the dimer is shown in Fig. 6. The calculated values of binding free energy of the SARS-CoV-2 M pro dimer are − 137.81 ± 5.64 kcal/mol at neutral pH and − 70.28 ± 4.73 kcal/mol at acidic pH. The standard error is relatively small with respect to ΔG bind , implying that the simulation length associated with binding free-energy calculations is sufficiently enough. A comparison of ΔG bind values suggested that the SARS-CoV-2 M pro dimer at a low pH value was markedly weaker than that at neutral pH. This indicated that the interactions of dimer tend to become weaker, resulting in a decrease in dimer stability under acidic conditions. It should be noted that the decreased dimer stability was also be observed in SARS-CoV (− 110.33 ± 4.83 kcal/mol in neutral and − 68.97 ± 3.96 kcal/mol in acid) and Bat CoV-HKU4 (− 84.38 ± 4.80 kcal/mol in neutral and − 60.87 ± 3.71 kcal/mol in acid). From the results, the M pro exhibits the same pH-dependence of protein stability within the coronavirus family.

Dimer and domain interaction energy
To assess the contributions and their changes in the interactions within the M pro dimer at different pH values, we computed the nonbonded interaction energy based on the molecular mechanic calculation using NamdEnergy Plugin. The nonbonded interaction can play a part in stabilizing protein since the structure stability is a result of a net balance of interactions within protein and between protein and solvent. We quantified the energetic contribution of particular interactions at the level of subunit-subunit and subunit-domain. The calculated energy of the nonbonded interaction within the protein is composed of electrostatic contributions and van der Waals interactions and cooperativity between the two nonbonded energies is critical to the stability. Fig. 7 shows intersubunit and interdomain nonbonded interaction energy together with electrostatic and van der Waals interaction energies. The graphs showed that the electrostatic interaction is a major contribution of the interaction energy. The van der Waals energy contribution appears to be relatively small and nearly equivalent to all systems. From a relative comparison, the electrostatic interactions play the dominant role in the dimer stability of the enzyme. For SARS-CoV-2 and SARS-CoV, the interaction energy between the two monomers of M pro was attractive at neutral pH while a low pH value caused the intersubunit interaction to become more repulsive due to strong charge-charge repulsion (the positive energy values in Fig. 7A). We found that an increase of the electrostatic repulsion involved disruption of salt-bridge interactions in the protein. Apparently, the formation of a salt bridge strongly depends on the protonation state between two charged residues and thus on pH. A change in the protonation could generate large electrostatic effect to eliminate the salt-bridges and therefore destabilize the intersubunit and interdomain interactions of the protein. We also found that at neutral pH, the attractive interaction between subunit and domain I and III has a relatively large contribution to the net interaction energy (Fig. 7B-7D). The subunit-domain II interaction is relatively weak for SARS-CoV-2 and SARS-CoV but such interaction plays a significant role in the dimer stability of the Bat CoV-HKU4 enzyme.
Further, per-residue interaction energy analysis of SARS-CoV-2 M pro was performed to elucidate the contribution of important amino acids involved in the stability of the protein dimer. Per residue interaction energy (E res− subunit Nonb ) is computed based on the non-bonded molecular mechanic (MM) energy between all atoms of a given residue and all other atoms of another subunit. The per-residue energy plot showed that the N-terminal and DM III residues make significant contribution to the interaction energy of the dimer (Supplementary Fig. S2). From the energy plot, strong attractive (highly negative energy) interactions were found in residues Arg4, Lys12, Glu14, Glu288 and Glu290 at neutral pH (Fig. S2A). However, interactions of these residues became very weak attractive (near zero energy) or repulsive (highly positive energy) in acidic solution (Fig. S2B). It appears that most affected residues are charged residues which are located within 10 Å at the dimer interface. Particularly, Arg4, Lys12, Glu14 and Glu290 contribute significantly over other amino acids on the M pro subunit toward the overall energy of dimer interactions.

Intersubunit hydrogen bond
Since the enzyme is functional in the dimeric form, the hydrogen bond of residue pairs between subunits plays a major role in the stability of the protein dimer. Analysis of intersubunit hydrogen bonds revealed residues involving domain interactions within the M pro dimer. The presence and absence of Intersubunit hydrogen bonding upon a pH change from are mapped on the protein domains as shown in Fig. 8, together with residue pairs involving in the interactions (Table 3 and Fig. 9). From the simulation of SARS-CoV-2, we found intersubunitinterdomain hydrogen bonds, involving residue pairs Ser1-Phe140', Ser1-His172', Ala7-Val125' and Ser1-Glu166' for NT-DM II interaction (Fig. 9B), Arg4-Glu290' for NT-DM III interaction (Fig. 9C), Ser10-Ser10', Glu14-Ser10' and Glu14-Gly11' for DM I-DM I interaction (Fig. 9D), and Ser139-Gln299' for DM II-DM III interaction (Fig. 9E). The intersubunit and interdomain hydrogen bonds observed for SARS-CoV and Bat CoV-HKU4 M pro were slightly different. These observed hydrogen bonds which are characterized as strong and moderate based on the high percent occupancy calculated from the trajectory data, contribute to the stability of the dimer. Strong hydrogen bond is defined as the percent hydrogen bond occupancy greater 80% while moderate and weak hydrogen bond are characterized by 40-80% and < 40% of the occupancy, respectively [40,45]. In the acidic pH, some residues show a loss of hydrogen bonds as well as contacts between subunits. These pairs of residues are: Ser1-His172' and Ser1-Glu166' for NT-DM II interaction, Arg4-Glu290' for NT-DM II interaction, Glu14-Ser10' and Glu14-Lys12' for DM I-DM I interaction, and Ser139-Glu299' for DM II-DM III interaction ( Table 3). As can be seen in Table 3 and Fig. 9B-9E, a disruption of most of these hydrogen bonds caused the enzyme being less stable in dimeric state at the acidic pH. As was evidenced by the larger separation distance between subunits and the loose compactness of dimer, a decline in the dimer stability must be associated with the loss of these hydrogen bonds. This loss has also substantially contributed to the loss in binding free energy (ΔG bind ) and the interaction energy. Overall, our simulation results provide strong support for the decreased stability of the M pro dimer at low pH condition.
It becomes clear that neutralization of the negatively charged residues of the protein in the acidic pH makes the overall interactions more repulsive. The repulsive interaction energy is mainly driven by electrostatic effects, eliminates salt-bridges in the protein dimer and hence destabilizes the dimer structure. However, the calculated binding free energy at acidic pH shown in Fig. 6 is still negative. This could be attributed to the solvation free energy that play the dominant contribution to binding free energy of the dimer.
An earlier study of pH-dependence SARS CoV-2 main protease showed structural stability during 250 ns MD simulation over the pH range 4 to 7 while the enzyme became unstable at pH 3 [46]. The main difference between this work and ours is the charge state of the proteins used in the simulations. In our protein model, acidic residues (Asp and Glu) have a neutral form at which the pH value was not explicitly given to our simulations in an acidic solution. This is because there is an argument for obtaining an accurate protonation state of the protein at a specific pH value using pK a calculations [35]. On the other hand,  Barazorda-Ccahuana et al. used the protein with a different charge state at different pH depending on the prediction of pK a values of ionizable groups. This difference affects significantly to protein conformations which were generated from different initial protonation state of the protein. It should be noted that our protein has a closer value of the total charge to their protein at pH 3 where both proteins show similar instability. In addition, we analyzed the change in conformation around the S1 pocket which is implicated in substrate binding and protein dimerization. The Cα-atom distance of three pairs of residues including Cys145-His41, Met49-Gln189 and His172-Phe140, was used to estimate a change in the pocket volume of the protein. We found that the separation distance of residues Met49-Gln189 and His172-Phe140 increases remarkably (~5 Å) when the pH changes from neutral to acidic conditions ( Supplementary Fig. S3). This result is consistent with an increased volume of the S1 pocket at pH 3 proposed by Barazorda-Ccahuana [46].
Previous theoretical and experimental studies have shown that the interaction involving the NT residue is important in dimer stability [3,43,[47][48][49][50][51]. Especially, the strong hydrogen bond for the salt-bridge interaction between Arg4 of the NT and Glu290' of the DM III found in both SARS-CoV-2 and SARS-CoV at neutral pH was not detected at acidic condition (Table 3 and Fig. 9). Hsu et al. proposed that the Arg4 in NT and residues in DM III are critical for stabilizing the M pro dimer by maintaining a correct conformation of the active site [52]. Lowering the pH led to a decrease in the M pro dimerization and activity [29]. Our simulations indicated that the Arg4-Glu290 hydrogen bond was not only involved in the observed decreased stability of M pro dimer but intersubunit hydrogen-bonding residues of NT-DM II, DM I-DM I and DM II-DM III interactions should also be taken into consideration. The disruption of hydrogen bonds formed by these residues can interrupt enzyme dimerization, and thus this offers an alternative drug design approach to interfere the viral replication.
The discovery of potential repurposed antiviral agents has brought viral main protease, M pro , to be an appealing therapeutic target for the treatment of COVID-19. Several available drugs targeting viral proteases such as HIV-1 and hepatitis C virus (HCV) protease inhibitors are   currently being studied to test the safety and efficacy of the drug against this pandemic. As shown by x-ray crystallographic and biochemical studies, the FDA-approved anti-HCV drugs (boceprevir, narlaprevir and telaprevir) have been shown to effectively inhibit the enzyme by binding to the active site [53]. Since the enzyme is functional as dimer, a potential alternative mode of inhibition is to design inhibitors that bind efficiently to the charged residues critical for the dimer interaction. Our study showed that the interface-forming residues located in the N-terminus and DM III are capable of either disrupting the dimer or preventing dimer formation. The characteristics of subunit interface from this study can potentially be a novel site for drugs to interrupt the dimer formation as well as to reduce the stability of the M pro .

Conclusions
Stability of the enzyme dimer was investigated in neutral and acidic solutions using MD simulations. The dimer stability was assessed by examining the loss of residue contacts, enhanced separation between the domains of protein monomers, reduced binding free energy and the interaction energy of the dimer in acidic medium in comparison to neutral environment. The structural fluctuation was larger in acidic pH than that with the neutral pH. Based on our MD data, the enzyme maintained its dimeric structure better in neutral pH condition with a substantial loss of its stability with acidic pH. Separation distance between domains (I, II, III) of one protein monomer and that of the second in the dimer was found to increase due to acidity. We have identified the residues that can stabilize the dimer conformation via hydrogen bonds of domains. The interactions of these residues may be helpful for rational drug design process against COVID-19.

Notes
The authors declare no competing financial interest.

Declaration of Competing Interest
There are no conflicts of interest to declare.

Data availability
Data will be made available on request.