1. Introduction

Ion channels are proteins that allow the passage of ions, especially Na⁺, K⁺, and Ca⁺⁺ across membranes. In this discussion, we will give most examples from a particular potassium channel. However, we will be primarily concerned with the question of the best means of obtaining information from calculations, which we will apply much more generally. In essence, there are two approaches to calculations on proteins, and channels in particular. The calculations can be done using classical or quantum physics, and we will argue that quantum physics is necessary. For one thing, there are certain quantum phenomena that are essentially impossible to simulate using classical physics. Of these, tunneling is almost certainly relevant to ion channels. Furthermore, certain energy terms have no classical analogues. While most of the energy has classical analogues (kinetic energy and Coulomb terms in the potential energy, for example), exchange and correlation energies are strictly quantum phenomena. Not all interactions scale with the same distance dependence; in short-range interactions (i.e., at roughly bonding distances, where quantum effects dominate), the dependence on interatomic distance does not scale in a simple manner. Coulomb interactions, on the other hand, scale as 1/r, and hence are much longer range. At interatomic distances greater than about 4 Å, Coulomb interactions dominate, and a classical approximation can be accurate; bonding interactions occur at somewhat shorter distances; those interactions that are inherently quantum are much more difficult to model classically, and switches in bonding, especially hydrogen bonding, are particularly difficult. For ion channels, therefore, where the interactions concern chemical bonding, and hydrogen bonding, it is necessary to use quantum calculations to be sure that all interactions are correctly modeled, with the correct distance dependence. We believe that proton tunneling, and proton transport, are critical. Proton transport is known to happen, but in gating models that have become standard, this is not considered relevant to channel function. Our calculations suggest that it is fundamental to channel function. The attempt to produce classical force fields that are sufficiently accurate to simulate channels, for example (other protein systems have similar requirements) has led to a great deal of effort ^[1][2][3]; this has succeeded only to a limited extent. This said, molecular dynamics (MD) calculations are much more commonly used (there is also a less commonly used classical physics alternative, Monte Carlo simulations, which do not do time dependence, but give a Boltzmann distribution of states). Simulations have been reviewed repeatedly (e.g., for ion channel interactions and other protein cases ^{[4][5][6][7][8]}).

There is also the possibility of doing QM/MM calculations (QM = quantum mechanics, MM = molecular mechanics, which is classical); this allows the key section of the protein to be done properly as a QM calculation, while the background is approximated much better than a vacuum; a pure QM calculation of a section of protein has a vacuum outside, albeit at a greater distance than the boundary between the QM and MM sections of the QM/MM calculations ^[9][10][11]. This is not the only multiscale method; the development of multiscale methods is ongoing, from many laboratories. The major advantage is that much more of the protein can be included, albeit at the MM level, while the main section of the calculation is done accurately at the QM level; the entire calculation is much faster than doing a large section of protein entirely at the QM level. The principal disadvantage is that one must create a boundary between the QM and the MM sections. Sometimes this can be done without causing major disruption. However, covalent bonds must usually be broken to create the boundary. This said, the smaller QM section makes the boundary closer to the key section. If the boundary introduces an extra dipole, or otherwise affects the local field at the main QM section, it can be worse than a vacuum with a somewhat more distant boundary. However, QM/MM is often considered as a useful compromise between the requirements for computer resources and the requirements for QM accuracy in calculation. For QM calculations that are intended to trace a path, whether of a K+ or H+, only differences are relevant, so that the boundary effects are largely subtracted out. This said, it would be better to include more, and it is not possible with limited computer resources.

In addition to splitting the calculation into QM and MM sections, other multiscale simulations are possible, in which different levels of computation can be applied to different parts of the system. For example, an outer part can be coarse-grained, in which groups of atoms are combined into a single item that can be treated with the same computer resources as a single atom. There are other possibilities, and a number of these have been explored ^{[12][13][14][15]}.

This is a very limited summary of the reported calculations on ion channels—only a tiny fraction of the literature can be cited here. However, this sample suggests what alternate possibilities have been applied to ion channels and gives a suggestion of what alternatives exist for calculations on ion channels. We chose the full quantum mechanical approach, and we will outline our reasons below. However, as exascale computing becomes available, less approximation will be necessary, so that quantum calculations will become more practical. Although exascale computing is not yet generally available, it makes quantum calculations become of greater significance as the method of choice.

2. Our Calculations

Our calculations have other limitations. For one, the calculation is at 0 K. When comparing energy, this is not terrible, as the comparison is between two systems of the same size and approximately the same conformation. Therefore, the vibrations should be similar enough that the difference in energy at around 300 K would not normally be significantly different from the difference in ground state energy alone, which is what we calculate. There is another danger, that the system settles into a local minimum far from the global minimum. In one sense, we depend on this happening: we start the system with a proton or ion in a particular location, and do not expect that proton or ion to travel a considerable distance to find a different minimum. Instead, the optimization may, for example, move a proton to the other side of a hydrogen bond, and rotate a side chain as a consequence; if the calculation continued to find a global minimum every calculation would give the same result, after much too long a time to be practical. However, there exists a chance of convergence to a nearby very shallow, hence unimportant, minimum. Thus far, we have found no evidence of this, as the algorithm appears to move the system out of any minimum it comes close to, until it settles into convergence in a minimum too deep to leave. Nevertheless, it would be useful to confirm this, and we cannot at this point do so.

Gaussian09 ^[16], which is the computational suite we used the most, allows for the calculation of thermodynamic functions at finite temperature, in the harmonic approximation. However, the determination of the normal coordinates demands too much memory, using present resources, for systems above about 300 atoms. There is an additional question in doing such a calculation: normal coordinate analysis is a reasonable approximation for systems that are not too large. However, for a 1000 atom system, there should be many modes with anharmonic behavior. It would seem more likely that the appropriate approximation would be something more like a solid, with long wavelength modes in addition to more local, and more nearly harmonic modes of higher frequency. By analogy with solids, the latter would be more like Einstein modes. At this point, we do not have the means to examine this question, nor to determine whether it is significant or not. One thing that gives us confidence is that the modes should be very similar for two systems with the same structure, albeit with a proton shift or a water rotation. As a consequence, we continue to have confidence in comparing energy. This is supported by the fact that completely independent calculations differ by a small amount in reasonable directions; that is, each energy calculation starts from essentially the same position, with a proton shift or other small change; each optimization is independent. If the results were extremely sensitive to small perturbations, this would appear as large random changes in energy, something we do not see.