Development and test of highly accurate endpoint free energy methods. 2: Prediction of logarithm of n-octanol-water partition coefficient (logP) for druglike molecules using MM-PBSA method.

The logarithm of n-octanol-water partition coefficient (logP) is frequently used as an indicator of lipophilicity in drug discovery, which has substantial impacts on the absorption, distribution, metabolism, excretion, and toxicity of a drug candidate. Considering that the experimental measurement of the property is costly and time-consuming, it is of great importance to develop reliable prediction models for logP. In this study, we developed a transfer free energy-based logP prediction model-FElogP. FElogP is based on the simple principle that logP is determined by the free energy change of transferring a molecule from water to n-octanol. The underlying physical method to calculate transfer free energy is the molecular mechanics-Poisson Boltzmann surface area (MM-PBSA), thus this method is named as free energy-based logP (FElogP). The superiority of FElogP model was validated by a large set of 707 structurally diverse molecules in the ZINC database for which the measurement was of high quality. Encouragingly, FElogP outperformed several commonly-used QSPR or machine learning-based logP models, as well as some continuum solvation model-based methods. The root-mean-square error (RMSE) and Pearson correlation coefficient (R) between the predicted and measured values are 0.91 log units and 0.71, respectively, while the runner-up, the logP model implemented in OpenBabel had an RMSE of 1.13 log units and R of 0.67. Given the fact that FElogP was not parameterized against experimental logP directly, its excellent performance is likely to be expanded to arbitrary organic molecules covered by the general AMBER force fields.

Machine learning on ligand-residue interaction profiles to significantly improve binding affinity prediction.

Structure-based virtual screenings (SBVSs) play an important role in drug discovery projects. However, it is still a challenge to accurately predict the binding affinity of an arbitrary molecule binds to a drug target and prioritize top ligands from an SBVS. In this study, we developed a novel method, using ligand-residue interaction profiles (IPs) to construct machine learning (ML)-based prediction models, to significantly improve the screening performance in SBVSs. Such a kind of the prediction model is called an IP scoring function (IP-SF). We systematically investigated how to improve the performance of IP-SFs from many perspectives, including the sampling methods before interaction energy calculation and different ML algorithms. Using six drug targets with each having hundreds of known ligands, we conducted a critical evaluation on the developed IP-SFs. The IP-SFs employing a gradient boosting decision tree (GBDT) algorithm in conjunction with the MIN + GB simulation protocol achieved the best overall performance. Its scoring power, ranking power and screening power significantly outperformed the Glide SF. First, compared with Glide, the average values of mean absolute error and root mean square error of GBDT/MIN + GB decreased about 38 and 36%, respectively. Second, the mean values of squared correlation coefficient and predictive index increased about 225 and 73%, respectively. Third, more encouragingly, the average value of the areas under the curve of receiver operating characteristic for six targets by GBDT, 0.87, is significantly better than that by Glide, which is only 0.71. Thus, we expected IP-SFs to have broad and promising applications in SBVSs.

In silico binding profile characterization of SARS-CoV-2 spike protein and its mutants bound to human ACE2 receptor.

Severe acute respiratory syndrome coronavirus (SARS-CoV-2), a novel coronavirus, has brought an unprecedented pandemic to the world and affected over 64 million people. The virus infects human using its spike glycoprotein mediated by a crucial area, receptor-binding domain (RBD), to bind to the human ACE2 (hACE2) receptor. Mutations on RBD have been observed in different countries and classified into nine types: A435S, D364Y, G476S, N354D/D364Y, R408I, V341I, V367F, V483A and W436R. Employing molecular dynamics (MD) simulation, we investigated dynamics and structures of the complexes of the prototype and mutant types of SARS-CoV-2 spike RBDs and hACE2. We then probed binding free energies of the prototype and mutant types of RBD with hACE2 protein by using an end-point molecular mechanics Poisson Boltzmann surface area (MM-PBSA) method. According to the result of MM-PBSA binding free energy calculations, we found that V367F and N354D/D364Y mutant types showed enhanced binding affinities with hACE2 compared to the prototype. Our computational protocols were validated by the successful prediction of relative binding free energies between prototype and three mutants: N354D/D364Y, V367F and W436R. Thus, this study provides a reliable computational protocol to fast assess the existing and emerging RBD mutations. More importantly, the binding hotspots identified by using the molecular mechanics generalized Born surface area (MM-GBSA) free energy decomposition approach can guide the rational design of small molecule drugs or vaccines free of drug resistance, to interfere with or eradicate spike-hACE2 binding.

Phosphorylation at Ser289 Enhances the Oligomerization of Tau Repeat R2.

In tauopathies such as Alzheimer's disease (AD), aberrant phosphorylation causes the dissociation of tau proteins from microtubules. The dissociated tau then aggregates into sequent forms from soluble oligomers to paired helical filaments and insoluble neurofibrillary tangles (NFTs). NFTs is a hallmark of AD, while oligomers are found to be the most toxic form of the tau aggregates. Therefore, understanding tau oligomerization with regard to abnormal phosphorylation is important for the therapeutic development of AD. In this study, we investigated the impact of phosphorylated Ser289, one of the 40 aberrant phosphorylation sites of full-length tau proteins, on monomeric and dimeric structures of tau repeat R2 peptides. We carried out intensive replica exchange molecular dynamics simulation with a total simulation time of up to 0.1 ms. Our result showed that the phosphorylation significantly affected the structures of both the monomer and the dimer. For the monomer, the phosphorylation enhanced ordered-disordered structural transition and intramolecular interaction, leading to more compactness of the phosphorylated R2 compared to the wild-type one. As to the dimer, the phosphorylation increased intermolecular interaction and ß-sheet formation, which can accelerate the oligomerization of R2 peptides. This result suggests that the phosphorylation at Ser289 is likely to promote tau aggregation. We also observed a phosphorylated Ser289-Na+-phosphorylated Ser289 bridge in the phosphorylated R2 dimer, suggesting an important role of cation ions in tau aggregation. Our findings suggest that phosphorylation at Ser289 should be taken into account in the inhibitor screening of tau oligomerization.

Stable Cavitation Interferes with Aß16-22 Oligomerization.

Ultrasound and microbubbles are used for many medical applications nowadays. Scanning ultrasound can remove amyloid-ß (Aß) aggregates in the mouse brain and restores memory in an Alzheimer's disease mouse model. In vitro studies showed that amyloid fibrils are fragmented due to the ultrasound-induced bubble inertial cavitation, and ultrasonic pulses accelerate the depolymerization of Aß fibrils into monomers at 1 µM of concentration. Under applied ultrasound, microbubbles can be in a stable oscillating state or unstable inertial cavitation state. The latter occurs when ultrasound causes a dramatic change of bubble sizes above a certain acoustic pressure. We have developed and implemented a nonequilibrium molecular dynamics simulation algorithm to the AMBER package, to facilitate the investigation of the molecular mechanism of Aß oligomerization under stable cavitation. Our results indicated that stable cavitation not only inhibited oligomeric formation, but also prevented the formation of ß-rich oligomers. The network analysis of state transitions revealed that stable cavitation altered the oligomerization pathways of Aß16-22 peptides. Our simulation tool may be applied to optimize the experimental conditions to achieve the best therapeutical effect.

Elastic moduli of normal and cancer cell membranes revealed by molecular dynamics simulations.

Recent studies indicate that there are mechanical differences between normal cells and cancer cells. Because the cell membrane takes part in a variety of vital processes, we test the hypothesis of whether or not two fundamental alterations in the cell membrane, i.e., the overexpression of phosphatidylserine lipids in the outer leaflet and a reduction in cholesterol concentration, could cause the softening in cancer cells. Adopting ten models of normal and cancer cell membranes, we carry out 1 µs all-atom molecular dynamics simulations to compare the structural properties and elasticity properties of two membrane types. We find that the overexpression of the phosphatidylserine lipids in the outer leaflet does not significantly alter the area per lipid, the membrane thickness, the lipid order parameters and the elasticity moduli of the cancer membranes. However, a reduction in the cholesterol concentration leads to clear changes in those quantities, especially decreases in the bending, tilt and twist moduli. This implies that the reduction of cholesterol concentration in the cancer membranes could contribute to the softening of cancer cells.

A multiple-step in silico screening protocol to identify allosteric inhibitors of Spike-hACE2 binding.

While the COVID-19 pandemic continues to worsen, effective medicines that target the life cycle of SARS-CoV-2 are still under development. As more highly infective and dangerous variants of the coronavirus emerge, the protective power of vaccines will decrease or vanish. Thus, the development of drugs, which are free of drug resistance is direly needed. The aim of this study is to identify allosteric binding modulators from a large compound library to inhibit the binding between the Spike protein of the SARS-CoV-2 virus and human angiotensin-converting enzyme 2 (hACE2). The binding of the Spike protein to hACE2 is the first step of the infection of host cells by the coronavirus. We first built a compound library containing 77 448 antiviral compounds. Molecular docking was then conducted to preliminarily screen compounds which can potently bind to the Spike protein at two allosteric binding sites. Next, molecular dynamics simulations were performed to accurately calculate the binding affinity between the spike protein and an identified compound from docking screening and to investigate whether the compound can interfere with the binding between the Spike protein and hACE2. We successfully identified two possible drug binding sites on the Spike protein and discovered a series of antiviral compounds which can weaken the interaction between the Spike protein and hACE2 receptor through conformational changes of the key Spike residues at the Spike-hACE2 binding interface induced by the binding of the ligand at the allosteric binding site. We also applied our screening protocol to another compound library which consists of 3407 compounds for which the inhibitory activities of Spike/hACE2 binding were measured. Encouragingly, in vitro data supports that the identified compounds can inhibit the Spike-ACE2 binding. Thus, we developed a promising computational protocol to discover allosteric inhibitors of the binding of the Spike protein of SARS-CoV-2 to the hACE2 receptor, and several promising allosteric modulators were discovered.

Molecular dynamics simulation of cancer cell membrane perforated by shockwave induced bubble collapse.

It has been widely accepted that cancer cells are softer than their normal counterparts. This motivates us to propose, as a proof-of-concept, a method for the efficient delivery of therapeutic agents into cancer cells, while normal cells are less affected. The basic idea of this method is to use a water jet generated by the collapse of the bubble under shockwaves to perforate pores in the cell membrane. Given a combination of shockwave and bubble parameters, the cancer membrane is more susceptible to bending, stretching, and perforating than the normal membrane because the bending modulus of the cancer cell membrane is smaller than that of the normal cell membrane. Therefore, the therapeutic agent delivery into cancer cells is easier than in normal cells. Adopting two well-studied models of the normal and cancer membranes, we perform shockwave induced bubble collapse molecular dynamics simulations to investigate the difference in the response of two membranes over a range of shockwave impulse 15-30 mPa s and bubble diameter 4-10 nm. The simulation shows that the presence of bubbles is essential for generating a water jet, which is required for perforation; otherwise, pores are not formed. Given a set of shockwave impulse and bubble parameters, the pore area in the cancer membrane is always larger than that in the normal membrane. However, a too strong shockwave and/or too large bubble results in too fast disruption of membranes, and pore areas are similar between two membrane types. The pore closure time in the cancer membrane is slower than that in the normal membrane. The implications of our results for applications in real cells are discussed in some details. Our simulation may be useful for encouraging future experimental work on novel approaches for cancer treatment.

Molecular Mechanism of Ultrasound-Induced Structural Defects in Liposomes: A Nonequilibrium Molecular Dynamics Simulation Study.

The use of ultrasound in combination with liposomes is a promising approach to improve drug delivery. To achieve an optimal drug release rate, it is important to understand how ultrasound induces pathways on the liposome surface where drugs can be released from the liposome. To this end, we carry out large-scale ultrasound-induced molecular dynamics simulations for three single lipid component liposomes formed from the commonly used phospholipids: 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC), 1,2-dipalmitoylphosphatidylcholine (DPPC), or phosphatidylcholine (POPC). The results show that ultrasound induces the detachment of two leaflets of the DOPC surface, suggesting that the drug release pathway may be through the low lipid packing areas on the stretched surface. In contrast, ultrasound induces pore formation on the surface of DPPC and DOPC, where drugs could escape from the liposomes. While the leaflet detachment and transient pore formation are the mechanisms of DOPC and DPPC, respectively, in both liquid-ordered and liquid-disordered phases, the leaflet detachment mechanism is switched to the transient pore formation mechanism on going from the liquid-ordered phase to the liquid-disordered phase in the POPC liposome. By adding 30% mol cholesterol, the leaflet detachment mechanism is observed in all liposomes. We found that the molecular origin that determines a mechanism is the competition between the intraleaflet and interleaflet interacting energy of lipids. The connection to experimental and theoretical modeling is discussed in some detail.

Molecular mechanism of ultrasound interaction with a blood brain barrier model.

The brain is strictly protected by the blood brain barrier preventing the crossing of therapeutics to treat brain diseases. The high and low intensity focused ultrasound methods have been used to temporarily open the blood brain barrier, facilitating the transport of drugs. The methods are very promising because the opening is transient, localized, and noninvasive. However, the molecular mechanism of the opening is unknown, and this limits the development and application of these methods. With this in mind, we carry out a molecular dynamics simulation study to understand the interaction of ultrasound with the cell membrane and the tight junction. Our minimal blood brain barrier model is composed of two lipid bilayers, mimicking two portions of neighboring cells, connected together by a tight junction formed by a pair of two cis-dimers of the claudin-5 protein. Using an experimental ultrasound frequency of 50 MHz, simulations show that at low intensities, ultrasound does not impact the structure of the cell membranes and tight junction, implying that the direct interaction of ultrasound with the blood brain barrier is not responsible for the experimentally observed opening. At high intensities, the ultrasound pulls the monolayers of individual cell membrane lipid bilayers apart, creating air compartments inside the bilayers. This reduces the free energy barrier for the translocation of drugs across the lipid bilayer and enhances drug permeability. At very high intensities, the two monolayers are largely separated, resulting in cell damage and implying that the blood brain barrier is primarily opened at the experimentally observed damaged areas.

E-motif formed by extrahelical cytosine bases in DNA homoduplexes of trinucleotide and hexanucleotide repeats.

Atypical DNA secondary structures play an important role in expandable trinucleotide repeat (TR) and hexanucleotide repeat (HR) diseases. The cytosine mismatches in C-rich homoduplexes and hairpin stems are weakly bonded; experiments show that for certain sequences these may flip out of the helix core, forming an unusual structure termed an 'e-motif'. We have performed molecular dynamics simulations of C-rich TR and HR DNA homoduplexes in order to characterize the conformations, stability and dynamics of formation of the e-motif, where the mismatched cytosines symmetrically flip out in the minor groove, pointing their base moieties towards the 5'-direction in each strand. TRs have two non-equivalent reading frames, (GCC)n and (CCG)n; while HRs have three: (CCCGGC)n, (CGGCCC)n, (CCCCGG)n. We define three types of pseudo basepair steps related to the mismatches and show that the e-motif is only stable in (GCC)n and (CCCGGC)n homoduplexes due to the favorable stacking of pseudo GpC steps (whose nature depends on whether TRs or HRs are involved) and the formation of hydrogen bonds between the mismatched cytosine at position i and the cytosine (TRs) or guanine (HRs) at position i - 2 along the same strand. We also characterize the extended e-motif, where all mismatched cytosines are extruded, their extra-helical stacking additionally stabilizing the homoduplexes.

New application of in silico methods in identifying mechanisms of action and key components of anti-cancer herbal formulation YIV-906 (PHY906).

YIV-906 (formally PHY906, KD018) is a four-herb formulation that is currently being developed to improve the therapeutic index and ameliorate the side effects of many chemotherapeutic drugs including sorafenib, irinotecan, and capecitabine. However, as a promising anti-cancer adjuvant, the molecular mechanism of action of YIV-906 remains unrevealed due to its multi-component and multi-target features. Since YIV-906 has been shown to induce apoptosis and autophagy in cancer cells through modulating the negative regulators of ERK1/2, namely DUSPs, it is of great interest to elucidate the key components that cause the therapeutic effect of YIV-906. In this work, we investigated the mechanism of YIV-906 inhibiting DUSPs, using a broad spectrum of molecular modelling techniques, including molecular docking, molecular dynamics (MD) simulations, and binding free energy calculations. In total, MD simulations and binding free energy calculations were performed for 99 DUSP-ligand complexes. We found that some herbal components or their metabolites could inhibit DUSPs. Based on the docking scores and binding free energies, the sulfation and glucuronidation metabolites of the S ingredient in YIV-906 play a leading role in inhibiting DUSPs, although several original herbal chemicals with carboxyl groups from the P and Z ingredients also make contributions to this inhibitory effect. It is not a surprise that the electrostatic interaction plays the dominant role in the ligand binding process, given the fact that several charged residues reside in the binding pockets of DUSPs. Our MD simulation results demonstrate that the sulfate moieties and carboxyl moieties of the advantageous ligands from YIV-906 can occupy the enzymes' catalytic sites, mimicking the endogenous phosphate substrates of DUSPs. As such, the ligand binding can inhibit the association of DUSPs and ERK1/2, which in turn reduces the dephosphorylation of ERK1/2 and causes cell cycle arrest in the tumor. Our modelling study provides useful insights into the rational design of highly potent anti-cancer drugs targeting DUSPs. Finally, we have demonstrated that multi-scale molecular modelling techniques are able to elucidate molecular mechanisms involving complex molecular systems.

Nonequilibrium atomistic molecular dynamics simulation of tubular nanomotor propelled by bubble propulsion.

We develop a molecular nanoscaled model for tubular motors propelled by bubble propulsion. The motor is modeled by a carbon nanotube, and the bubble is represented by a particle interacting with water by a time-dependent potential. Effects of liquid viscosity, fuel concentration, geometry, and size of the tube on the performance of the motor are effectively encoded into two parameters: time scales of the bubble expansion and bubble formation. Our results are qualitatively consistent with experimental data of much larger motors. Simulations suggest that (i) the displacement of the tube is optimized if two time scales are as short as possible, (ii) the compromise between the performance and fuel consumption is achieved if the bubble formation time is shorter than the velocity correlation time of the tube, (iii) the motor efficiency is higher with slow expansion, short formation of the bubble than fast growth but long formation time, and (iv) the tube is propelled by strong forces on the order of mN, reaching high speeds up to ∼60 m/s. Our simulation may be useful for refining and encouraging future experimental work on nanomotors having the size of a few nanometers. The tiny size and high speed motors could have great potential applications in real life.

Interaction mechanism between the focused ultrasound and lipid membrane at the molecular level.

Focused ultrasound (FUS) has a wide range of medical applications. Nowadays, the diagnostic and therapeutic ultrasound procedures are routinely used; effects of ultrasound on biological systems at the molecular level are, however, not fully understood. Experimental results on the interaction of the cell membrane, a simplest but important system component, with ultrasound are controversial. Molecular dynamics (MD) simulations could provide valuable insights, but there is no single study on the mechanism of the FUS induced structural changes in cell membranes. With this in mind, we develop a simple method to include FUS into a standard MD simulation. Adopting the 1,2-dioleoyl-sn-glycero-3-phosphocholine lipid membrane as a representative model described by the MARTINI coarse-grained force field, and using experimental values of the ultrasound frequency and intensity, we show that the heat and bubble cavitation are not the primary direct mechanisms that cause structural changes in the membrane. The spatial pressure gradients between the focused and free regions and between the parallel and perpendicular directions to the membrane are the origin of the mechanism. These gradients force lipids to move out of the focused region, forming a lipid flow along the membrane diagonal. Lipids in the free region move in the opposite direction due to the conservation of the total momentum. These opposite motions create wrinkles along the membrane diagonal at low FUS intensities and tear up the membrane at high FUS intensities. Once the membrane is torn up, it is not easy to reform. The implication of our findings in the FUS-induced drug delivery is discussed in some detail.

Structure and Dynamics of DNA and RNA Double Helices of CAG and GAC Trinucleotide Repeats.

CAG trinucleotide repeats are known to cause 10 late-onset progressive neurodegenerative disorders as the repeats expand beyond a threshold, whereas GAC repeats are associated with skeletal dysplasias and expand from the normal five to a maximum of seven repeats. The TR secondary structure is believed to play a role in CAG expansions. We have carried out free energy and molecular dynamics studies to determine the preferred conformations of the A-A noncanonical pairs in (CAG)n and (GAC)n trinucleotide repeats (n = 1, 4) and the consequent changes in the overall structure of the RNA and DNA duplexes. We find that the global free energy minimum corresponds to A-A pairs stacked inside the core of the helix with anti-anti conformations in RNA and (high-anti)-(high-anti) conformations in DNA. The next minimum corresponds to anti-syn conformations, whereas syn-syn conformations are higher in energy. Transition rates of the A-A conformations are higher for RNA than DNA. Mechanisms for these various transitions are identified. Additional structural and dynamical aspects of the helical conformations are explored, with a focus on contrasting CAG and GAC duplexes. The neutralizing ion distribution around the noncanonical pairs is described.

Comparative melting and healing of B-DNA and Z-DNA by an infrared laser pulse.

We explore the use of a fast laser melting simulation approach combined with atomistic molecular dynamics simulations in order to determine the melting and healing responses of B-DNA and Z-DNA dodecamers with the same d(5'-CGCGCGCGCGCG-3')2 sequence. The frequency of the laser pulse is specifically tuned to disrupt Watson-Crick hydrogen bonds, thus inducing melting of the DNA duplexes. Subsequently, the structures relax and partially refold, depending on the field strength. In addition to the inherent interest of the nonequilibrium melting process, we propose that fast melting by an infrared laser pulse could be used as a technique for a fast comparison of relative stabilities of same-sequence oligonucleotides with different secondary structures with full atomistic detail of the structures and solvent. This could be particularly useful for nonstandard secondary structures involving non-canonical base pairs, mismatches, etc.

Investigation of the Structure of Full-Length Tau Proteins with Coarse-Grained and All-Atom Molecular Dynamics Simulations.

Tau proteins not only have many important biological functions but also are associated with several neurodegenerative diseases, such as Parkinson's disease and Alzheimer's disease (AD). However, it is still a challenge to identify the atomic structure of full-length tau proteins due to their lengthy and disordered characteristics and the factor that there are no crystal structures of full-length tau proteins available. We performed multi- and large-scale molecular dynamics simulations of the full-length tau monomer (the 2N4R isoform and 441 residues) in aqueous solution under biological conditions with coarse-grained and all-atom force fields. The obtained atomic structures produced radii of gyration and chemical shifts that are in excellent agreement with those of experiment. The generated monomer structure ensemble would be very useful for further studying the oligomerization mechanism and discovering tau oligomerization inhibitors, which are important events in AD drug development.

Phosphorylation of Tau R2 Repeat Destabilizes Its Binding to Microtubules: A Molecular Dynamics Simulation Study.

Phosphorylation, the most popular post-translational modification of tau protein, plays an important role in regulating tau physiological functions. However, aberrant phosphorylation attenuates the binding affinity of tau to a microtubule (MT), resulting in MT destabilization followed by accumulation of neurofibrillary tangles in the brain. There are in total 85 potential phosphorylation sites in a full-length tau protein, and about half of them are abnormal as they occur in tau of Alzheimer's disease (AD) brain only. In this work, we investigated the impact of abnormal Ser289, Ser293, and Ser289/Ser293 phosphorylation on tau R2-MT binding and the conformation of tau R2 using molecular dynamics simulation. We found that the phosphorylation significantly affected R2-MT interaction and reduced the binding affinity of tau R2 peptides to MTs. Free energy decomposition analysis suggested that the post-translational modified residues themselves made a significant contribution to destabilize tau repeat R2-MT binding. Therefore, the phosphorylation may attenuate the binding affinity of tau to MTs. Additionally, the phosphorylation also enhanced helix-coil transition of monomeric R2 peptides, which may result in the acceleration of tau aggregation. Since these phosphorylated sites have not been examined in previous experimental studies, our finding through all-atom molecular dynamics simulations and free energy analysis can inspire experimental scientists to investigate the impact of the phosphorylation on MT binding and aggregation of full-length tau and the pathological roles of the phosphorylation at those sites in AD development through in vitro/in vivo assays.

Comparison and summary of in silico prediction tools for CYP450-mediated drug metabolism.

The cytochrome P450 (CYP450) enzyme system is responsible for the metabolism of more than two-thirds of xenobiotics. This review summarizes reports of a series of in silico tools for CYP450 enzyme-drug interaction predictions, including the prediction of sites of metabolism (SOM) of a drug and the identification of inhibitor/substrates for CYP subtypes. We also evaluated four prediction tools to identify CYP inhibitors utilizing 52 of the most frequently prescribed drugs. ADMET Predictor and CYPlebrity demonstrated the best performance. We hope that this review provides guidance for choosing appropriate enzyme prediction tools from a variety of in silico platforms to meet individual needs. Such predictions are useful for medicinal chemists to prioritize their designed compounds for further drug discovery.

Unpolarized laser method for infrared spectrum calculation of amide I CO bonds in proteins using molecular dynamics simulation.

The investigation of the strong infrared (IR)-active amide I modes of peptides and proteins has received considerable attention because a wealth of detailed information on hydrogen bonding, dipole-dipole interactions, and the conformations of the peptide backbone can be derived from the amide I bands. The interpretation of experimental spectra typically requires substantial theoretical support, such as direct ab-initio molecular dynamics simulation or mixed quantum-classical description. However, considering the difficulties associated with these theoretical methods and their applications are limited in small peptides, it is highly desirable to develop a simple yet efficient approach for simulating the amide I modes of any large proteins in solution. In this work, we proposed a comprehensive computational method that extends the well-established molecular dynamics (MD) simulation method to include an unpolarized IR laser for exciting the CO bonds of proteins. We showed the amide I frequency corresponding to the frequency of the laser pulse which resonated with the CO bond vibration. At this frequency, the protein energy and the CO bond length fluctuation were maximized. Overall, the amide I bands of various single proteins and amyloids agreed well with experimental data. The method has been implemented into the AMBER simulation package, making it widely available to the scientific community. Additionally, the application of the method to simulate the transient amide I bands of amyloid fibrils during the IR laser-induced disassembly process was discussed in details.