The role of acidic amino acid residues in the structural stability of snake cardiotoxins

Specific Amino Acid Residues in the Three Loops of Snake Cytotoxins Determine Their Membrane Activity and Provide a Rationale for a New Classification of These Toxins

Cytotoxins (CTs) are three-finger membrane-active toxins present mainly in cobra venom. Our analysis of the available CT amino acid sequences, literature data on their membrane activity, and conformational equilibria in aqueous solution and detergent micelles allowed us to identify specific amino acid residues which interfere with CT incorporation into membranes. They include Pro9, Ser28, and Asn/Asp45 within the N-terminal, central, and C-terminal loops, respectively. There is a hierarchy in the effect of these residues on membrane activity: Pro9 > Ser28 > Asn/Asp45. Taking into account all the possible combinations of special residues, we propose to divide CTs into eight groups. Group 1 includes toxins containing all of the above residues. Their representatives demonstrated the lowest membrane activity. Group 8 combines CTs that lack these residues. For the toxins from this group, the greatest membrane activity was observed. We predict that when solely membrane activity determines the cytotoxic effects, the activity of CTs from a group with a higher number should exceed that of CTs from a group with a lower number. This classification is supported by the available data on the cytotoxicity and membranotropic properties of CTs. We hypothesize that the special amino acid residues within the loops of the CT molecule may indicate their involvement in the interaction with non-lipid targets.

Best Practices in Constant pH MD Simulations: Accuracy and Sampling

Various approaches have been proposed to include the effect of pH in molecular dynamics (MD) simulations. Among these, the λ-dynamics approach proposed by Brooks and co-workers [Kong, X.; Brooks III, C. L. J. Chem. Phys.1996, 105, 2414−2423] can be performed with little computational overhead and hfor each typeence be used to routinely perform MD simulations at microsecond time scales, as shown in the accompanying paper [Aho, N. et al. J. Chem. Theory Comput.2022, DOI: 10.1021/acs.jctc.2c00516]. At such time scales, however, the accuracy of the molecular mechanics force field and the parametrization becomes critical. Here, we address these issues and provide the community with guidelines on how to set up and perform long time scale constant pH MD simulations. We found that barriers associated with the torsions of side chains in the CHARMM36m force field are too high for reaching convergence in constant pH MD simulations on microsecond time scales. To avoid the high computational cost of extending the sampling, we propose small modifications to the force field to selectively reduce the torsional barriers. We demonstrate that with such modifications we obtain converged distributions of both protonation and torsional degrees of freedom and hence consistent pK_a estimates, while the sampling of the overall configurational space accessible to proteins is unaffected as compared to normal MD simulations. We also show that the results of constant pH MD depend on the accuracy of the correction potentials. While these potentials are typically obtained by fitting a low-order polynomial to calculated free energy profiles, we find that higher order fits are essential to provide accurate and consistent results. By resolving problems in accuracy and sampling, the work described in this and the accompanying paper paves the way to the widespread application of constant pH MD beyond pK_a prediction.

Extensive Assessment of Various Computational Methods for Aspartate's pKa Shift

A series of computational methods for pK_a shift prediction are extensively tested on a set of benchmark protein systems, aiming at identifying pitfalls and evaluating their performance on high variants. Including 19 ASP residues in 10 protein systems, the benchmark set consists of both residues with highly shifted pK_a values as well as those varying little from the reference value, with an experimental RMS free energy differences of 2.49 kcal/mol with respect to blocked amino acid, namely the RMS pK_a shift being 1.82 pK_a units. The constant pH molecular dynamics (MD), alchemical methods, PROPKA3.1, and multiconformation continuum electrostatics give RMSDs of 1.52, 2.58, 1.37, and 3.52 pK_a units, respectively, on the benchmark set. The empirical scoring method is the most accurate one with extremely low computational cost, and the pH-dependent model is also able to provide accurate results, while the accuracy of MD sampling incorporating alchemical free energy simulation is prohibited by convergence achievement and the performance of conformational search incorporating multiconformation continuum electrostatics is bad. Former research works did not define statistical uncertainty with care and yielded the questionable conclusion that alchemical methods perform well in most benchmarks. In this work the traditional alchemical methods are thoroughly tested for high variants. We also performed the first application of nonequilibrium alchemical methods to the pK_a cases.

Enhancing constant-pH simulation in explicit solvent with a two-dimensional replica exchange method

We present a new method for enhanced sampling for constant-pH simulations in explicit water based on a two-dimensional (2D) replica exchange scheme. The new method is a significant extension of our previously developed constant-pH simulation method, which is based on enveloping distribution sampling (EDS) coupled with a one-dimensional (1D) Hamiltonian exchange method (HREM). EDS constructs a hybrid Hamiltonian from multiple discrete end state Hamiltonians that, in this case, represent different protonation states of the system. The ruggedness and heights of the hybrid Hamiltonian's energy barriers can be tuned by the smoothness parameter. Within the context of the 1D EDS-HREM method, exchanges are performed between replicas with different smoothness parameters, allowing frequent protonation-state transitions and sampling of conformations that are favored by the end-state Hamiltonians. In this work, the 1D method is extended to 2D with an additional dimension, external pH. Within the context of the 2D method (2D EDS-HREM), exchanges are performed on a lattice of Hamiltonians with different pH conditions and smoothness parameters. We demonstrate that both the 1D and 2D methods exactly reproduce the thermodynamic properties of the semigrand canonical (SGC) ensemble of a system at a given pH. We have tested our new 2D method on aspartic acid, glutamic acid, lysine, a four residue peptide (sequence KAAE), and snake cardiotoxin. In all cases, the 2D method converges faster and without loss of precision; the only limitation is a loss of flexibility in how CPU time is employed. The results for snake cardiotoxin demonstrate that the 2D method enhances protonation-state transitions, samples a wider conformational space with the same amount of computational resources, and converges significantly faster overall than the original 1D method.

Constant pH Molecular Dynamics in Explicit Solvent with Enveloping Distribution Sampling and Hamiltonian Exchange

We present a new computational approach for constant pH simulations in explicit solvent based on the combination of the enveloping distribution sampling (EDS) and Hamiltonian replica exchange (HREX) methods. Unlike constant pH methods based on variable and continuous charge models, our method is based on discrete protonation states. EDS generates a hybrid Hamiltonian of different protonation states. A smoothness parameter s is used to control the heights of energy barriers of the hybrid-state energy landscape. A small s value facilitates state transitions by lowering energy barriers. Replica exchange between EDS potentials with different s values allows us to readily obtain a thermodynamically accurate ensemble of multiple protonation states with frequent state transitions. The analysis is performed with an ensemble obtained from an EDS Hamiltonian without smoothing, s = ∞, which strictly follows the minimum energy surface of the end states. The accuracy and efficiency of this method is tested on aspartic acid, lysine, and glutamic acid, which have two protonation states, a histidine with three states, a four-residue peptide with four states, and snake cardiotoxin with eight states. The pKa values estimated with the EDS-HREX method agree well with the experimental pKa values. The mean absolute errors of small benchmark systems range from 0.03 to 0.17 pKa units, and those of three titratable groups of snake cardiotoxin range from 0.2 to 1.6 pKa units. This study demonstrates that EDS-HREX is a potent theoretical framework, which gives the correct description of multiple protonation states and good calculated pKa values.

Fusogenicity of Naja naja atra cardiotoxin-like basic protein on sphingomyelin vesicles containing oxidized phosphatidylcholine and cholesterol

This study investigated the effect of oxidized phosphatidylcholine (oxPC) and cholesterol (Chol) on Naja naja atra cardiotoxin-like basic protein (CLBP)-induced fusion and leakage in sphingomyelin (SM) vesicles. Compared with those on PC/SM/Chol vesicles, CLBP showed a lower activity to induce membrane permeability but a higher fusogenicity on oxPC/SM/Chol vesicles. A reduction in inner-leaflet fusion elucidated that CLBP fusogenicity was not in parallel to its membrane-leakage activity on oxPC/SM/Chol vesicles. The lipid domain formed by Chol and SM supported CLBP fusogenicity on oxPC/SM/Chol vesicles, while oxPC altered the interacted mode of CLBP with oxPC/SM/Chol vesicles as evidenced by Fourier transform infrared spectra analyses and colorimetric phospholipid/polydiacetylene membrane assay. Although CLBP showed similar binding affinity with PC/SM/Chol and oxPC/SM/Chol vesicles, the binding capability of CLBP with PC/SM/Chol and oxPC/SM/Chol vesicles was affected differently by NaCl. This emphasized that CLBP adopted different membrane interaction modes upon binding with PC/SM/Chol and oxPC/SM/Chol vesicles. CLBP induced fusion in vesicles containing oxPC bearing the aldehyde group, and aldehyde scavenger methoxyamine abrogated the CLBP ability to induce oxPC/SM/Chol fusion. Taken together, our data indicate that Chol and oxPC bearing aldehyde group alter the CLBP membrane-binding mode, leading to fusogenicity promotion while reducing the membrane-damaging activity of CLBP.

pH replica-exchange method based on discrete protonation states

We propose a new algorithm for obtaining proton titration curves of ionizable residues. The algorithm is a pH replica-exchange method (PHREM), which is based on the constant pH algorithm of Mongan et al. (J Comput Chem 2004;25:2038-2048). In the original replica-exchange method, simulations of different replicas are performed at different temperatures, and the temperatures are exchanged between the replicas. In our PHREM, simulations of different replicas are performed at different pH values, and the pHs are exchanged between the replicas. The PHREM was applied to a blocked amino acid and to two protein systems (snake cardiotoxin and turkey ovomucoid third domain), in conjunction with a generalized Born implicit solvent. The performance and accuracy of this algorithm and the original constant pH method (PHMD) were compared. For a single set of simulations at different pHs, the use of PHREM yields more accurate Hill coefficients of titratable residues. By performing multiple sets of constant pH simulations started with different initial states, the accuracy of predicted pK(a) values and Hill coefficients obtained with PHREM and PHMD methods becomes comparable. However, the PHREM algorithm exhibits better samplings of the protonation states of titratable residues and less scatter of the titration points and thus better precision of measured pK(a) values and Hill coefficients. In addition, PHREM exhibits faster convergence of individual simulations than the original constant pH algorithm.

A fast and accurate method for predicting pKa of residues in proteins

Predicting the pH-activities of residues in proteins is an important problem in enzyme engineering and protein design. A novel predictor called 'Pred-pK(a)' was developed based on the physicochemical properties of amino acids and protein 3D structure. The Pred-pK(a) approach considers the influence of all other residues of the protein to predict the pK(a) value of an ionizable residue. An empirical equation was formulated, in which the pK(a) value was a distance-dependent function of physicochemical parameters of 20 amino acid types, describing their electrostatic and van der Waals interaction, as well as the effects of hydrogen bonds and solvation. Two sets of coefficients, {a(alpha)} and {b(l)}, were used in the predictor: {a(alpha)} is the weight factors of 20 amino acid types and {b(l)} is the weight factors of physicochemical properties of amino acids. An iterative double least square procedure was proposed to solve the two sets of weight factors alternately and iteratively in a training set. The two coefficient sets {a(alpha)} and {b(l)} thus obtained were used to predict the pK(a) values of residues in a protein. The average predictive error is +/-0.6 pH in less than a minute in common personal computer.

Electrostatic contributions to the stabilities of native proteins and amyloid complexes

The ability to predict electrostatic contributions to protein stability from structure has been a long-standing goal of experimentalists and theorists. With recent advances in NMR spectroscopy, it is possible to determine pK(a) values of all ionizable residues for at least small proteins, and to use the pK(a) shift between the folded and unfolded states to calculate the thermodynamic contribution from a change in charge to the change in free energy of unfolding. Results for globular proteins and for α-helical coiled coils show that electrostatic contributions to stability are typically small on an individual basis, particularly for surface-exposed residues. We discuss why NMR often suggests smaller electrostatic contributions to stability than X-ray crystallography or site-directed mutagenesis, and discuss the type of information needed to improve structure-based modeling of electrostatic forces. Large pK(a) shifts from random coil values are observed for proteins bound to negatively charged sodium dodecyl sulfate micelles. The results suggest that electrostatic interactions between proteins and charges on the surfaces of membrane lipid bilayers could be a major driving force in stabilizing the structures of peripheral membrane proteins. Finally, we discuss how changes in ionization states affect amyloid-β fibril formation and suggest that electrostatic repulsion may be a common destabilizing force in amyloid fibrils.

Benchmarking pK(a) prediction

Background

pK_a values are a measure of the protonation of ionizable groups in proteins. Ionizable groups are involved in intra-protein, protein-solvent and protein-ligand interactions as well as solubility, protein folding and catalytic activity. The pK_a shift of a group from its intrinsic value is determined by the perturbation of the residue by the environment and can be calculated from three-dimensional structural data.

Results

Here we use a large dataset of experimentally-determined pK_as to analyse the performance of different prediction techniques. Our work provides a benchmark of available software implementations: MCCE, MEAD, PROPKA and UHBD. Combinatorial and regression analysis is also used in an attempt to find a consensus approach towards pK_a prediction. The tendency of individual programs to over- or underpredict the pK_a value is related to the underlying methodology of the individual programs.

Conclusion

Overall, PROPKA is more accurate than the other three programs. Key to developing accurate predictive software will be a complete sampling of conformations accessible to protein structures.

Chemically accurate protein structures: validation of protein NMR structures by comparison of measured and predicted pKa values

A new method is presented for evaluating the quality of protein structures obtained by NMR. This method exploits the dependence between measurable chemical properties of a protein, namely pKa values of acidic residues, and protein structure. The accurate and fast empirical computational method employed by the PROPKA program ( http://www.propka.chem.uiowa.edu) allows the user to test the ability of a given structure to reproduce known pKa values, which in turn can be used as a criterion for the selection of more accurate structures. We demonstrate the feasibility of this novel idea for a series of proteins for which both NMR and X-ray structures, as well as pKa values of all ionizable residues, have been determined. For the 17 NMR ensembles used in this study, this criterion is shown effective in the elimination of a large number of NMR structure ensemble members.

Cobra cardiotoxin and phospholipase A2 as GAG-binding toxins: on the path from structure to cardiotoxicity and inflammation

Glycosaminoglycans (GAGs) represent the sulfated carbohydrate moieties of proteoglycans which occur abundantly in tissues of the cardiovascular system. Many proteins bind specifically to GAGs and perform an important role in inflammation, cell proliferation, and blood coagulation processes. Recently, in vitro GAG-binding studies of cardiotoxins (CTXs) and basic phospholipase A(2) (PLA(2)) from cobra venom established the toxins as two new families of GAG-binding proteins. In particular, discontinuous basic residues in beta-sheet CTXs may form a cationic cradle suitable for heparin binding, as in the case of fibronectin module III-13. The binding specificity of beta-sheet proteins to different GAGs can be further enhanced by involving other cationic clusters near the flexible loop of the molecule. Since the three-dimensional structures of many CTXs and PLA(2) are available, these two toxins may serve as models for the elucidation of the molecular recognition of GAG-binding proteins and also as polypeptide templates for further improvement of the binding specificity suitable for future biomedical application. Research along the line of GAG-guided toxicity of cobra venom components may help us to understand the functional role of GAGs and the action mechanism of cobra venom components in the cardiovascular system.

Structural basis of membrane-induced cardiotoxin A3 oligomerization

Cobra cardiotoxins (CTXs) have previously been shown to induce membrane fusion of vesicles formed by phospholipids such as cardiolipin or sphingomyelin. CTX can also form a pore in membrane bilayers containing a anionic lipid such as phosphatidylserine or phosphatidylglycerol. Herein, we show that the interaction of CTX with negatively charged lipids causes CTX dimerization, an important intermediate for the eventual oligomerization of CTX during the CTX-induced fusion and pore formation process. The structural basis of the lipid-induced oligomerization of CTX A3, a major CTX from Naja atra, is then illustrated by the crystal structure of CTX A3 in complex with SDS; SDS likely mimics anionic lipids of the membrane under micelle conditions at 1.9-A resolution. The crystal packing reveals distinct SDS-free and SDS-rich regions; in the latter two types of interconnecting CTX A3 dimers, D1 and D2, and several SDS molecules can be identified to stabilize D1 and D2 by simultaneously interacting with residues at each dimer interface. When the three CTXSDS complexes in the asymmetric unit are overlaid, the orientation of CTX A3 monomers relative to the SDS molecules in the crystal is strikingly similar to that of the toxin with respect to model membranes as determined by NMR and Fourier transform infrared methods. These results not only illustrate how lipid-induced CTX dimer formation may be transformed into oligomers either as inverted micelles of fusion intermediates or as membrane pore of anionic lipid bilayers but also underscore a potential role for SDS in x-ray diffraction study of protein-membrane interactions in the future.

Empirical relationships between protein structure and carboxyl pKa values in proteins

Relationships between protein structure and ionization of carboxyl groups were investigated in 24 proteins of known structure and for which 115 aspartate and 97 glutamate pK(a) values are known. Mean pK(a) values for aspartates and glutamates are < or = 3.4 (+/-1.0) and 4.1 (+/-0.8), respectively. For aspartates, mean pK(a) values are 3.9 (+/-1.0) and 3.1 (+/-0.9) in acidic (pI < 5) and basic (pI > 8) proteins, respectively, while mean pK(a) values for glutamates are approximately 4.2 for acidic and basic proteins. Burial of carboxyl groups leads to dispersion in pK(a) values: pK(a) values for solvent-exposed groups show narrow distributions while values for buried groups range from < 2 to 6.7. Calculated electrostatic potentials at the carboxyl groups show modest correlations with experimental pK(a) values and these correlations are not improved by including simple surface-area-based terms to account for the effects of desolvation. Mean aspartate pK(a) values decrease with increasing numbers of hydrogen bonds but this is not observed at glutamates. Only 10 pK(a) values are > 5.5 and most are found in active sites or ligand-binding sites. These carboxyl groups are buried and usually accept no more than one hydrogen bond. Aspartates and glutamates at the N-termini of helices have mean pK(a) values of 2.8 (+/-0.5) and 3.4 (+/-0.6), respectively, about 0.6 units less than the overall mean values.

Effect of D57N mutation on membrane activity and molecular unfolding of cobra cardiotoxin

Cobra cardiotoxins (CTXs) are able to adopt a three-fingered beta-strand structure with continuous hydrophobic patch that is capable of interacting with zwitterionic phospholipid bilayer. In addition to the four disulfide bonds that form the rigid core of CTXs, Asp57 near the C-terminus interacts electrostatically with Lys2 near the N-terminus (Chiang et al. 1996. Biochemistry. 35:9177-9186). We indicate herein, using circular dichroism and the time-resolved polarized tryptophan fluorescence measurement, that Asp57 to Asn57 (D57N) mutation perturbs the structure of CTX molecules at neutral pH. The structural stability of the D57N mutant was found to be lower, as evidenced by the reduced effective concentration of the 2,2,2-trifluoethanol (TFE)-induced beta-sheet to alpha-helix transition. Interestingly, the single mutation also allows a greater degree of molecular unfolding, because the rotational correlation time of the TFE-induced unfolding intermediate is larger for the D57N mutant. It is suggested that the electrostatic interaction between N- and C-termini also contributes to the formation of the functionally important continuous hydrophobic stretch on the distant end of CTX molecules, because both the binding to anilinonaphthalene fluorescent probe and the interaction with phospholipid bilayer were also reduced for D57N mutant. The result emphasizes the importance of the hydrophobic amino acid residues near the tip of loop 3 as a continuous part of the three-fingered beta-strand CTX molecule and indicates how a distant electrostatic interaction might be involved. It is also implicated that electrostatic interaction plays a role in expanding the radius of gyration of the folding/unfolding intermediate of proteins.

Snake venom cardiotoxins-structure, dynamics, function and folding

Snake cardiotoxins are highly basic (pI > 10) small molecular weight (approximately 6.5 kDa), all beta-sheet proteins. They exhibit a broad spectrum of interesting biological activities. The secondary structural elements in these toxins include antiparallel double and triple stranded beta-sheets. The three dimensional structures of these toxins reveal an unique asymmetric distribution of the hydrophobic and hydrophilic amino acids. The 3D structures of closely related snake venom toxins such as neurotoxins and cardiotoxin-like basic proteins (CLBP) fail to show similar pattern(s) in the distribution of polar and nonpolar residues. Recently, many novel biological activities have been reported for cardiotoxins. However, to-date, there is no clear structure-function correlation(s) available for snake venom cardiotoxins. The aim of this comprehensive review is to summarize and critically evaluate the progress in research on the structure, dynamics, function and folding aspects of snake venom cardiotoxins.

Analysis of binding of cobra cardiotoxins to heparin reveals a new beta-sheet heparin-binding structural motif

Heparin and heparan sulfate have recently been shown to bind to snake cardiotoxin (CTX) and to potentiate its penetration into phospholipid monolayer under physiological ionic conditions. Herein we analyze the heparin-binding domain of CTX using 10 CTXs from Taiwan and African cobra venom. We also performed computer modeling to obtain more information of the binding at molecular level. The results provide a molecular model for interaction of CTX-heparin complex where the cationic belt of the conserved residues on the concave surface of three finger beta-sheet polypeptides initiates ionic interaction with heparin-like molecules followed by specific binding of Lys residues near the tip of loop 2 of CTX. The dissociation constants of CTXs differ by as much as 4 orders of magnitude, ranging from approximately 140 microM for toxin gamma to approximately 20 nM for CTX M3, depending on the presence of Lys residues near the tip of loop 2. High affinity heparin binding becomes possible due to the presence of Arg-28, Lys-33, or the so-called consensus heparin binding sequence of XKKXXXKRX near the tip of the loop. The well defined three-finger loop structure of CTX provides an interesting template for the design of high affinity heparin-binding polypeptides with beta-sheet structure. The finding that several cobra CTXs and phospholipase A2 bind to heparin with different affinity may provide information on the synergistic action of the two venom proteins.

Heparin and heparan sulfate bind to snake cardiotoxin. Sulfated oligosaccharides as a potential target for cardiotoxin action

Cardiotoxins (CTXs) from cobra venom show cytotoxicity toward several cell types. They cause systolic heart arrest and severe tissue necrosis. Their interaction with phospholipids is established but by itself fails to explain the specificity of these toxins; other component(s) of membrane must, therefore, intervene to direct them toward their target. We herein show, for the first time, that sulfated glycosaminoglycans, heparin, heparan sulfate (HS), chondroitin sulfate (CS), and dermatan sulfate (DS), interact with CTX A3, a major component of Taiwan cobra venom, by use of affinity chromatography, circular dichroism, absorbance, and fluorescence intensity and anisotropy measurements. The relative strength of binding, determined by the NaCl concentration required to dissociate the CTX-glycosaminoglycan complex, varied as follows: heparin > DS > CS > HS. In physiological buffer (8 mM Na2HPO4, 2.7 mM KCl, 1.8 mM KH2PO4, 138 mM NaCl, pH 7.4), however, only heparin and HS bound to CTX, with respective dissociation constants of 1.4 and 16 microM, while CS and DS failed to exhibit well defined binding behavior, as indicated by fluorescence measurements. We estimate that CTX makes 3-4 ionic contacts with heparin based on a salt-dependent binding study and that approximately 40% of binding free energy is derived from purely electrostatic interactions under physiological conditions. Sulfated pentasaccharide may be sufficient to bind to CTX. We also found that heparin accentuates the penetration of CTX into phospholipid membranes as analyzed by Langmuir monolayer measurement. In view of these results we propose that heparin-like moieties of the cell surface may modulate the action of CTX.