Computational Analysis of Viral Protease Inhibitors and Drug Resistance

Introduction

Viral proteases represent one of the most extensively validated targets for antiviral intervention across multiple viral families. These enzymes catalyze the proteolytic processing of viral polyproteins into functional subunits, a step essential for viral replication and maturation. The development of protease inhibitors (PIs) has been a cornerstone of antiviral therapy, yet the emergence of drug resistance mutations (DRMs) within protease active sites poses a persistent challenge to therapeutic efficacy. Computational approaches, including molecular docking, molecular dynamics (MD) simulations, quantum mechanical calculations, and machine learning, have become indispensable tools for characterizing inhibitor binding mechanisms and predicting resistance profiles [1, 2, 3]. This review examines the computational analysis of viral protease inhibitors and drug resistance, with emphasis on active site rotamer modifications, structural overlays of resistant mutants, and the biophysical principles underlying inhibitor escape.

Structural Biology of Viral Proteases

Viral proteases belong to several mechanistic classes, including aspartic proteases (e.g., HIV-1 protease), serine proteases (e.g., HCV NS3/4A, dengue NS2B-NS3), and cysteine proteases (e.g., SARS-CoV-2 main protease Mpro and papain-like protease PLpro) [4, 5, 6]. Despite divergent evolutionary origins, these enzymes share a conserved functional architecture centered on a catalytic dyad or triad positioned within a substrate-binding cleft. The active site geometry, including the arrangement of backbone atoms and side chain rotamers, determines substrate specificity and inhibitor susceptibility [7, 8].

HIV-1 protease is a homodimeric aspartic protease with a C2-symmetric active site. Each monomer contributes one catalytic aspartate residue (Asp25 and Asp25') to form the catalytic center. The flap regions (residues 45-55) undergo conformational rearrangements upon substrate or inhibitor binding, transitioning from an open to a closed state [7, 9]. The active site cavity accommodates peptidomimetic inhibitors that mimic the natural cleavage site of the Gag-Pol polyprotein. Residue-level affinity decomposition using quantum electron density approaches has revealed that hydrogen bonding networks and hydrophobic packing within the S1/S1' and S2/S2' subsites are critical determinants of binding affinity [7].

SARS-CoV-2 Mpro (also termed 3CLpro) is a cysteine protease with a chymotrypsin-like fold. The catalytic dyad consists of Cys145 and His41. The substrate-binding cleft is divided into subsites S1 through S4, which accommodate the P1 through P4 side chains of the substrate or inhibitor [5, 10]. The S1 subsite exhibits a strong preference for glutamine, a feature exploited by covalent inhibitors such as nirmatrelvir and ensitrelvir [5, 11]. The S2 subsite accommodates hydrophobic residues, while the S4 subsite is more solvent-exposed [10, 12].

HCV NS3/4A protease is a serine protease with a chymotrypsin-like fold. The catalytic triad comprises Ser139, His57, and Asp81. The NS4A cofactor is essential for structural stability and catalytic activity [13]. The substrate-binding groove recognizes a consensus sequence with a cysteine or threonine at the P1 position. Inhibitors such as simeprevir and grazoprevir occupy the active site through extensive hydrogen bonding and hydrophobic contacts [13].

Dengue virus NS2B-NS3 protease is a serine protease that requires the NS2B cofactor for proper folding and activity. The active site contains the catalytic triad Ser135, His51, and Asp75 [9, 14, 15]. The substrate-binding pocket is relatively shallow and electrostatically polarized, presenting challenges for inhibitor design [9, 14].

Mechanisms of Protease Inhibitor Binding

Protease inhibitors can be classified as covalent or noncovalent, reversible or irreversible, and peptidomimetic or non-peptidomimetic. Covalent inhibitors form a stable chemical bond with the catalytic nucleophile, typically a cysteine or serine residue. Noncovalent inhibitors rely on hydrogen bonding, van der Waals interactions, and desolvation effects to achieve high affinity [3, 16, 17].

For SARS-CoV-2 Mpro, covalent inhibitors such as nirmatrelvir and ensitrelvir contain an electrophilic warhead (e.g., nitrile or aldehyde) that reacts with the thiol group of Cys145 [5, 11]. The binding mechanism involves initial noncovalent recognition followed by covalent bond formation. The stereochemistry of the warhead and the adjacent beta-lactone moiety in certain inhibitor classes critically influences reactivity and selectivity [1]. Deep learning approaches have identified novel covalent inhibitor fragments by screening virtual libraries against the Mpro active site, demonstrating that target-specific neural networks can predict warhead reactivity and binding pose simultaneously [3].

Noncovalent inhibition of SARS-CoV-2 PLpro has been investigated through multiscale computational studies combining classical MD simulations with quantum mechanics/molecular mechanics (QM/MM) calculations [16]. These studies revealed that noncovalent inhibitors occupy the BL2 loop region and the S4 subsite, inducing conformational changes that stabilize the enzyme in an inactive state [16]. Similarly, benzopyran derivatives have been evaluated as inhibitors of Chikungunya virus nsP2 and nsP4 proteases through integrated docking, MD, and density functional theory (DFT) studies, highlighting the importance of electrostatic complementarity and solvation effects [2].

For HIV-1 protease, peptidomimetic inhibitors such as darunavir and ritonavir bind in an extended conformation within the active site cavity. The central hydroxyl group forms hydrogen bonds with the catalytic aspartates, while the P1/P1' and P2/P2' moieties occupy the corresponding subsites [7, 10]. Co-amorphous systems containing ritonavir and coformers have been studied using molecular simulations to understand the physical stability of drug formulations, though the focus here is on the molecular interactions at the protease active site [10].

Drug Resistance Mutations: Structural and Computational Perspectives

Drug resistance mutations arise from selective pressure exerted by inhibitor exposure. These mutations can be classified as active site mutations, which directly alter inhibitor binding contacts, or non-active site mutations, which affect conformational dynamics, dimer stability, or flap flexibility [18, 29, 35]. Computational methods are essential for predicting the impact of mutations on inhibitor binding and for forecasting resistance evolution.

HIV-1 Protease Resistance

HIV-1 protease resistance is well-characterized, with over 50 DRMs identified across the protease coding region. Major mutations include D30N, I50V, V82A, I84V, and L90M [18, 29]. These mutations reduce inhibitor affinity while maintaining sufficient catalytic activity for viral replication. The D30N mutation, selected by nelfinavir, disrupts a hydrogen bond network between the inhibitor and the backbone of Asp30 [18]. The I50V mutation, associated with amprenavir resistance, alters the conformation of the flap region, reducing the closed-state stability [18].

Machine learning models have been benchmarked for predicting HIV-1 protease inhibitor resistance, with performance heavily dependent on data set construction and feature representation [29]. Sequence-based features, structural features (e.g., residue depth, solvent accessibility), and physicochemical descriptors all contribute to predictive accuracy. Random forest and gradient boosting methods have shown superior performance compared to logistic regression or support vector machines when trained on large phenotypic resistance data sets [29].

Forecasting drug-resistant HIV protease evolution requires integrating fitness landscapes with mutational pathways. Aggarwal and Periwal developed a computational framework that combines sequence co-variation analysis with structural stability calculations to predict the most probable resistance trajectories under inhibitor pressure [18]. This approach identified that resistance mutations often follow a specific order, with primary mutations reducing inhibitor binding and compensatory mutations restoring catalytic efficiency or structural stability [18].

The HMM-SA structural alphabet method has been applied to HIV-2 protease, revealing that resistance mutations induce local structural rearrangements detectable through hidden Markov model analysis of backbone torsion angles [35]. This method captures subtle conformational shifts that may not be apparent from simple sequence alignment or static crystal structures [35].

SARS-CoV-2 Mpro Resistance

The emergence of SARS-CoV-2 Mpro inhibitors has been accompanied by the identification of resistance mutations in cell culture and clinical isolates. The T21I and E166A mutations confer differential resistance to simnotrelvir, bofutrelvir, and ensitrelvir [5]. The E166A mutation eliminates a key hydrogen bond between the inhibitor and the Glu166 side chain, which is critical for stabilizing the inhibitor in the active site. The T21I mutation, located near the S2 subsite, alters the hydrophobic packing and reduces inhibitor binding affinity through steric effects [5].

Cross-resistance patterns between nirmatrelvir and ensitrelvir have been characterized, revealing asymmetrical resistance profiles [11]. Mutations such as S144A and S144G confer resistance to nirmatrelvir but remain sensitive to ensitrelvir, while mutations at position E166 confer resistance to both inhibitors [11]. Structural analysis shows that the S144A mutation disrupts a hydrogen bond with the nitrile warhead of nirmatrelvir, whereas ensitrelvir's larger P2 moiety compensates through additional hydrophobic contacts [11].

Multi-scale modeling of the WU-04 inhibitor has identified mutation-induced resistance mechanisms involving both direct contact changes and long-range conformational effects [12]. The L50F mutation, located distal to the active site, alters the dynamics of the S3-S4 loop, indirectly affecting inhibitor binding through an allosteric mechanism [12]. This finding underscores the importance of considering non-active site mutations in resistance prediction.

HCV NS3/4A Protease Resistance

HCV NS3/4A protease resistance mutations include R155K, A156T, D168V, and V170A [13]. These mutations reduce inhibitor binding by altering the electrostatic environment of the active site or by disrupting specific hydrogen bonds. The R155K mutation eliminates a salt bridge with the inhibitor's P2 moiety, while the D168V mutation introduces a bulky side chain that sterically clashes with the inhibitor [13]. Computational exploration of flavonoids as HCV NS3/4A inhibitors has identified natural product scaffolds that may circumvent common resistance mutations by engaging alternative binding interactions [13].

Dengue NS2B-NS3 Protease Resistance

Resistance to dengue NS2B-NS3 protease inhibitors is less well-characterized clinically, but computational studies have identified potential resistance hotspots. Molecular dynamics-driven optimization of triterpenoid, amidinium, and flavonoid inhibitors has revealed that mutations at positions 129 and 131 in the NS3 domain can reduce inhibitor binding by altering the conformation of the oxyanion hole [9]. Repurposing of FDA-approved drugs has identified several compounds with activity against dengue NS2B-NS3 protease, though resistance profiles remain to be fully characterized [15].

Active Site Rotamer Modifications and Structural Overlay Analysis

A critical aspect of computational resistance analysis involves examining rotameric changes in active site residues upon mutation. Rotamer libraries, derived from statistical analysis of protein crystal structures, define the preferred side chain dihedral angles for each amino acid. Resistance mutations often introduce side chains with different rotamer preferences, leading to altered steric and electrostatic environments within the binding pocket [7, 31].

For HIV-1 protease, the V82A mutation replaces a valine with a smaller alanine side chain. This change reduces hydrophobic contacts with the inhibitor's P1' moiety and creates additional void volume within the active site. Structural overlay of the wild-type and V82A mutant crystal structures reveals that the loss of the valine gamma-methyl groups results in a 1.5-2.0 Angstrom shift in the inhibitor binding pose, reducing hydrogen bond distances and increasing conformational strain [7, 18].

For SARS-CoV-2 Mpro, the E166A mutation removes the carboxylate side chain of glutamate, which normally forms a bidentate hydrogen bond with the inhibitor's P1 moiety. Rotamer analysis shows that the alanine side chain adopts a single conformation with no hydrogen bonding capacity. Structural overlay of wild-type and E166A mutant complexes demonstrates that the inhibitor repositions to compensate for the lost interaction, but the resulting pose has reduced binding affinity [5, 11].

The following table summarizes key resistance mutations and their structural consequences across viral proteases.

Computational Workflow for Resistance Analysis

The computational analysis of viral protease inhibitors and drug resistance typically follows a multi-step workflow. The following Mermaid diagram illustrates a representative pipeline.

flowchart TD
    A[Target Protease Structure], > B[Inhibitor Library Preparation]
    A, > C[Mutation Library Generation]
    B, > D[Molecular Docking]
    C, > D
    D, > E[Binding Pose Clustering]
    E, > F[Binding Free Energy Calculation]
    F, > G[MM-PBSA / MM-GBSA]
    F, > H[QM/MM Refinement]
    G, > I[Resistance Prediction]
    H, > I
    I, > J[MD Simulation Validation]
    J, > K[Rotamer Analysis]
    K, > L[Structural Overlay & Visualization]
    L, > M[Resistance Profile Summary]

The workflow begins with preparation of the target protease structure, typically obtained from X-ray crystallography or cryo-electron microscopy. Inhibitor libraries are generated through virtual screening or rational design. Mutation libraries are constructed based on known DRMs or predicted resistance hotspots [18, 19]. Molecular docking places inhibitors into the active site of both wild-type and mutant proteases. Binding poses are clustered and ranked using scoring functions [2, 8, 14].

Binding free energy calculations, using methods such as molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) or generalized Born surface area (MM-GBSA), provide quantitative estimates of affinity changes upon mutation [9, 12, 31]. QM/MM refinement captures electronic effects, such as charge redistribution and polarization, that are not accounted for by classical force fields [1, 7, 16].

MD simulations validate the stability of predicted binding modes and capture conformational dynamics. Rotamer analysis of MD trajectories identifies preferred side chain conformations and their population shifts upon mutation [31, 35]. Structural overlay of representative snapshots from wild-type and mutant simulations visualizes the spatial rearrangements underlying resistance.

Machine Learning and Data-Driven Approaches

Machine learning has become integral to resistance prediction and inhibitor design. Quantitative structure-activity relationship (QSAR) models correlate molecular descriptors with inhibitory activity. Integrated machine learning-driven QSAR combined with systems biology approaches has been applied to identify potential SARS-CoV-2 3CLpro inhibitors, incorporating both compound features and target pathway information [20].

Support vector regression with multiple kernel functions and particle swarm optimization has been used for QSAR modeling of cathepsin L inhibitors as SARS-CoV-2 therapeutics [30]. The enhanced SVR approach outperformed traditional regression methods in predicting pIC50 values, demonstrating the utility of optimized kernel functions for capturing nonlinear structure-activity relationships [30].

Deep learning methods, including convolutional neural networks and graph neural networks, have been applied to predict inhibitor binding affinity and resistance. Target-specific deep learning identified novel covalent inhibitor fragments for SARS-CoV-2 Mpro by training on a library of electrophilic compounds and their reactivity profiles [3]. The model learned to recognize favorable warhead geometries and complementary binding site features without explicit feature engineering [3].

Benchmarking studies have emphasized the importance of data set construction for machine learning-based resistance prediction. Class imbalance, sequence redundancy, and feature representation all significantly affect model performance [29]. Proper cross-validation strategies, including temporal or clade-based splitting, are necessary to avoid overestimating predictive accuracy [29].

Novel Computational Pipelines and Broad-Spectrum Inhibitor Design

Novel computational pipelines have been developed to identify target sites for broad-spectrum antiviral drugs. Sears et al. described a pipeline that integrates sequence conservation analysis, structural pocket identification, and dynamics-based druggability assessment to prioritize conserved binding sites across viral proteases [19]. This approach identified a conserved pocket in the S2 subsite of coronavirus Mpro enzymes that is amenable to broad-spectrum inhibitor design [19].

Structure-guided design of mutant peptide inhibitors has been explored for SARS-CoV-2 Mpro. Bhagat et al. designed peptide inhibitors based on the natural substrate sequence and evaluated their binding to wild-type and mutant proteases using triplicate MD simulations [31]. The designed peptides incorporated mutations at non-cleavage positions to enhance binding affinity and reduce susceptibility to resistance mutations [31].

Ligand-based virtual screening has identified potential SARS-CoV-2 Mpro inhibitors from large compound libraries. Kaur and Goyal used pharmacophore modeling and 3D-QSAR to screen millions of compounds, identifying several hits with novel scaffolds that bind to the active site with high predicted affinity [32]. These scaffolds may provide starting points for developing inhibitors with reduced resistance liability [32].

Integration with Experimental Validation

Computational predictions require experimental validation to confirm inhibitor activity and resistance profiles. High-throughput in vitro screening combined with in silico analysis has been used for Zika virus inhibitor identification, demonstrating the synergy between computational and experimental approaches [21]. Similarly, marine actinomycetes extracts have been screened for SARS-CoV-2 3CLpro inhibitory activity using integrated phylogeny-based metabolomics with functional screening and bioinformatic analysis [22].

In vitro validation of resistance mutations typically involves site-directed mutagenesis, recombinant protease expression, and enzymatic assays with inhibitor titration. The IC50 values for wild-type and mutant proteases are compared to calculate fold-change resistance [5, 11]. Computational predictions of resistance should be interpreted in the context of these experimental measurements, as discrepancies may arise from force field inaccuracies, incomplete sampling, or unaccounted allosteric effects [12, 11].

Conclusion

Computational analysis of viral protease inhibitors and drug resistance has matured into a sophisticated discipline that integrates structural biology, biophysics, and data science. Molecular docking, MD simulations, QM/MM calculations, and machine learning methods each contribute unique insights into inhibitor binding mechanisms and resistance pathways. Active site rotamer modifications, visualized through structural overlay of resistant mutants, provide a mechanistic understanding of how single amino acid substitutions alter inhibitor affinity. The continued development of computational methods, combined with experimental validation, will be essential for designing next-generation inhibitors with improved resistance profiles and for predicting the evolutionary trajectories of drug-resistant viruses.

References

[1] Świderek K, Moliner V. The pivotal role of β-lactone stereochemistry in the development of SARS-CoV-2 M(pro) inhibitors. Protein Sci. 2026. https://pubmed.ncbi.nlm.nih.gov/42316787/

[2] Dansana J, Acharya SS, Meher BR. Elucidating structure-activity relationships of benzopyran derivatives as inhibitors of Chikungunya virus nsP2 and nsP4 via integrated docking, molecular dynamics, and DFT studies. J Mol Graph Model. 2026. https://pubmed.ncbi.nlm.nih.gov/42275761/

[3] Zhou W, D Oliviera A, Dai X, et al. A Novel Covalent Inhibitor Fragment for the SARS-CoV-2 Main Protease Identified by Target-Specific Deep Learning. ACS Chem Biol. 2026. https://pubmed.ncbi.nlm.nih.gov/42066065/

[4] Khan DA, Sheikh MA, Imon RR, et al. Boesenbergia rotunda-Derived Phytochemicals as Potent Inhibitors of SARS-CoV-2 Papain-Like Protease (PLpro): Insights From Molecular Docking and Dynamic Simulation. Scientifica (Cairo). 2026. https://pubmed.ncbi.nlm.nih.gov/42064585/

[5] Chen L, Su H, Shang W, et al. SARS-CoV-2 3CLpro mutations T21I and E166A confer differential resistance to simnotrelvir, bofutrelvir, and ensitrelvir. J Virol. 2026. https://pubmed.ncbi.nlm.nih.gov/42037412/

[6] Anderson KJ, Lee MS, Sevillano N, et al. Structural Basis of Serine Protease Inhibition by Antibodies from Biased Fab Phage-Display Libraries. bioRxiv. 2026. https://pubmed.ncbi.nlm.nih.gov/41959119/

[7] Gutiérrez-Flores J, Padilla-Bernal G, Sánchez-Juárez C, et al. Residue-Level Affinity Decomposition via Quantum Electron Density: A Multivariable Framework Applied to HIV-1 Protease Inhibitors. J Phys Chem B. 2026. https://pubmed.ncbi.nlm.nih.gov/41923454/

[8] Malik A, Noreen S, Ijaz B. Inhibitory potential of Morus alba leaf extract and its phytoconstituent against SARS-CoV-2 main protease: An integrative in silico and in vitro analysis. J Mol Graph Model. 2026. https://pubmed.ncbi.nlm.nih.gov/41875707/

[9] Bantun F, Chaturvedi N, Jalal NA, et al. Molecular dynamics-driven optimization of triterpenoid, amidinium, and flavonoid inhibitors targeting dengue NS2B-NS3 protease. Comput Biol Chem. 2026. https://pubmed.ncbi.nlm.nih.gov/41855844/

[10] Budiman A, Paujiah ASU, Fakih TM, et al. Evaluation of co-amorphous systems containing ritonavir and coformers: An experimental and molecular simulation study. J Pharm Sci. 2026. https://pubmed.ncbi.nlm.nih.gov/41771489/

[11] Zhou Y, Gammeltoft KA, Tjørnelund-Sjursen HD, et al. SARS-CoV-2 Mpro inhibitor ensitrelvir: asymmetrical cross-resistance

[12] Liu M, Zhao D, Duan H, et al. Mutation-Induced Resistance of SARS-CoV-2 M(pro) to WU-04 Revealed by Multi-Scale Modeling. Int J Mol Sci. 2026. https://pubmed.ncbi.nlm.nih.gov/41596645/

[13] Mahmoud E, Rehan M. Computational Exploration of Flavonoids as HCV NS3/4A Protease Inhibitors: Advancing Antiviral Therapies to Mitigate Liver Cancer Risk. Curr Pharm Des. 2025. https://pubmed.ncbi.nlm.nih.gov/41088918/

[14] Uddin MAR, Paul AC, Islam MS, et al. Efficacy of phytochemicals derived from Artocarpus heterophyllus (Jackfruit) as inhibitors against NS2B/NS3 protease of dengue virus: an in-silico investigation. Sci Rep. 2026. https://pubmed.ncbi.nlm.nih.gov/41652071/

[15] Adiba, Uzma A, Babar MM, et al. Repurposing of FDA-approved drugs to identify potential ns2b/ns3 protease inhibitors against dengue virus: In silico and in vitro evaluation. Acta Trop. 2026. https://pubmed.ncbi.nlm.nih.gov/41571221/

[16] da Silva Ribeiro FV, da Penha Valente RP, Kruger HG, et al. Mechanistic insights into the noncovalent inhibition of SARS-CoV-2 PLpro: a multiscale computational study. J Comput Aided Mol Des. 2026. https://pubmed.ncbi.nlm.nih.gov/41639476/

[17] Sá Magalhães Serafim M, Kronenberger T, Francisco KR, et al. Discovery of benzyl carbamate inhibitors of coronavirus M(pro) enzymes from a legacy collection of cysteine protease inhibitors. J Enzyme Inhib Med Chem. 2025. https://pubmed.ncbi.nlm.nih.gov/41246830/

[18] Aggarwal M, Periwal V. Forecasting drug resistant HIV protease evolution. PLoS Comput Biol. 2026. https://pubmed.ncbi.nlm.nih.gov/41592116/

[19] Sears JD, Popov KI, Sylvester PA, et al. Novel computational pipeline to identify target sites for broad spectrum antiviral drugs. Antiviral Res. 2026. https://pubmed.ncbi.nlm.nih.gov/41338303/

[20] Manaithiya A, Bhowmik R, Ray R, et al. Integrated machine learning-driven QSAR modelling and systems biology approach for the identification of potential SARS-CoV-2 3CLpro inhibitors. SAR QSAR Environ Res. 2025. https://pubmed.ncbi.nlm.nih.gov/41255262/

[21] Batool S, Chokkakula S, Kim BK, et al. High-throughput in vitro screening and in silico analysis for Zika virus inhibitor identification. Sci Rep. 2025. https://pubmed.ncbi.nlm.nih.gov/41310171/

[22] Dalisay DS, Mateo JC, Teodosio JJR, et al. Identification of SARS-CoV-2 3CLpro inhibitors from marine actinomycetes through integrated phylogeny-based metabolomics with functional screening and bioinformatic analysis. Sci Rep. 2025. https://pubmed.ncbi.nlm.nih.gov/41285897/

[23] Sardar P, Ranjit B, Biswas S, et al. Modeling NRTIs and PIs class drug therapy on the dynamics of HIV infection with real patient data analysis and optimized control strategy. Sci Rep. 2026. https://pubmed.ncbi.nlm.nih.gov/41501180/

[24] Mishra S, Kumar SB, Kumari A, et al. Computational insights into a protease inhibitor from Streptomyces globosus VITSMAB-2 molecular docking and dynamics simulations against SARS-CoV-2 main protease. Sci Rep. 2025. https://pubmed.ncbi.nlm.nih.gov/41390850/

[25] Al-Thiabat MG, Alhawarri MB. Computational Insights Into New Portulaca oleracea Alkaloids as Potential SARS-CoV-2 PL(pro) Inhibitors. Chem Biodivers. 2026. https://pubmed.ncbi.nlm.nih.gov/41388231/

[26] Behrens GMN, Assoumou L, Liegeon G, et al. Integrase versus protease inhibitor therapy in advanced HIV disease (LAPTOP): a multicountry, randomised, open-label, non-inferiority trial. Lancet Infect Dis. 2026. https://pubmed.ncbi.nlm.nih.gov/41344354/

[27] Oladipo SD, Olalekan SO, Obakachi VA, et al. Thiuram disulfides as prospective inhibitors of 3-chymotrypsin-like cysteine protease: an in-silico approach targeting SARS-CoV-2. In Silico Pharmacol. 2025. https://pubmed.ncbi.nlm.nih.gov/41322711/