A little walk from physical to biological complexity: protein folding and stability

TL;DR

This paper explores how physics principles help understand protein folding and stability, and how this knowledge can be used to study diseases.

Contribution

The paper introduces thermodynamics-based tools for predicting protein stability and resistance to amino acid changes.

Findings

Statistical potentials derived from experimental structures predict protein stability.

These methods help understand protein aggregation in neuronal diseases.

Inverse Boltzmann law is used to develop knowledge-driven bioinformatics tools.

Abstract

As an example of topic where biology and physics meet, we present the issue of protein folding and stability, and the development of thermodynamics-based bioinformatics tools that predict the stability and thermal resistance of proteins and the change of these quantities upon amino acid substitutions. These methods are based on knowledge-driven statistical potentials, derived from experimental protein structures using the inverse Boltzmann law. We also describe an application of these predictors, which contributed to the understanding of the mechanisms of aggregation of a particular protein known to cause a neuronal disease.

Figures

(a) Part of a protein chain. Proteins are sequences of about 50 to more than 10,000 covalently bound amino acids. These are of 20 types: the 20 amino acids that are encoded in the genome. They share a common part (called main chain) and differ by their side chains which carry molecular groups with different physicochemical properties (charged, aromatic, aliphatic, etc). They are represented by uppercase letters. (b) Native (folded) structure of triose phosphate isomerase (Protein DataBank (PDB) [11] code: 7tim). Helices are represented in blue, β-strands in red and loops in magenta.

Stability curves of proteins belonging to the adenylate kinase family and identified by the PDB [11] code 1aky (from Saccharomyces cerevisiae), 1s3g (from Sporosarcina globispora), 1ank (from Escherichia coli) and 1zip (from Geobacillus stearothermophilus). Temperatures are measured in K and folding free energies in kJ/mol. (a) Experimental stability curves; (b) Stability curves predicted by SCooP [13].

Examples of distance potentials obtained from sets of mesostable proteins (blue curves) and thermostable proteins (red curves). The folding free energy (in kcal/mol) is given as a function of inter-residue distance d (in Â). (a) ΔW(R, d, T) potential based on the propensity of the amino acid of type R (arginine) to be at distance d from any other amino acid. (b) ΔW(V, V, d, T) potential based on the propensity of two amino acids of type V (valines) to be separated by a distance d.

Schematic representation of the ANNs used for the parameter identifications. (a) PoPMuSiC and HoTMuSiC: 2-layer ANN (perceptron) with sigmoid activation functions and input neurons encoding T-independent potentials, volume terms and an independent term; (b) Tm-HoTMuSiC: 3-layer ANN, consisting of 3 perceptrons with sigmoid weights; the neurons of the first (second) perceptron encode the high-T (low-T) potentials, and those of the third perceptron the volume and independent terms. The outputs of these three perceptrons are the inputs of another perceptron with polynomial weight functions.

Distribution of the ΔΔG0(Tr) values for all possible amino acid substitutions in a set of 15 protein structures predicted by PoPMuSIC [14, 15].

Tables

Equations

Full text

INTRODUCTION

I.

Molecular and cellular biology were traditionally purely experimental sciences, and it is only recently that theoretical biology – including bioinformatics, biomodeling and biomath-ematics – has started to develop. The need for informatics approaches has become essential with the advent of high-throughput technologies, which generate a huge amount of data in a single experiment, such as full genome sequences or the ensemble of transcribed RNA molecules in a cell population. These technologies are at the basis of omics science, defined as the collective study and characterization of biomolecules: proteomics is the large-scale study of proteins, genomics of genomes, transcriptomics of transcribed RNA, mutatomics of all mutations in proteins or DNA, and so on.

In parallel to omics science, (computational) systems biology has started to attract a lot of attention. It has become clear that the functioning of cells cannot be understood only through the study of individual macromolecules, such as proteins, RNA, and DNA. Indeed, the interactions between all biomolecules present inside and outside the cells are of prime importance, and they are indispensable to understand and model cellular life. In other words, the division of the cell into independent subsystems is a much too strong approximation. Obviously, this has the consequence of dramatically increasing the complexity of the systems to be studied.

An interesting feature of bioinformatics and biomodeling is that application is never far. Computational approaches are increasingly exploited to rationalize existing data and to develop knowledge-based predictors. These are unable to give a unique optimal solution – as the systems are far too complex –, but propose valuable candidate hypotheses, which need to be further analyzed and validated by experimental means. These targeted conjectures save a lot of time as they allow focusing on a few relevant possibilities. In the medical sector computational approches are more and more used in rational drug [1], antibody [2] and vaccine [3] design, or for identifying genome variants that cause human diseases [4, 5]. In the very near future, therapies will be personalized for each patient on the basis of the variant that appear in his genome. On the other hand, enzymes are used for their unique catalyti properties in many biotechnology sectors, and their optimization is an important challeng for setting up efficient and environment-friendly bioprocesses [6–8]. Finally, synthetic biology will certainly gain importance in the future, and modified or synthetic cells can be expecte to be used as little factories that produce specific compounds to be used for example a energy source (e.g. bioethanol, biohydrogen) [9, 10].

We focus in this paper on a specific subdomain of structural bioinformatics, i.e. proteii folding and stability. Proteins are chains of amino acids and a large class of them fold into a well-defined 3-dimensional (3D) structure in physiological conditions (Fig. 1). Generally when correctly folded, they fulfill their biological role (catalysis, molecular transport, sig nal transduction, genome regulation, etc). Protein folding is (usually) a reversible proces: basically ruled by classical physics, in which the 3D structure is stabilized by interaction between the different types of amino acids. The problem is thus well defined and can b tackled with physical methods, specifically through statistical mechanics. However, it i; extremely difficult to actually predict the 3D structure from the amino acid sequence, fo two reasons: the conformational space of a protein is enormous and the energy function that describe inter-residue interactions are not perfectly accurate.

(a) Part of a protein chain. Proteins are sequences of about 50 to more than 10,000 covalently bound amino acids. These are of 20 types: the 20 amino acids that are encoded in the genome. They share a common part (called main chain) and differ by their side chains which carry molecular groups with different physicochemical properties (charged, aromatic, aliphatic, etc). They are represented by uppercase letters. (b) Native (folded) structure of triose phosphate isomerase (Protein DataBank (PDB) [11] code: 7tim). Helices are represented in blue, β-strands in red and loops in magenta.

We present some bioinformatics tools that were developed in our group in view of estimating the thermodynamic and thermal stability of a protein, and the changes that occur upon substitution of one amino acid into another [12–16]. Such tools are very helpful in protein design, when aiming at rationally modifying or optimizing certain protein characteristics, whether in biopharmaceutical, biomedical or biotechnological applications, or simply when attempting to elucidate the complex relationships between protein sequence, structure, stability, interactions and function.

PROTEIN FOLDING AND STABILITY

II.

In general, proteins occur in two distinct structural states in their usual environment, that is the intra‐ or extracellular medium, or the cellular membrane. These are the folded and unfolded state, which correspond to two free energy minima. The folded, biologically active, state has a low entropy but favorable enthalpy, while the unfolded state has a large entropy but unfavorable enthalpy. The (un)folding transition is an equilibrium process, where the folded structure is preferred for a range of conditions including the physiological conditions (e.g. the living temperature and pH of the host organism), and the unfolded structure is preferred for the other conditions. The transition is usually reversible, at least if aggregation or post-translational modifications are overlooked and if large multidomain proteins are not considered. Note that this is true for the class of proteins called structured proteins. Another class consists of proteins that are intrinsically unstructured or only fold in the presence of some biomolecular partner; we will not consider such proteins here.

Proteins play crucial roles in and outside the cells – they do most of the actual work required for cellular life. Mutations in proteins define polymorphisms – which are at the basis of the phenotypic variability between individuals. However, though the large majority of mutations have a neutral effect, some significantly affect the protein structure, stability, interactions or function, and cause diseases.

The stability of a protein structure is defined by its standard folding free energy: [eqn] which depends on all environmental parameters such as the solvent, pH, pressure and tem-perature. If we only consider the dependence on the temperature, the stability curve of a protein [17] has a inverted bell-shape (Fig. 2) and is described by the Gibbs-Helmholtz equation, which is valid for proteins that fold according to a two-state transition: [eqn] where Tm is the melting temperature of the protein, i.e. the temperature at which ΔG^0^ vanishes, [eqn] is the standard folding enthalpy measured at Tm and [eqn] the folding heat capacity assumed to be T-independent. Note that Tm is here the denaturation temperature at high T; there is also a denaturation temperature at low T which we do not consider here as water usually freezes before the protein denatures.

Stability curves of proteins belonging to the adenylate kinase family and identified by the PDB [11] code 1aky (from Saccharomyces cerevisiae), 1s3g (from Sporosarcina globispora), 1ank (from Escherichia coli) and 1zip (from Geobacillus stearothermophilus). Temperatures are measured in K and folding free energies in kJ/mol. (a) Experimental stability curves; (b) Stability curves predicted by SCooP [13].

Two complementary descriptions of protein stability are commonly used: the thermodynamic stability is identified by ΔG^0^(Tr) at room temperature (Tr), whereas the thermal stability is defined by the melting temperature Tm.

Thermal and thermodynamic stabilities are not always correlated, which is obvious from the Gibbs-Helmholtz equation (2.2) when T = Tr. The molecular reason of the lack of correlation between the two stabilities is that the 20 amino acids are of different physicochemical types, and that their energetic contributions have different temperature dependences.

Upon substitution of one amino acid into another, or of several amino acids, the so-called wild-type protein is transformed into a mutant protein, and this has an impact on all its biophysical properties. The folding free energy gets modified as: ΔΔG^0^ = ΔG^0 mutant^ – ΔG^0 wild^, the melting temperature as: [eqn], and similarly for the changes in folding heat capacity [eqn] and enthalpy [eqn]. The changes in thermodynamic and thermal stabilities upon mutations are more – but still imperfectly – correlated than the stabilities themselves, as can be deduced from the relation: [eqn]

The linear correlation between ΔΔG^0^(Tr) and ΔTm values was found to be equal to -0.7 on about 450 experimentally characterized mutations [18].

STATISTICAL POTENTIALS

III.

To estimate the energy of a protein, two types of energy functions are commonly used. The semi-empirical force fields assume an analytic form for the different types of interactions with parameters that are transferred from experiments on small molecules or from quantum chemistry calculations [19]. The statistical potentials do not assume any analytical form, and are based on a coarse-grained representation of protein structure, which is a computational advantage if they have to be applied on an omics scale or in other time-consuming calculations. They have moreover the advantage of considering implicitly the effect of the solvent (i.e. water for soluble proteins). They are obtained by deriving a potential of mean force (PMF) from frequencies of associations of structure and sequence elements in a dataset of known 3D protein structures, using the Boltzmann law [20, 21]. The simplest PMF can be written as: [eqn] where c and s are structure and sequence elements respectively, F represent the relative frequencies of c and/or s, and k is the Boltzmann constant. Sequence elements are single amino acids, amino acid pairs, triplets, etc, and structure elements are inter-residue distances, main chain torsion angle domains, solvent accessibility, etc. More complex potentials can be constructed by considering more than two structure elements and/or sequence elements. Examples of such potentials are given in Fig. 3.

Examples of distance potentials obtained from sets of mesostable proteins (blue curves) and thermostable proteins (red curves). The folding free energy (in kcal/mol) is given as a function of inter-residue distance d (in Â). (a) ΔW(R, d, T) potential based on the propensity of the amino acid of type R (arginine) to be at distance d from any other amino acid. (b) ΔW(V, V, d, T) potential based on the propensity of two amino acids of type V (valines) to be separated by a distance d.

The derivation of these potentials is based on several assumptions. All the conformations accessible to a given protein are known to follow the Boltzmann law, but here it is assumed that the sequence-structure elements observed in the lowest free energy conformation of the whole set of different proteins follow the same law. Another strong assumption is that the unfolded state may be modeled by conformations that are sequence-independent.

A property of these potentials that can be either inconvenient or interesting – according to the issue – is that they depend on the dataset from which they are derived. Such a dependence can lead to unduly biased potentials if the dataset is noisy, redundant, or unbalanced. But it is also possible to take advantage of this dependence. For example, by dividing the protein structure set into a subset of thermostable proteins (Tm ≥ 65°C, average T̄_m_ ≅ 80°C, labeled by Δ) and a subset of mesostable proteins (Tm < 65°C, T̄m ≅ 50°C, labeled by ∇), and deriving the potentials separately from both subsets, we obtain T-dependent potentials [22]: [eqn] with T ≅ T̄m. The potentials derived from mesostable (thermostable) proteins yield energy estimations at low (high) T. We would like to emphasize that this is currently the only way of obtaining T-dependent free energy estimations. These potentials allowed us to get insight into the T-dependences of specific interactions. For example, as shown in Fig. 3, the relative propensity of the positively charged residue R (arginine) to be at a distance d from any other amino acid is strongly T-dependent, whereas the relative weight of the effective interactions between two hydrophobic residues V (valines) is T-independent. Note also that, in principle, one could define (overlapping) sets of protein structures with increasing Tm-values, to get almost continuous T-dependence. However, this is currently impossible due to the limited number of protein structures of experimental determined Tm, which do not ensure reliable statistics.

PREDICTION OF THE STABILITY CURVE OF PROTEINS

IV.

Using the standard and the T-dependent statistical potentials described above, we de-veloped a method called SCooP to predict the Gibbs-Helmholtz ΔG^0^(T) curve of a target protein with known 3D structure [12, 13]. SCooP is thus able to predict all the thermodynamic quantities that characterize the folding transition, namely the melting temperature Tm as well as the folding heat capacity [eqn], enthalpy [eqn] and free energy ΔG^0^(Tr) at room temperature.

These thermodynamic descriptors were obtained from specific combinations of the statistical potentials described in Eqs (3.1, 3.2) with other protein characteristics such as their number of residues N and their total surface area ASA. In particular, the functional form for the folding heat capacity [eqn] was expressed as: [eqn] where αv and βv are parameters to be optimized and the T-dependent potentials [eqn] are defined in Eq.(3.2)^1^. Analogous expressions with slightly different combinations of features were used for Tm and [eqn] [13]. These functional forms were chosen on the basis of biophysical insights and empirical evidence. For example, we considered the first two terms in the right hand side of expression (4.1) on the basis of the experimental observation that ΔCp, N and ASA are quite well correlated. The last terms in (4.1) come from the ΔGp definition : [eqn] and the approximation consisting of expressing the folding entalphy ΔH^0^(T) as a specific combination of the statistical potentials at temperature T weighted with parameters that depend on the number of residues [eqn]. Since the continuous T-dependence of the potentials is not known, we replaced the derivative by the finite difference between thermo‐ and mesostable potentials: [eqn] where ΔT is a constant that was integrated in the parameters.

All the parameters that appear in the functional form of [eqn] were optimized by minimizing the mean square error between the predicted and the experimental quantities of ℓ ≈ 220 proteins with known 3D structure and stability curve: [eqn]. The same procedure was used to identify the parameters appearing in the functional forms of Tm and [eqn] (not shown). The predicted values of [eqn], Tm and [eqn] were then used to draw the full stability curve using the Gibbs-Helmholtz equation (2.2). The SCooP predictor is freely available on the website http://babylone.ulb.ac.be/SCooP.

One of the strengths of the SCooP method is the computation speed of the full T-dependence of the folding free energy, which allows large-scale analyses of the entire structurome. Another of its strengths is its good performance: the linear correlation coefficients between the experimental and the predicted values of Tm, [eqn] and [eqn] for the dataset of about 220 proteins are all three equal to 0.7, as computed with a leave-one-out cross validation procedure. The scores increase up to 0.8 after the exclusion of 10% outliers. Examples of stability curves predicted by SCooP for proteins belonging to the adenylate kinase family are given in Fig. 2b. As seen from the comparison with Fig. 2a in which the experimentally determined stability curves are plotted, the SCooP predictions reproduce rather well the full stability curves.

PREDICTION OF THE STABILITY CHANGES UPON MUTATIONS

V.

Thermodynamic stability

A.

It can be expected to be easier to predict changes in stability upon amino acid substitutions than the actual stability of proteins. Indeed, the large majority of single-site mutations do not affect the structure of the protein, but simply modify its stability, solubility, affinity for biomolecular partners or function. The PoPMuSiC program [14, 15], which is available on http://www.dezyme.com, predicts changes in the thermodynamic stability of a protein on the basis of its 3D structure and statistical potentials. The ΔΔG^0^(Tr) of a point mutation is expressed as the sum of 16 terms: [eqn] where the 13 ΔΔWv terms are changes in folding free energy upon mutation computed with 13 different statistical potentials ΔW defined in Eq.(3.1)^2^. The coefficients αv are taken as sigmoid functions of the solvent accessibility of the mutated residue: [eqn] where ωv, rv, bv and ϕv ∈ ℝ are parameters to be optimized. This A-dependence is justified by the fact that the weight of the interactions vary continuously from the protein core to the surface. The ΔV± are volume terms defined by ΔV± = θ(±ΔV)||ΔV||, where θ(V) is the Heaviside function. They are associated with the difference in volume of the mutant and wild-type amino acids and provide a description of the impact of the creation of a cavity or the accommodation of stress inside the protein structure.

The optimal (ωv, rv, bv, ϕv) parameter values were determined on the basis a dataset of ℓ = 2648 single-site mutations in 131 different proteins with experimentally determined 3D structure and ΔΔG^0^(Tr). For that purpose, an artificial feedforward neural network (ANN) was used, as shown in Fig. 4.a, with the mean square error between the experimental and predicted folding free energy changes upon mutation as cost function: [eqn]. The parameter identification was particularly complex due to the size and roughness of the parameter space, which makes the search for the global minimum complicated and yields a high probability of getting stuck in local minima. Therefore, the initial values of the weights were chosen randomly, the initialization and training processes were repeated 30 times, and the solution reaching the lowest σ-value was chosen.

Schematic representation of the ANNs used for the parameter identifications. (a) PoPMuSiC and HoTMuSiC: 2-layer ANN (perceptron) with sigmoid activation functions and input neurons encoding T-independent potentials, volume terms and an independent term; (b) Tm-HoTMuSiC: 3-layer ANN, consisting of 3 perceptrons with sigmoid weights; the neurons of the first (second) perceptron encode the high-T (low-T) potentials, and those of the third perceptron the volume and independent terms. The outputs of these three perceptrons are the inputs of another perceptron with polynomial weight functions.

The performance of PoPMuSiC was computed in 5-fold cross validation. The linear correlation coefficient between experimental and predicted ΔΔG^0^ (Tr) values is equal to r = 0.6 and increases to 0.8 with 10% outliers suppressed; the mean square error σ is equal to 1.1 kcal/mol and decreases to 0.9 kcal/mol without 10% outliers. This score is comparable or better than that of other approaches. It must be emphasized that PoPMuSiC is extremely fast, and predicts the ΔΔG^0^ (Tr) of all possible point mutations in a medium-size protein in a few minutes. It can thus be applied on a large, structuromics, scale. Of course this high speed comes at the expense of some accuracy.

The distribution of ΔΔG^0^ (Tr) values computed by PoPMuSiC for all possible amino acid substitutions (about 10^5^ mutations in total) in a set of 15 protein structures is shown in Fig. 5 [23]. The large majority of the mutations have positive ΔΔG^0^ (Tr) values and are thus destabilizing with our conventions. This is of course expected, as protein sequences are relatively well optimized for the stability of the native fold. Note, however, that natural proteins are never perfectly optimal, because some residues are optimized for function rather than for stability, but also, probably, because natural evolution has not had the time to reach optimality.

Distribution of the ΔΔG0(Tr) values for all possible amino acid substitutions in a set of 15 protein structures predicted by PoPMuSIC [14, 15].

Thermal stability

B.

As thermodynamic and thermal stabilities are far from perfectly correlated, as clear from Eq.(2.3), it is necessary to develop methods that predict directly ΔTm values, instead of approximating them on the basis of ΔΔG^0^(Tr) predictions. This is the goal of HoTMuSiC [16], which is available on http://www.dezyme.com. The first version of this program uses the relation: [eqn] where N is the number of residues in the protein, a, c ∈ ∈, and the coefficients αv(A) are sigmoid functions of A. This equation is inspired by Eq.(2.3) with all but the first term on the right hand side overlooked, and is thus very similar to Eq.(5.1) for ΔΔG^0^(Tr). The 9 ΔΔWv terms are computed from the first 9 of the 13 statistical potentials used in Eq.(5.1). The optimal parameter values appearing in this equation were identified on the basis of ℓ = 1531 single-site mutations with known 3D structure and ΔTm, the cost function [eqn], and a feedforward ANN (Fig. 4a).

The second HoTMuSiC version is inspired by the full Eq.(2.3), and uses the T-dependent statistical potentials of Eq.(3.2); it requires as input, in addition to the protein structure, the melting temperature [eqn] of the wild-type protein: [eqn]

The 5 [eqn] terms are folding free energy changes computed from the dataset of thermostable (mesostable) proteins using Eq.(3.2)^3^. β^Δ^, β^∇^ and β^V^ are polynomial functions of [eqn] and the number of residues N. Their functional form guessed from Eq.(2.3) is: [eqn] with I = Δ, ∇ or V. The dependence on the number of residues N comes from the enthalpy factor ΔHm in Eq.(2.3), as these two quantities show a good correlation.

To identify the parameters of this second method, a 3-layer ANN is used which is shown schematically in Fig. 4b. The input layer consists of three sets of neurons that encode the mesostable potentials, the thermostable potentials, and the volume and independent terms, respectively. These three sets of neurons are the inputs of three perceptrons, whose outputs are the three neurons of the hidden layer. These are the input neurons of yet another perceptron, which yields a ΔTm-prediction as output. This predicted ΔTm is finally averaged with the ΔTm predicted by the first HoTMuSiC method (Eq.(5.3)). The initialization and identification procedures of all parameters and the cross validation procedure are the same as for the first method.

The root mean square error between experimental and predicted ΔTm-values is equal to σ = 4.6°C for HoTMuSiC and σ = 4.5°C for Tm-HoTMuSiC; the associated Pearson correlation coefficients r are both equal to 0.6. When 10% outliers are excluded, σ decreases to 3.3 and 3.2°C and r rises to 0.7.

Finally note that, because of the two volume terms and the energy-independent term in Eqs (5.1,5.3,5.4), the symmetry: [eqn] is explicitly broken, and similarly for ΔΔG^0^(Tr) [23]. Due to these symmetry-breaking terms and the fact that the large majority of experimentally characterized mutations in the training dataset are destabilizing, the predictions of ΔΔG^0^(Tr) and ΔTm are on the average more accurate for destabilizing mutations, at the price of a decrease in the prediction performance of stabilizing mutations. But actually, the stabilizing mutations are the most interesting ones in most protein design and optimization applications. To solve this problem and achieve higher prediction scores for stabilizing mutations, we are currently developing new PoPMuSiC and HoTMuSiC versions that preserve the symmetry of Eq.(5.6).

APPLICATION TO A CONFORMATIONAL DISEASE PROTEIN

VI.

To illustrate the power of the above described approaches, we describe one of their suc-cessful applications in the context of conformational diseases [24]. Such disorders are caused by the misfolding of specific proteins, which become inactive or toxic [25, 26]. Often these proteins form soluble or insoluble aggregates. Creutzfeldt-Jakob and Alzheimer are well known examples of this class of diseases. Misfolded conformations correspond to free energy minima that differ from those of the folded and unfolded states. They may be thermodynamically or kinetically favored under certain conditions or for certain protein variants. The precise mechanism by which misfolding is initiated and propagates to other proteins is only partially known.

Spinocerebellar ataxia type 3 is a less known example of conformational disease [27, 28]. It is characterized by progressive neuronal dysfunction and the presence of neuronal inclusions containing fibrillar aggregates of ataxin-3 proteins. Ataxin-3 contains a globular, structured domain, called Josephin, and a flexible polyglutamine (polyQ) domain whose repeat-length modulates pathogenicity. It has been suggested that the fibrillogenesis pathway of ataxin-3 starts with a polyQ-independent step mediated by Josephin domain interactions, followed by a polyQ-dependent step. To test the involvement of the Josephin domain in ataxin-3 fibrillogenesis, the isolated Josephin domain was in turn stabilized and destabilized, and the consequences on the aggregation tendencies of this domain taken individually and of the full ataxin-3 protein were analyzed. Two ataxin-3 variants were considered, a non-pathological variant with a 15-residue-long polyQ tract (ataxin-3(Q15)), and a pathological variant with a 64-glutamine tract (ataxin-3(Q64)) [24].

To modify the stability of the Josephin domain, the PoPMuSiC software was used. Two mutations were selected, which were predicted as the most stabilizing and are located suffi-ciently far from the active site to avoid affecting the protein’s function: R103G and S81A^4^. Furthermore, four destabilizing mutations were chosen in silico, but only one resulted in soluble protein expression: L169H. The Tm and aggregation propensities of the wild-type and of these 3 protein mutants, as well as of the double mutant R103G/S81A, were experimentally measured. As shown in Table I, the mutations predicted to stabilize the Josephin domain actually do so. They also tend to stabilize the complete ataxin-3 protein. Moreover, these mutants significantly decrease the aggregation propensities: the aggregation midpoint is increased both for the Josephin domain and the two ataxin-3 variants Q15 and Q64. Moreover, the mutation predicted to destabilize the Josephin domain actually does so and moreover drastically decreases the aggregation midpoint.

Table I:: Predicted stabilizing and destabilizing mutations in the Josephin domain and experimental characterization of the thermal denaturation and aggregation midpoints (Tm and Agg) of the Josephin domain and two ataxin-3 variants of different polyQ length (Q15 and Q64). No results are shown for the L169H mutant in the two ataxin-3 variants because soluble expression was not achieved.

These results show that changing the thermodynamic stability of the Josephin domain modulates ataxin-3 fibrillogenesis, and support the hypothesis that the first stage of ataxin-3 fibrillogenesis is caused by interactions involving the non-polyQ containing Josephin domain. They also show that the use of in silico predictors such as PoPMuSiC is very helpful to detect very rapidly relevant candidate mutations, which can be exploited to optimize a protein system or to get insights into biological or biomedical issues.