Molecular insight into the specific enzymatic properties of TREX1 revealing the diverse functions in processing RNA and DNA/RNA hybrids.

In various autoimmune diseases, dysfunctional TREX1 (Three prime Repair Exonuclease 1) leads to accumulation of endogenous single-stranded DNA (ssDNA), double-stranded DNA (dsDNA) and DNA/RNA hybrids in the cytoplasm and triggers immune activation through the cGAS-STING pathway. Although inhibition of TREX1 could be a useful strategy for cancer immunotherapy, profiling cellular functions in terms of its potential substrates is a key step. Particularly important is the functionality of processing DNA/RNA hybrids and RNA substrates. The exonuclease activity measurements conducted here establish that TREX1 can digest both ssRNA and DNA/RNA hybrids but not dsRNA. The newly solved structures of TREX1-RNA product and TREX1-nucleotide complexes show that 2'-OH does not impose steric hindrance or specific interactions for the recognition of RNA. Through all-atom molecular dynamics simulations, we illustrate that the 2'-OH-mediated intra-chain hydrogen bonding in RNA would affect the binding with TREX1 and thereby reduce the exonuclease activity. This notion of higher conformational rigidity in RNA leading TREX1 to exhibit weaker catalytic cleavage is further validated by the binding affinity measurements with various synthetic DNA-RNA junctions. The results of this work thus provide new insights into the mechanism by which TREX1 processes RNA and DNA/RNA hybrids and contribute to the molecular-level understanding of the complex cellular functions of TREX1 as an exonuclease.

Structural basis for water modulating RNA duplex formation in the CUG repeats of myotonic dystrophy type 1.

Secondary structures formed by expanded CUG RNA are involved in the pathobiology of myotonic dystrophy type 1. Understanding the molecular basis of toxic RNA structures can provide insights into the mechanism of disease pathogenesis and accelerate the drug discovery process. Here, we report the crystal structure of CUG repeat RNA containing three U-U mismatches between C-G and G-C base pairs. The CUG RNA crystallizes as an A-form duplex, with the first and third U-U mismatches adopting a water-mediated asymmetric mirror isoform geometry. We found for the first time that a symmetric, water-bridged U-H2O-U mismatch is well tolerated within the CUG RNA duplex, which was previously suspected but not observed. The new water-bridged U-U mismatch resulted in high base-pair opening and single-sided cross-strand stacking interactions, which in turn dominate the CUG RNA structure. Furthermore, we performed molecular dynamics simulations that complemented the structural findings and proposed that the first and third U-U mismatches are interchangeable conformations, while the central water-bridged U-U mismatch represents an intermediate state that modulates the RNA duplex conformation. Collectively, the new structural features provided in this work are important for understanding the recognition of U-U mismatches in CUG repeats by external ligands such as proteins or small molecules.

Edge weights in a protein elastic network reorganize collective motions and render long-range sensitivity responses.

The effects of inter-residue interactions on protein collective motions are analyzed by comparing two elastic network models (ENM)-structural contact ENM (SC-ENM) and molecular dynamics (MD)-ENM-with the edge weights computed from an all-atom MD trajectory by structure-mechanics statistical learning. A theoretical framework is devised to decompose the eigenvalues of ENM Hessian into contributions from individual springs and to compute the sensitivities of positional fluctuations and covariances to spring constant variation. Our linear perturbation approach quantifies the response mechanisms as softness modulation and orientation shift. All contacts of Cα positions in SC-ENM have an identical spring constant by fitting the profile of root-of-mean-squared-fluctuation calculated from an all-atom MD simulation, and the same trajectory data are also used to compute the specific spring constant of each contact as an MD-ENM edge weight. We illustrate that the soft-mode reorganization can be understood in terms of gaining weights along the structural contacts of low elastic strengths and loosing magnitude along those of high rigidities. With the diverse mechanical strengths encoded in protein dynamics, MD-ENM is found to have more pronounced long-range couplings and sensitivity responses with orientation shift identified as a key player in driving the specific residues to have high sensitivities. Furthermore, the responses of perturbing the springs of different residues are found to have asymmetry in the action-reaction relationship. In understanding the mutation effects on protein functional properties, such as long-range communications, our results point in the directions of collective motions as a major effector.

Structural basis for overhang excision and terminal unwinding of DNA duplexes by TREX1.

Three prime repair exonuclease 1 (TREX1) is an essential exonuclease in mammalian cells, and numerous in vivo and in vitro data evidenced its participation in immunity regulation and in genotoxicity remediation. In these very complicated cellular functions, the molecular mechanisms by which duplex DNA substrates are processed are mostly elusive because of the lack of structure information. Here, we report multiple crystal structures of TREX1 complexed with various substrates to provide the structure basis for overhang excision and terminal unwinding of DNA duplexes. The substrates were designed to mimic the intermediate structural DNAs involved in various repair pathways. The results showed that the Leu24-Pro25-Ser26 cluster of TREX1 served to cap the nonscissile 5'-end of the DNA for precise removal of the short 3'-overhang in L- and Y-structural DNA or to wedge into the double-stranded region for further digestion along the duplex. Biochemical assays were also conducted to demonstrate that TREX1 can indeed degrade double-stranded DNA (dsDNA) to a full extent. Overall, this study provided unprecedented knowledge at the molecular level on the enzymatic substrate processing involved in prevention of immune activation and in responses to genotoxic stresses. For example, Arg128, whose mutation in TREX1 was linked to a disease state, were shown to exhibit consistent interaction patterns with the nonscissile strand in all of the structures we solved. Such structure basis is expected to play an indispensable role in elucidating the functional activities of TREX1 at the cellular level and in vivo.

Basal leakage in oscillation: Coupled transcriptional and translational control using feed-forward loops.

The circadian clock is a complex system that plays many important roles in most organisms. Previously, many mathematical models have been used to sharpen our understanding of the Arabidopsis clock, which brought to light the roles of each transcriptional and post-translational regulations. However, the presence of both regulations, instead of either transcription or post-translation, raised curiosity of whether the combination of these two regulations is important for the clock's system. In this study, we built a series of simplified oscillators with different regulations to study the importance of post-translational regulation (specifically, 26S proteasome degradation) in the clock system. We found that a simple transcriptional-based oscillator can already generate sustained oscillation, but the oscillation can be easily destroyed in the presence of transcriptional leakage. Coupling post-translational control with transcriptional-based oscillator in a feed-forward loop will greatly improve the robustness of the oscillator in the presence of basal leakage. Using these general models, we were able to replicate the increased variability observed in the E3 ligase mutant for both plant and mammalian clocks. With this insight, we also predict a plausible regulator of several E3 ligase genes in the plant's clock. Thus, our results provide insights into and the plausible importance in coupling transcription and post-translation controls in the clock system.

Compound Molecular Logic in Accessing the Active Site of Mycobacterium tuberculosis Protein Tyrosine Phosphatase B.

Protein tyrosine phosphatase B (PtpB) from Mycobacterium tuberculosis (Mtb) extends the bacteria's survival in hosts and hence is a potential target for Mtb-specific drugs. To study how Mtb-specific sequence insertions in PtpB may regulate access to its active site through large-amplitude conformational changes, we performed free-energy calculations using an all-atom explicit solvent model. Corroborated by biochemical assays, the results show that PtpB's active site is controlled via an "either/or" compound conformational gating mechanism, an unexpected discovery that Mtb has evolved to bestow a single enzyme with such intricate logical operations. In addition to providing unprecedented insights for its active-site surroundings, the findings also suggest new ways of inactivating PtpB.

Importance of Internal Porosity for Glucan Adsorption in Mesoporous Carbon Materials.

To better understand the adsorption of long-chain poly(1 → 4)-ß-D-glucans on carbon surfaces as well as interactions responsible for this adsorption, we use a comparative study involving mesoporous carbon-silica composite materials that have been etched to varying degrees and all-atom molecular dynamics simulations. The materials synthesized as part of this etching study consist of an as-synthesized composite material (MCN-MSN), MCN-MSN-0.5 (composite materials consisting of 50% carbon by mass), MCN-MSN-0.3 (composite materials consisting of 70% carbon by mass), and MCN, in which silica etching was conducted using an aqueous ethanolic solution of either NaOH or HF. Data for the adsorption of long-chain glucans to these materials from concentrated aqueous HCl (37 wt %) solution demonstrate a direct relationship between the amount of ß-glu adsorption and the magnitude of exposed carbon mesopore surface area, which systematically increases and is also accompanied by an increase in the mesopore size during silica etching. This demonstrates ß-glu adsorption as occurring on internal carbon mesopores rather than exclusively on the external carbon surface. These experimental data on adsorption were corroborated by molecular dynamics (MD) simulations of ß-glu adsorption to a graphene bilayer separated by a distance of 3.2 nm, chosen to correspond to the carbon mesopore diameter of the experimental system. Simulation results using a variety of ß-glu solvent systems demonstrate the rapid adsorption of a ß-glu strand on the graphitic carbon surface via axial coupling and are consistent with experimentally observed trends in fast adsorption kinetics. Solvent-mediated effects such as small-scale hydrophobicity and preferential interactions with ions are shown to play important roles in modulating glucan adsorption to carbon surfaces, whereas experimental data on hydrophobically modified silica demonstrate that hydrophobicity in and of itself is insufficient to cause ß-glu adsorption from concentrated aqueous HCl solution.

Systems-level modeling with molecular resolution elucidates the rate-limiting mechanisms of cellulose decomposition by cellobiohydrolases.

Interprotein and enzyme-substrate couplings in interfacial biocatalysis induce spatial correlations beyond the capabilities of classical mass-action principles in modeling reaction kinetics. To understand the impact of spatial constraints on enzyme kinetics, we developed a computational scheme to simulate the reaction network of enzymes with the structures of individual proteins and substrate molecules explicitly resolved in the three-dimensional space. This methodology was applied to elucidate the rate-limiting mechanisms of crystalline cellulose decomposition by cellobiohydrolases. We illustrate that the primary bottlenecks are slow complexation of glucan chains into the enzyme active site and excessive enzyme jamming along the crowded substrate. Jamming could be alleviated by increasing the decomplexation rate constant but at the expense of reduced processivity. We demonstrate that enhancing the apparent reaction rate required a subtle balance between accelerating the complexation driving force and simultaneously avoiding enzyme jamming. Via a spatiotemporal systems analysis, we developed a unified mechanistic framework that delineates the experimental conditions under which different sets of rate-limiting behaviors emerge. We found that optimization of the complexation-exchange kinetics is critical for overcoming the barriers imposed by interfacial confinement and accelerating the apparent rate of enzymatic cellulose decomposition.

Dynamic mesoscale model of dipolar fluids via fluctuating hydrodynamics.

Fluctuating hydrodynamics (FHD) is a general framework of mesoscopic modeling and simulation based on conservational laws and constitutive equations of linear and nonlinear responses. However, explicit representation of electrical forces in FHD has yet to appear. In this work, we devised an Ansatz for the dynamics of dipole moment densities that is linked with the Poisson equation of the electrical potential ϕ in coupling to the other equations of FHD. The resulting ϕ-FHD equations then serve as a platform for integrating the essential forces, including electrostatics in addition to hydrodynamics, pressure-volume equation of state, surface tension, and solvent-particle interactions that govern the emergent behaviors of molecular systems at an intermediate scale. This unique merit of ϕ-FHD is illustrated by showing that the water dielectric function and ion hydration free energies in homogeneous and heterogenous systems can be captured accurately via the mesoscopic simulation. Furthermore, we show that the field variables of ϕ-FHD can be mapped from the trajectory of an all-atom molecular dynamics simulation such that model development and parametrization can be based on the information obtained at a finer-grained scale. With the aforementioned multiscale capabilities and a spatial resolution as high as 5 Å, the ϕ-FHD equations represent a useful semi-explicit solvent model for the modeling and simulation of complex systems, such as biomolecular machines and nanofluidics.

Trajectory Statistical Learning of the Potential Mean of Force and Diffusion Coefficient from Molecular Dynamics Simulations.

Central to studying the conformational changes of a complex protein is understanding the dynamics and energetics involved. Phenomenologically, structural dynamics can be formulated using an overdamped Langevin model along an observable, e.g., the distance between two residues in the protein. The Langevin model is specified by the deterministic force (the potential of mean force, PMF) and stochastic force (characterized by the diffusion coefficient, D). It is therefore of great interest to be able to extract both PMF and D from an observable time series but under the same computational framework. Here, we approach this challenge in molecular dynamics (MD) simulations by treating it as a missing-data Bayesian estimation problem. An important distinction in our methodology is that the entire MD trajectory, as opposed to the individual data elements, is used as the statistical variable in Bayesian imputation. This idea is implemented through an eigen-decomposition procedure for a time-symmetrized Fokker-Planck equation, followed by maximizing the likelihood for parameter estimation. The mathematical expressions for the functional derivatives used in learning PMF and D also provide new physical insights for the manner by which the information on both the deterministic and stochastic forces is encoded in the dynamics data. An all-atom MD simulation of a nontrivial biomolecule case is used to illustrate the application of this approach. We show that, interestingly, the results of trajectory statistical learning can motivate new order parameters for an improved description of the kinetic bottlenecks in conformational changes. Complementing purely data-driven or black-box methods, this work underscores the advantages of physics-based machine learning in gaining chemical insights from quantitative parameter estimation.

Subdomain dynamics enable chemical chain reactions in non-ribosomal peptide synthetases.

Many peptide-derived natural products are produced by non-ribosomal peptide synthetases (NRPSs) in an assembly-line fashion. Each amino acid is coupled to a designated peptidyl carrier protein (PCP) through two distinct reactions catalysed sequentially by the single active site of the adenylation domain (A-domain). Accumulating evidence suggests that large-amplitude structural changes occur in different NRPS states; yet how these molecular machines orchestrate such biochemical sequences has remained elusive. Here, using single-molecule Förster resonance energy transfer, we show that the A-domain of gramicidin S synthetase I adopts structurally extended and functionally obligatory conformations for alternating between adenylation and thioester-formation structures during enzymatic cycles. Complementary biochemical, computational and small-angle X-ray scattering studies reveal interconversion among these three conformations as intrinsic and hierarchical where intra-A-domain organizations propagate to remodel inter-A-PCP didomain configurations during catalysis. The tight kinetic coupling between structural transitions and enzymatic transformations is quantified, and how the gramicidin S synthetase I A-domain utilizes its inherent conformational dynamics to drive directional biosynthesis with a flexibly linked PCP domain is revealed.

Fisher information metric for the Langevin equation and least informative models of continuous stochastic dynamics.

The evaluation of the Fisher information matrix for the probability density of trajectories generated by the over-damped Langevin dynamics at equilibrium is presented. The framework we developed is general and applicable to any arbitrary potential of mean force where the parameter set is now the full space dependent function. Leveraging an innovative Hermitian form of the corresponding Fokker-Planck equation allows for an eigenbasis decomposition of the time propagation probability density. This formulation motivates the use of the square root of the equilibrium probability density as the basis for evaluating the Fisher information of trajectories with the essential advantage that the Fisher information matrix in the specified parameter space is constant. This outcome greatly eases the calculation of information content in the parameter space via a line integral. In the continuum limit, a simple analytical form can be derived to explicitly reveal the physical origin of the information content in equilibrium trajectories. This methodology also allows deduction of least informative dynamics models from known or available observables that are either dynamical or static in nature. The minimum information optimization of dynamics is performed for a set of different constraints to illustrate the generality of the proposed methodology.

Statistical learning of protein elastic network from positional covariance matrix.

Positional fluctuation and covariance during protein dynamics are key observables for understanding the molecular origin of biological functions. A frequently employed potential energy function for describing protein structural variation at the coarse-gained level is elastic network model (ENM). A long-standing issue in biomolecular simulation is thus the parametrization of ENM spring constants from the components of positional covariance matrix (PCM). Based on sensitivity analysis of PCM, the direct-coupling statistics of each spring, which is a specific combination of position fluctuation and covariance, is found to exhibit prominent signal of parameter dependence. This finding provides the basis for devising the objective function and the scheme of running through the effective one-dimensional optimization of every spring by self-consistent iteration. Formal derivation of the positional covariance statistical learning (PCSL) method also motivates the necessary data regularization for stable calculations. Robust convergence of PCSL is achieved in taking an all-atom molecular dynamics trajectory or an ensemble of homologous structures as input data. The PCSL framework can also be generalized with mixed objective functions to capture specific property such as the residue flexibility profile. Such physical chemistry-based statistical learning thus provides a useful platform for integrating the mechanical information encoded in various experimental or computational data.

Mechanical codes of chemical-scale specificity in DNA motifs.

In gene transcription, certain sequences of double-stranded (ds)DNA play a vital role in nucleosome positioning and expression initiation. That dsDNA is deformed to various extents in these processes leads us to ask: Could the genomic DNA also have sequence specificity in its chemical-scale mechanical properties? We approach this question using statistical machine learning to determine the rigidity between DNA chemical moieties. What emerges for the polyA, polyG, TpA, and CpG sequences studied here is a unique trigram that contains the quantitative mechanical strengths between bases and along the backbone. In a way, such a sequence-dependent trigram could be viewed as a DNA mechanical code. Interestingly, we discover a compensatory competition between the axial base-stacking interaction and the transverse base-pairing interaction, and such a reciprocal relationship constitutes the most discriminating feature of the mechanical code. Our results also provide chemical-scale understanding for experimental observables. For example, the long polyA persistence length is shown to have strong base stacking while its complement (polyAc) exhibits high backbone rigidity. The mechanical code concept enables a direct reading of the physical interactions encoded in the sequence which, with further development, is expected to shed new light on DNA allostery and DNA-binding drugs.

"Fluctuograms" reveal the intermittent intra-protein communication in subtilisin Carlsberg and correlate mechanical coupling with co-evolution.

The mechanism of intra-protein communication and allosteric coupling is key to understanding the structure-property relationship of protein function. For subtilisin Carlsberg, the Ca²+-binding loop is distal to substrate-binding and active sites, yet the serine protease function depends on Ca²+ binding. The atomic molecular dynamics (MD) simulations of apo and Ca²+-bound subtilisin show similar structures and there is no direct evidence that subtilisin has alternative conformations. To model the intra-protein communication due to Ca²+ binding, we transform the sequential segments of an atomic MD trajectory into separate elastic network models to represent anharmonicity and nonlinearity effectively as the temporal and spatial variation of the mechanical coupling network. In analogy to the spectrogram of sound waves, this transformation is termed the "fluctuogram" of protein dynamics. We illustrate that the Ca²+-bound and apo states of subtilisin have different fluctuograms and that intra-protein communication proceeds intermittently both in space and in time. We found that residues with large mechanical coupling variation due to Ca²+ binding correlate with the reported mutation sites selected by directed evolution for improving the stability of subtilisin and its activity in a non-aqueous environment. Furthermore, we utilize the fluctuograms calculated from MD to capture the highly correlated residues in a multiple sequence alignment. We show that in addition to the magnitude, the variance of coupling strength is also an indicative property for the sequence correlation observed in a statistical coupling analysis. The results of this work illustrate that the mechanical coupling networks calculated from atomic details can be used to correlate with functionally important mutation sites and co-evolution.

Entropy of cellulose dissolution in water and in the ionic liquid 1-butyl-3-methylimidazolim chloride.

The entropic driving forces of cellulose dissolution in water and in the ionic liquid 1-butyl-3-methylimidazolium chloride (BmimCl) are investigated via molecular dynamics simulations and the two-phase thermodynamic model. An atomistic model of cellulose was simulated at a dissociated state and a microfibril state to represent dissolution. The calculated values of entropy and internal energy changes between the two states inform the interplay of energetic and entropic driving forces in cellulose dissolution. In both water and BmimCl, we found that the entropy associated with the solvent degrees of freedom (DOF) decreases upon cellulose dissolution. However, solvent entropy reduction in BmimCl is much smaller than that in water and counteracts the entropy gain from the solute DOF to a much lesser extent. Solvent entropy reduction in water also plays a major role in making the free energy change of cellulose dissolution unfavorable at room temperature. In BmimCl, the interaction energies between solvent molecules and glucan chains and the total entropy change both contribute favorably to the dissolution free energy of cellulose. Calculations at different temperatures in the two solvents indicate that changes in total internal energy are a good indicator of the sign of the free energy change of cellulose dissolution.

Fluctuating hydrodynamics for multiscale modeling and simulation: energy and heat transfer in molecular fluids.

This work illustrates that fluctuating hydrodynamics (FHD) simulations can be used to capture the thermodynamic and hydrodynamic responses of molecular fluids at the nanoscale, including those associated with energy and heat transfer. Using all-atom molecular dynamics (MD) trajectories as the reference data, the atomistic coordinates of each snapshot are mapped onto mass, momentum, and energy density fields on Eulerian grids to generate a corresponding field trajectory. The molecular length-scale associated with finite molecule size is explicitly imposed during this coarse-graining by requiring that the variances of density fields scale inversely with the grid volume. From the fluctuations of field variables, the response functions and transport coefficients encoded in the all-atom MD trajectory are computed. By using the extracted fluid properties in FHD simulations, we show that the fluctuations and relaxation of hydrodynamic fields quantitatively match with those observed in the reference all-atom MD trajectory, hence establishing compatibility between the atomistic and field representations. We also show that inclusion of energy transfer in the FHD equations can more accurately capture the thermodynamic and hydrodynamic responses of molecular fluids. The results indicate that the proposed MD-to-FHD mapping with explicit consideration of finite molecule size provides a robust framework for coarse-graining the solution phase of complex molecular systems.

Structural distributions from single-molecule measurements as a tool for molecular mechanics.

A mechanical view provides an attractive alternative for predicting the behavior of complex systems since it circumvents the resource-intensive requirements of atomistic models; however, it remains extremely challenging to characterize the mechanical responses of a system at the molecular level. Here, the structural distribution is proposed to be an effective means to extracting the molecular mechanical properties. End-to-end distance distributions for a series of short poly-L-proline peptides with the sequence P(n)CG(3)K-biotin (n = 8, 12, 15 and 24) were used to experimentally illustrate this new approach. High-resolution single-molecule Förster-type resonance energy transfer (FRET) experiments were carried out and the conformation-resolving power was characterized and discussed in the context of the conventional constant-time binning procedure for FRET data analysis. It was shown that the commonly adopted theoretical polymer models-including the worm-like chain, the freely jointed chain, and the self-avoiding chain-could not be distinguished by the averaged end-to-end distances, but could be ruled out using the molecular details gained by conformational distribution analysis because similar polymers of different sizes could respond to external forces differently. Specifically, by fitting the molecular conformational distribution to a semi-flexible polymer model, the effective persistence lengths for the series of short poly-L-proline peptides were found to be size-dependent with values of ~190 Å, ~67 Å, ~51 Å, and ~76 Å for n = 8, 12, 15, and 24, respectively. A comprehensive computational modeling was carried out to gain further insights for this surprising discovery. It was found that P(8) exists as the extended all-trans isomaer whereas P(12) and P(15) predominantly contained one proline residue in the cis conformation. P(24) exists as a mixture of one-cis (75%) and two-cis (25%) isomers where each isomer contributes to an experimentally resolvable conformational mode. This work demonstrates the resolving power of the distribution-based approach, and the capacity of integrating high-resolution single-molecule FRET experiments with molecular modeling to reveal detailed structural information about the conformation of molecules on the length scales relevant to the study of biological molecules.

Structure-mechanics statistical learning uncovers mechanical relay in proteins.

A protein's adaptive response to its substrates is one of the key questions driving molecular physics and physical chemistry. This work employs the recently developed structure-mechanics statistical learning method to establish a mechanical perspective. Specifically, by mapping all-atom molecular dynamics simulations onto the spring parameters of a backbone-side-chain elastic network model, the chemical moiety specific force constants (or mechanical rigidity) are used to assemble the rigidity graph, which is the matrix of inter-residue coupling strength. Using the S1A protease and the PDZ3 signaling domain as examples, chains of spatially contiguous residues are found to exhibit prominent changes in their mechanical rigidity upon substrate binding or dissociation. Such a mechanical-relay picture thus provides a mechanistic underpinning for conformational changes, long-range communication, and inter-domain allostery in both proteins, where the responsive mechanical hotspots are mostly residues having important biological functions or significant mutation sensitivity.

Methotrexate inhibition of SARS-CoV-2 entry, infection and inflammation revealed by bioinformatics approach and a hamster model.

Background: Drug repurposing is a fast and effective way to develop drugs for an emerging disease such as COVID-19. The main challenges of effective drug repurposing are the discoveries of the right therapeutic targets and the right drugs for combating the disease. Methods: Here, we present a systematic repurposing approach, combining Homopharma and hierarchal systems biology networks (HiSBiN), to predict 327 therapeutic targets and 21,233 drug-target interactions of 1,592 FDA drugs for COVID-19. Among these multi-target drugs, eight candidates (along with pimozide and valsartan) were tested and methotrexate was identified to affect 14 therapeutic targets suppressing SARS-CoV-2 entry, viral replication, and COVID-19 pathologies. Through the use of in vitro (EC50 = 0.4 µM) and in vivo models, we show that methotrexate is able to inhibit COVID-19 via multiple mechanisms. Results: Our in vitro studies illustrate that methotrexate can suppress SARS-CoV-2 entry and replication by targeting furin and DHFR of the host, respectively. Additionally, methotrexate inhibits all four SARS-CoV-2 variants of concern. In a Syrian hamster model for COVID-19, methotrexate reduced virus replication, inflammation in the infected lungs. By analysis of transcriptomic analysis of collected samples from hamster lung, we uncovered that neutrophil infiltration and the pathways of innate immune response, adaptive immune response and thrombosis are modulated in the treated animals. Conclusions: We demonstrate that this systematic repurposing approach is potentially useful to identify pharmaceutical targets, multi-target drugs and regulated pathways for a complex disease. Our findings indicate that methotrexate is established as a promising drug against SARS-CoV-2 variants and can be used to treat lung damage and inflammation in COVID-19, warranting future evaluation in clinical trials.