PLA2-like proteins myotoxic mechanism: a dynamic model description

Phospholipase A2-like (PLA2-like) proteins contribute to the development of muscle necrosis in Viperidae snake bites and are not efficiently neutralized by current antivenom treatments. The toxic mechanisms of PLA2-like proteins are devoid of catalytic activity and not yet fully understood, although structural and functional experiments suggest a dimeric assembly and that the C-terminal residues are essential to myotoxicity. Herein, we characterized the functional mechanism of bothropic PLA2-like structures related to global and local measurements using the available models in the Protein Data Bank and normal mode molecular dynamics (NM-MD). Those measurements include: (i) new geometric descriptions between their monomers, based on Euler angles; (ii) characterizations of canonical and non-canonical conformations of the C-terminal residues; (iii) accessibility of the hydrophobic channel; (iv) inspection of ligands; and (v) distance of clustered residues to toxin interface of interaction. Thus, we described the allosteric activation of PLA2-like proteins and hypothesized that the natural movement between monomers, calculated from NM-MD, is related to their membrane disruption mechanism, which is important for future studies of the inhibition process. These methods and strategies can be applied to other proteins to help understand their mechanisms of action.


Results and Discussion
Tertiary structure variability: local measurement using the Cα and Cβ distances. PLA 2 -like proteins are small (molecular weight of approximately 14 kDa), with seven disulfide bridges and the following secondary structural elements: a N-terminal helix, putative calcium binding loop residues, two long antiparallel α-helices at the core of the protein that make up the hydrophobic channels, two antiparallel β-sheets known as β wings, one short 3 10 helix and C-terminal loop region. Among the PLA 2 -like models, two specific regions have greater variation, with root-mean-square fluctuations (RMSF) of Cα values (calculation description is provided in Material and Methods -geometrical description between identical monomers) above 1.5 Å (black line in Fig. 1). Flexible regions are usually related to functions, as dynamic structural alterations are part of substrate binding and enzymatic actions or interactions 29 . First, the labile region at residues 86 and 87 (with 2.2 and 1.6 Å variations, respectively) was previously reported as a region that oscillates between two metastable conformations 25 . A second region, which contains the most variable residues, is in the C-terminal region of the toxin (black line in Fig. 1), where two residues (121 and 125) were suggested to be responsible for toxin-induced membrane disruption and named MDiS 14 .
MDiS is composed of residues 121 and 125, which are hydrophobic and exposed to the solvent in a short 3 10 helix, together with residues 122 and 123 24 . This secondary structure is stabilized by one hydrogen bond between the main chain atoms of residues 121 and 125 24 . K122 is the only conserved C-terminal residue among PLA 2 -like proteins, although the hydrophobic character of residues 121 and 125 is also conserved 14 . MDiS is characterized by a Cβ distance between residues 121 and 125 (dMDiS) below 5.5 Å. In active state structures, both monomers are similar with this canonical conformation of MDiS (Tables 1 and 2). In inactive state structures, monomers are asymmetrical, with one in a non-canonical conformation and a dMDiS higher than 9 Å, as these residues are in an extended loop structure (Tables 1 and 2) 24 . Inact is the only inactive state structure with a non-canonical conformation, with a dMDiS of 8.5 Å; MDiS is in an intermediate state 24 . Characterization and accessibility of the hydrophobic channel of PLA 2 -like toxins. Bothropic PLA 2 -like proteins have a cavity that is surrounded by positive and hydrophobic residues that participate in the binding of fatty acids 1 . The entrance of fatty acids into the hydrophobic channel of the toxin is hypothesized to be essential to induce dimer reorientation from the inactive to active state, a key step of the myotoxic mechanism of Figure 1. Graph of the RMSF in Å by residues of available bothropic PLA 2 -like crystallographic structures. The RMSF values of Cα and centroid of 7 consecutive Cα are coloured in black and grey, respectively. Secondary structure is shown below residue numbers with helices and sheets represented by cylinders and arrows, respectively (extracted from PDB). action 14,22 . To better characterize the allosteric activation induced by the entrance of hydrophobic molecules, we evaluated the accessibility of the hydrophobic channel through the iFace of the protein using CAVER software.
All canonical monomers contain a hydrophobic channel that is accessible through C-terminal residues. The PLA 2 -like structures in the inactive state only contain the accessible hydrophobic channel in monomers that are in the canonical conformation with a volume ~266 Å 3 (Table 2). Curiously, the hydrophobic channels of the canonical monomers of BthTX-I/PEG400 and apo1 BthTX-I are closed by the side chain of H120, whereas they are open in the available models of other bothropic PLA 2 -like proteins 24 .
PLA 2 -like structures in the active state have symmetrical monomers, and the calculated tunnel is occupied by ligands, such as fatty acids, PEG molecules or inhibitors (Supplementary Table S1). Most of the calculated hydrophobic channels have a volume of ~400 Å 3 , with the exception of PrTX-II/n-tridecanoic acid, whose volumes are ~223 Å 3 ( Table 2). MjTX-II has an insertion at residue 120 compared to the other bothropic toxins (BaspTX-II, BbTX-II, BnIV, BnSP-VII, BthTX-I, PrTX-I and PrTX-II); thus, instead of the hydrophobic channel exiting through the middle region of the N-terminal helix, it is shifted to the beginning of the N-terminal helix 15 . Despite these differences, all of the calculated hydrophobic channels structurally connect the C-terminal region of one monomer to the N-terminal region and interior of the two parallel α-helices of the other monomer. The N-terminal region with its hydrophobic residues may participate in toxin activity 30,31 and may be related to membrane interactions, similar to the MDiS.
Global measurements: orientation and translations between almost identical objects. The orientation between the monomers of available PLA 2 -like models are remarkably different and are related to the symmetrical and asymmetrical hydrophobic channel accessibility of active and inactive states structures, respectively. Previous efforts attempted to characterize the different orientations of models by one 25 or two angles descriptions 22 . Herein, we update this angle description with a more intuitive, general and informative methodology. Since these structures are homodimers, we characterize the monomers orientation by describing the necessary rotation in Euler angles to superpose one monomer onto the other.
Han et al. used a similar approach to study the interdomain geometry between two domains of homologous proteins 32 . They described domain motion with a single translation and single rotation. Hayward wrote the DYNDOM software (fizz.cmp.uea.ac.uk/dyndom) to determine the hinge axes and measure closure or twist motions in a single angle description by using two conformations of the same protein 33  available crystallographic structures, our approach differs, as we distinguish PLA 2 -like models using replicable and independent metrics based on Euler angles (Fig. 2). In our angle description and inspired by the DNA angular base pair parameters 34 , we called the angles tilt, twist and roll for the X, Y and Z axes, respectively ( Fig. 2b-f and the calculation description is provided in Material and Methods -geometrical description between identical monomers). The X-axis is the centre of the longest α-helix (H3), the Y-axis is the vector that connects the X-axis to the other antiparallel α-helix (H2), and the Z-axis (orange axis in Fig. 2b) is the normal vector of the plane XY. Among the evaluated dimers, twist (Ts) is the angle with the lowest divergence and a value close to 25° in all samples (Tables 1 and 2). The roll angles (R) separate inactive and active states. The former is composed of the apo BnSP-VII, apo PrTX-I, apo1&2 BthTX-I, apo BbTX-II, BthTX-I/Zn 2+ and BthTX-I/PEG400 (inact) structures, which have a roll of ~140° and tilt (Ti) of ~55° (Tables 1 and 2). All of these models are apo structures (Supplementary Table S1), with exception of BthTX-I/Zn 2+ and apo PrTX-I, which have charged molecules or ions bound. Moreover, inact has a small PEG molecule inside the hydrophobic channel that induces dMDiS shortening (8.5 Å) 24 .
The structures in the active state have a roll angle greater than 160° (Tables 1 and 2). They are complexed with inhibitor molecule(s) or natural compound(s) in the hydrophobic channel (Supplementary Table S1), which agrees with the accessibility calculation of the previous sub item. Most of these structures have abundant negative molecules surrounding the protein, being the R34 and PLA 2 -like iFace [considered the side chains of residues 16, 17, 20, 115 and 118 (MDoS)] the most common region. This environment, which is rich in hydrophobic and negative molecules, may mimic a membrane interaction 23 . Roll angles of the active state close to 180° are closer to a symmetrical orientation, which may be the optimal exposure of same residues of both monomers in the protein iFace.
Evaluation of the flexibility of PLA 2 -like protein in active and inactive states. Analyzing the local and global measurements of available bothropic PLA 2 -like proteins, we identified two different states based on their roll angles: i) an inactive state (roll ~ 140°) that has an asymmetrical conformation with one non-canonical monomer and one hydrophobic channel accessible and ii) an active state (roll > 160°) with symmetrical canonical monomers and with both hydrophobic channels filled by hydrophobic molecules. It was hypothesized that for these toxins, the first state is converted to the second by an allosteric activation induced by the entrance of a hydrophobic molecule into the accessible hydrophobic channel 24 .
To address the allosteric activation hypothesis in silico, we evaluated the flexibility via NM-MD analysis of the BthTX-I models, as both inactive and active states were structurally determined. We chose inact and acti1 as representatives of the former and latter groups, respectively. In both cases, we evaluated flexibility in absence and presence of fatty acids, whereas PEG molecules available in the model were replaced by myristic acid. PLA 2 -like activation. All of the models from inact and acti1 calculated from a ± 6.0 Å displacement of NM 7-9 had satisfactory energies, as the measured potential energies were smaller than 0.5 kcal/mol ( Supplementary  Fig. S2) and within order of magnitude of evaluated structures of other studies 35,36 . To observe the transition from one state to another, we compared the inact and acti1 models by superposing them on the corresponding NM-MD structures (Supplementary Table S3). Inact-FA/NM9/-3.1 Å (inact structure complexed to fatty acid minimized and displaced -3.1 Å among NM9) is the most similar to acti1, with a RMSD of 1.9 Å, and as a reference, the model itself prior to minimization has a minimum RMSD of 1.3 Å (Fig. 3a) Table S3).
To confirm that lower amplitude NMs of apo-inact would not reach a better representation of acti1 model, the NM10-21 were also calculated, but with no better similarity than 2.8 Å. Therefore, the fatty acid in the canonical monomer of BthTX-I was fundamental for the toxin reach its active state (described by comparison with acti1) since the lowest RMSD obtained from inact-FA was almost 1 Å higher than apo-inact. Moreover, inact-FA/NM9 from −0.4 to −3.1 Å are the calculated structures that better describes the BthTX-I activation; such conformation changes will be further described as activation (Fig. 3a, Supplementary Fig. S4 and Supplementary Video S5). For acti1-FA, the minimized structure on NM8/−2.7 Å is the most similar to inact, with a RMSD of 1.5 Å, and as a reference, the model itself prior to minimization has a minimum RMSD of 1.3 Å (Fig. 3b). The other calculated NMs of acti1-FA have a minimum RMSD, compared to inact, above 2.3 Å (Supplementary Table S3). For apo-acti1, the structure NM8/−2.1 Å is the most similar to inact with a RMSD of 1.7 Å, and as a reference, the model itself prior to minimization had a minimum RMSD of 1.8 Å (Supplementary Table S3). The presence of a fatty acid in both of canonical monomers of BthTX-I seems not essential for toxin acquiring its inactive state, as both apo-acti1 and acti1-FA NM8 reached similar representations of inact. Therefore, acti1-FA/NM8 from −2.7 to 0 Å are the calculated structures that better describes BthTX-I inactivation; such conformation changes will be further described as inactivation (Fig. 3b and Supplementary Video S6).
To observe whether the activation and inactivation processes were complementary, we compared them by superposing inact-FA/NM9/±6 Å and acti1-FA/NM8/±6 Å and then plotting the RMSD in a heat map (Fig. 3c),  (in a). N-and C-terminal residues compose this iFace, which includes a membrane docking site represented by green sticks in (a) that interacts with sulphates. In (b) the PLA 2 -like monomer is represented by a grey cartoon with Cc7 represented as spheres. The axes are established by Cc7/98 as the origin and coloured black; by Cc7/108 as the X-axis and coloured magenta; by a perpendicular vector to the X-axis, in the plane composed by X-axis and Cc7/48 sphere, as the Y-axis and coloured dark blue; and the normal vector to the XY plane coloured vermilion. In (c-f) the two antiparallel helices are illustrated to indicate different orientations of monomer A (light blue) corresponding to monomer B; the calculated angles are shown in the box. Tilt is coloured in vermilion, twist in dark blue and roll in dark blue. In (c-f) monomer B, coloured in black, is close to the orientation of BnIV/myristic acid structure (PDB id: 3MLM), which is established as a reference to show a difference of 45° in each of the axes of movement. The axes in (d-f) are shown as illustrations of the rotation that describes monomer B in transparent black relative to the other monomer, which is not transparent.
with darker colours indicating similar structures. These two processes demonstrate that there is better agreement between inact-FA/NM9/−3.7 Å and acti1-FA/NM8/0 Å to inact-FA/NM9/0 Å and acti1-FA/NM8/−3.7 Å that are 38 consecutive superpositions with an average RMSD of 2.0 Å (shown by the dashed white transparent line in Fig. 3c). The angles from these paired structures correspond, as shown in Fig. 3d, which demonstrates the activation process as a movement that increases the roll angle and decreases the tilt angle. Therefore, our results suggest that the inactive state of PLA 2 -like proteins may transition to the active state through physically plausible movements in the presence of a fatty acid in the hydrophobic channel. PLA 2 -like flexibility among active state models. PLA 2 -like structures may change to an active state from an inactive state after fatty acid interaction, which is described by an increase in the roll angles. Active state models may be further separated by tilt angles, with at least one model with both hydrophobic channels occupied by a fatty acid: acti1 (BthTX-I), with ~30°; acti2 (MjTX-II), with ~50°; and acti3 (PrTX-II), with ~80° (Table 2). Are these angles exclusive to a single toxin or are they a consequence of crystal packing? To address this question in silico, we evaluated their flexibility with NM-MD. Only the fatty acids in the toxin hydrophobic channels were maintained in the model.
All of the models from acti1, acti2, and acti3 containing fatty acids and calculated from a ± 6.0 Å displacement of NM 7-9 had satisfactory energies, as they had potential energies of less than 0.5 kcal/mol ( Supplementary  Fig. S2) and within order of magnitude of evaluated structures of other studies 35,36 . Among the calculated NM, NM7 of acti1-3 were the movements that best described the conversion from one conformation to the other, as observed by RMSD comparisons (Fig. 4a-c). For acti1/NM7, the minimum RMSD for the acti2 and acti3 models were 1.6 and 2.0 Å, respectively, which corresponded to the structures acti1/NM7/0.6 Å and acti1/ NM7/5.4 Å, respectively; as a reference, the model itself prior to minimization had a minimum RMSD of 1.4 Å (Supplementary Table S3). For acti2/NM7, the minimum RMSD for acti1 and acti3 were 2.0 and 2.0 Å, Figure 3. Activation and inactivation descriptions using RMSD and proposed angles. In (a) inact/NM9/0 Å and -3Å structures resemble inact (orange line) and acti1 (blue line), respectively, describing the BthTX-I activation process. In (b) acti1/NM8/−3Å and 0 Å resemble inact (orange line) and acti1 (blue line), respectively, describing the BthTX-I inactivation process. In (c) the RMSD heat map demonstrates darker and lighter colours corresponding to high and low similarities, respectively. The structures from the dashed white transparent line from the white circle, which starts in act/NM8/−3.7 Å and inact/NM9/0 Å and ends in acti1/NM8/0 Å and inact/ NM9/−3.7 Å, are similar and describe complementarity between activation and inactivation. In (d) the pairs of structures delimitated by the dashed white transparent line in (c) have similar angles, and activation is characterized by increase of roll, stability of twist and decrease of tilt angles. In (d) roll, twist, and tilt are coloured magenta, dark blue and vermilion, respectively, and models of acti1 and inact are represented by X and ◆, respectively. respectively, which corresponded to the structures acti2/NM7/3.0 Å and acti2/NM7/−3.8 Å, respectively; as a reference, the model itself prior to minimization had a minimum RMSD of 1.5 Å. For acti3/NM7, the minimum RMSD for acti1 and acti2 were 3.3 Å and 3.1 Å, respectively, which both corresponded to the structure acti3/NM7/−4.5 Å; as a reference, the model itself prior to minimization had a minimum RMSD of 1.0 Å. The other structures from calculated NMs 8 and 9 of acti1, acti2, and acti3 had minimum RMSD values compared to the independent models higher than the obtained by NM7 (Supplementary Table S3). Therefore, NM7 better describes the movements of conversion between different active state models (Fig. 5a,b and Supplementary Video S7).
The NM7 of these structures were correlated, as observed by the white dashed lines with a circle on the edges in the RMSD heat maps in Fig. 4d-f. These white dashed lines are the lowest 60 consecutive superpositions: 6 Å of highly correlated pairs of structures. The average RMSD for the superposition of these white dashed lines was 2.2, 2.6, and 2.8 Å for acti1/acti2, acti1/acti3, and acti2/acti3, respectively. Each NM7 had a different shift with respect to the other NM7, which was 2.1, 5.0, and 3.0 Å for acti1/acti2 (Fig. 4d), acti1/acti3 (Fig. 4e), and acti2/acti3 (Fig. 4f), respectively. The orientations between these paired structures, which were obtained by applying these shifts, corresponded (Fig. 4g and Supplementary Video S7). Therefore, our results suggest that the active state of PLA 2 -like structures with different tilt angles may be interchangeable by the described NM7s, which suggests that their mechanism of action is probably similar. PLA 2 -like membrane disruption description. As the X-axis (measured by roll) is perpendicular to the plane of the protein iFace, the movements described by the other two orthogonal rotations, such as tilt, will change the orientations of the residues in the iFace. Therefore, the experimental flexibility observed in the tilt of the active state models and calculated NM7 models (Fig. 5) could be related to the disruption of the membrane in a clawing movement, as the C-terminal residues (MDiS) would move toward the plane of the iFace. Such a movement could pull phospholipids out of the membrane, disrupting its integrity, as hypothesized by Da Silva-Giotto et al. 25 . We calculated, along the direction of the tilt angle, a reduction in NM7 of active state PLA 2 -like proteins with fatty acids, tilt angles and MDiS approximations to the iFace ( Fig. 5c and Supplementary Video S8). This approximation was evaluated by calculating the distance of the centroids of 121 Cβ and 125 Cβ to the plane of sulphates (dMDiSiFace). We calculated this distance using the plane equation with three sulphur atoms from the plane of sulphates in BnIV/myristic acid (shown in Fig. 2) after the dimers were superposed. The reduction of the tilt angles was correlated with dMDoSiFace and dMDiSiFace reductions (Fig. 5c). Therefore, the tilt variation among the crystal and NM-MD structures supports the hypothesis that reduced tilt angles correlate the MDoS and MDiS approximations to the toxin iFace.
Action mechanism. Based on these global and local analyses, we can better describe the various states of bothropic PLA 2 -like protein structures (Fig. 6). The inactive state of bothropic PLA 2 -like protein has monomers in an asymmetrical conformation. One monomer (A) is in a non-canonical conformation that does not have a MDiS cluster (dMDiS > 9 Å), and the hydrophobic channel is inaccessible (as observed by grey triangles of the orange structures in the dMDiS and the hydrophobic channel accessibility graphs shown in Fig. 6). The geometric relationship of the monomers is far from a symmetrical orientation, as the roll is distant to 180° (observed from angles of the orange structures in the monomer-monomer angle graph shown in Fig. 6). Both monomers of these structures lack hydrophobic molecules inside the hydrophobic channel.
Due to the entrance of a hydrophobic molecule in the accessible toxin channel, the monomers are reoriented, with an increased roll angle towards 180°, to a more symmetrical relationship, leading to toxin activation. Active state models feature both monomers as being symmetrical with MDiS (dMDiS < 6 Å), and both hydrophobic channels are accessible (observed on the light blue structures in the MDiS distance and hydrophobic channel accessibility graphs in Fig. 6) and occupied by hydrophobic molecules. Despite these similarities, they differ according to their tilt angles in three primary categories, which are interchangeable according to NM calculations. Our data support that these differences are not related to a specific protein but are exclusive to a minimum energy level reached in the crystal packing, and different PLA 2 -like toxins may achieve the conformations of one another in solution. The reduction of the tilt angles precludes the MDiS approximation to the iFace, supporting this clawing movement as a membrane disruption mechanism (Supplementary Video S8).
In conclusion, we observed new structural features of PLA 2 -like toxins by measuring the MDiS distance, hydrophobic channel accessibility and geometric orientation between monomers from the available crystal and NM-MD structures. Our proposed geometrical description, with Euler angles, has proven to be useful for generating a detailed description of the dimeric arrangement of different PLA 2 -like toxin models, as we identified an inactive state and three active states. We also used this general description to interpret interchain movements according to MD studies after validating the model pairing by RMSD calculations. Moreover, this methodology can be used as a primary tool to characterize geometrically homodimeric protein models by quickly separating structures into groups and using the angles to intuitively interpret global flexibility, describing the mechanism of action.  15 , and MjTX-II/suramin+PEG4k (4YV5) 44 . The structures BthTX-I/BPB (3I03), apo BthTX-I (3I3I), apo BnSP-VI (1PC9), apo BaspTX-II (1CLP), and BbTX-II (4DCF) were not included because in the first two, their asymmetric unit (ASU) content is composed of a monomer, and in the others, the resolution was lower than or equal to 2.5 Å.

Methods
Bioinformatic tools. The software CAVER 3.0 (command line version) 45 was used to characterize the toxin hydrophobic channel of the PLA 2 -like crystallographic models. Hydrophobic channels are referred to as tunnels by this algorithm. To improve the convergence of the different tunnels, the probe and shell radii were increased to 1.8 and 4, respectively. The calculated tunnels were manually inspected in PyMOL (The PyMOL Molecular Graphics System, Version 1.7 Schrödinger, LLC.), and their volumes were calculated using the diameter and length of the tunnel. When PEGs or fatty acids were present in the model, the tunnel whose path filled this molecule was chosen. PyMOL and VMD (http://www.ks.uiuc.edu/Research/vmd/) 46 was used to create cartoon and stick images.
Geometrical description between identical monomers. The geometric orientation between the monomers was measured by superposing one on the other and converting the matrix rotation into Euler angles . Correspondence between NM7 of acti1-3 shown with the RMSD and proposed angles. In A, acti1/ NM7/0 Å, 0.6 Å, and 5.4 Å structures resemble acti1, acti2, and acti3 crystallographic models, respectively. In (b) acti2/NM7/−3.8 Å, 0 Å, and 3.0 Å structures resemble acti1, acti2, and acti3 models, respectively. In (c) acti3/NM7/−4.5 Å, −4.5 Å, and 0 Å structures resemble acti1, acti2, and acti3 models, respectively. In (a-c) acti3, acti2, and acti1 are coloured in light grey, dark grey, and black, respectively. In (c-e) the RMSD heat map demonstrates darker and lighter colours as high and low similarities, respectively. The dashed white lines with circles on their edges correspond to the 60 lowest consecutive superpositions; therefore, 6 Å of highly correlated pair of structures. In (d) the NM7s of acti1 and acti2 are similar, with a shift of 2.1 Å, as shown by dashed white arrow. In (e) the NM7s of acti1 and acti3 are similar, with a shift of 5.0 Å, as shown by dashed white arrow. In (f) the NM7s of acti2 and acti3 are similar, with a shift of −3.0 Å, as shown by the dashed white arrow. In (g) the angles of the acti1-3/NM7 structures are plotted following the shift of (d-f). In (g) roll, twist, and tilt are coloured magenta, dark blue and vermilion, respectively, and the acti1, acti2, and acti3 structures are represented by X, │, and ─, respectively.
SCIENtIFIC REPORtS | 7: 15514 | DOI:10.1038/s41598-017-15614-z (following the sequence 3,2,1 using the formulas from Section 8.11 of Diebel's publication 47 ). This is the most intuitive 3D rotation description as it uses each of the 3 axes once, such as the Tait-Bryan angles (row, pitch and yaw) used in aeronautics 48 .
Prior to the superposition calculations, the models were placed into immutable and standard orthogonal axes that relate structure to function. As the PLA 2 -like protein iFace forms a plane that is mostly composed of the N-and C-terminal regions of the protein, which are perpendicular to the two antiparallel α-helices (H2 and H3) (Fig. 2a), these two helices were used as references for two of the axes. To obtain the first axis in the direction and centre of the largest α-helix, the centroid of 7 consecutive carbon α (named herein as Cc7 and shown as grey spheres in Fig. 2b) was calculated. Approximately seven residues complete two α-helix turns. To facilitate the placement of different models in a common set of orthogonal axes, three Cc7s with small fluctuations were selected from the two longest helices, Cc7/48, Cc7/98 and Cc7/104, which had RMSF of 0.13, 0.08, and 0.08 Å, respectively (grey line in Fig. 1). The RMSF of Cα (black line in Fig. 1) was calculated after monomers were separated and superposed to a template (monomer A of BthTX-I/PEG4k). The RMSF of Cc7 was obtained using previous superposed structures. Cc7/98 was chosen as the origin, and Cc7/98 to Cc7/104 was used as the X-axis (magenta axis in Fig. 2b). The Y-axis was constructed as the perpendicular vector to the X-axis inside the plane composed of the X-axis and coordinate Cc7/48 (blue axis in Fig. 2b). One of the monomers was chosen as a reference to calculate the rotation matrix that was used to superpose it onto the other monomer.
The software lsqkab 49 was used to generate the required distances to calculate the RMSF and obtain the matrix rotation of the superpositions. The Ti, Ts and R angles were obtained from the rotation matrix. For the translation calculation, the distance between the centre of gravity of each monomer was obtained using the function gravity of mass available in PyMOL.
Normal mode analysis. The LF NMs movements of the bothropic PLA 2 -like toxin containing fatty acids in the hydrophobic channel were calculated by the CHARMM v.36b1 program 50 and CHARMM36 force field 51 . The inactive state model inact and active state models acti1-3 were selected. As some of these models have PEGs inside the toxin hydrophobic channels, the ligands in acti1 and inact canonical monomers were replaced by myristic acids from the BnIV/myristic acid model after monomer-monomer Cα superposition. Since inact only has one canonical monomer, only one fatty acid was introduced in the structure. Inact and acti1 models absent of fatty acids were included in NM calculations. The topology and parameter files were generated by the CHARMM-GUI server (www.charmm-gui.org) 52 by employing an additional energy minimization. The conjugated gradient algorithm was applied with harmonic constraints that were progressively decreased from 250 to 5 kcal/mol·Å 2 , with 100 steps of minimization at each decrease. The constraints were removed to carry out an additional 10,000 conjugated gradient steps, followed by 300,000 steps of the basis Newton-Raphson algorithm. The final minimized structure was used to calculate the three lowest frequencies NMs (7-9) using the VIBRAN module of CHARMM for the structures.
The structures were generated along each NM based on a short MD simulation at a low temperature (30 K) using the VMOD facility of CHARMM followed by a minimization afterwards. A maximum displacement range of 6.0 Å was established for each direction of the NM with a 0.1 Å projection step based on the values of the mass-weighted root mean square. For each of these steps, a constant harmonic force was applied over the Cα atoms (increasing from 1,000 until 10,000 kcal/mol·Å 2 ), and a short MD simulation (1 ps) was carried out after each constant value, totalling 10 ps of simulation. The final structures were obtained by an additional 1,000 steps of conjugated gradient energy minimization maintaining the restraints. The miscellaneous mean field potential facility of CHARMM of the final structures was used for energy validation.
To relate NM to function, the distance of the essential cluster of residues to the PLA 2 -like iFace was calculated. This was accomplished by extracting the plane equation from the coordinates of 3 sulphurs of sulphates Figure 6. Distance between Cβs between residues 121 and 125, hydrophobic channel accessibility with the tunnel volume calculation, and monomer-monomer angles for all available bothropic PLA 2 -like toxins. The points refer to the structure shown in the abscissa of the third graph. In the first two graphs, monomers A and B are represented as ▴ and ⦁, respectively. In the monomer-monomer angle chart, the roll, twist, and tilt angles are coloured in magenta, vermilion and dark blue, respectively. The structures are coloured according to state, with inactive in orange and active in light blue. Toxins in an inactive state feature asymmetrical monomers, as one is non-canonical monomer with a MDiS distance greater than 8 Å and has hydrophobic channel closed (tunnel with volume 0) colored in grey in first two graphs, and its orientation relative to the other monomer is not symmetrical, as the numbers are far from 0 and 180°. Toxins in the active state possess both canonical monomers with MDiS distances close to 5 Å and both hydrophobic channels open, and the monomer-monomer orientation tends to a symmetrical relationship, as the angles tend to be either 0 or 180°. The abbreviations are related to polyethylene glycol (PEG), suramin (sur), α-tocopherol (VitE), and fatty acid (FA).
interacting with the protein iFace and calculating the distance from the chosen residues. To compare different NM motions in respect to the different states of PLA 2 -like proteins crystallographic structures, RMSD was calculated using these and the structures from NM analysis. The RMSD was calculated with lsqkab program. To illustrate NM motions, the principal component was calculated using the quasi-routine of the module VIBRAN of CHARMM to obtain the covariance matrix o Cα atomic displacements (according to Floquet et al. 53 ).