PyMC: a modern, and comprehensive probabilistic programming framework in Python.

PyMC is a probabilistic programming library for Python that provides tools for constructing and fitting Bayesian models. It offers an intuitive, readable syntax that is close to the natural syntax statisticians use to describe models. PyMC leverages the symbolic computation library PyTensor, allowing it to be compiled into a variety of computational backends, such as C, JAX, and Numba, which in turn offer access to different computational architectures including CPU, GPU, and TPU. Being a general modeling framework, PyMC supports a variety of models including generalized hierarchical linear regression and classification, time series, ordinary differential equations (ODEs), and non-parametric models such as Gaussian processes (GPs). We demonstrate PyMC's versatility and ease of use with examples spanning a range of common statistical models. Additionally, we discuss the positive role of PyMC in the development of the open-source ecosystem for probabilistic programming.

Exploring the quality of protein structural models from a Bayesian perspective.

We explore how ideas and practices common in Bayesian modeling can be applied to help assess the quality of 3D protein structural models. The basic premise of our approach is that the evaluation of a Bayesian statistical model's fit may reveal aspects of the quality of a structure when the fitted data is related to protein structural properties. Therefore, we fit a Bayesian hierarchical linear regression model to experimental and theoretical 13 Cα chemical shifts. Then, we propose two complementary approaches for the evaluation of such fitting: (a) in terms of the expected differences between experimental and posterior predicted values; (b) in terms of the leave-one-out cross-validation point-wise predictive accuracy. Finally, we present visualizations that can help interpret these evaluations. The analyses presented in this article are aimed to aid in detecting problematic residues in protein structures. The code developed for this work is available on: https://github.com/BIOS-IMASL/Hierarchical-Bayes-NMR-Validation.

The Marginal Stability of Proteins: How the Jiggling and Wiggling of Atoms is Connected to Neutral Evolution.

Here we propose that the upper bound marginal stability of proteins is a universal property that includes macro-molecular complexes and is not affected by molecular changes such as mutations and post-translational modifications. We theorize that its existence is a consequence of Afinsen's thermodynamic hypothesis rather than a result of an evolutionary process. This result enables us to conjecture that neutral evolution should also be, with respect to protein stability, a universal phenomenon.

Assessing the One-Bond Cα-H Spin-Spin Coupling Constants in Proteins: Pros and Cons of Different Approaches.

In the present work, we explore three different approaches for the computation of the one-bond spin-spin coupling constants (SSCC) 1JCαH in proteins: density functional theory (DFT) calculations, a Karplus-like equation, and Gaussian process regression. The main motivation of this work is to select the best method for fast and accurate computation of the 1JCαH SSCC, for its use in everyday applications in protein structure validation, refinement, and/or determination. Our initial results showed a poor agreement between the DFT-computed and observed 1JCαH SSCC values. Further analysis leads us to the understanding that the model chosen for the DFT computations is inappropriate and that more complex models will require a higher, if not prohibitively, computational cost. Finally, we show that the Karplus-like equation and Gaussian Process regression provide faster and more accurate results than DFT-based calculations.

Classification of RNA backbone conformations into rotamers using 13C' chemical shifts: exploring how far we can go.

The conformational space of the ribose-phosphate backbone is very complex as it is defined in terms of six torsional angles. To help delimit the RNA backbone conformational preferences, 46 rotamers have been defined in terms of these torsional angles. In the present work, we use the ribose experimental and theoretical 13C' chemical shifts data and machine learning methods to classify RNA backbone conformations into rotamers and families of rotamers. We show to what extent the experimental 13C' chemical shifts can be used to identify rotamers and discuss some problem with the theoretical computations of 13C' chemical shifts.

Nature representation in South American protected areas: country contrasts and conservation priorities.

BACKGROUND: South America faces strong environmental pressures as a result of agriculture and infrastructure expansion and also of demographic growth, demanding immediate action to preserve natural assets by establishing protected areas. Currently, 7.1% of the (sub)continent is under strict conservation categories (I to IV, IUCN), but the spatial distribution of these 1.3 × 106 km2 is poorly understood. We evaluated the representation of nature within the networks of protected areas, map conservation priorities and assess demographic, economic or geopolitical causes of existing protection patterns. METHODS: We characterized nature representation by looking at two components: the extent and the equality of protection. The first refers to the fraction of territory under protection, while the second refers to the homogeneity in the distribution along natural conditions of this protected fraction. We characterized natural conditions by either 113 biogeographical units (specifically, ecoregions) or a series of limited and significant climatic, topographic and edaphic traits. We analyzed representation every ten years since 1960 at national and continental levels. In the physical approach, histograms allowed us to map the degree of conservation priorities. Finally, we ranked the importance of different economic or geopolitical variables driving the observed distributions with a random forest technique. RESULTS: Nature representation varied across countries in spite of its priority in conservation agendas. In Brazil, Peru and Argentina there are still natural conditions with no formal protection, while in Bolivia and Venezuela, protected areas incorporate the natural diversity in a more balanced manner. As protected networks have increased their extent, so did their equality across and within countries over time. Our maps revealed as top continental priorities the southern temperate, subhumid and fertile lowland environments, and other country-specific areas. Protection extent was generally driven by a low population density and isolation, while other variables like distance to frontiers, were relevant only locally (e.g., in Argentina). DISCUSSION: Our description of the spatial distribution of protected areas can help societies and governments to improve the allocation of conservation efforts. We identified the main limitations that future conservation efforts will face, as protection was generally driven by the opportunities provided by low population density and isolation. From a methodological perspective, the physical approach reveals new properties of protection and provides tools to explore nature representation at different spatial, temporal and conceptual levels, complementing the traditional ones based on biodiversity or biogeographical attributes.

Outline of an experimental design aimed to detect a protein A mirror image in solution.

There is abundant theoretical evidence indicating that a mirror image of Protein A may occur during the protein folding process. However, as to whether such mirror image exists in solution is an unsolved issue. Here we provide outline of an experimental design aimed to detect the mirror image of Protein A in solution. The proposal is based on computational simulations indicating that the use of a mutant of protein A, namely Q10H, could be used to detect the mirror image conformation in solution. Our results indicate that the native conformation of the protein A should have a pKa, for the Q10H mutant, at ≈6.2, while the mirror-image conformation should have a pKa close to ≈7.3. Naturally, if all the population is in the native state for the Q10H mutant, the pKa should be ≈6.2, while, if all are in the mirror-image state, it would be ≈7.3, and, if it is a mixture, the pKa should be largerthan 6.2, presumably in proportion to the mirror population. In addition, evidence is provided indicating the tautomeric distribution of H10 must also change between the native and mirror conformations. Although this may not be completely relevant for the purpose of determining whether the protein A mirror image exists in solution, it could provide valuable information to validate the pKa findings. We hope this proposal will foster experimental work on this problem either by direct application of our proposed experimental design or serving as inspiration and motivation for other experiments.

Opportunities drive the global distribution of protected areas.

BACKGROUND: Protected areas, regarded today as a cornerstone of nature conservation, result from an array of multiple motivations and opportunities. We explored at global and regional levels the current distribution of protected areas along biophysical, human, and biological gradients, and assessed to what extent protection has pursued (i) a balanced representation of biophysical environments, (ii) a set of preferred conditions (biological, spiritual, economic, or geopolitical), or (iii) existing opportunities for conservation regardless of any representation or preference criteria. METHODS: We used histograms to describe the distribution of terrestrial protected areas along biophysical, human, and biological independent gradients and linear and non-linear regression and correlation analyses to describe the sign, shape, and strength of the relationships. We used a random forest analysis to rank the importance of different variables related to conservation preferences and opportunity drivers, and an evenness metric to quantify representativeness. RESULTS: We find that protection at a global level is primarily driven by the opportunities provided by isolation and a low population density (variable importance = 34.6 and 19.9, respectively). Preferences play a secondary role, with a bias towards tourism attractiveness and proximity to international borders (variable importance = 12.7 and 3.4, respectively). Opportunities shape protection strongly in "North America & Australia-NZ" and "Latin America & Caribbean," while the importance of the representativeness of biophysical environments is higher in "Sub-Saharan Africa" (1.3 times the average of other regions). DISCUSSION: Environmental representativeness and biodiversity protection are top priorities in land conservation agendas. However, our results suggest that they have been minor players driving current protection at both global and regional levels. Attempts to increase their relevance will necessarily have to recognize the predominant opportunistic nature that the establishment of protected areas has had until present times.

Detection of methylation, acetylation and glycosylation of protein residues by monitoring (13)C chemical-shift changes: A quantum-chemical study.

Post-translational modifications of proteins expand the diversity of the proteome by several orders of magnitude and have a profound effect on several biological processes. Their detection by experimental methods is not free of limitations such as the amount of sample needed or the use of destructive procedures to obtain the sample. Certainly, new approaches are needed and, therefore, we explore here the feasibility of using (13)C chemical shifts of different nuclei to detect methylation, acetylation and glycosylation of protein residues by monitoring the deviation of the (13)C chemical shifts from the expected (mean) experimental value of the non-modified residue. As a proof-of-concept, we used (13)C chemical shifts, computed at the DFT-level of theory, to test this hypothesis. Moreover, as a validation test of this approach, we compare our theoretical computations of the (13)Cε chemical-shift values against existing experimental data, obtained from NMR spectroscopy, for methylated and acetylated lysine residues with good agreement within ∼1 ppm. Then, further use of this approach to select the most suitable (13)C-nucleus, with which to determine other modifications commonly seen, such as methylation of arginine and glycosylation of serine, asparagine and threonine, shows encouraging results.

Azahar: a PyMOL plugin for construction, visualization and analysis of glycan molecules.

Glycans are key molecules in many physiological and pathological processes. As with other molecules, like proteins, visualization of the 3D structures of glycans adds valuable information for understanding their biological function. Hence, here we introduce Azahar, a computing environment for the creation, visualization and analysis of glycan molecules. Azahar is implemented in Python and works as a plugin for the well known PyMOL package (Schrodinger in The PyMOL molecular graphics system, version 1.3r1, 2010). Besides the already available visualization and analysis options provided by PyMOL, Azahar includes 3 cartoon-like representations and tools for 3D structure caracterization such as a comformational search using a Monte Carlo with minimization routine and also tools to analyse single glycans or trajectories/ensembles including the calculation of radius of gyration, Ramachandran plots and hydrogen bonds. Azahar is freely available to download from http://www.pymolwiki.org/index.php/Azahar and the source code is available at https://github.com/agustinaarroyuelo/Azahar .

Factors affecting the computation of the 13C shielding in disaccharides.

Knowledge of the three-dimensional structures of glycans and glycoproteins is useful for a full understanding of molecular processes in which glycans are involved, such as antigen-recognition and virus infection, to name a few. Among the ubiquitous nuclei in glycan molecules, the (13)C nucleus is an attractive candidate for computation of theoretical chemical shifts at the quantum chemical level of theory to validate and determine glycan structures. For this purpose, it is important to determine, first, which carbons can be used as probes to sense conformational changes and, second, all factors that affect the computation of the shielding, at the density functional theory (DFT) level of theory, of those carbons. To answer such questions, we performed a series of analyses on low-energy conformations, obtained by sampling the glycosidic torsional angles (ϕ, ψ) every 10°, of 12 disaccharides. Our results provide evidence that: (i) the carbons that participate in the glycosidic linkage are the most sensitive probes with which to sense conformational changes of disaccharides; (ii) the rotation of the hydroxyl groups closest to the glycosidic linkage significantly affects the computation of the shieldings of the carbons that participate in the glycosidic linkage; (iii) it is not possible to obtain the shieldings of one disaccharide from the computed values of a different disaccharide or from those disaccharides that differ in the anomeric state; and (iv) a proper basis set distribution, a functional, and a step size, with which to sample the conformational space, are necessary to compute shieldings accurately and rapidly.

Accounting for a mirror-image conformation as a subtle effect in protein folding.

By using local (free-energy profiles along the amino acid sequence and (13)C(α) chemical shifts) and global (principal component) analyses to examine the molecular dynamics of protein-folding trajectories, generated with the coarse-grained united-residue force field, for the B domain of staphylococcal protein A, we are able to (i) provide the main reason for formation of the mirror-image conformation of this protein, namely, a slow formation of the second loop and part of the third helix (Asp29-Asn35), caused by the presence of multiple local conformational states in this portion of the protein; (ii) show that formation of the mirror-image topology is a subtle effect resulting from local interactions; (iii) provide a mechanism for how protein A overcomes the barrier between the metastable mirror-image state and the native state; and (iv) offer a plausible reason to explain why protein A does not remain in the metastable mirror-image state even though the mirror-image and native conformations are at least energetically compatible.

Are accurate computations of the 13C' shielding feasible at the DFT level of theory?

The goal of this study is twofold. First, to investigate the relative influence of the main structural factors affecting the computation of the (13)C' shielding, namely, the conformation of the residue itself and the next nearest-neighbor effects. Second, to determine whether calculation of the (13)C' shielding at the density functional level of theory (DFT), with an accuracy similar to that of the (13)C(α) shielding, is feasible with the existing computational resources. The DFT calculations, carried out for a large number of possible conformations of the tripeptide Ac-GXY-NMe, with different combinations of X and Y residues, enable us to conclude that the accurate computation of the (13)C' shielding for a given residue X depends on the: (i) (ϕ,ψ) backbone torsional angles of X; (ii) side-chain conformation of X; (iii) (ϕ,ψ) torsional angles of Y; and (iv) identity of residue Y. Consequently, DFT-based quantum mechanical calculations of the (13)C' shielding, with all these factors taken into account, are two orders of magnitude more CPU demanding than the computation, with similar accuracy, of the (13)C(α) shielding. Despite not considering the effect of the possible hydrogen bond interaction of the carbonyl oxygen, this work contributes to our general understanding of the main structural factors affecting the accurate computation of the (13)C' shielding in proteins and may spur significant progress in effort to develop new validation methods for protein structures.

Physics-based method to validate and repair flaws in protein structures.

A method that makes use of information provided by the combination of (13)C(α) and (13)C(ß) chemical shifts, computed at the density functional level of theory, enables one to (i) validate, at the residue level, conformations of proteins and detect backbone or side-chain flaws by taking into account an ensemble average of chemical shifts over all of the conformations used to represent a protein, with a sensitivity of ∼90%; and (ii) provide a set of (χ1/χ2) torsional angles that leads to optimal agreement between the observed and computed (13)C(α) and (13)C(ß) chemical shifts. The method has been incorporated into the CheShift-2 protein validation Web server. To test the reliability of the provided set of (χ1/χ2) torsional angles, the side chains of all reported conformations of five NMR-determined protein models were refined by a simple routine, without using NOE-based distance restraints. The refinement of each of these five proteins leads to optimal agreement between the observed and computed (13)C(α) and (13)C(ß) chemical shifts for ∼94% of the flaws, on average, without introducing a significantly large number of violations of the NOE-based distance restraints for a distance range ≤ 0.5 , in which the largest number of distance violations occurs. The results of this work suggest that use of the provided set of (χ1/χ2) torsional angles together with other observables, such as NOEs, should lead to a fast and accurate refinement of the side-chain conformations of protein models.

CheShift-2: graphic validation of protein structures.

UNLABELLED: The differences between observed and predicted (13)C(α) chemical shifts can be used as a sensitive probe with which to detect possible local flaws in protein structures. For this reason, we previously introduced CheShift, a Web server for protein structure validation. Now, we present CheShift-2 in which a graphical user interface is implemented to render such local flaws easily visible. A series of applications to 15 ensembles of conformations illustrate the ability of CheShift-2 to locate the main structural flaws rapidly and accurately on a per-residue basis. Since accuracy plays a central role in CheShift predictions, the treatment of histidine (His) is investigated here by exploring which form of His should be used in CheShift-2. AVAILABILITY: CheShift-2 is free of charge for academic use and can be accessed from www.cheshift.com

Toxoplasma gondii Sis1-like J-domain protein is a cytosolic chaperone associated to HSP90/HSP70 complex.

Toxoplasma gondii is an obligate intracellular protozoan parasite in which 36 predicted Hsp40 family members were identified by searching the T. gondii genome. The predicted protein sequence from the gene ID TGME49_065310 showed an amino acid sequence and domain structure similar to Saccharomyces cerevisiae Sis1. TgSis1 did not show differences in its expression profile during alkaline stress by microarray analysis. Furthermore, TgSis1 showed to be a cytosolic Hsp40 which co-immunoprecipitated with T. gondii Hsp70 and Hsp90. Structural modeling of the TgSis1 peptide binding fragment revealed structural and electrostatic properties different from the experimental model of human Sis1-like protein (Hdj1). Based on these differences; we propose that TgSis1 may be a potentially attractive drug target for developing a novel anti-T. gondii therapy.

In silico study of the inhibition of DNA polymerase by a novel catalpol derivative.

In this work, a novel catalpol derivative (6,10,2',6'-tetraacetyl-O-catalpol), which was previously obtained by our group and shown experimentally to inhibit a type of Taq DNA polymerase, was studied in silico. Studies of the interaction of 6,10,2',6'-tetraacetyl-O-catalpol with the Klentaq fragment of the Taq DNA polymerase I from Thermus aquaticus helped to elucidate the mechanism of inhibition of the enzyme, and offered valuable information that can be used to propose substrate structural modifications aimed at increasing the binding affinity. Classical and semi-empirical methods were used to characterize the conformational preferences of this organic compound in solution. Using docking simulations, the most probable binding mode was found, and the stabilities of the docked solutions were tested in a series of molecular dynamics experiments. Results indicated that the mechanism of inhibition may be competitive, which agrees with previous binding experiments done with 6,10,2',6'-tetraacetyl-O-catalpol.

Analysis of 13Calpha and 13Cbeta chemical shifts of cysteine and cystine residues in proteins: a quantum chemical approach.

Cysteines possess a unique property among the 20 naturally occurring amino acids: it can be present in proteins in either the reduced or oxidized form, and can regulate the activity of some proteins. Consequently, to augment our previous treatment of the other types of residues, the 13Calpha and 13Cbeta chemical shifts of 837 cysteines in disulfide-bonded cystine from a set of seven non-redundant proteins, determined by X-ray crystallography and NMR spectroscopy, were computed at the DFT level of theory. Our results indicate that the errors between observed and computed 13Calpha chemical shifts of such oxidized cysteines can be attributed to several effects such as: (a) the quality of the NMR-determined models, as evaluated by the conformational-average (ca) rmsd value; (b) the existence of high B-factor or crystal-packing effects for the X-ray-determined structures; (c) the dynamics of the disulfide bonds in solution; and (d) the differences in the experimental conditions under which the observed 13Calpha chemical shifts and the protein models were determined by either X-ray crystallography or NMR-spectroscopy. These quantum-chemical-based calculations indicate the existence of two, almost non-overlapped, basins for the oxidized and reduced -SH 13Cbeta, but not for the 13Calpha, chemical shifts, in good agreement with the observation of 375 13Calpha and 337 13Cbeta resonances from 132 proteins by Sharma and Rajarathnam (2000). Overall, our results indicate that explicit consideration of the disulfide bonds is a necessary condition for an accurate prediction of 13Calpha and 13Cbeta chemical shifts of cysteines in cystines.

Quantum-mechanics-derived 13Calpha chemical shift server (CheShift) for protein structure validation.

A server (CheShift) has been developed to predict (13)C(alpha) chemical shifts of protein structures. It is based on the generation of 696,916 conformations as a function of the phi, psi, omega, chi1 and chi2 torsional angles for all 20 naturally occurring amino acids. Their (13)C(alpha) chemical shifts were computed at the DFT level of theory with a small basis set and extrapolated, with an empirically-determined linear regression formula, to reproduce the values obtained with a larger basis set. Analysis of the accuracy and sensitivity of the CheShift predictions, in terms of both the correlation coefficient R and the conformational-averaged rmsd between the observed and predicted (13)C(alpha) chemical shifts, was carried out for 3 sets of conformations: (i) 36 x-ray-derived protein structures solved at 2.3 A or better resolution, for which sets of (13)C(alpha) chemical shifts were available; (ii) 15 pairs of x-ray and NMR-derived sets of protein conformations; and (iii) a set of decoys for 3 proteins showing an rmsd with respect to the x-ray structure from which they were derived of up to 3 A. Comparative analysis carried out with 4 popular servers, namely SHIFTS, SHIFTX, SPARTA, and PROSHIFT, for these 3 sets of conformations demonstrated that CheShift is the most sensitive server with which to detect subtle differences between protein models and, hence, to validate protein structures determined by either x-ray or NMR methods, if the observed (13)C(alpha) chemical shifts are available. CheShift is available as a web server.

What can we learn by computing 13Calpha chemical shifts for X-ray protein models?

The room-temperature X-ray structures of ubiquitin (PDB code 1ubq) and of the RNA-binding domain of nonstructural protein 1 of influenza A virus (PDB code 1ail) solved at 1.8 and 1.9 A resolution, respectively, were used to investigate whether a set of conformations rather than a single X-ray structure provides better agreement with both the X-ray data and the observed 13Calpha chemical shifts in solution. For this purpose, a set of new conformations for each of these proteins was generated by fitting them to the experimental X-ray data deposited in the PDB. For each of the generated structures, which show R and Rfree factors similar to those of the deposited X-ray structure, the 13Calpha chemical shifts of all residues in the sequence were computed at the DFT level of theory. The sets of conformations were then evaluated by their ability to reproduce the observed 13Calpha chemical shifts by using the conformational average root-mean-square-deviation (ca-r.m.s.d.). For ubiquitin, the computed set of conformations is a better representation of the observed 13Calpha chemical shifts in terms of the ca-r.m.s.d. than a single X-ray-derived structure. However, for the RNA-binding domain of nonstructural protein 1 of influenza A virus, consideration of an ensemble of conformations does not improve the agreement with the observed 13Calpha chemical shifts. Whether an ensemble of conformations rather than any single structure is a more accurate representation of a protein structure in the crystal as well as of the observed 13Calpha chemical shifts is determined by the dispersion of coordinates, in terms of the all-atom r.m.s.d. among the generated models; these generated models satisfy the experimental X-ray data with accuracy as good as the PDB structure. Therefore, generation of an ensemble is a necessary step to determine whether or not a single structure is sufficient for an accurate representation of both experimental X-ray data and observed 13Calpha chemical shifts in solution.