Transmembrane Protease Serine 2 Proteolytic Cleavage of the SARS-CoV-2 Spike Protein: A Mechanistic Quantum Mechanics/Molecular Mechanics Study to Inspire the Design of New Drugs To Fight the COVID-19 Pandemic

Despite the development of vaccines against the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) virus, there is an urgent need for efficient drugs to treat infected patients. An attractive drug target is the human transmembrane protease serine 2 (TMPRSS2) because of its vital role in the viral infection mechanism of SARS-CoV-2 by activation of the virus spike protein (S protein). Having in mind that the information derived from quantum mechanics/molecular mechanics (QM/MM) studies could be an important tool in the design of transition-state (TS) analogue inhibitors, we resorted to adiabatic QM/MM calculations to determine the mechanism of the first step (acylation) of proteolytic cleavage of the S protein with atomistic details. Acylation occurred in two stages: (i) proton transfer from Ser441 to His296 concerted with the nucleophilic attack of Ser441 to the substrate’s P1-Arg and (ii) proton transfer from His296 to the P1′-Ser residue concerted with the cleavage of the ArgP1-SerP1′ peptide bond, with a Gibbs activation energy of 17.1 and 15.8 kcal mol–1, relative to the reactant. An oxyanion hole composed of two hydrogen bonds stabilized the rate-limiting TS by 8 kcal mol–1. An analysis of the TMPRSS2 interactions with the high-energy, short-lived tetrahedral intermediate highlighted the limitations of current clinical inhibitors and pointed out specific ways to develop higher-affinity TS analogue inhibitors. The results support the development of more efficient drugs against SARS-CoV-2 using a human target, free from resistance development.


MOLECULAR DYNAMICS SIMULATION ANALYSIS
In order to characterize the enzyme:substrate complex, we have performed a molecular dynamics (MD) simulation study. This simulation was done after two minimization steps and one water relaxation step, which were performed to obtain a starting enzyme:substrate complex for the QM/MM study. Our MD protocol consisted of four additional steps: 1) system heating; 2) system relaxation; 3) equilibration; and 4) production.
The system heating step was employed during 200 ps, with a starting temperature of 0.0 K and at a constant volume. The temperature was increased every 0.5 ps until the system reached 300 K.
Then, the temperature was maintained with a weak-coupling algorithm 2 . During this step, the protein and the octapeptide residues were restrained with a restrain weight of 5.0 kcal/(mol•Å 2 ) (this was also the case for the following system relaxation step).
The system relaxation step consisted of a 300 ps simulation with a reference temperature of 300 K. A Langevin dynamics protocol with a collision frequency of 1.0 ps -1 was used 2 . Concerning the system pressure, an isotropic position scaling with a pressure relaxation time of 2.0 ps and a Berendsen barostat 3 were employed to maintain a 1 bar pressure.

5
The last two steps did not have any restraints applied. The equilibration step consisted in a simulation of 2 ns, with the same reference temperature and scaling method as the one employed in the system relaxation step, except for the collision frequency, which in this step was 2.0 ps -1 (this was also the case for the production step).
The production step consisted in a simulation of 100 ns. In every 5000 integration steps the coordinates were saved. Like the system relaxation step, the pressure of the system was maintained using an isotropic position scaling with a pressure relaxation time of 2.0 ps. However, in this case we employed a Monte Carlo barostat 2 .
To evaluate the structure and stability of our system during the simulation, we measured the root mean square deviation (RMSd) of i) all protein atoms, ii) the three catalytic residues atoms, iii) the atoms included in the QM region of the QM/MM study, and iv) all substrate atoms. These measures were done using the minimized structure as a reference (i.e., the starting structure used in the QM/MM study). In addition, the root mean square fluctuation (RMSf) of all substrate's Cα, and the key distances for the reaction were also analyzed. In Figure S2 (black curve), we observed a relatively good stability of the TMPRSS2 protein throughout the simulated time, with an average RMSd of 1.60 Å. Concerning the three catalytic residues ( Figure S2, red), they remained stable throughout the entire simulation, having only observed larger fluctuations (but still below 1 Å) in discrete periods of the 100 ns production (e.g. 7 35-48 ns). This was also observed in the atoms that were included in the QM region of the QM/MM study ( Figure S2, blue).
In addition, we have also evaluated key distances for the reaction during the production stage.
These included the nucleophilic attack and proton transfer distances, the oxyanion distances, and the hydrogen bond distances between the side chains of His296 and Asp345 ( Figure S3). A) The distance between Ser441 Oγ and Arg P1 C1 (nucleophilic attack) is represented in black, meanwhile the distance between Ser441 H1 and His296 Nε (proton transfer) is 8 represented in red. B) The distance between Arg P1 O1 and Gly H is represented in blue, whilst the distance between Arg P1 O1 and Ser441 H is represented in green. C) The distance between Asp345 Oδ1 and His296 Hδ is represented in orange, meanwhile the distance between Asp345 Oδ2 and His296 Hδ is represented in violet.
In Figure S3.A, it is possible to observe that during most of the simulation both interatomic distances remained stable. However, and in agreement with the observation made in Figure S2, between 35 ns and 48 ns there was a higher destabilization of these two distances, with the proton transfer distance showing a larger increase than the nucleophilic attack one.
In Figure S3.B it is possible to observe that the oxyanion hole interactions remained stable during the MD simulation. However, it appears that the interatomic distance between Arg P1 O1 and Ser441 H suffers a slight fluctuation throughout the 100 ns simulation.
In Figure S3.C one can observe that the interatomic distance between both Asp345's sidechain oxygen atoms and the His296's Hδ remained stable during the simulation, meaning that the interaction between these catalytic residues was stable. This interaction is important to maintain His296 with the required pK a for the proton transfer.

9
Considering both the RMSd analysis and the variation of the key interatomic distances for the reaction mechanism of TMPRSS2, we observed that the active site was stable during most of the simulation. Still, between 35 and 48 ns of the 100 ns production stage, we observed a destabilization of the active site, with pronounced fluctuations in the nucleophilic attack and proton transfer distances (involving the catalytic Ser441).
To understand the reason for the deviation shown on the RMSd graphics and the interatomic distances, five structures from different periods of the simulation, were retrieved. These structures were compared with the minimized structure used in the QM/MM study ( Figure S4).
Concerning the structure retrieved at 27.40 ns, in agreement with the observations made during the RMSd and the interatomic distance analysis, it is possible to see that the position of the catalytic residues and the oxyanion hole Gly was very similar. This allowed the enzyme to maintain these residues in the required conformation for initiating the catalytic mechanism. Contrary to this structure, in the structure retrieved at 47.86 ns, we observed that Ser441 suffer a significant change in its orientation. In this structure, the Ser cat 's sidechain rotated, which led to a significant increase between the proton transfer participating atoms. This rotation also led to an increase between the nucleophilic attack participating atoms. Due to this change in orientation, the hydrogen bond interaction between Ser441's H1 and His296's Nε was broken, resulting in a slightly different conformation of His296. However, the hydrogen bond interactions between His296's Hδ and Asp345 sidechain oxygen atoms were preserved. After this period, the graphical analysis showed that the system evolved back to the catalytically competent conformation. In fact, the structure retrieved at 49.95 ns showed that Ser441's sidechain returned to the required orientation, allowing the re-establishment of the hydrogen bond interaction with His296.
Contrary to the protein, the RMSd of all substrate's residues showed that the octapeptide was more flexible, especially in the last 15 ns of the simulation ( Figure S3, green). To explore the influence of this flexibility on the QM region, two structures were retrieved: i) one at 88.79 ns, and ii) another at the end of the simulation. Concerning the 88.79 ns structure, we observed that the position of most of the atoms included in the QM region was very similar to the minimized structure, with the exception being Ser441's Oγ. However, the last structure of the simulation showed that this atom evolved back to the required orientation for the reaction. In fact, if we 13 analyze the nucleophilic and proton transfer distances in the last 15 ns of the simulation, we do not observe a high destabilization of these two distances. This means that the large fluctuation of the substrate probably occurred beyond the P1-P1' region. To explore this effect, we analyzed the RMSf of the substrate's Cα ( Figure S5). In Figure S5, one can observe that the residues near the N-terminal and C-terminal of the peptide (and especially the C-terminal), are significantly more flexible than the residues near the middle of the octapeptide. In fact, Arg P1 is the most rigid residue in the substrate, with next being Ser P1' . So, the high RMSd variation for the ligand during the 100 ns simulations, does not originate in residues of the substrate that are included in the QM region. To confirm this observation, several 14 structures of the substrate were retrieved from the simulation, and subsequently aligned to compare the conformation of all octapeptide's residues ( Figure S6).
The structural comparison of the evolution of the substrate throughout the simulation confirmed the observation made with the RMSf analysis. In fact, one can observe that the residues near both the N-terminal and C-terminal regions are more flexible throughout the simulation, than the ones involved in the reaction. Both, Arg P1 and Ser P1' are significantly more rigid, and consequently more stable than Pro P4 and Glu P4' .