Approaches of Computational Biophysics and Chemistry in Molecular Biology

The following sections are included:

• Preface
• About the Editor
• Contents

Bonded (or valence) interactions, which directly determine the local structures of the molecules, are fundamental parts of molecular mechanics force fields (FFs). Most popular classical FFs adopt the simple harmonic models for bond stretching and angle bending and ignore cross-coupling effects among the valence terms. This may lead to less accurate vibrational properties and configurations in molecular dynamics (MD) simulations. AMOEBA models utilize an MM3(MM4)-style bonded interaction model, in which the vibrational anharmonicity, the coupling effects among different energy terms, and the out-of-plane bending for sp2-hybridized atoms are considered. In this work, we report the development of bonded interaction parameters for a wide range of chemistry based on quantum mechanics (QM). About 270 atomic types defined by SMARTS strings were used to model the valence interactions. Our results indicate that the resulting valence parameters produce accurate vibrational frequencies (RMSD from QM is less than ∼ 36:6 cm⁻¹) over a large set of molecules with diverse functional groups (445 molecules). By contrast, the harmonic models usually give an RMS error greater than 60 cm⁻¹. Meanwhile, this model accurately reflects the potential energy surface of the out-of-plane bending. Our model can generally be applied to the AMOEBA family and any MM3 (MM4)-based molecular mechanics FFs.

The correct reproduction of conformational substates of amino acids was tested for the CHARMM Drude polarizable force field. This was achieved by evaluating the reorganization energies for all low lying energy minima occurring in all 15 neutral blocked amino acids on a quantum-mechanical (QM) energy surface at the MP2/cc-pVDZ level. The results indicate that the bonded parameters of the N-acetyl (ACE) and N-Methylamide (CT3) blocking groups lead to significant discrepancies. A reparametrization of five bond angles significantly improved the agreement with the QM energy surface. The corrected Drude force field exhibits almost the same average reorganization energies relative to the MP2 energy surface as the AM1 and PM3 semi-empirical methods.

Molecular dynamics (MD) simulations play a crucial role in modeling biomolecular systems in which the electrostatic interactions are critical in dictating the structural and dynamical properties. Thus, the treatment of the electrostatic interactions defined in the underlying force field (FF) strongly affects the simulation accuracy. Most FFs use fixed partial atomic charges to include electrostatic interactions, and therefore lack the electronic polarization response, representing an intrinsic limitation. To address this limitation, polarizable FFs have been developed that treat atomic polarizabilities explicitly. Here we present the application of the allatom polarizable (Drude) and non-polarizable (CHARMM) nucleic acid FFs in RNA hairpin systems to investigate the impact of polarization on structural properties, dipole moment distributions, and cation interactions. Results show that the presence of polarizability in the FF significantly improves the stabilization of RNA hairpin structure. As expected, the distributions of dipole moments show more fluctuations when simulated using the polarizable FF, with the variation in dipoles contributing to the stabilization of the structures of the loop regions of the RNAs. Contact map analyses between the bases and cations show that the variation of the ion distribution around the entire hairpin is larger for the polarizable FF and the cations occupy the outer hydration shell to a greater extent. The presented results indicate the importance of the explicit treatment of electronic polarizability in molecular simulations of RNA, including in non-canonical regions.

Receptor–ligand interactions are involved in various biological processes, therefore understanding the binding mechanism and ability to predict the binding mode are essential for many biological investigations. While many computational methods exist to predict the 3D structure of the corresponding complex provided the knowledge of the monomers, here we use the newly developed DelPhiForce steered Molecular Dynamics (DFMD) approach to model the binding of barstar to barnase to demonstrate that first-principles methods are also capable of modeling the binding. Essential component of DFMD approach is enhancing the role of long-range electrostatic interactions to provide guiding force of the monomers toward their correct binding orientation and position. Thus, it is demonstrated that the DFMD can successfully dock barstar to barnase even if the initial positions and orientations of both are completely different from the correct ones. Thus, the electrostatics provides orientational guidance along with pulling force to deliver the ligand in close proximity to the receptor.

Intrinsically disordered proteins (IDPs) are a class of proteins without stable three-dimensional structures under physiological conditions. IDPs exhibit high dynamic nature and could be described by structural ensembles. As one of the most widely used tools, molecular dynamics (MD) simulation could provide the atomic descriptions of the structural ensemble of IDPs. However, the accuracy of the MD simulation largely depends on the accuracy of the force field. In this paper, we compared the structural ensembles of the activation domain 1 (AD1) in p53 tumor suppressor obtained from the widely used force fields, AMBER99SB-ILDN, CHARMM27, CHARMM36m with different water models. The results show that CHARMM36m generates more extended conformations than other force fields, while CHARMM27 prefers to sample the α-helical structure. Moreover, the chemical shifts obtained by CHARMM36m are the closest to the experimental measurements. These results indicate that the CHARMM36m force field performs best in characterizing the structure properties of p53 AD1. Water models are also critical to describe the structural ensemble of IDPs.

TIP4P water model can obtain more extended conformations and produce more local helical conformations than the TIP3P model in our simulation. In addition, we also compare the chemical shifts predicted by different chemical shift predicting programs with experimental measurements, the results show that SHIFTX2 obtains the best performance in the chemical shifts prediction.

p53 is a transcription factor with intrinsically disordered regions that plays an essential role in many cellular processes. As a tumor suppressor, the dysfunction of p53 causes various cancers. p53 can be activated by binding with cofactors in the cell due to stresses or DNA damages. The N-terminal transactivation domain (TAD) of p53 can regulate cell apoptosis by interacting with its binding partners, such as the transcriptional adaptor zinc-binding 2 (Taz2) domain of p300, and cofactors, i.e. cyclic-AMP response element-binding protein (CREB)-binding protein (CBP). The experimentally solved structure of p300 Taz2 and p53 TAD2 has provided insights into the interactions that potentially lead to the structural changes of p53 TAD2 and stabilize the complex. To explore the structural changes as well as the residues that lead to such changes from an isolated state to a bound state in p53 TAD2, we used all-atom molecular dynamics (MD) to simulate two different systems: (1) the p300 Taz2–p53 TAD2 complex and (2) isolated p53 TAD2 (residues 35–59). Although still largely unstructured, residues across the p53 TAD2 contribute significantly to stabilizing the binding between p300 Taz2 and p53 TAD2. Our results suggest that the binding affinity of the p300 Taz2–p53 TAD2 complex originates from hydrophobic and electrostatic interactions. The results are in agreement with previous reports and experimental data. By comparing the two simulated systems, our results not only demonstrate the structural changes of p53 TAD2 after binding with Taz2 but also identify the key residues leading to such changes. We also identify the critical residues that can provide insight into the interaction network between p300 Taz2 and p53 TAD2.

A common approach to computing protein pKas uses a continuum dielectric model in which the protein is a low dielectric medium with embedded atomic point charges, the solvent is a high dielectric medium with a Boltzmann distribution of ionic charges, and the pKa is related to the electrostatic free energy which is obtained by solving the Poisson– Boltzmann equation. Starting from the model pKa for a titrating residue, the method obtains the intrinsic pKa and then computes the protonation probability for a given pH including site–site interactions. This approach assumes that acid dissociation does not affect protein conformation aside from adding or deleting charges at titratable sites. In this work, we demonstrate our treecode-accelerated boundary integral (TABI) solver for the relevant electrostatic calculations. The pKa computing procedure is enclosed in a convenient Python wrapper which is publicly available at the corresponding author’s website. Predicted results are compared with experimental pKas for several proteins. Among ongoing efforts to improve protein pKa calculations, the advantage of TABI is that it reduces the numerical errors in the electrostatic calculations so that attention can be focused on modeling assumptions.

Water molecules play a key role in all biochemical processes. They help define the shape of proteins, and they are reactant or product in many reactions and are released as ligands are bound. They facilitate the transfer of protons through transmembrane proton channel, pump and transporter proteins. Continuum electrostatics (CE) force fields used by program Multiconformation CE (MCCE) capture electrostatic interactions in biomolecules with an implicit solvent, which captures the averaged solvent water equilibrium properties. Hybrid CE methods can use explicit water molecules within the protein surrounded by implicit solvent. These hybrid methods permit the study of explicit hydrogen bond networks within the protein and allow analysis of processes such as proton transfer reactions. Yet hybrid CE methods have not been rigorously tested. Here, we present an explicit treatment of water molecules in the Gramicidin A (gA) channel using MCCE and compare the resulting distributions of water molecules and key hydration features against those obtained with explicit solvent Molecular Dynamics (MD) simulations with the nonpolarizable CHARMM36 and polarizable Drude force fields. CHARMM36 leads to an aligned water wire in the channel characterized by a large absolute net water dipole moment; the MCCE and Drude analysis lead to a small net dipole moment as the water molecules change orientation within the channel. The correct orientation is not as yet known, so these calculations identify an open question.

Apolipophorin-III (ApoLp-III) is required for stabilization of molecular shuttles of lipid fuels in insects and is found to contribute to the insect immune reaction. Rearrangement of its five α-helices enables ApoLp-III to reversibly associate with lipids. We investigate computationally the conformational changes of ApoLp-III and the pH-dependence of the binding free energy of ApoLp-III association with a lipid disk. A dominant binding mode along with several minor, low population, modes of the ApoLp-III binding to a lipid disk was identified. The pH-dependence of the binding energy for ApoLp-III with the lipid disk is predicted to be significant, with the pH-optimum at pH= 8.0. The calculations suggest that there are no direct interactions between the lipid head groups and titratable residues of ApoLp-III. In the physiological pH range from 6.0 to 9.0, the binding free energy of ApoLp-III with the lipid disk decreases significantly with respect to its optimal value at pH 8.0 (at pH = 6.0, it is 1.02 kcal/mol and at pH = 9.0 it is 0.23 kcal/mol less favorable than at the optimal pH = 8.0), indicating that the pH is an important regulator of ApoLp-III lipid disk association.

Cell trajectory inference is very important to understand the details of tissue cell development, state differentiation and gene dynamic regulation. However, due to the high noise and heterogeneity of the single-cell data, it is challenging to infer cell trajectory in complex biological processes. Here, we proposed a new trajectory inference method, called Metric Learning Bhattacharyya Kernel Feature Decomposition (MLBKFD). In MLBKFD, a statistical model was used to infer cell trajectory by calculating the Markov transition matrix and Bhattacharyya kernel matrix between cells. Before that, to expedite the matrix calculation in the statistical model, a deep feedforward neural network was used to perform dimensionality reduction on single-cell data. The MLBKFD was evaluated on four typical datasets as well as seven recent human fetal lung datasets. Comparisons with the two outstanding methods (i.e., DTFLOW and MARGARET) demonstrate that the MLBKFD is capable of accurately inferring cell development and differentiation trajectories from single-cell data with different sizes and sources. Notably, MLBKFD exhibits nearly twice the speed of DTFLOW while maintaining high precision, particularly when dealing with large datasets. MLBKFD provides accurate and efficient trajectory inference, empowering researchers to gain deeper insights into the complex dynamics of cell development and differentiation.

Genetically identical cell populations growing in the same environment show a large degree of heterogeneity in gene expression profiles, resulting in significant phenotypic consequences. One of the important sources of this heterogeneity is transcriptional bursting. In recent years, single cell transcriptome sequencing (scRNA-seq) data have been widely used to infer the kinetics of transcriptional bursting. Moreover, the increasing evidence showed that genes are jointly regulated by several competitive pathways when activated by external signals, and there is molecular memory in the process of gene state switching. Therefore, a gene expression model considering both gene activation pathway and molecular memory can better reflect the nature of transcriptional bursting. In this paper, we used an approximate Bayesian computation algorithm called ABC-PRC combined with a gene expression model that considered crosstalk and memory, to infer the kinetics of transcriptional bursting. We analyze the internal relationship between transcriptional bursting and gene expression using scRNA-seq data obtained from mouse embryonic stem cells and mouse embryonic fibroblasts. Both model analysis and data simulations demonstrate that the bimodal coefficient increases with the increase of burst size, and the noise intensity decreases with the increase of burst frequency. Additionally, we show that some previously proposed models can be simplified as special cases of the gene expression model and inference algorithm proposed here under certain circumstances.

Enzymes are protein molecules that play a crucial role in various biological processes in living organisms. They function as catalysts in biological reactions such as digestion, metabolism, DNA replication and other physiological processes. Furthermore, enzymes are widely used in food production, pharmaceuticals and biofuel production. In these industries, they accelerate desired chemical reactions as biocatalysts. Therefore, applying computational methods and data-driven algorithms to predict enzyme properties is essential. Over the past decade, deep learning has made remarkable advancements in science and technology. Deep learning is a subset of machine learning algorithms that rely on artificial neural networks. These algorithms can be employed for supervised, semi-supervised and unsupervised learning. Here, to provide an update on the current literature, we provide an overview of various deep learning algorithms and recent advancements in their application to enzyme science. These applications can generally be categorized into diverse subjects: function prediction, enzyme kinetic parameters prediction, enzyme-substrate identification, condition optimization, thermophilic property prediction, enzyme catalytic site prediction and enzyme design. In conclusion, we discuss the convergence of enzyme science and deep learning, highlighting the potential opportunities and challenges.

Hereditary spastic paraplegias (HSPs) are a genetically heterogeneous collection of neurodegenerative disorders categorized by progressive lower-limb spasticity and frailty. The complex HSP forms are characterized by various neurological features including progressive spastic weakness, urinary sphincter dysfunction, extra pyramidal signs and intellectual disability (ID). The kinesin superfamily proteins (KIFs) are microtubule-dependent molecular motors involved in intracellular transport. Kinesins directionally transport membrane vesicles, protein complexes, and mRNAs along neurites, thus playing important roles in neuronal development and function. Recent genetic studies have identified kinesin mutations in patients with HSPs. In this study, we used the computational approaches to Vidhyanand Mahase and Adebiyi Sobitan contributed equally to this work. investigate the 40 missense mutations associated with HSP and ID in KIF1A and KIF5A. We performed homology modeling to construct the structures of kinesin–microtubule binding domain and kinesin–tubulin complex. We applied structure-based energy calculation methods to determine the effects of missense mutations on protein stability and protein–protein interaction. The results revealed that the most of disease-causing mutations could change the folding free energy of kinesin motor domain and the binding free energy of kinesin–tubulin complex. We found that E253K associated with ID in KIF1A decrease the protein stability of kinesin motor domains. We showed that the HSP mutations located in kinesin–tubulin complex interface, such as K253N and R280C in KIF5A, can destabilize the kinesin–tubulin complex. The computational analysis provides useful information for understanding the roles of kinesin mutations in the development of ID and HSPs.

Background: The multiple endocrine neoplasia type 1 (MEN1) gene located on chromosome 11q13 encodes menin protein. Previously reported mutations were thought to result in loss of function of menin protein and that they are associated with multiple endocrine neoplasia 1 disorder. However, recently menin has also been characterized as an oncosuppressor protein and it was suggested that mutations in it are associated with various other tumors. Studies indicate that the menin protein stimulates the estrogen receptor (ER) that in turn increases the predisposition for inherited breast cancer.

Methods: Here, we used our supervised in-house combinatory in-silico predictor method to investigate the impact of unclassified missense mutations in MEN1 gene found in breast cancer tissue. We also examined the biophysical and biochemical properties to predict the effects of these missense variants on the menin protein stability and interactions. The results are compared with the effects of known pathogenic mutations in menin causing neoplasia.

Results: Our analysis indicates that some of the variants found in breast cancer tissue show similar pattern of destabilizing the menin protein and its interactions as the pathogenic variants associated with neoplasia. Taking together with the results of our in-silico consensus predictor, we classify missense mutations in menin protein found in breast cancer tissue into pathogenic and benign, and thus, suggesting as an indicator for early detection of elevated breast cancer risk.

Rett Syndrome is a rare genetic disorder exclusively seen in girls. Approximately 95% of RTT cases is caused by mutations in the MeCP2 gene which codes for Methyl-CpG-binding protein 2 (MeCP2). In this review, first, a brief introductory review of Rett Syndrome, MeCP2 protein structure and function, mutation types and frequencies, and phenotype– genotype relationships were provided. After that, the current knowledge on the wild-type and mutant MeCP2 protein structure and dynamics as well as its binding to DNA is reviewed. The review particularly focuses on computational (such as molecular dynamics) and experimental (such as electrophoretic mobility shift assays) studies on the MeCP2 binding to different types of DNA as well as the computational and experimental (such as circular dichroism) studies on the stability changes upon mutations. In the end, a brief opinion on future outlook for further computational studies is provided.