key: cord-0550826-b5bui141 authors: Boldrini, Alberto; Terruzzi, Luca; Spagnolli, Giovanni; Astolfi, Andrea; Massignan, Tania; Lolli, Graziano; Barreca, Maria Letizia; Biasini, Emiliano; Faccioli, Pietro; Pieri, Lidia title: Identification of a Druggable Intermediate along the Folding Pathway of the SARS-CoV-2 Receptor ACE2 date: 2020-04-28 journal: nan DOI: nan sha: 42b5ab0401796d011f8a5d9016960abc7803bef1 doc_id: 550826 cord_uid: b5bui141 The steep climbing of victims caused by the new coronavirus disease 2019 (COVID-19) throughout the planet is sparking an unprecedented effort to identify effective therapeutic regimens to tackle the pandemic. The SARS-CoV-2 virus is known to gain entry into various cell types through the binding of one of its surface proteins (spike) to the host Angiotensin Converting Enzyme 2 (ACE2). Thus, the spike-ACE2 interaction represents a major target for vaccines and antiviral drugs. We recently described a novel method to pharmacologically down-regulate the expression of target proteins at the post-translational level. This technology builds on computational advancements in the simulation of folding mechanisms to rationally block protein expression by targeting folding intermediates, hence hampering the folding process. Here, we report the all-atom simulation of the entire sequence of events underlying the folding pathway of ACE2. Our data reveal the existence of a folding intermediate showing druggable pockets hidden in the native conformation. Both pockets were targeted by a virtual screening repurposing campaign aimed at quickly identifying drugs capable to decrease the expression of ACE2. Importantly, among the different virtual hits, we identified mefloquine, a quinoline-derivative belonging to a class of antimalaria agents (e.g. chloroquine and hydroxychloroquine) recently described for their effects on ACE2 maturation. Thus, our results suggest that these drugs could act against SARS-CoV-2 by altering the folding pathway of its receptor ACE2. Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is the causative agent of the viral pneumonia outbreak named Coronavirus Disease 2019 (COVID-19), a pandemic started at the end of 2019 in the Hubei province of China and later spread to the rest of the planet [1] . The disease shows unprecedented transmission, morbidity and mortality rates [2, 3] . With vaccines still far from being available, effective therapeutics to control the disease are in urgent need. SARS-CoV-2 belongs to the broad family of Coronaviridae, positive-sense, single-stranded RNA viruses of humans and animals causing pathologies ranging from a common cold to severe respiratory diseases, such as the Severe Acute Respiratory Syndrome (SARS) and the Middle East Respiratory Syndrome (MERS) [4] . Similarly to other coronaviruses, SARS-CoV-2 has four main structural proteins, known as the nucleocapsid (N), membrane (M), envelope (E) and spike (S) [5] . The latter promotes the first adhesion of the virus to host cells through direct interaction with extracellular domains of the angiotensin-converting enzyme 2 (ACE2) [6] [7] [8] . Subsequently, another protein of the host, the transmembrane protease serine 2 (TMPRSS2), cleaves the spike protein, exposing a fusion peptide that directly promotes virus entry into the cell [9] . ACE2 belongs to the protein family of angiotensin-converting enzymes [10] . The protein is reported to be expressed in cells of the kidney, intestine, arteries, heart, and lung, although discrepancies about its expression profiles in different organs recently emerged [11, 12] . The main function of ACE2 is connected to its ability to lower blood pressure by catalyzing the conversion of angiotensin II into the vasodilator angiotensin [13] . The protein is a promising drug target for treating cardiovascular diseases [14, 15] . Consistent with the mode of entry of SARS-CoV-2 into human cells, blocking the initial interaction of Spike with ACE2 should effectively counteract virus replication. This objective could be achieved by different strategies, including antibodies against Spike, compounds interfering with the Spike-ACE2 interaction, or molecules capable of lowering the exposure of ACE2 at the cell surface [16] . The latter strategy could overcome possible evolutionary variations of the virus conferring resistance to treatments directed against viral proteins. However, the function of ACE2 on vascular regulation has been shown to play protective roles against virus-induced lung injury, and the complete silencing of the protein would likely result in severe secondary effects [6, 17] . For this reason, any pharmacological approach aimed at targeting ACE2 should be able to down-regulate its expression rather than abrogating it. We have recently described an entirely novel method for selectively reducing the level of target proteins, called Pharmacological Protein Inactivation by Folding Intermediate Targeting (PPI-FIT) [18] . The technology is based on the concept of targeting folding intermediates of proteins rather than native conformations and is made possible by computational algorithms allowing the allatom reconstruction of folding pathways [18, 19] . The rationale underlying the PPI-FIT method is that stabilizing a folding intermediate of a protein with small ligands should promote its degradation by the cellular quality control machinery, which could recognize such artificially stabilized intermediates as improperly folded species. Here, we report the full atomistic reconstruction of the folding pathway of ACE2, which predicts the existence of a folding intermediate presenting two unique druggable pockets hidden into the native state. These pockets were employed as novel target sites for virtual screening campaigns aimed at repurposing compounds capable of modulating the expression of ACE2. To obtain the all-atom reconstruction of the folding pathway of apo-ACE2, we employed the so-called Bias Functional approach (see Methods) [20] . This enhanced path sampling scheme requires to provide the protein native conformation, which was retrieved from PDB 1R42 (Figure 1 ). This structure includes the atomic coordinates of the catalytic domain of the enzyme (residues 19-615). The zinc ion in the active site, as well as the glycans, were both removed, due to uncertainties about their binding kinetics. The Charmm36m forcefield [21] was employed throughout the simulations. The native structure of ACE2 was employed to generate 32 denatured conformations by a combination of unfolding ratchet-and-pawl Molecular Dynamics (rMD) and high temperature MD simulations (see Methods). For each conformation, 40 trajectories were generated by folding rMD, for a total of 1280 trajectories. All the non-productive trajectories were discarded while folding trajectories of the remaining sets (N = 18) were ranked based on the value of their Bias Functional (Equation 5 in Methods). For each set, we selected the least biased trajectory (LBT), corresponding to the folding pathway with the highest probability to occur in the absence of biasing forces [20] . To gain information about the folding mechanism, we considered a specific collective variable R defined as a normalized linear combination of the fraction of native contacts (Q) and Root Mean Squared Deviation (RMSD) of atomic positions from the native structure. According to our definition, R = 0 corresponds to the fully denatured state of the protein, while R = 1 to the native configuration. The 18 resulting least biased trajectories were employed to generate the transition path energy function G(R) defined in Methods. This quantity provides a lower bound to reaction limiting free-energy barriers that are overcome along the folding process. Therefore, it can be used to infer the existence of folding intermediates. Since the bias functional approach cannot be used to explore the dynamics within the reactant state, we restricted our analysis to the transition region R > 0.55, where significant tertiary-structure content begins to appear. This reactive section of the pathway contains two non-native metastable states, referred to as "early" and "late" intermediates ( Figure 2 ). The analysis of protein conformations along the trajectories led to the description of the complete folding mechanism for ACE2 ( Figure 3 ). The first event along the ACE2 folding pathway is the formation of three distinct regions (foldons) encompassing: (i) the two N-terminal α-helices (residues 19-82); (ii) the globular domain (residues 295-425) containing the zinc coordinating amino acids; (iii) the central core (residues 149-284 and 430-587), followed by the docking of the C-terminal tail (residues 589-615). These regions then sequentially fold onto each other, giving rise to the native state. A k-mean clustering was performed to obtain representative protein structures from the two non-native free energy wells (i.e. early and late intermediates) appearing along the folding pathways (Supp. Figure 1 ). Four cluster centers per well were selected as representative conformations ( Figure 4 ). The early intermediate is characterized by the lack of docking of the three foldons, which may also appear as not completely structured in some conformations. Conversely, the late ACE2 folding intermediate differs from the native conformation mainly for the topology of the two N-terminal α-helices, which are not docked to the rest of the structure in all the conformations. To account for minor variations among the structures belonging to the same cluster, we selected an additional protein conformation for each cluster of the late intermediate (Supp. Figure 2 ). This process yielded a final number of 8 protein conformations, which were then used for the subsequent pocket detection step. These structures were subjected to 50 ns of MD simulations at T = 310 K with position restrain on Cα, to explore the conformational variability of the residues side-chains. Then, for each MD simulation, 200 structures, equispaced in time, were extracted and analyzed with the SiteMap [22] and DogSiteScorer [23] software. Pockets were then ranked based on a series of druggability descriptors (Supp. Table 1 ). This analysis identified two predicted druggable pockets, exclusively present in the ACE2 late intermediate state but not in the native conformation ( Figure 5 ). In particular, pocket 1 results from a local structural variation into the core region (residues 468-498), while pocket 2 results from the missing docking of the C-terminal tail onto the core region. To determine the reliability of each druggable site, we measured the number of folding pathways presenting the relevant structural variation underlying each pocket ( Figure 6 ). The least biased trajectories were projected on two graphs plotting the RMSD of each relevant region (residues 468-498 or C-terminal tail) against the RMSD of the corresponding docking site. These analyses revealed that the pocket 1 is present in a single trajectory, while the pocket is predicted to appear in 9 different trajectories. The identification of potential ACE2 folding intermediate ligands was pursued by employing a drug repositioning strategy. We built a unique collection of 9187 compounds by combining libraries of drugs approved by the U.S. Food and Drug Administration (FDA) and molecules at different stages of currently ongoing clinical trials (see Material and Methods). The chemical collection was screened against the two identified pockets by following a consensus virtual screening workflow ( Figure 7A ). Two different docking software, Glide [22] and LeadIT [24] , were employed in parallel to predict the binding affinity of each compound to the ACE2 folding intermediate pockets.Two different docking software, Glide [22] and LeadIT [24] , were employed in parallel to predict the binding affinity of each compound to the ACE2 folding intermediate pockets. Only compounds showing promising predicted affinity (i.e. Glide ds ≤ -6 kcal/mol; LeadIT HYDE aff ≤ 50 µM) in both docking protocols were submitted to a third docking round based on AutoDock [25] . This process identified two consensus sets (AD LBE ≤ -6 kcal/mol, AD NiC ≥ 25), including 145 compounds for pocket 1 and 238 for pocket 2. The top scoring compounds from Glide (Glide ds ≤ -9 kcal/mol) and LeadIT (HYDE aff ≤ 5 µM) were also added to these sets. Finally, a visual inspection of binding mode and chemical similarity annotation for each ligand allowed the selection of 14 virtual hits for pocket 1 and additional 21 for pocket 2 (Supp. Table 2 ). Collectively, these results predicted 35 potential ligands for the ACE2 folding intermediate. Interestingly, among the selected hits, we identified mefloquine ( Figure 7B and 7C) , a quinoline analogue belonging to a class of anti-malaria agents recently reported to exert anti-SARS-CoV-2 effects in vitro. Of note, retrospective analysis of 7 additional members of the same chemical class (shown in Supp. Table 2) identified hydroxychloroquine, a drug currently in clinical trials for COVID-19, as a putative ligand for pocket 2 of the ACE2 folding intermediate. Multiple pieces of evidence indicate that down-regulating the expression of ACE2 in the early stages of the SARS-CoV2 infection should effectively inhibit virus replication. However, selectively decreasing the expression of a target protein could be a difficult task. RNA silencing or CRISPR-based strategies represent valid options, but their use could be limited by delivery issues [26, 27] . These problems could be overcome by emerging pharmacological technologies like the proteolysis targeting chimeras (PROTACs), which build on the principle of designing bi-functional compounds capable of interacting with the target protein with one side and engaging the E3 ubiquitin ligase with the other, leading to the degradation of the polypeptide by the proteasome [28] . Similarly, the PPI-FIT method capitalizes on the cellular quality control machinery to promote the degradation of the target polypeptide, although it does not require the development of bi-functional molecules. PPI-FIT-derived compounds aim at stimulating the removal of the target protein by directly blocking its folding pathway [18] . In this manuscript, we described the first step along with the application of the PPI-FIT paradigm to ACE2 by reporting the all-atom reconstruction of the folding pathway of the protein. Our analyses predict the existence of a folding intermediate showing two unique druggable pockets not present in the native ACE2 conformation. In order to respond to the urgent need for an effective therapy against SARS-CoV-2, we targeted both pockets by a in silico virtual screening approach aimed at repurposing compounds currently in clinical trials or already approved by the FDA. Interestingly, among the predicted ligands for pocket 2 we found the quinoline-derivative mefloquine. This drug belongs to a class of antimalaria agents recently described for their potential effect against SARS-CoV-2 [29, 30] . While the precise mechanism by which chloroquine and its more active derivative hydroxychloroquine inhibit virus replication is not known, reports suggest that the compounds may act by reducing the glycosylation of ACE2 [31] . This evidence is highly consistent with a PPI-FIT-like effect of these drugs, which could alter the correct maturation of ACE2 glycosidic moieties by acting on the folding pathway of the protein. Based on our previous experience with the PPI-FIT technology, we anticipate that the possible down-regulating effects of our in silico predicted compounds targeting the ACE2 folding intermediate will be easily tested in cell-based assays by mean of standard biochemical techniques like western blotting. The simulations were performed on the high-performance-computing facilities of the Italian Institute for Nuclear Physics (INFN) and on the computer cluster of Sibylla Biotech. A modified version of Gromacs 2019 [32] was employed to run the calculations. Data analysis was carried out using Python and its libraries: NumPy, SciPy, Matplotlib and MDAnalysis. The structure of ACE2 was retrieved from PDB 1R42. The file contains the catalytic domain (residues 19-615) solved by X-ray crystallography at 2.2 Å of resolution. Protein topology was generated using Charmm36m after removal of the catalytic zinc ion and sugar moieties. Cysteines were treated in the reduced state. Water molecules were modeled using the Charmm-modified TIP3P. The ACE2 structure was positioned in a cubic box with minimum wall distance equal to 60 Å. The system was filled with water molecules, neutralized with 28 Na + ions and brought to a final 150 mM NaCl concentration. This setup was followed by an energy minimization using the steepest descent algorithm. An NVT equilibration at 800 K was carried out using the V-rescale thermostat. Then, so-called ratchet and pawl (rMD) molecular dynamics simulations [33, 34] were employed to obtain the denatured states. In this method, an external biasing force is added, acting along a single collective variable, denoted as z(X): where C ij (X) is an entry of the instantaneous contact map and C ij (X R ) is an entry of the reference contact map, both defined as: where x i and x j are the coordinates of the i th and j th atoms of the protein; while r 0 is the reference contact-distance that is set to 7.5 Å. In applications of rMD to protein denaturation, all the entries of the reference contact map were set to 0. In applications of rMD to protein folding simulations, the reference contact map C ij (X R ) is calculated from Eq. (2) using the protein's native conformation, obtained by performing energy minimization starting from the PDB structure. In any rMD simulation, the system evolves according to plain MD as long as the value of z(x) spontaneously decreases during time. Instead, system fluctuations leading to an increase in z(X) result in the activation of the biasing force defined as: where z m (t) is the minimum value assumed by z(X) up to time t, and the index i indicates the atom on which the force is acting. To obtain the denatured protein conformations, 32 rMD trajectories of 1 ns each with a force constant k R = 8 · 10 −4 kJ/mol. These simulations were followed by 4 ns of plain MD at 800 K in the NVT ensemble. For each plain MD unfolding trajectory, the protein conformation minimizing the following function was selected: where n is the total number of atoms, m is the number of ignored subsequent atoms and r 0 is a parameter, set to 40 Å. This metrics is a norm of the contact maps, where atoms with distant indexes weight more in the calculation. Conformations with low value of S(X) are characterized by high degree of denaturation, vice-versa, conformations with high value of S(X) are less denatured. Selected initial conditions (N = 32) were positioned in a cubic box with minimum wall-distance of 10 Å. Each system was filled with water molecules, neutralized with 28 Na + ions and brought to a final 150 mM NaCl concentration. Energy minimization was subsequently performed using the steepest descent algorithm. Then, each system was subjected to 500 ps of NVT equilibration at 350 K (using the V-rescale thermostat [35] ), followed by 500 ps of NPT equilibration at 350 K and 1 Bar (using the V-rescale thermostat and the Parrinello-Rahman barostat). During equilibration, position restraints with force constant 1000 kJ·mol −1 nm −2 where introduced on heavy atoms. For each initial condition, 40 rMD trajectories were generated, each one consisting in 5 · 10 6 steps with integration time-step of 2 fs. In this set of rMD simulations, reference contact map entering Eq. 1 was calculated using the native structure of ACE2. Cutoff for Coulomb interactions was set to 12 Å, cutoff for Van der Waals interactions was set to 12 Å with force-switch having 10 Å radius. Long range electrostatics was treated with particle mesh Ewald. Trajectories reaching a configuration with an RMSD lower than 4 Å were considered productive folding events. For each set of folding trajectories (starting from the same initial condition) the pathway with the highest probability to occur in absence of external bias was identified by selecting the one minimizing the Bias Functional [20] , defined as: Where m i and γ i are the mass and the friction coefficient of the i th atom, F rM D i is the force acting on it, while t is the trajectory folding time. The collective variables Q and RMSD were linearly combined and normalized to obtain a one-dimensional reaction coordinate, R: From the frequency histogram of the least biased trajectories we estimated the probability P (R) to observe a given value of the collective variable. We refer to the quantity G(R) = −ln[P (R)] as to the "transition path energy". In a biased dynamics, G(R) yields a lower-bound to the rate limiting freeenergy barriers, thus provides a useful tool to identify metastable states. The standard deviation σ was estimated through a jackknifing procedure, in particular, the free energy profile was computed for 18 times by leaving a different trajectory out in each calculation. For each bin, the standard deviation σ was calculated as: where n is the total number of leas biased trajectories,θ i is the mean value of the jackknife replicate (where the ith trajectory was removed) andθ (·) is the mean of the jackknife replicates. Protein conformations were sampled from the two identified energy wells, in particular, the earlyintermediate was defined with 0.59 < R < 0.68, while the late-intermediate was defined with 0.797 < R < 0.865 (Supp. Figure 1) . The two sets of structures were then subjected to k-mean clustering (k = 4), using the Frobenius norm of the contact maps as a distance metrics. Then, for each cluster, a single representative conformation was selected by choosing the structure minimizing its distance from the relative cluster center. The 4 cluster centers of the late intermediate and 4 additional conformations (each one extracted as a cluster variant) were positioned in a cubic box with minimum wall-distance of 10 Å, that was filled with water molecules. The system was neutralized with 28 Na + ions and brought to a final 150 mM NaCl concentration. Energy minimization was subsequently performed, using the steepest descent algorithm. Then, each system was equilibrated for 500 ps in the NVT ensemble at 310 K (using the V-rescale thermostat), followed by 500 ps of NPT equilibration at 310 K and 1 Bar (using the V-rescale thermostat and the Parrinello-Rahman barostat). During equilibration, position restraints with force constant 1000 kJ·mol −1 nm −2 where introduced on heavy atoms. Subsequently, 50 ns of MD were performed for each structure by retaining the positional restraints on the Cα of the protein backbone. For each trajectory, 200 frames, equally spaced in time, were extracted. Each frame was analyzed by means of SiteMap [36] and DogSiteScorer [23] software in order to identify druggable pockets. We considered as druggable a pocket falling within the following thresholds: volume ≥ 300 Å 3 ; depth ≥ 10 Å; balance ≥ 1.0; exposure ≤ 0.5; enclosure ≥ 0.70; SiteScore ≥ 0.8; DScore ≥ 0.90; DrugScore ≥ 0.5; SimpleScore ≥ 0.5 (Supp. Figure 7A ). The FDA-approved drug library included compounds from the following chemical collections: Selleck (accession date 21/03/2020, 2684 compounds), Prestwick (accession date 21/03/2020, 1520 compounds) and eDrug-3D (accession date 21/03/2020, 1930 compounds). Compounds from the DrugBank collection (DrugBank Release Version 5.1.5; 1784 molecules), which includes approved, experimental and investigational drugs, were added. The four chemical libraries were merged using an in-house developed KNIME workflow to obtain a unique library of non-redundant entries (total of 9187 compounds). The final collection was prepared with the LigPrep tool [37] . Ionization/tautomeric states were generated at pH range 7-8 using Epik [38] . Furthermore, at most 32 stereoisomers per ligand and three lowest energy conformations per ligand ring were produced. Where not defined, all the chiral form of each stereocenter was produced. In total, 12771 docking clients were generated. For the Glide docking [22] , the N-and C-terminal residues of the ACE2 intermediate were capped with acetyl (ACE) and N-methylamine (NMA) groups, respectively, using the Schrödinger Protein Preparation Wizard. The capped protein structure was used to generate the receptor grid, with no scaling of Van der Waals radii for non-polar receptor atoms. The docking space was centered on the centroid of the residues defining the pockets (x: 75.4, y: 58.2, z: 67.6 for pocket 1; x: 50.8, y: 57.0, z: 69.1 for pocket 2). The docking space was defined as a 35 Å 3 cubic box, while the diameter midpoint of docked ligands was required to remain within a smaller, nested 15 Å 3 cubic box. Docking experiments were performed in the Glide standard precision mode using a 0.8 factor to scale the Van der Waals radii of the ligand atoms with partial atomic charge less than 0.15. For the BioSolveIT docking, LeadIT (version 3.2.0) was used for protein preparation and docking parameters definition. The binding site was defined on the basis of the residues composing the identified druggable pocket. The residue protonation states, as well as the tautomeric forms, were automatically assessed in LeadIT using the ProToss method, that generates the most probable hydrogen positions on the basis of an optimal hydrogen bonding network using an empirical scoring function. The virtual screening workflow was developed by using the KNIME analytic platform and the BioSolveIT KNIME nodes. Specifically, the workflow was organized as follows: (i) the Compute LeadIT Docking node was selected to perform the docking simulations of the ∼11·10 3 docking clients by using the FlexX algorithm [39] . Ten poses for each ligand were produced; (ii) the Assess Affinity with HYDE in SeeSAR node generated refined binding free energy and HYDE [40] predicted activity (HYDE aff ) for each ligand pose using the HYDE rescoring function; and (iii) for each ligand, the pose with the lowest HYDE aff was extracted. For AutoDock docking, ligands and receptor structures were converted to the AD4 format using AutoDockTools and the Gasteiger-Marsili empirical atomic partial charges were assigned. The dimensions of the grid were 60 X 60 X 60 points, with grid points separated by a 0.375 Å. The grid was centered on the centroids of pockets 1 and 2. The Lamarckian genetic algorithm was used, and for each compound, the docking simulation was composed of 100 runs. Clustering of docked conformations was performed on the basis of their RMSD (tolerance = 2.0 Å) and the results were ranked based on the estimated free energy of binding. The obtained poses were filtered based on the Glide docking score (Glide ds ), the HYDE aff and AutoDock ligand binding energy (AD LBE ) as well as on the number of clusters (AD NiC ). The consensus virtual screening workflow was applied for both pocket 1 and pocket 2 ( Figure 7A ). The best scored structures for the three software were visually inspected and selected based on their binding mode and predicted interactions. In addition, the top scoring compounds for Glide docking and SeeSAR rescoring [41] were also evaluated. In particular, the Glide set was generated following these criteria: Glide ds ≤ -9 kcal/mol, giving rise to 97 molecules for pocket 1 and 30 for pocket 2. The SeeSAR set was instead defined by applying the following filters: (i) HYDE aff ≤ 5 µM; (ii) Ligand efficiency (LE) and lipophilic ligand efficiency (LLE) category from 1 to 4 (where 1 corresponds to favorable LE/LLE, 5 corresponds to unfavorable LE/LLE); (iii) Good ligand conformer quality, judged on the basis of torsional quality of the rescored pose rotamers and intramolecular clashes; (iv) no protein-ligand intermolecular clashes, giving rise to 89 molecules for pocket 1 and 205 molecules for pocket 2. Atomic resolution structures of the folding intermediates together with the corresponding druggable pockets are freely available upon request. Supp. Figure 1 : Supp. Figure 2 : Selected structure variants for each cluster. Representative conformations of the additional structures selected for each cluster. Supp. The COVID-19 epidemic The epidemiology and pathogenesis of coronavirus disease (COVID-19) outbreak Early Transmission Dynamics in Wuhan, China, of Novel Coronavirus-Infected Pneumonia Human Coronaviruses: A Review of Virus-Host Interactions A new coronavirus associated with human respiratory disease in China Characterization of the receptorbinding domain (RBD) of 2019 novel coronavirus: implication for development of RBD protein as a viral attachment inhibitor and vaccine Structural basis of receptor recognition by SARS-CoV-2 Structure of the SARS-CoV-2 spike receptorbinding domain bound to the ACE2 receptor SARS-CoV-2 Cell Entry Depends on ACE2 and TMPRSS2 and Is Blocked by a Clinically Proven Protease Inhibitor ACE2, a new regulator of the renin-angiotensin system Angiotensin converting enzyme 2 is primarily epithelial and is developmentally regulated in the mouse lung SARS-CoV-2 receptor ACE2 and TMPRSS2 are primarily expressed in bronchial transient secretory cells The angiotensin-converting enzyme gene family: genomics and pharmacology Angiotensin-converting enzyme 2 is a key modulator of the renin-angiotensin system in cardiovascular and renal disease Targeting the ACE2-Ang-(1-7) pathway in cardiac fibroblasts to treat cardiac remodeling and heart failure A Review of SARS-CoV-2 and the Ongoing Clinical Trials Angiotensin-converting enzyme 2 (ACE2) as a SARS-CoV-2 receptor: molecular mechanisms and potential therapeutic target Pharmacological Protein Inactivation by Targeting Folding Intermediates. BioRXiv. 2020 Full atomistic model of prion structure and conversion Variational scheme to compute protein reaction pathways using atomistic force fields with explicit solvent CHARMM36m: an improved force field for folded and intrinsically disordered proteins DoGSiteScorer: a web server for automatic binding site prediction, analysis and druggability assessment LeadIT version 3.2.0. wwwbiosolveitde/LeadIT AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility Small interfering RNA from the lab discovery to patients' recovery CRISPR-Based Therapeutic Genome Editing: Strategies and In Vivo Delivery by AAV Vectors Small-Molecule PROTACS: New Approaches to Protein Degradation Identification of antiviral drug candidates against SARS-CoV-2 from FDA-approved drugs Hydroxychloroquine, a less toxic derivative of chloroquine, is effective in inhibiting SARS-CoV-2 infection in vitro New insights on the antiviral effects of chloroquine against coronavirus: what to expect for COVID-19? GROMACS: fast, flexible, and free Forced unfolding of fibronectin type 3 modules: an analysis by biased molecular dynamics simulations Hierarchy of folding and unfolding events of protein G, CI2, and ACBP from explicit-solvent simulations Canonical sampling through velocity rescaling Identifying and characterizing binding sites and assessing druggability Molecular designing and in silico evaluation of darunavir derivatives as anticancer agents Epik: a software program for pK( a ) prediction and protonation state generation for drug-like molecules A fast flexible docking method using an incremental construction algorithm Towards an integrated description of hydrogen bonding and dehydration: decreasing false positives in virtual screening with the HYDE scoring function SeeSAR Version 9.2. wwwbiosolveitde/SeeSAR The table shows the different hits emerging from the virtual screening as well as the mefloquine analogues analyzed subsequently The computational reconstruction of ACE2 folding pathway was performed thanks to the strong support of the high performance computational infrastructures made available by Italian National Institute of Nuclear Physics (INFN). This work was supported by a grant from Fondazione Telethon (Italy, TCP14009). GS is a recipient of a fellowship from Fondazione Telethon. EB is an Assistant Telethon Scientist at the Dulbecco Telethon Institute.