{'fraction_non_alphanumeric': 0.03126169973792587, 'fraction_numerical': 0.03172968925496069, 'mean_word_length': 5.147583429228999, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Identification of novel TMPRSS2 inhibitors for COVID-19 using e-pharmacophore modelling, molecular docking, molecular dynamics and quantum mechanics studies

SARS coronavirus 2 (SARS-CoV-2) has spread rapidly around the world and continues to have a massive global health effect, contributing to an infectious respiratory illness called coronavirus infection-19 (COVID-19). TMPRSS2 is an emerging molecular target that plays a role in the early stages of SARS-CoV-2 infection; hence, inhibiting its activity might be a target for COVID-19. This study aims to use many computational approaches to provide compounds that could be optimized into clinical candidates. As there is no experimentally derived protein information, initially we develop the TMPRSS2 model. Then, we generate a pharmacophore model from TMPRSS2 active site consequently, and the developed models were employed for the screening of one million molecules from the Enamine database. The created model was then screened using e-pharmacophore-based screening, molecular docking, free energy estimation and molecular dynamic simulation. Also, ADMET prediction and density functional theory calculations were performed. Three potential molecules (Z126202570, Z46489368, and Z422255982) exhibited promising stable binding interactions with the target. In conclusion, these findings empower further in vitro and clinical assessment for these compounds as novel anti-COVID19 agents.


Introduction
The current COVID-19 coronavirus disease pandemic is a worldwide health concern associated with SARS-CoV-2 illness [1,2]. As of August 7th, 2021, the pandemic has grown to over 200 million illnesses and over 4 million deaths [3]. This epidemic has impacted negatively on national healthcare systems and rained down the global economy, necessitating the development of efficient therapies to halt its spread and control its severity [4].
Coronaviruses can affect both animals and humans, producing a wide range of illnesses [5]. MERS coronavirus, SARS coronavirus and the 2019 coronavirus that originated from Wuhan in Dec 2019 are the three contagious viral human coronaviruses (hCoVs) that are being identified to date [6,7]. However, from WHO statistics, the incidence of SAR-CoV-2 transmission in people has far outpaced that of SARScoronavirus and MERS-coronavirus [8].
Transmembrane spike glycoproteins (S proteins) extending from the viral surface to give the coronavirus its crown-like ("corona") appearance [9]. The alpha, beta, gamma, and delta are by far the four classes of coronaviruses [10]. Bats are thought to be the source of SARS-CoV-2, a beta coronavirus that can be transmitted without causing illness. The infectious agent of is a virus that has 4 glycoproteins: spike, envelope, membrane, and nucleocapsid [11,12].
The virus's interaction with the host influences the incidence and progression of COVID-19 illness [13]. For viral proliferation, all corona viruses have unique genes that encode protein in downstream ORF1 regions. Spike glycoprotein controls the virus adhesion and entry [14]. The virus's transmembrane spike glycoprotein attaches to the ACE2 protein to accelerate virus penetration into the host cell [15][16][17].
Various viruses, including SARS coronavirus, Ebola virus, influenza virus and MERS coronavirus, utilize host cell proteases to activate their envelope glycoproteins, according to prior [18][19][20]. Transmembrane protease/serine subfamily member 2 cleaves and activates the SARS-CoV-2 spike protein, which is critical for membrane fusion and host cell penetration [21][22][23]. It's worth mentioning that when the corona spike protein binds to the ACE2 receptor, the host cell's surrounding serine protease TMPRSS-2 receptor facilitates in the virion's entry via priming the S-2 domain of the spike protein [10,22,[24][25][26]. Since TMPRSS2 mRNA is abundant in a variety of tissues including the prostate, kidney, small intestine, lung, and stomach and others, it is linked to the progression and severity of the illness [1,18,27,28].As a result, TMPRSS2 processing is amongst the most important stages in activating the SARS-CoV-2 S protein's membrane activity [29]. Because of the extreme role of TMPRSS2, it might be a promising therapeutic target in the struggle against COVID-19 infection [30,31].
Intensive attempts by several groups globally are in full turn to find novel agents, but this may take a considerable time since current drugs have to progress through long preclinical and clinical safety testing [32]. Structure-based drug discovery is one of the most and newest used techniques to investigating a medication's action [33]. When experimental approaches fail or there is insufficient data about the target, homology modeling is a successful in silico strategy for developing a tertiary structure of the target protein to estimate structural information for docking studies [34]. The number of compounds that a medicinal chemist can synthesize can be substantially decreased by using pharmacophore model creation. As a result, drug research and development time and cost can be decreased [35]. It offers the estimation of activity in a wide range of structurally varied compounds and could also assist in the discovery of novel molecules with enhanced activity [36,37]. The pharmacophore model shows the key chemical characteristics that are responsible for activity in a logical manner. Pharmacophore modeling has been one of the most fundamental and promising approaches for novel drug development in recent years [38].
In present study, we aim to find novel TMPRSS2 inhibitors through using in silico computational approaches namely homology modeling to construct TMPRSS2 protein model, pharmacophore mapping and screening, docking to calculate their binding affinities and their free binding energies. Molecular dynamics and density functional theory calculations were also performed to validate the results.

Method
All the in silico studies were carried out by Maestro v 12.8 (Schrodinger 2020) as depicted in Fig. 1.

Homology modeling search and model preparation
TMPRSS2 sequence was downloaded from UniProt (accession number: O15393) [39]. A similarity search was done using BLAST server [40] in order to find the best template for modeling of TMPRSS2. Human plasma kallikrein was the best template for homology modeling based on sequence identity of 42.21% and sequence similarity of 47%. Accordingly, human plasma kallikrein 3D structure (PDB ID: 5TJX) [41] was obtained from protein databank and used as a template in the Prime module of Schrodinger suite [42]. Alignment was performed using single template alignment (STA) and knowledge-based methods of the protein predication wizard were used to build a model using Prime of Schrodinger. Then, using the protein preparation wizard tool [43], the protein structure model was preprocessed to correct the order of all bonds in the structure and through Epik, proper ligands ionization and tautomeric state were confirmed in the specified pH range. Hydrogen atoms were added to the structure; bonds were broken to metals and correct the formal charge on the metal and the neighboring atoms to treat the bonds as ionic, which is required for use with force fields, which treat metal compounds as ionic rather than covalent. Then add zero-order bonds between the metal and its ligands, so that it is still considered part of the same molecule. It also finds sulfur atoms that are within 3.2 Å of each other and adds bonds between them. Also, the model was refined by optimizing hydrogen bond order using PROPKA in PH = 7 and restrained minimization was performed in OPLS3e force field to decrease side-chain energies.
The Ramachandran plot obtained by Prime module of Schrodinger was used for further verification in order to estimate model validity for the optimized homology model [44].

Active site prediction
The SiteMap tool was used to predict receptor sites for TMPRSS2 protein. This software uses grid points on the surface of the protein that could be suitable for ligands to bind to the receptor to provide data regarding binding site qualities. The final part of the process is the evaluation, which includes rating each site by computing different characteristics: two characteristics that provide distinct measurements of the site's availability to the solvent are exposure and enclosure, Donor/acceptor character, H donor/acceptor areas; contact, SiteScore and Ability of the site to absorb and bind to small molecules (druggability).

Ligand library preparation and grid generation
The chemical structures of million compounds from Enamine database [45] were downloaded and prepared using LigPrep [43]. LigPrep assures that compound structure is optimized, that structurally undesirable ligands are excluded, and that a low-energy 3-dimensional structure with optimal chiralities is developed for the input ligand. Briefly, the ionization states of ligands were defined at pH 7.0 ± 2.0, and the salt was removed using Epik. For each ligand, a maximum of 32 conformations was allowed to be generated. The energies of each ligand were minimized by OPLS3e. The Grid Generation panel was used to create the grid. The grid was constructed at the active site's centroid point for consistency.

E-pharmacophore generation
Using the option Generate pharmacophore from Receptor Cavity in the Phase module of Schrodinger [46], an energy-optimized pharmacophore hypothesis was developed using receptor centroid of the amino acids of the active site predicted by SiteMap. Before generating the e-pharmacophore, Phase first recruits the Glide XP tool to dock the already bound ligand inside the binding pocket of the protein under study and selects the top ranked pose for auto-pharmacophore generation, taking into account the following six chemical features: aromatic ring (R), hydrogen bond acceptor (A), hydrophobic group (H), hydrogen bond donor (D), negative ionic feature (N) and positive ionic feature (P).

E-pharmacophore based database screening
For the e-pharmacophore screening against TMPRSS2, the prepared conformations of the Enamine database were utilized. The compounds had to match a minimum of five features on the created e-pharmacophore model in order to get TMPRSS2 inhibitors with the appropriate chemical characteristics.

Docking of the generated compounds and MM-GBSA
The molecular docking makes it possible to predict the strength and affinity of the generated compounds for the target and also to determine the compound's binding site residues in the target protein and their interactions. The resulted conformations were docked intoTMPRSS2 active site via HTVS mode of glide [42], and then the top-ranked poses were subjected to SP mode of glide and further to XP mode. The results were compared against approved TMPRSS2 inhibitors ambroxol acefyllinate and nafamostat mesylate.
By using the default settings of the Prime-MMGBSA function [42], MM-GBSA estimations have also been done for the top-ranked XP docked poses and compared with known TMPRSS2 inhibitors ambroxol acefyllinate and nafamostat mesylate.

Molecular dynamics
The stability of binding of the top compounds in complex with TMPRSS2 was studied using by molecular dynamics (MD). It helps assess patterns of interaction between proteins and ligands and timedependent protein-ligand interaction associations and conformational dynamics of complex systems. Using Academic Desmond v6.5 by D.E. Shawn Research [47], the protein-ligand complexes were minimized and relaxed under NPT using Desmond's default method. Simulation time was 50 nanosecond, the temperature (300 K) and pressure (1.01325 bar) were kept constant throughout the process by Nose-Hoover thermostat and Martyna-Tobias-Klein barostat. The OPLS3e force field was used in MD simulations. 5208 frames for each system during the simulations were collected. Lastly, the interaction analysis was performed.

Density functional theory
Using Jaguar [47] the top molecules were subjected to density functional theory calculations using optimization panel. After structures were included in the workspace, the theory was adjusted to be DFT which allows employing a variety of functional to describe exchange and correlation for either open or closed-shell systems. Utilizing a basis set of 6-31G** level, complete geometry optimization was done by using functional B3LYP. The geometries were optimized 100 times. The accuracy for SCF calculations is quick, which Use mixed pseudo-spectral grids with loose cutoffs (iacc = 3). The Poisson Boltzmann solver was used to calculate energy in a water body to model physiological circumstances. Molecular electrostatic properties were analyzed.

Homology modeling
The TMPRSS2-spike structure is currently unavailable, its 3-dimensional structure with 234 amino acids was modeled utilizing human plasma kallikrein as a template to better explain the substrate binding mechanism to the TMPRSS2 receptor. With the TMPRSS2 sequence, the plasma kallikrein template has the highest level of similarity and identity. TMPRSS2 model was subjected to model refinement and energy minimization ( Fig. 2A). The Ramachandran plot (Fig. 2B) shows the quality of the projected models, since most residues are in favorable and permitted areas of the plot.

Active site prediction
Characterization of the binding site and its physicochemical characteristics is essential not only in the creation of new compounds, but also in the identification and optimization of molecules with known binding profiles. The template structures derived from the PDB database did not contain any co-crystallized ligands. As a result, the binding sites were obtained using SiteMap (Table 1 and Fig. 3). Site1 was selected as

E-pharmacophore model and pharmacophore-based database screening
Based on amino acid residues of the selected active site, a sevenfeatured pharmacophore hypothesis consisting of one donor (D), one hydrophobic region (H) and five aromatic rings (R) were obtained (Fig. 4).
Validated e-pharmacophore models of TMPRSS2 were used to screen over 4 million prepared molecules from the enamine database. From the total numbers of compounds and conformers present in database. 31340 compounds match the screened hypothesis. Following that, these hits went through a three-tier Glide docking method for docking-based ligand screening.

Docking of the generated compounds and MM-GBSA
In computational chemistry, docking has now become a common method for estimating the affinity of molecules to their targets. The 31059 compounds that obtained from phase screen were narrowed down to only 10 top-ranked molecules with docking score − 9.68 to − 8.017 kcal/mol, through successive docking by Glide HTVS, SP and XP in order to determine their binding interactions and affinity to TMPRSS2 compared to docking score of − 5.42 and − 6.46 kcal/mol for references nafamostat and ambroxol as shown in Table 2. In addition, the top 10 compounds were investigated for MMGBSA binding free energy calculations. They showed binding energies better than the reference compounds, ranging from − 63.55 to -89-75 kcal/mol. The interactions of top 3 molecules were selected for further analysis (Fig. 5). Z126202570 molecule forms four hydrogen bonds, one pication interaction, one pi-pi interaction and hydrophobic interactions with the active pocket of TMPRSS2. Amino acids involved in the hydrogen bond formation are GLU299, GLY464, ASP435 and PO4501. One pi-cation interaction was observed between molecule and HIS296 residue and the one pi-pi interaction with TRP461. The hydrophobic interactions were formed with CYS437, CYC465, TRP461, TYR337, VAL473, ILE420, LEU419, MET424 and TYR33 residues.
Binding affinity between active pocket amino acids of TMPRSS2 and Z46489368 was maintained with two hydrogen bonds, one pi-pi interaction and hydrophobic interactions. SER460 and PO4501 made the hydrogen bonding, the one pi-pi interaction with TRP461 and the hydrophobic interaction were observed with the CYS437, CYC465,  TRP461, VAL473, ILE420, LEU419 and MET424 residues. Many interactions were formed between TMPRSS2 active pocket and Z422255982. Of these interactions, three hydrogen bonds with SER460, LYS300 and PO4501 residues, one is a pi-pi interaction with TRP461, one is a pi-cation interaction with HIS296 residue and the remaining were hydrophobic interactions with CYS437, CYC465, TRP461 and TYR337 residues.

Molecular dynamics
Docking methods are often quick and imprecise, but they lack protein flexibility, which might compromise the precision of the resulting ligand-protein complexes. As a result, more precise molecular dynamic modeling approaches may provide a better docking complement. Here, we examined the various molecular simulation results including RMSD, Table 2 Docking score, MM-GBSA results, interactions and fitness score of the best ten compounds with TMPRSS2.  RMSF and PLC. RMSD represents the dynamic stability and conformation changes occur in Cα backbone of the TMPRSS-inhibitor complexes during the simulation. The average RMSD values of the TMPRSS2 protein and TMPRSS2-inhibitor complexes are shown in Table 3 and Fig. 6. The Cα backbone RMSD organization patterns of the TMPRSS2 and its inhibitors complexes showed no significant changes, indicating their stability. The RMSF was used to evaluate the flexibility of the ligands and proteins. Throughout the simulation, the coordinates of the Cα atom fluctuate from their typical location for example the average Cα fluctuation in Z126202570 was from 0.66 to 1.04 Å. The average RMSF values of the TMPRSS2 protein and ligands are shown in Table 3 and Fig. 7. The lower RMSF and RMSD (≤2.35 Å) average for both the TMPRSS2 and its inhibitors indicate how reliable our predicated model and the affinity of these promising inhibitors.
While in Z422255982-TMPRSS2 complex (Fig. 8C), only two types of interactions were formed. GLU299 maintains strong hydrogen bond interaction for 80% of the stimulation time. Weak water bridge interactions were shown by HIS296, CYS297, and LYS340.

Density functional theory
The molecular orbitals provide details on electrical and optical characteristics, as well as quantum chemistry, and are utilized to understand the different types of conjugated system processes. The reactivity of compounds is also linked to molecular orbital energy and orbital coefficients according to molecular orbital theory. The B3LYP/6-311** basis set was used to calculate the graphical representations of MOs and the solvation energy for the top three compounds as shown in Table 4 and Fig. 9.
The HOMO orbitals in Z126202570 are located on the N8 that interacts with SER460, aromatic ring and O1 that interact with SER441 and O10 that interacts with PO4501.The LUMO orbitals are found on S19, N17, N32, C16, C18, C20, C18 of 3H,4H-thieno[2,3-d]pyrimidin-4- The Molecular Electrostatic Potential (MEP) is a crucial way of understanding electrostatic allocation potential and visualizing the reactive site in chemical reactions involving both electrophilic and nucleophilic interactions. The ESP surface map provides a significant visual method to comprehend about reactive sites. MEP maps for the top three compounds were computed at the B3LYP/6-311** level to address potential sites for nucleophilic and electrophilic interactions and represented in Fig. 10.
The MEP map of Z126202570 range from − 1.087 to 155.8 kcal/mol. The highest negative (red) regions were spread over O31 of 3H,4Hthieno [2,3-d]pyrimidin-4-one moiety and O10 of the propan-2-one moiety that interact with PO4501. While the highest positive (purple)   regions are located on N14 of the piperidine moiety that interacts with GLU299 and HIS296. The MEP map of Z46489368 range from − 53.9 to 64.8 kcal/mol. The red regions are spread over 2,5-dimethylpyridine moiety that interact with TRP461 and on O20 of N-(2-methylprop-2-en-1-yl)acetamide moiety. While the purple regions are presented on N14 of ethyl(methyl) amine moiety.
The MEP map of Z422255982 range from − 4.5 to 125.4. The red regions are located on O2 of the benzenolate moiety that interacts with LYS300 and on N10 of the 1H-pyrazol-1-ide moiety. While the purple regions are found on N33 of the imidazo[1,2-a]pyridin-1-ium moiety that interacts with HIS296 and TRP461.

Conclusion
The TMPRSS2 protein is an important target of SARS-CoV-2 viral replication owing to its proteolytic activity. In the current study, a good homology model for TMPRSS2 has been built. E-pharmacophore model was developed and screened against one million compounds from enamine database, followed by molecular docking and MM-GBSA binding energy calculations to validate their affinity to the target. The top three molecules (Z126202570, Z46489368, and Z422255982) were also subjected to molecular dynamics simulations to confirm the stability of interactions obtained from docking analysis as well as density functional theory calculation. These results indicated that these compounds might be potential hits against TMPRSS2 and if experimentally confirmed may serve as future anti-covid agents.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.