Combined molecular dynamics and continuum solvent studies of the pre-pore Cry4Aa trimer suggest its stability in solution and how it may form pore

Cry4Aa toxin is one of the highly specific mosquito-larvicidal proteins produced by the bacterium Bacillus thuringiensis subspecies israelensis. It is thought to form pores in the larval midgut membrane that cause membrane leakage and subsequent insect death. Therefore, Cry4Aa and other Cry toxins have been used as efficient and safe bacterial insecticides to control the disease-carrying mosquitoes such as Aedes, Anopheles, and Culex. However, we still do not clearly understand how Cry toxins kill mosquito-larvae at molecular details. Recent electron crystallographic images of Cry4Ba toxin, another toxin closely related to Cry4Aa toxin, have suggested that the protein forms trimer in aqueous solution and in lipid monolayer. Moreover, the unit cell of X-ray crystal structure of Cry4Ba toxin has been shown to be trimeric. In this study, we constructed the first full-atom structural model of Cry4Aa trimer using the trimeric unit cell structure of Cry4Ba toxin as a template and then used the methods of molecular dynamics (MD) and molecular mechanics combined with Poisson-Boltzmann and surface area (MM-PBSA) to show that the trimeric structure of Cry4Aa toxin is stable in 150 mM KCl solution on 10 ns timescale. The results reveal that Cry4Aa toxins use polar amino acid residues on α-helices 3, 4, and 6 to form trimer and suggest that the proteins form trimer to reduce their non-polar interactions with surrounding water. Based on the obtained trimeric structure of Cry4Aa toxins, we propose that pore formation of Cry toxins may involve a 90°-hairpin rotation during the insertion of their three α4-α5 hairpins into the membrane. This process may be mediated by water and ions. PACS Codes: 87.15.ap, 87.15.bk, 87.14.ep


Introduction
Cry4Aa is a highly specific mosquito-larvicidal protein produced by the Gram-positive bacterium Bacillus thuringiensis subspecies israelensis (Bti). It is one of several toxins produced by the bacterium in the form of crystalline inclusions during its sporulation phase [1]. These cytoplasmic inclusions consist of at least four types of mosquito-larvicidal proteins (Cry4Aa, Cry4Ba, Cry11Aa and Cyt1Aa) which can be classified into two families known as crystal (Cry) and cytolytic (Cyt) δ-endotoxins [2,3]. In particular, the Bti Cry toxins have been used as efficient and safe bacterial insecticides to control the disease-carrying mosqui-toes such as Aedes, Anopheles, and Culex [1]. However, we still do not clearly understand how Cry toxins kill mosquito-larvae at molecular details.
Pore formation is thought to be one of the molecular mechanisms for larvicidal activity of the Bti Cry toxins [3]. It can be considered as a series of important steps, including proteolytic activation, toxin oligomerization, and membrane insertion. The activation occurs when the toxin inclusions are solubilized by the highly alkaline solution in the midgut lumen. This releases soluble protoxins which are then digested by larval midgut proteases to form 65-kDa active toxins. During the oligomerizaion step, the monomeric active toxins have been suggested to bind together to form pre-pore trimers in solution [4] prior to binding to specific receptors located on the apical brush-border membrane of midgut epithelium. Subsequent conformational changes allow the insertion of their pore-forming portions into the cell membrane to form pores [5]. These pores are permeable to ions and small solutes causing membrane leakage that leads to osmotic lysis of midgut cells and finally insect death.
Structural data from X-ray crystallography of six different Bt Cry toxins, Cry1Aa [6], Cry2Aa [7], Cry3Aa [8], Cry3Ba [9], Cry4Aa [10], and Cry4Ba [11], reveal high degree of overall structural similarity. These proteins contain three distinct domains. The N-terminal domain I is a bundle of seven α-helices in which the α5 is surrounded by other helices. This domain has been shown to be responsible for pre-pore oligomerization [12][13][14], membrane insertion and pore formation [15,16]. Domain II is composed of three antiparallel β-sheets with highly varying amino acid sequence at the surface-exposed loops. This domain has been shown to participate and determine the specificity of the toxin for target insect larvae [17,18]. The C-terminal domain III contains a sandwich of two antiparallel β-sheets. This domain has been implicated in insect specificity [19][20][21] or ion regulation [22].
Recent structural studies of Cry toxins have allowed us to construct a full-atom model of pre-pore Cry4Aa trimer. Several evidence has suggested that Cry toxins may form pre-pore in solution before inserting into the membrane [4,[12][13][14]. Electron crystallographic studies have also shown that the pre-pore structure of Cry4Ba toxin, another toxin closely related to Cry4Aa toxin, is trimeric in solution [4] and in lipid monolayer [4,5]. Moreover, a unit-cell structure of Cry4Ba toxin solved by X-ray crystallography has been shown to have a three-fold symmetry [11] suggesting that the trimeric coordinate of Cry4Ba toxin could be used as a full-atom model of a Cry4Ba pre-pore trimer in MD studies. However, the X-ray structure displays a missing part comprising helices α1 and α2 of domain I [11]. One possible alternative is to use the monomeric structure of Cry4Aa toxin which has been fully solved with all parts including α1 and α2 [10]. In addition, structure and sequence alignments between Cry4Aa and Cry4Ba showed high similarity in sequences (55% similarity) [23] as well as structures (RMSD 2.5 Å). Therefore, we decided to build a full-atom model of pre-pore Cry4Aa trimer.
Is this pre-pore trimeric structure of Cry toxins stable in an aqueous solution on a molecular dynamics (MD) simulation timescale? To address this question, we set up atomic-scale MD simulations and used molecular mechanics combined with Poisson-Boltzmann and surface area calculations (MM-PBSA) to investigate the structural stability of the pre-pore Cry4Aa trimer in 150 mM KCl solution on 10 ns timescale by comparing it energetically with a 10-ns MD simulation of Cry4Aa monomer in the same solution. This allowed us to observe how amino acid residues at the inter-subunit inter-faces of the pre-pore trimer adjusted during course of the simulations as well as to study their key interactions. A conceivable notion for membrane insertion and pore formation by this toxin was also proposed. The results of this work would be paving the way for further MD simulations and mutagenic studies to gain more insights into the toxicity mechanism of the Cry toxins, particularly at the toxin oligomerization and pore formation steps.

Preparation of pre-pore Cry4Aa trimer and Cry4Aa monomer in KCl solution
The crystal structure of the 65-kDa activated Cry4Aa toxin was used to construct a structural model of pre-pore Cry4Aa trimer with its atomic coordinates taken from the Protein Data Bank (PDB) with PDB code 2C9K [10]. The structure consists of 598 amino acid residues, 252 water molecules and one methyl-2,4-pentanediol molecule which was then removed. The amino acid residues start from residues 68 to 679 with 14 missing residues between residues 235 and 244 (10 residues); and between residues 485 and 488 (4 residues). These missing residues correspond to two loops connecting α5 and α6; and β9 and β10, respectively. The two missing loops were modeled and incorporated into the structure using Loopy program [24,25] with the first loop having Arg 235 replaced by a glutamine to reproduce the mutant used in the X-ray crystallographic experiment [10].
Since there is no trimeric structure of Cry4Aa toxin available, a structural model of pre-pore Cry4Aa trimer was constructed by performing structure alignment using the incomplete trimeric structure of Cry4Ba from PDB code 1W99 as a template [11]. The Cry4Aa trimeric structure was obtained by aligning α3, α4 and α5 of the Cry4Aa monomeric structure with that of the Cry4Ba trimeric structure using DaliLite v.3 program [26]. Only α3, α4 and α5 were used for structure alignment because they are at the contact interfaces between Cry4Aa monomers.
A full-atom system of Cry4Aa trimer in a 150 mM KCl solution box was constructed using several molecular builder tools in VMD program [27]. A protein structure file (PSF) of Cry4Aa trimer was generated using Auto PSF Builder Tool with CHARMM27 force field [28]. All ionizable amino acid residues such as Arg, Lys, Asp and Glu were assigned to be in their charged states corresponding to an alkaline condition with pH 9.0. Then the Cry4Aa trimer was solvated in a water box with 12 to 15 Å buffering distance. Finally, 94 potassium (K + ) and 88 chloride (Cl -) ions were placed randomly into the system to make 150 mM KCl solution. All K + and Clions were placed at distances more than 5 Å away from the proteins. The final system of the modeled Cry4Aa trimer has 223,439 atoms consisting of three 65-kDa Cry4Aa proteins (3 × 9,771 atoms), 64,648 water molecules, 94 K + ions, and 88 Clions, and has a dimension of 150 × 150 × 105 Å 3 .
Similarly, an all-atom system of Cry4Aa monomer in a 150 mM KCl solution box was also constructed for energy comparison with the trimer system. A Cry4Aa monomer was solvated in a water box of size 92 × 102 × 78 Å 3 . Then 32 potassium and 30 chloride ions were added randomly around the protein. The final system has 76,712 atoms (about one third of the trimer system) and consists of one Cry4Aa toxin (9,771 atoms), 22,293 water molecules, 32 K + ions and 30 Clions.

MD simulations of pre-pore Cry4Aa trimer and Cry4Aa monomer in KCl solution
MD simulations of the modeled Cry4Aa trimer in KCl solution were performed using NAMD 2.6 program [29], CHARMM27 force field [28], and TIP3P model for water [30] on 96 CPU-cores of a 2.3 GHz AMD Opteron computer cluster at 2 days per ns. The sys-tem was energy minimized and then equilibrated for 10 ns at temperature 298 K and pressure 1 atm with all heavy atoms of the proteins under harmonic restraints before equilibration using a force constant of 2 kcal/mol/Å 2 . The simulations were performed using periodic boundary conditions. Temperature and pressure were controlled using Langevin temperature control with damping constant of 5 ps -1 and Langevin piston control [31] with an oscillation period of 200 fs and a damping time of 100 fs. An integration time step of 1.0 fs was used to allow a multiple time-stepping algorithm [32,33] to be employed such that long-range electrostatic interactions were computed every 4 time steps using the particle mesh Ewald method [34,35] with a grid spacing of ~1 Å for the reciprocal sum. van der Waals interactions were calculated every step using a 13 Å cutoff and a switching function. MD simulations of the Cry4Aa monomer system were performed in the same way as that of the trimer using the same simulation conditions on 32 CPU-cores of a 2.3 GHz AMD Opteron computer cluster at about 2 days per ns. From the simulations, both Cry4Aa trimer and monomer systems reached equilibrium within the first 5 ns and were very stable during the 10-ns equilibration.

Binding free energy calculations of the pre-pore Cry4Aa trimer in continuum solvent
To show that the pre-pore Cry4Aa trimer in solution is more stable than three isolated Cry4Aa monomers in solution, the method of molecular mechanics combined with Poisson-Boltzmann and surface area (MM-PBSA) calculations [36,37] was used to calculate the binding free energies of three Cry4Aa monomers in Cry4Aa trimer. Fifty snapshots of protein structures were taken at 100 ps interval from the last 5 ns of the two MD trajectories of Cry4Aa trimer and Cry4Aa monomer. For each snapshot, protein coordinate and structure files in AMBER file format were prepared for MM-PBSA calculations using AMBER force fields [38]. Then MM-PBSA calculation was performed on each snapshot using AMBER 9 [39] and then energies obtained were averaged over the fifty snapshots.
The binding free energies ΔG bind of the pre-pore Cry4Aa trimer were calculated as the difference between the free energy of Cry4Aa trimer ( ) and that of the three isolated Cry4Aa monomers ( ) according to a simple thermodynamics cycle: where the free energies and were defined as the sum of the molecular mechanical energies H MM , the solvation free energies G sol and the entropic energies - The molecular mechanical energy H MM was computed by summing five contributions, i.e. three bonded interactions due to bond stretching, bond bending and bond torsioning, and two non-bonded interactions due to electrostatic and van der Waals interactions, according to the following AMBER potential energy function [38] G tot The first term corresponds to bond stretching, where k b is a bond force constant and b b 0 is a bond extension from equilibrium. The second term is due to bond bending, where k θ is an anglular force constant and θθ 0 is an angular deviation from equilibrium.
The third term is for twisting of dihedral angles around bonds where k φ is a dihedral force constant, n is the multiplicity of the function, φ is a dihedral angle and δ is a phase shift. The fourth term accounts for the van der Waals (vdW) potential calculated using 12-6 Lennard-Jones potential. The fifth term corresponds to electrostatic potential between all (i, j) pairs of atoms in the system, where ε is the dielectric constant of the medium.
The solvation free energy G sol was calculated as the sum of its polar and non-polar proteins complexes seem to suggest that the vibrational entropic energy change due to complex formation is of order a few kcal/mol [40,41], so the vibrational entropy of the Cry4Aa trimer was approximated to be about three times that of Cry4Aa monomer.

Results and discussion
The results from the MD simulations suggest that the trimeric structure of the pre-pore Cry4Aa toxin is stable in aqueous solution during 10 ns as shown in Figure 1a and 1b.
Overall, the structure of Cry4Aa trimer did not change much from the initial structure with all three monomers still intact as can be seen from Figure 1c and 1d. This suggests that the interactions among the three Cry4Aa monomers and their surrounding water may be strong enough to bind them together in the solution on 10-ns timescale. So we investigated the structural changes and the binding between monomers by viewing the system, and calculating radii of gyration (R gyr ), root mean square deviations (RMSD), root mean square fluctuations (RMSF) and binding free energy. Finally, we proposed how Cry toxins may form pore.

Radius of gyration and RMSD of C α atoms versus time
To measure the overall change of the trimeric structure of the pre-pore Cry4Aa, the radius of gyration (R gyr ) and RMSD values of all C α atoms in the modeled trimer and its three monomers were computed. The R gyr plot in Figure 2a shows that the trimer has expanded radially from its center of mass by about 1.4 Å while the three monomers have expanded by about 0.5 Å from their centers of mass. The R gyr values of each monomer is about 1 Å less than that of the trimer because Cry4Aa trimer has a propeller-like shape which causes a small increase in R gyr of each monomer to result in a larger increase in R gyr of the trimer. The RMSD plot in Figure 2b also supports the 1-Å radial expansion of Cry4Aa trimer as it can be seen that the RMSD value for the trimer is about 1 Å larger than that of the monomer. Most of the changes in R gyr and RMSD values occur during the first 5 ns and after that the values become steady indicating that after 5 ns the system approaches equilibrium.

Structural fluctuation of C α atoms by residue
To investigate the structural flexibility of the pre-pore Cry4Aa trimer, RMSF value was calculated for each monomer on Cry4Aa trimer as shown in Figure 3. and α4, α4 and α5, and α5 and α6 in domain I. The loop connecting α4 and α5, which has been shown to be critically involved in forming pore [15], has low fluctuation, indicating that this loop is more rigid than other loops in Cry4Aa trimer. This rigidity has been

Amino acid residues at the inter-subunit interfaces
To indentify amino acid residues that participate in binding between two Cry4Aa monomers in the pre-pore trimeric structure, we counted amino acid residues to see how often they were within a distance of 3.5 Å from the inter-subunit interfaces of each Cry4Aa monomer from 1,000 frames of the 10-ns MD trajectory (10 ps per frame). The results are presented as a histogram shown in Figure 4. It can be seen that most of amino acid residues at the inter-subunit interfaces are mainly residues on α3 and α4 and partially on α2, the loop connecting α5 and α6, and the N-terminal of α6. From Figure 4, we have classified the amino acid residues that are within 3.5 Å distance from the inter-subunit interfaces of three Cry4Aa monomers into two groups. Group 1 consists of amino acid residues that are within the inter-subunit interfaces for almost 100% of the time or more with frequency more than 2,980 or more than ~1,000 from each interface (The frequencies of more than 3,000 come from some amino acid residues being at the interfaces near the center of the trimer, so they are counted more than once per interface). Group 2 consists of residues that are within the inter-subunit interfaces for more than 70% of the time with frequency more than 2,100 but less than 2,980. These two groups of amino acid are located on α2, α3, α4 and α6, and their vicinity as listed in Table 1. Residues in Group 1 are mainly at the center of the inter-subunit interface of Cry4Aa monomer while residues in Group 2 are on the peripheral of the inter-subunit interface as shown in Figure 5. It can be seen clearly that most of the amino acid residues at the inter-subunit interfaces of Cry4Aa monomers are polar residues (e.g. Tyr 179 , His 180 , Asn 183 and Ser 191 ) with a few charged (Arg 143 , Glu 187 and Lys 139 ) and non-polar residues (e.g. Val 147 , Leu 176 , and Pro 186 ).
By investigating the interfaces in details, it can be shown that the binding between two Cry4Aa monomers occurs specifically at the interfaces, for example, between α2 and α3 of one monomer (M3; in cyan and green), and α4, and α6 of adjacent monomer (M1; in blue and purple) as shown in Figure 6. The amino acid residues located at such interfaces are those listed earlier in Table 1 which dominantly are polar residues. In particular, Figure 6d shows that the polar residues from Group 1, such as Lys 139 , Asn 142 , Arg 143 , Asn 146 on α3, and Tyr 179 , His 180 , Gln 182 , Asn 183 , Glu 187 , Asn 190 , Ser 191 on α4, and Thr 245 on α6, are key residues that interlock with each other and form one of three inter-subunit interfaces within the Cry4Aa trimer.

Binding free energy of the pre-pore Cry4Aa trimer
From Table 2, the negative binding free energy (-100 kcal/mol) of the pre-pore Cry4Aa trimer in a continuum solvent obtained from the method of molecular mechanics combined with Poisson-Boltzmann and surface area (MM-PBSA) [36,37] suggests that Cry4Aa trimer is more stable than three separated Cry4Aa monomers in solution. The binding free energy is contributed negatively by the non-polar interacions between the proteins and water (-446 kcal/mol), the van der Waals interactions within the proteins (-10 kcal/mol), and possibly the vibrational entropy of the proteins. It is contributed positively by the polar interactions between the proteins and water (43 kcal/mol), the electrostatic interactions within the proteins (143 kcal/mol), the conformational changes of the proteins (105 kcal/mol) as well as the entropies due to translations and rotations of the proteins (34 and 36 kcal/mol). The large negative binding free energy contribution of Figure 4 Frequency distribution of amino acid residues at inter-subunit interfaces of the pre-pore Cry4Aa trimer. Frequencies that amino acid residues appear within a distance of 3.5 Å from the inter-subunit interfaces of three Cry4Aa monomers in the pre-pore trimer during 10-ns MD simulations. The dotted line at frequency of 2,980 is used to indicate amino acid residues that are in contact 100% of the simulation time while the dotted line at 2,100 is for 70%.
-446 kcal/mol due to the non-polar interactions between the proteins and water possibly came from the formation of the inter-subunit interfaces between the three Cry4Aa monomers during their trimer formation. This removed water molecules from the intersubunit interfaces of the three monomers and therefore decreased their water accessible surface area.
On the other hand, the large positive binding free energy of 143 kcal/mol caused by electrostatic interactions within the proteins possibly arose from the net charge of -2 on each Cry4Aa monomer which may cause the increase in the free energy after trimer formation. Another 105 kcal/mol contribution of binding free energy probably resulted from the conformational stress within the proteins due to the inter-subunit interfaces of Cry4Aa trimer. Moreover, the binding free energy of 43 kcal/mol due to polar-solvation energy was likely contributed by the removal of water molecules from the inter-subunit interfaces as this decreased the electrostatic attractions between the proteins and water. However, further investigation using a combination of MD simulations of Cry4Aa mutants and MM-PSBA may be needed to identify key residues at the interfaces that mainly contribute to the negative binding free energy.

Figure 5
Amino acid residues at the inter-subunit interface of a Cry4Aa monomer within the pre-pore trimer. (a) side view with Group 1's residues rendered as van der Waals (vdW) spheres and Group 2's residues rendered as sticks. Green is for polar and white is for non-polar. Blue is for positively charged and red is for negatively charged. One of Cry4Aa is depicted as surface and colored in spectrum from red to white to blue by residue number. The other two Cry4Aa monomers M2 and M3 are rendered as lines and colored red and green. (b) top view of (a), (c) Group 1's residues rendered as sticks, and (d) Group 2's residues depicted as sticks with Group 1's residues depicted as vdW spheres. a b c d

Proposal of a pore-forming mechanism for Cry toxins
Why Cry toxins utilize hydrophilic core to form trimer and reduce their non-polar interactions with water? It has been shown that globular proteins can utilize a various means of protein interaction to form oligomers such as dimers [43]. One is to have a hydrophobic core at their interface while another is to have the opposite which is a hydrophilic core like Cry4Aa as observed here. Other different interaction may use a combination of both hydrophobic and hydrophilic patches at protein interface. In the case of Cry4Aa toxin, this may be related to their pore-forming function that requires to have polar environment for the pore interior. In this work, we have also observed a chloride ion being trapped by three Lys 139 in a large water cavity at the bottom of the interface within the Cry4Aa trimeric structure as shown in Figure 7. It is possible that, during pore formation, water molecules and ions are involved in breaking of some hydrogen bonds, and salt bridges, and overcoming steric interactions between helices in domain I of Cry toxins so that α4 and α5 can be inserted into the membrane bilayer and rearranged to form pore.
Based on several experimental evidences and the obtained full-atom structural model of Cry4Aa trimer, we would like to propose a mechanism that suggests how Cry toxins Two groups of amino acid residues located at the inter-subunit interface of Cry4Aa monomers M i and M i+1 in Cry4Aa trimer where i and i+1 run in cyclic of 1, 2, and 3. a The percentage of more than 100% was due to some residues were counted more than once because they were within two or more inter-subunit interfaces of Cry4Aa trimer. b Residues on the loop between α5 and α6. c Residue on the loop between α4 and α5.
may form pore. Recent electron crystallographic images of Cry4Ba toxin have shown that there is a pore at the center of Cry4Ba trimer [5]. These pores formed by Cry4Ba toxins can release calcein molecules trapped in lipid vesicles [4]. The pore diameters of Cry4Ba and Cry1C toxins [4,44] have been estimated to be about 20-26 Å which is about the diameter of an alpha helix (including sidechains). It has been shown that the pore of  Cry4Ba toxin is formed by the insertion of α4-α5 hairpins into the membrane [45,46]. From our obtained structure of Cry4Aa trimer shown in Figure 8a, we think that the prepore Cry4Aa trimer has to be oriented relative to the membrane such that all three domain II of Cry4Aa trimer can interact with the receptors on the membrane surface [17,18] as well as all three α4-α5 hairpins are oriented in their inserting positions with their loops pointing toward the membrane. However, from the same figure if the three α4-α5 hairpins are inserted directly into the membrane by a distance of ~40 Å they cannot form pore and some reorientations of the three hairpins are required for all six helices to form a barrel. Therefore, during insertion the three hairpins in Figure 8a and 8b have to be rotated by about 90° clockwise so that all six helices can form a barrel that has a channel with the size of an alpha helix (about 15-20 Å) at the center as shown in Figure  8c and 8d. The opening and closing of the pore may be controlled by the tilting of the three α4-α5 helices relative to the normal to the membrane.

Conclusions
In this work, we have constructed the first all-atom structural model of the pre-pore Cry4Aa trimer based on the results from electron and X-ray crystallographic studies of Cry4Aa and Cry4Ba toxins. By using MD simulations, we have found that Cry4Aa toxin uses polar amino acid residues on helices α3, α4, and α6 to form trimer. A combination of MD and MM-PBSA calculations has suggested that Cry4Aa trimer is stable in aqueous solution mainly due to its formation of the inter-subunit interfaces which reduces the non-polar interactions between the proteins and their surrounding water. The results can later be investigated further by mutagenic studies to confirm the simulation results and help us understand more on the toxin oligomerization process. Moreover, we have proposed a pore-forming mechanism by suggesting that three α4-α5 hairpins are inserted into the membrane with each hairpin rotated by 90° clockwise (viewed from the top) to form pore during the insertion step of Cry toxins. This process may involve water molecules and ions.