Computational methods and theory for ion channel research

ABSTRACT

Ion channels are fundamental biological devices that act as gates in order to ensure selective ion transport across cellular membranes; their operation constitutes the molecular mechanism through which basic biological functions, such as nerve signal transmission and muscle contraction, are carried out. Here, we review recent results in the field of computational research on ion channels, covering theoretical advances, state-of-the-art simulation approaches, and frontline modeling techniques. We also report on few selected applications of continuum and atomistic methods to characterize the mechanisms of permeation, selectivity, and gating in biological and model channels.

GraphicalAbstract

KEYWORDS:

• Ion channels
• continuum models
• molecular dynamics
• rare events
• conductance
• selectivity
• gating
• machine learning
• biomimetic nanopores

1. Introduction

Ion channels are proteins delimiting pores that allow a regulated flow of water and ions across biological membranes. The first ion channel to be structurally characterized was, in 1998, KcsA, a potassium ion channel [1] (Figure 1, left panel). Ion channels are involved in key biological functions including control of homeostasis, muscle contraction, and the propagation of nerve signals [2]. The importance of ion channels cannot be overstated when considering that they represent the molecular basis of sensory perception and human thought. Their biomedical relevance is related both to their role in channelopathies [3] and to their role as drug targets for a wide range of diseases.

Figure 1. Left) Crystal structure of KcsA in the closed state (PDB ID: 1J95 [6],); two diagonally opposite monomers are shown for clarity, in cartoon representation. Potassium ions at the binding sites S₁-S₄ are represented as spheres using their van der Waals radius. Extracellular (EC) and intracellular (IC) boundaries of the lipid bilayer are also indicated. Protein sketch created using VMD [7]. Right) Cartoon illustrating the complex electrokinetic phenomenology occurring in a wide pore (biological or solid-state) under the effect of an external electric field. The walls of the pore can acquire surface charge, the field generates a current ($j^{\pm }$) of positive and negative ions. The unbalanced flow of ions, in turn, can result in a net transport of water molecules, $Q$, (electroosmosis). The arrow profile in the middle of the pore, represents the presence of the electroosmotic flow, or any other flow generated, e.g. by external gradients.

Ion channels are characterized by three main features [2]: (i) high permeability, that allows conduction rates close to the free diffusion limit; (ii) high selectivity, that allows discrimination between ions with the same or different charge; (iii) pore gating, that is, the highly regulated opening and closing of the pore in response to specific stimuli. Most research on ion channels focuses on the clarification of the molecular basis of these properties or how they can be modulated by ligands. In this review, we present current theoretical and numerical approaches to address this challenge.

The high conduction rates of ion channels (10⁷–10⁸ ions/second for KcsA [4]) are somehow counterintuitive, given the narrow pore connecting the two sides of the membrane. The high complexity of the permeation mechanism which allows for this extraordinary efficiency is exemplified by the case of voltage-gated potassium channels (see [5]). The Selectivity Filter (SF) of K⁺ channels is formed by the interface of four subunits where the carbonyl groups of five or six residues point toward the center of the pore defining four binding sites (labeled S₁ to S₄ from the outermost to the innermost, see Figure 1, left panel), which can be occupied by potassium ions.

Traditionally the anomalous scattering data from the bacterial K⁺ channels KcsA from Streptomyces lividans were interpreted [8] as a superposition of ion- and water-occupied states leading to a model of cotranslocations of ions and water. This model was further supported by numerous simulation studies [9–14]. However, this model has more recently been questioned [15,16], and an alternative view has been proposed, with ions passing through the selectivity filter in direct contact with each other without intervening water molecules. This interpretation directly challenged the classical view on permeation: the direct ion–ion interactions had been so far ruled out on account of an excessive electrostatic repulsion. While the issue is far from settled, this example simultaneously highlights the complexity of structural biology and electrophysiology data whose interpretation is far from trivial, and the power of computational approaches, which serve as a bridge between the two, providing molecular-level dynamical information.

Ion channel selectivity is another area of active research. It is again instructive to consider the case of KcsA that conducts K⁺ and Na⁺ ions with a ratio 150:1 [17]. Again, this property is surprising since K⁺ ions are larger than Na⁺ and selection cannot be explained by simple sieving based on ion size. This pattern was originally explained with a snug fit model [1,18] according to which the channel was potassium-selective because the K⁺ ion perfectly fitted into the SF, but Na⁺ was too small to favourably interact with the pore walls. This model, however, became untenable when it was discovered that the SF of this channel was capable of fluctuations larger that the difference in radius between the potassium and sodium ions [19]. The model was thus updated with a more sophisticated one also accounting for the electrostatic repulsion between the carbonyls of the SF [20,21] and limitations on the number and motion of coordinating carbonyls [22–24] as well as the interactions between multiple ions in the pore [25]. The emerging modern view of ion selectivity was reviewed in Ref [26]. In general, the size-exclusion model does not explain the counter-intuitive evidence that often large channels are selective to small ions while small channels are selective to large ions [27]. This pattern emerges as a result of the hydration shell that surrounds the ions. A small ion can cross a large pore keeping its hydration shell almost intact. Conversely, both small and large ions are massively dehydrated when crossing narrow pores but the energy penalty paid by a large ion is normally smaller than that of a small ion that interacts more tightly with the surrounding water. The redistribution of water around an ion crossing a nanopore is thus a topic of great relevance in bio- and nano-technology and can be studied computationally using, for example, Kirkwood approximation [28] (see Section 2.4).

The high selectivity of ion channels is seemingly at odds with their high permeation rates, leading to the so-called conductivity-selectivity paradox. The traditional argument [29] relies on the scheme that, if the channel is selective to a given ion, there should be a strong ion-binding site, with a deep energy well. However, the deeper the well, the slower the release of the ion, which leads to the prediction of small permeation rates in disagreement with experimental data. Many research efforts via theoretical and simulation-based approaches are currently underway to try and reconcile selectivity and permeability of ion channels [4,30,31]. A possible solution to the paradox relies on the presence of two or more neighbouring binding sites [32], or small energy steps that allow an ion to escape the deep energy well [33]. More research efforts are required in order to understand which solution of the paradox applies to each of the different families of ion channels.

Another key characteristics of ion channels is their gating capability. Traditionally [2] ion channels were classified as ligand-gated and voltage-gated, according to whether their opening is triggered by the binding of a ligand or by a change in the membrane potential. Further complex fine-tuning mechanisms of control, including allosteric modulators, can participate in the gating control. Recently, other gating mechanisms have been identified, based on channels sensitivity to pressure or temperature [34–37]. In general, how the information is transmitted to the pore gate and how gating can be influenced by a variety of other subtle factors remains unclear. Furthermore, there is growing interest to identify the physical means by which the pore blocks ion permeation: is steric occlusion the only one? For example, in the concept of hydrophobic gating [38], a pore with hydrophobic lining can be closed by the formation of a bubble that functionally occludes the pore even in the absence of steric block. This mechanism was originally observed in model nanopores [39,40] but it has been recently identified in an increasing number of biological channels [13,41–43]. Physically, it corresponds to a process of evaporation in hydrophobic nanoconfinement [44,45] which gives rise to bubble formation, see Sections 2.5 and 6.6. Hydrophobic gating implies that the analysis of the radius profile is no longer sufficient to classify a channel structure as open or closed. New criteria for the annotation of newly resolved channel structures are called for [46]. Both equilibrium and enhanced sampling simulations are proving to be efficient tools for the study of hydrophobic gating. A case study [47] of hydrophobic gating in artificial nanopores can be found in Section 6.6.

Even a well-established mechanism like voltage-gating is far from being fully understood. While the gating mechanism of domain-swapped K⁺ channels seems consolidated [48], the gating mechanism of non-domain-swapped channels is still the object of debate. Domain-swapped channels are characterized by a long linker helix (L45) that connects the Voltage Sensor domain (VSD) with the Pore domain (PD). It is currently accepted [48] that the displacement of the sensor helix S4 pushes down the L45 linker that acts as a mechanical lever on the pore helix S6, extending it and closing the pore. In non-domain-swapped channels, however, the linker element L45 is too short to act as a mechanical lever and the gating mechanism is still under investigation. Voltage-gated ion channels have recently been classified as allosteric machines [48] since they are characterized by a long-range coupling between the VSD and PD. This suggested the opportunity to use network theoretical approaches already used in the study of allosteric enzymes [49,50]. These techniques allow one to identify pathways of motion propagation from the VSD to the PD. Recent studies [51,52] applying this methodology identified two main pathways, the former reminiscent of the mechanism operating in domain-swapped channels, and the latter corresponding to a new, non-canonical mechanism important not only for activation/deactivation, but possibly also for inactivation [53]. A case study can be found in Section 6.5.

Ion channels have not only a biomedical relevance but they are also extremely important in nanotechnology and engineering. Currently, it is possible to embed ion channels or other nanopores in natural or artificial bilayers and use them as devices for single-molecule manipulation. One of the most successful applications is the use of $\alpha$-hemolysin for DNA sequencing [54] which is currently commercially exploited. Ion channels have also been engineered to respond to specific stimuli [55–57]. Unfortunately, ion channels are extremely fragile molecules that tend to unfold and lose their properties when placed outside their natural biological environment. This has motivated massive research efforts to develop artificial solid-state nanopores [58,59]. Unfortunately, so far synthetic nanopores are much less performing than their biological counterparts; in particular, they typically lack two distinguishing properties of ion channels, that is, selectivity and gating. The goal of much research in the field is thus to design biomimetic artificial nanopores [60–62].

As an example, a promising material to realise biomimetic nanopores is represented by graphene [63]. The lithography and beam irradiation techniques allow the creation of nanoscale pores in graphene mono-layers, which opens the way to a wide range of applications, from wastewater treatment to desalination [64,65]. In order to use graphene as a molecular sieve, it is necessary to enforce a uniformly small diameter of the pores. Small pores with a diameter up to 5.5 Å are naturally selective due to steric exclusion [66], but their fabrication is still a challenge [67]. Larger pores are easier to manufacture but they are not selective. However, even large nanopores can be made selective through an appropriate functionalization [62]. The quest for the ideal functionalization simultaneously allowing high selectivity and high permeation rates of the selected solute brings back to biological models. In a landmark study [68], Corry and coworkers created graphene nanopores functionalized with four carbonyl groups arranged so as to reproduce the geometry of the SF of the K⁺ channel KcsA: the synthetic pores indeed were potassium selective. The performance of biological channels, however, does not only depend on the arrangement of the SF groups. Indeed, graphene nanopores with four carboxylate groups, engineered to mimic the SF of the Na⁺ channel NavAb, resulted in pores not selective to sodium [68]. However, a voltage-tunable ion selectivity could be enforced using only three carboxylate groups.

Computational modelling has always been an integral part of the research on ion channels. For instance, Hodgkin and Huxley [69] described the onset and propagation of the action potential modeling the membrane as a non-linear electric circuit. It is notable that Hodgkin and Huxley, who won the 1963 Nobel Prize in Physiology or Medicine for this work, did not know about the existence of ion channels in the same way as Mendel did not know about the existence of genes when he discovered his laws of transmission and segregation of inherited traits. The discovery of specific proteins acting as pathways for ionic currents across the membranes dates back to the patch clamp studies by Sakmann and Neher [70] who were awarded the 1991 Nobel Prize in Physiology or Medicine. The first structural resolution of an ion channel, however, had to wait until McKinnon (2003 Nobel laureate in Chemistry) determined the X-ray structure of the KcsA channel [1]. We are currently living in a particularly favourable historical period for the computational study of ion channels. On the one hand, cryo-electron microscopy is providing a large number of high-resolution structures of ion channels [71]. On the other hand, hardware improvements are providing the scientific community with increasing computing power [72]: as of November 2021 the computing speed of Fugaku, the fastest supercomputer, is ca. 442 petaflops [73], that is, more than 10¹⁷ floating point operations per second. Finally, researchers can now rely on a rich toolbox of advanced computing techniques that are the topic of this review together with theoretical advancements in the field.

Notwithstanding the recent advancements, the characterization of ion channels still remains a theoretical and computational challenge; accordingly, the choice of the approach must be performed with care based on the task at hand. The calculation of electrostatic potential profiles or ion flows highly benefits from continuum approaches where protein, water, and membrane are modeled as continuous dielectric media while ions are described by a continuous density distribution. This approach, traditionally relying on the Poisson-Boltzmann equation [74] and on Poisson-Nernst-Planck (PNP) equation [75], is computationally convenient and extremely flexible. Indeed, combination of PNP equation with the Navier-Stokes equation allows to include hydrodynamic effects yielding the system of electrokinetic equations [76] presented in Section 2. The continuum approach, however, does not account for fluctuations, ion correlations, polarization effects, and finite-size effects. When these effects become relevant, atomistic approaches are the tool of choice, see Section 6. In particular, Molecular Dynamics (MD) simulations can be considered as a computational microscope [77] with a time resolution of the order of the femto-second and a space resolution of the order of the hydrogen atom. The main limitation of MD is that, due to the constraint to use a femto-second time-scale, the simulations are normally limited to a few hundreds of nanoseconds (although the growing use of GPUs [78] or dedicated hardware like the Anton supercomputer [79] allow for longer simulations). This limitation is particularly serious because many phenomena of interest in biology can be classified as rare events [80], see Section 4.2. Rare events are processes associated with the overcoming of an energy barrier higher than the thermal energy $k_{B}T$. As a result, the system lingers for a very long time in the pre-barrier minimum before a fluctuation endows it with sufficient kinetic energy to overcome the barrier. Rare events are thus characterized by long waiting times, but when they finally occur, the barrier crossing is very rapid and the sampling of the region near the barrier is very poor.

In order to overcome this limitation, one is forced to develop tools that enhance the sampling of MD simulations or bias them in a controllable way (see section 4.2), to focus on unlikely and yet important regions of the free-energy landscape. Methods like Replica Exchange Molecular Dynamics [81] do this by overcoming free-energy barriers thermally through exchanges with high-temperature replica simulations. These can suffer, however, from needing a very large number of replica simulations. Alternatively, biases can be applied to an appropriate set of collective variables, or CVs, to characterize those regions and to define the corresponding free-energy landscape. Choosing appropriate CVs, is however, a non-trivial task. A clever approach is to create hybrid algorithms where the strengths of Replica Exchange and CV-based algorithms (e.g. metadynamics) can be combined [82,83]. As an alternative, if the goal is not the reconstruction of the whole free-energy landscape but just the calculation of the free energy profile along the transition of interest, it is possible to use methods like Transition Path Sampling [84] or the String Method [85]. Finally, Machine Learning is disclosing the opportunity to train a deep neural network to automatically choose good CVs [86], for instance, those along the direction of maximal auto-correlation [87,88], see Section 7.

This Review ensues from the talks and the scientific discussions held in February 2021 at the conference ‘Frontiers in ion channels and nanopores: theory, experiments, and simulation’ [89]; while the overview of the approaches attempts to be general, the selection of applications reflects its origin. The Review has the following main organization. Section 2 is devoted to continuum methods and Section 3 to their applications; these Sections introduce a general physical framework, which is useful to understand the fundamental transport phenomena occurring in ion channels but also in synthetic nanopores; while the continuum approach is highly succinct and informative, it requires some formalism to be introduced. We propose a number of selected applications referring to the use of electrokinetic equations, the methods to compute hydration patterns, and the use of density functional theory. Section 4 of the review covers atomistic simulation approaches of ion channels which provide a detailed picture of the same phenomena and whose treatment benefits from a less formal introduction: after an overview, we introduce the different families of methods for the study of rare events and we review a hybrid continuum/atomistic approach: Brownian Dynamics. In Section 6, we discuss few selected applications of atomistic methods to characterize selectivity, permeation, and gating in biological and model channels. In Section 7, we discuss a number of forefront techniques that are still in a development stage but have a potential to revolutionize the field. In particular, we focused on polarizable force fields, hybrid enhanced sampling algorithms, multi-resolution approaches, and machine learning. In the final section, we draw the conclusions of the work.

2. Continuum methods

The transport of solutes, including ions, across narrow pores is not only important for ion channel research, but is also at the core of many nanotechnology applications in the field of nanopore sensing [90], nanofiltration [91], and nanoporous energy materials [92]. Two main computational approaches are useful to address this general problem: atomistic/coarse-grained simulations and continuum methods. While the former provides a detailed (atomistic) or approximated (coarse-grained) description of chemical interactions occurring at relatively short-time scales, continuum theories, which neglect the granular nature of matter, allow to explore transport mechanisms on longer time scales, moreover yielding, in some cases, analytical estimates of the relevant observables characterising transport. In the following, we provide a rather general continuum framework to study the said transport phenomena and we specialise them for ion channel research.

2.1. Electrokinetic equations

In the continuum approach, the complexity of the confined electrolyte transport within artificial or biological nanopores, such as that in Figure 1, right panel, is drastically reduced by using local densities and fields associated to each component of the system shown in Figure 1, left panel. For instance, the solvent (e.g. water) is generally treated as a homogeneous medium with a constant electric permittivity $ϵ$, the ion species are described by their local concentrations, and channel walls contribute with a surface charge density. Within this framework, the transport of ions is determined by electrokinetic equations [93,94], that for systems with axial symmetry, read [94]

(1a) $\frac{\mathrm{\partial }{\rho }^{\pm }}{\mathrm{\partial }t}=-\mathrm{\nabla }\cdot {\mathit{j}}^{\pm }$(1a)

(1b) $\varrho \left(\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+\mathit{v}\cdot \mathrm{\nabla }\mathit{v}\right)=-\mathrm{\nabla }{P}_{tot}+\eta {\mathrm{\nabla }}^2\mathit{v}+{\mathit{F}}_{tot}$(1b)

(1c) $\mathrm{\nabla }\cdot \mathit{v}=0.$(1c)

Equation (1a)(1a) $\frac{\mathrm{\partial }{\rho }^{\pm }}{\mathrm{\partial }t}=-\mathrm{\nabla }\cdot {\mathit{j}}^{\pm }$(1a) is the continuity equation specifying that the variation of the density, ${\rho }^{\pm }$, of positive/negative ions in a volume is related to their fluxes across the surface of that volume. Equation (1b)(1b) $\varrho \left(\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+\mathit{v}\cdot \mathrm{\nabla }\mathit{v}\right)=-\mathrm{\nabla }{P}_{tot}+\eta {\mathrm{\nabla }}^2\mathit{v}+{\mathit{F}}_{tot}$(1b) is the Navier-Stokes equation governing the evolution of the velocity profile of a solution of density $\varrho$ and viscosity $\eta$, driven by external forces ${\mathit{F}}_{tot}$ and, possibly, pressure gradients $\mathrm{\nabla }{P}_{tot}$. Finally, Equation (1c)(1c) $\mathrm{\nabla }\cdot \mathit{v}=0.$(1c) defines the incompressibility condition for the solution, an approximation usually well satisfied by common salt solutions. The physical nature of the ion transport is assigned by specifying the constitutive equation of the ion fluxes in Equation (1a)(1a) $\frac{\mathrm{\partial }{\rho }^{\pm }}{\mathrm{\partial }t}=-\mathrm{\nabla }\cdot {\mathit{j}}^{\pm }$(1a) . A traditional approach employs the so-called Nernst-Planck current

(2) ${\mathit{j}}^{\pm }(\mathit{r},x,t)={\rho }^{\pm }(\mathit{r},x,t)[\mathit{v}(\mathit{r},x,t)-D\beta \mathrm{\nabla }{\mu }^{\pm }(\mathit{r},x,t)],$(2)

where $x$ is the coordinate along the longitudinal axis of the channel (transport direction) and $\mathit{r}=(y,z)$ indicates the radial coordinate, see Figure 1, right panel. Equation (2)(2) ${\mathit{j}}^{\pm }(\mathit{r},x,t)={\rho }^{\pm }(\mathit{r},x,t)[\mathit{v}(\mathit{r},x,t)-D\beta \mathrm{\nabla }{\mu }^{\pm }(\mathit{r},x,t)],$(2) expresses the ion currents ${\mathit{j}}^{\pm }$ in terms of advection by the velocity $\mathbf{v}$, diffusion, and electromigration due to the gradient of the electrochemical potential, which includes ion concentration and charge effects

(3) ${\mu }^{\pm }(\mathit{r},x,t)=k_{B}Tln[{\rho }^{\pm }(\mathit{r},x,t)]\pm ze\psi (\mathit{r},x,t),$(3)

where $D$ is the diffusivity of the ions, $z$ denotes their valence, $e$ is the elementary charge, and $\beta \equiv 1/k_{B}T$, with $k_B$ the Boltzmann constant and $T$ the absolute temperature. The electrostatic potential $\psi (\mathit{r},x)$ in the channel is obtained by solving the Poisson equation

(4) $ϵ{\mathrm{\nabla }}^2\psi \left(\mathit{r},x,t\right)=-ze\left[{\rho }^{+}\left(\mathit{r},x\right)-{\rho }^{-}\left(\mathit{r},x\right)\right]$(4)

relating $\psi (\mathit{r},x,t)$ to the concentrations ${\rho }^{\pm }(\mathit{r},x,t)$ of two symmetric ion species; $ϵ$ is the permittivity of the solvent taken as an homogeneous medium. The typical size of ion channels leads to vanishing small values of the Reynolds number¹ and hence the Navier-Stokes Equation (1b)(1b) $\varrho \left(\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+\mathit{v}\cdot \mathrm{\nabla }\mathit{v}\right)=-\mathrm{\nabla }{P}_{tot}+\eta {\mathrm{\nabla }}^2\mathit{v}+{\mathit{F}}_{tot}$(1b) reduces to the (steady-state) Stokes equation:

(5) $\eta {\mathrm{\nabla }}^2\mathit{v}(\mathit{r},x)=-{\mathit{F}}_{tot}(\mathit{r},x)+\mathrm{\nabla }{P}_{tot}(\mathit{r},x).$(5)

Equation (4)(4) $ϵ{\mathrm{\nabla }}^2\psi \left(\mathit{r},x,t\right)=-ze\left[{\rho }^{+}\left(\mathit{r},x\right)-{\rho }^{-}\left(\mathit{r},x\right)\right]$(4) and (5(5) $\eta {\mathrm{\nabla }}^2\mathit{v}(\mathit{r},x)=-{\mathit{F}}_{tot}(\mathit{r},x)+\mathrm{\nabla }{P}_{tot}(\mathit{r},x).$(5) ) are subject to standard channel impermeability and no-slip boundary conditions for the flow ($\mathit{v}=0$ at the channel walls) as well as $\mathit{n}\cdot \mathrm{\nabla }\psi =\sigma$ or $\psi =\zeta$ at the walls, which account for the dielectric (with surface charge density $\sigma$) or conductive (with $\zeta$ potential) nature of the walls [76], respectively, with $\mathit{n}$ the normal to the wall surface.

Analytical insight into the steady-state solution of Equation (1a)(1a) $\frac{\mathrm{\partial }{\rho }^{\pm }}{\mathrm{\partial }t}=-\mathrm{\nabla }\cdot {\mathit{j}}^{\pm }$(1a) and (4(4) $ϵ{\mathrm{\nabla }}^2\psi \left(\mathit{r},x,t\right)=-ze\left[{\rho }^{+}\left(\mathit{r},x\right)-{\rho }^{-}\left(\mathit{r},x\right)\right]$(4) ) can be obtained by assuming that, along the radial direction, the charge densities attain their equilibrium profile

(6) ${\rho }^{\pm }(\mathit{r},x)={\bar{\rho }}^{\pm }(x)e^{\mp ze\psi (\mathit{r},x)}$(6)

and by linearising the charge density [76]

(7) $q(\mathit{r},x)\simeq {\kappa }_D^2\psi (\mathit{r},x)$(7)

where ${\kappa }_D=1/{\lambda }_D$ is the inverse of the Debye length

(8) ${\lambda }_D=\sqrt{\frac{ϵ}{2\beta {\left(ze\right)}^2{\rho }_0}}$(8)

which measures the length scale at which the re-distribution of mobile ions screens an electric field, with ${\rho }_0$ the bulk concentration of the salt. In the case of smoothly-varying channel profile, a linear-response theory can be applied. In this regime, the total pressure varies solely along the pore axis, such that its gradient is $\mathrm{\nabla }{P}_{tot}(\mathit{r},x)=\mathrm{d}P(x)/\mathrm{d}x+\mathrm{\Delta }P/L$, where $\mathrm{d}P/\mathrm{d}x$ is the $x$ component of the geometrically induced local pressure gradient which is determined by the boundary conditions and by fluid incompressibility, while $\mathrm{\Delta }P$ indicates the pressure drop across the channel length $L$. Moreover, within such a regime, ${\mathit{F}}_{tot}(\mathit{r},x)=-zeq(\mathit{r},x)\mathrm{\partial }\psi (\mathit{r},x)/\mathrm{\partial }x$ is the electrostatic driving force.

Within the linear response approach, the fluxes are proportional to the generalized forces and the transport coefficient can be derived by the relaxation of equilibrium fluctuations, see Ref [94]. Within these approximations, the transport along the channel axis is captured by the following fluxes

(9) $j^{+}-j^{-}=J_q\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{c}\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t}$(9)

(10) $j^{+}+j^{-}=J_c\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{t}\mathrm{e}\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t}$(10)

(11) $Q\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t},$(11)

which can be induced by different forcing. Clearly, an electric current $J_q$ can be generated by an external electric field while a solvent current $Q$ by a pressure drop. However, also cross phenomena may appear where, e.g. an electric current can be prompted by applying a pressure drop on the solvent. In order to capture these non-trivial couplings it is insightful to write down the linear system of equations governing the relevant currents in the following way:

(12) $\left(\begin{array}{c}{J}_q\\ {J{ }^{\mathrm{\prime }}}_c\\ Q\end{array}\right)=\left(\begin{array}{ccc}{\mathcal{L}}_{11}& {\mathcal{L}}_{12}& {\mathcal{L}}_{13}\\ {\mathcal{L}}_{12}& {\mathcal{L}}_{22}& 0\\ {\mathcal{L}}_{13}& 0& {\mathcal{L}}_{33}\end{array}\right)\left(\begin{array}{c}ze\mathrm{\Delta }V\\ \mathrm{\Delta }\bar{\mu }\\ \mathrm{\Delta }P\end{array}\right)\frac{1}{L}$(12)

where ${\mathcal{L}}_{ij}$ are the coefficients of the Onsager matrix [95] associated with the drops across the channel of electrostatic potential $\mathrm{\Delta }V$, chemical potential $\mathrm{\Delta }\bar{\mu }$, and pressure $\mathrm{\Delta }P$ and where we have introduced $J_c^{\mathrm{\prime }}=J_c-Q$ in order to make the matrix symmetric. The symmetry of the Onsager matrix, that relies on the microscopic time reversibility of the underlying Hamiltonian dynamics [94], has a crucial practical outcome: the cross-terms, such as ${\mathcal{L}}_{12}$, can be computed in two ways: either measuring $J_q$ upon applying a chemical potential drop $\mathrm{\Delta }\bar{\mu }$ or by measuring the excess solute flow $J{{ }^{\mathrm{\prime }}}_c$ upon applying an electrostatic potential, $\mathrm{\Delta }V$. This approach is convenient in numerical simulations where applying chemical potential gradients requires to simulate large reservoirs (hence raising the computational time) whereas electrostatic forcing can be simulated with (computationally more convenient) periodic boundary conditions. Finally, we remark that ${\mathcal{L}}_{23}=0$ in Equation (12)(12) $\left(\begin{array}{c}{J}_q\\ {J{ }^{\mathrm{\prime }}}_c\\ Q\end{array}\right)=\left(\begin{array}{ccc}{\mathcal{L}}_{11}& {\mathcal{L}}_{12}& {\mathcal{L}}_{13}\\ {\mathcal{L}}_{12}& {\mathcal{L}}_{22}& 0\\ {\mathcal{L}}_{13}& 0& {\mathcal{L}}_{33}\end{array}\right)\left(\begin{array}{c}ze\mathrm{\Delta }V\\ \mathrm{\Delta }\bar{\mu }\\ \mathrm{\Delta }P\end{array}\right)\frac{1}{L}$(12) is due to the Debye-Hückel assumption [76] and that one expects ${\mathcal{L}}_{23}\ne 0$ for the case of the fully non-linear Poisson-Boltzmann equation.

In the case of smoothly varying channels, it is also possible to derive closed formulas for the coefficients elements of the Onsager matrix. Interestingly, the transport coefficients ${\mathcal{L}}_{ij}$ are particularly sensitive to the geometry and the conductive properties of the channel walls when the Debye length is comparable to the channel width. In this regime, one pair of off-diagonal Onsager matrix elements increases with the corrugation of the channel transport, in contrast to all other elements, which are either unaffected by or decrease with increasing corrugation [76].

2.2. Poisson-Nernst-Planck equations

A usual scenario is the one in which the flow of solvent is absent, $Q=\mathit{v}=0$. Indeed, in such a regime the (Navier-)Stokes equations are trivially fulfilled and Equation (1) simplify much reducing to the so-called Poisson-Nernst-Planck electrodiffusion equation in three-dimensional space (3d-PNP) governing the motion of positive and negative ions [96–98]:

${\mathit{j}}^{\pm }(\mathit{r},x)=-\beta D{\rho }^{\pm }(\mathit{r},x)\mathrm{\nabla }{\mu }^{\pm }(\mathit{r},x)$

(13) $=-D\left[\mathrm{\nabla }{\rho }^{\pm }(\mathit{r},x)\pm {\rho }^{\pm }(\mathit{r},x)\beta ze\mathrm{\nabla }\psi (\mathit{r},x)\right]$(13)

where $\psi$ is the solution of the Poisson Equation (4)(4) $ϵ{\mathrm{\nabla }}^2\psi \left(\mathit{r},x,t\right)=-ze\left[{\rho }^{+}\left(\mathit{r},x\right)-{\rho }^{-}\left(\mathit{r},x\right)\right]$(4) with the Boltzmann assumption, Equation (6)(6) ${\rho }^{\pm }(\mathit{r},x)={\bar{\rho }}^{\pm }(x)e^{\mp ze\psi (\mathit{r},x)}$(6) . At first, a one-dimensional reduced model was used to describe ion permeation (1d-PNP) [99], but eventually, treatments based on the full three-dimensional (3d-PNP) theory became possible [96,97]. In the absence of ion flux, the 3d-PNP theory reduces to the standard non-linear equilibrium Poisson-Boltzmann equation. As mentioned in the previous section, the 3d-PNP equation is based on a mean-field approximation that overlooks non-electrostatic ion–ion interactions (e.g. excluded volume) as well as polarization and ion correlation effects. This is why this equation is not well suited to treat highly charged molecules (like DNA) at high ion concentration or single file diffusion along a narrow pore. However, due to the low charge density of membrane systems, the PB equation is often used to compute the electrostatic potential profile along the pore of ion channels. The equation, allowing the comparison of relative stabilities of protonated and unprotonated residues, is also used for pKa calculations.

In practice, PNP is based on several simplifications: rigid channel structure, structureless dielectric solvent, and mean-field ion–ion interactions, which are of unknown validity in the context in which they are used. If one is to adopt a continuum electrodiffusion approach, such simplifications are necessary in order to have partial differential equations that can be solved numerically. The accuracy of the 3d-PNP approach for ion channels was discussed in a few studies [97,100,101], where 3d-PNP simulations were contrasted with Brownian dynamics. Results showed that 3d-PNP approach can predict the reversal potential for wide pores, but breaks down in narrow ion channels with radii smaller than the Debye length. Further extension of 3d-PNP to finite ion sizes [102], multiple ion species [103], and including dielectric repulsion terms [104] have been proposed, as well as a combination of one-dimensional (1d) PNP with classical density functional theory (see Section 2.5) [105]. However, depending on the situation, 3d-PNP may, or may not, be sufficiently accurate. Ultimately, the significance of that picture should not be expected to exceed that of the physical approximations upon which it is built.

2.3. Capture process of suspended particles

In many biologically and technologically relevant applications, such as drug delivery of molecular therapies, nanopore sensing, testing the impact of anesthetics on ion channels, etc., it is fundamental to characterize the capture of small and large molecules (analytes) from their bulk solution to narrow paths (channels and nanopores). An efficient capture, for instance, is crucial to obtain a correct behaviour of nanopore-sensing devices [106,107] which work detecting the ion-current variation produced by the passage of analytes into channels [108–111]. In this case, Equation (1) should be complemented with the equations governing the evolution of the density of analytes. Clearly, the density of analytes is coupled to that of the electrolyte and of the solute, hence the Onsager matrix will become a $4\times 4$ matrix including novel diagonal and cross terms.

If the analytes are highly diluted and weakly charged (or highly screened by the salt solution), their presence is assumed to mildly affect the electrostatics and the velocity profile, such that, they can be modelled as passive tracers. Accordingly, Equation (1) are solved and then the velocity profile and the electrostatic potential are used as external fields for the evolution of the density of analytes.

In the context of capture mechanism, Equation (1) are crucial to characterize and control the main electrokinetic effects that are known to drive the analytes to the capture regions [112–114]: electrophoresis (EP), dielectrophoresis (DEP), and electroosmosis (EO). EP is the motion of charged molecules relative to a solvent induced by a uniform electric field. DEP is the transport of neutral particles carrying a permanent or induced dipole in inhomogeneous electric fields, in which dipoles not only undergo a torque, but also a translational force. EO is the flow of an electroneutral solvent (e.g. water) through a membrane triggered by the motion of ions in electric fields; this flow can drag even large macromolecules [111]. The deep entanglement and competition of these three effects can also be seen from the structure of the Onsager matrix in Equation (14).

The continuum approach so far described allows deriving the dependence of the capture rates on the above effects. In the presence of an external electric field, the capture problem can be formulated in terms of a transport equation [115] for a diluted analyte concentration ${\rho }_p(r,t)$ that is a function of the distance $r$ from the center of the nanochannel entrance

(14a) $\frac{\mathrm{\partial }{\rho }_p(r,t)}{\mathrm{\partial }t}=-\mathrm{\nabla }\cdot \mathit{j}(r,t)$(14a)

(14b) $\mathit{j}(r,t)=-D_p\mathrm{\nabla }{\rho }_p(r,t)+\mathit{v}{\rho }_p(r,t)+\beta D_p\mathit{F}{\rho }_p(r,t)$(14b)

where the molecule flux $\mathit{j}$, pointing to the pore, has three components: i) diffusive, ii) advective due to the solvent motion, $\mathit{v}$, for example, induced by electroosmosis, and iii) phoretic due to external fields. $D_p$ denotes the molecule diffusivity. As compared to Equation (1), the presence of electrolyte solution manifests only through an EO flux and the dielectric constant.

It is important to notice that Equation (1) focus directly on the dynamics of the electrolytes and of the solution while in Equation (14a)(14a) $\frac{\mathrm{\partial }{\rho }_p(r,t)}{\mathrm{\partial }t}=-\mathrm{\nabla }\cdot \mathit{j}(r,t)$(14a) -(14b(14b) $\mathit{j}(r,t)=-D_p\mathrm{\nabla }{\rho }_p(r,t)+\mathit{v}{\rho }_p(r,t)+\beta D_p\mathit{F}{\rho }_p(r,t)$(14b) ) the electrolyte solution is approximated just as a background fluid, with an assigned dielectric constant, where macromolecules are transported by external fields and flows. This justifies the absence of Poisson and Navier-Stokes equations.

According to the hemi-spherical electrode approximation [116], the external voltage generates, far from the electrodes, a radial electric field of modulus $E(r)=-I/(2\pi r^2)$ directed towards the pore center. $I$ is the ion current due to the applied voltage. Thus, the problem Equation (14a)(14a) $\frac{\mathrm{\partial }{\rho }_p(r,t)}{\mathrm{\partial }t}=-\mathrm{\nabla }\cdot \mathit{j}(r,t)$(14a) -(14b(14b) $\mathit{j}(r,t)=-D_p\mathrm{\nabla }{\rho }_p(r,t)+\mathit{v}{\rho }_p(r,t)+\beta D_p\mathit{F}{\rho }_p(r,t)$(14b) ) acquires a radial symmetry with respect to the center of the pore entrance (origin).

The stationary solution ${\rho }_p(r)$ of Equation. (14a)(14a) $\frac{\mathrm{\partial }{\rho }_p(r,t)}{\mathrm{\partial }t}=-\mathrm{\nabla }\cdot \mathit{j}(r,t)$(14a) -(14b(14b) $\mathit{j}(r,t)=-D_p\mathrm{\nabla }{\rho }_p(r,t)+\mathit{v}{\rho }_p(r,t)+\beta D_p\mathit{F}{\rho }_p(r,t)$(14b) ) in the shell between $r=r_e$ (entrance radius) and $r\to \mathrm{\infty }$ is determined by the boundary conditions, ${\rho }_p(\mathrm{\infty })={\rho }_{p,\mathrm{\infty }}$, meaning that the molecule concentration is constant far from the channel, and $j(r_e)=-k{\rho }_p(r_e)$, indicating that the channel entrance is partially absorbing, that is, only a fraction of the incoming molecules actually enter.

Once ${\rho }_p(r)$ is known, the flux density $j(r)$ can be easily derived via Equation (14b)(14b) $\mathit{j}(r,t)=-D_p\mathrm{\nabla }{\rho }_p(r,t)+\mathit{v}{\rho }_p(r,t)+\beta D_p\mathit{F}{\rho }_p(r,t)$(14b) and the resulting capture frequency $\nu$ (number of molecules adsorbed per unit time) is nothing but $2\pi r_e^{2}j(r=r_e)$,

(15) $\nu =\frac{2\pi {\rho }_{p,\mathrm{\infty }}D}{\frac{D}{k{r}_e^2}{e}^{\varphi (r_e)}+{\int }_{r_e}^{\mathrm{\infty }}d\rho \frac{e^{\varphi (\rho )}}{\rho }}$(15)

where $\varphi (r)=\mu v(r)/D$ is the dimensionless effective potential, which accounts for the radial EP and DEP effects as well as for the contribution of the EO flow [115]. Equation (15)(15) $\nu =\frac{2\pi {\rho }_{p,\mathrm{\infty }}D}{\frac{D}{k{r}_e^2}{e}^{\varphi (r_e)}+{\int }_{r_e}^{\mathrm{\infty }}d\rho \frac{e^{\varphi (\rho )}}{\rho }}$(15) , written in the form $1/\nu =1/{\nu }_a+1/{\nu }_e$, shows that the capture frequency splits in two contributions, ${\nu }_a$, accounting for the approach to the pore, and ${\nu }_e$, associated with the molecules actually entering it. The entrance frequency reads ${\nu }_e=2\pi {\rho }_{p,\mathrm{\infty }}k{r}_e^{2}exp[-\varphi (r_e)]$, and depends on the parameter $k$, whose precise determination requires an atomistic description of the system. Nevertheless, in the same approximation of Equation (14a)(14a) $\frac{\mathrm{\partial }{\rho }_p(r,t)}{\mathrm{\partial }t}=-\mathrm{\nabla }\cdot \mathit{j}(r,t)$(14a) -(14b(14b) $\mathit{j}(r,t)=-D_p\mathrm{\nabla }{\rho }_p(r,t)+\mathit{v}{\rho }_p(r,t)+\beta D_p\mathit{F}{\rho }_p(r,t)$(14b) ), we can attempt an Eyring equation for the rate

(16) $k=\kappa \frac{k_{B}T}{h}exp\left(-\frac{\mathrm{\Delta }\mathcal{F}}{k_{B}T}\right),$(16)

where $\kappa$ is the transmission coefficient and $h$ the Planck constant. $\mathrm{\Delta }\mathcal{F}$ indicates the free-energy barrier, generally positive, which molecules need to cross to enter the pore, see Section 4 for methods to compute free-energy barriers of ion permeation and Section 6 for examples of atomistic calculations of free-energy barriers. High values of $\mathrm{\Delta }\mathcal{F}$ reduce exponentially the adsorption rate.

The simple expression (15) shows how continuum approaches are useful to derive estimates of the capture frequency that can be used to establish the order of magnitude of the capture rates and their dependence on the experimental conditions. Moreover, the Smoluchowski theory allows an immediate assessment of the competition among electroosmosis and various phoretic contributions involved in experiments [117]. Moreover, the approach allows to determine the Mean First Passage time [118] and the permeability [119] (see Section 3.2) in corrugated channels under the action of a chemical potential difference at the channel ends.

2.4. Kirkwood approximation: ionic hydration patterns in nanopores

Transport through nanopores is affected not only by ion-pore and ion–ion interactions, but also by the hydration barrier originating from the re-distribution of the water molecules when an ion traverses the nanopore. In addition, water distributions are also necessary to assess water-mediated ion-pore interactions. To address this issue Barabash et al. [28] extended an approach designed to characterize the DNA hydration [120,121] to analytically predict the ionic hydration patterns inside nanopores, using the radial distribution functions (RDFs) evaluated in the bulk electrolyte. Using a rigorous statistical mechanical correspondence between the potential of mean force (PMF) and RDFs, the density ${\rho }_w(r)$ of water-oxygen atoms can be rewritten via the multi-particle distribution functions [122]. By truncating the PMF decomposition at the 2-particle terms, one arrives at the Kirkwood approximation [123] involving only a product of the RDFs

(17) $\frac{{\rho }_w(r)}{{\rho }_{w,0}}=g_{iw}(|\mathit{r}-{\mathit{r}}^i|)\prod _{j=1}^{N_p}{g}_{pw}(|\mathit{r}-{\mathit{r}}_j^p|).$(17)

Here, $\mathit{r}$ is the position of interest, ${\mathit{r}}_j^p$ are the coordinates of the fixed atoms of the nanopore ($N_p$ in total), ${\mathit{r}}^i$ is the location of the ion, ${\rho }_{w,0}$ indicates the bulk water density, $g_{iw}$ and $g_{pw}$ are the ion-oxygen and nanopore-oxygen free bulk RDFs, respectively.

Equation (17)(17) $\frac{{\rho }_w(r)}{{\rho }_{w,0}}=g_{iw}(|\mathit{r}-{\mathit{r}}^i|)\prod _{j=1}^{N_p}{g}_{pw}(|\mathit{r}-{\mathit{r}}_j^p|).$(17) shows that the complex water patterns originate from the interference between the ionic hydration shells ($g_{iw}$) and the hydration cloud around the nanopore (${\prod }^{g_{pw}}$). Since the method requires the bulk RDFs $g(r)$ to be measured only once under given system structure and composition, a 10²–10⁴ speedup is achieved as compared to enhanced sampling all-atom MD simulations [28,120]. The truncation of the decomposition to 2-particle terms overlooks the interactions between the components of individual water molecules and their orientations, leading to overestimated densities near the pore walls. Nevertheless, the proposed analytical description generally agrees well with the results of molecular dynamics simulations for K⁺, Na⁺, and Cl ^– ions (see Figure 2) and predicts the locations of the trapped water molecules [28]. Hence, the method provides fundamental insights into the electrostatics, the dielectric properties, and the dynamics of ions in nanopores. A tentative analytical way of connecting the hydration patterns and the single-ion equilibrium PMF is proposed as well [28].

Figure 2. Three-dimensional ionic hydration patterns. Water distribution around a K⁺ ion (purple sphere) at 0.4 nm above the pore in a graphene lattice (black ball-and-stick representation), obtained from (a) MD simulations and (b) theory, Equation (17)(17) $\frac{{\rho }_w(r)}{{\rho }_{w,0}}=g_{iw}(|\mathit{r}-{\mathit{r}}^i|)\prod _{j=1}^{N_p}{g}_{pw}(|\mathit{r}-{\mathit{r}}_j^p|).$(17) . The layered structure of the hydration around the ion and the graphene lattice emerges due to the chosen isovalue of 1.15. A 5-point smoothing window has been applied to the original data in both panels. Reproduced without changes from Ref [28]. under the Creative Commons CC BY 4.0 license.

Moreover, the approach allows to study the effects on the hydration patterns imposed by external strain (stretching, skewing, bending, twisting) and chemical features (pore isomers, type and charge of the rim atoms, numbers of lattice layers, layer offset eclipse). Such flexibility should prove useful in designing the hydration pattern via the multi-parameter optimisation leading to pre-defined properties of the nanopores.

2.5. Classical density functional theory and ion channels

Classical density functional theory (C-DFT) [124] is a powerful approach, in the grand canonical ensemble, that allows, through a variational principle, to compute the structure and the energetics of a system in equilibrium and the evolution of the density distribution for systems out of equilibrium [125–127]. The starting point for C-DFT is the functional of the grand potential that, for a one-component system, takes the form

(18) $\mathrm{\Omega }[\rho ]=\mathcal{F}[\rho ]+\int d^{3}r\rho (\mathbf{r})\left(V_{ext}(\mathbf{r})-\mu \right),$(18)

where $\mathcal{F}={\mathcal{F}}_{id}+{\mathcal{F}}_{ex}$ is the intrinsic Helmholtz free-energy functional, which can be split into an exactly known ideal gas part ${\mathcal{F}}_{id}$ and an excess (over the ideal gas) contribution ${\mathcal{F}}_{ex}$, which contains all the information about the particle–particle interactions [124]. In Equation (18)(18) $\mathrm{\Omega }[\rho ]=\mathcal{F}[\rho ]+\int d^{3}r\rho (\mathbf{r})\left(V_{ext}(\mathbf{r})-\mu \right),$(18) , $V_{ext}(\mathbf{r})$ is the external and $\mu$ the chemical potential. For practically all systems of interest, ${\mathcal{F}}_{ex}$ is unfortunately known only approximately. The equilibrium density profile ${\rho }_0(\mathbf{r})$ can be calculated from the fact that $\mathrm{\Omega }[\rho ]$ is minimal in equilibrium [124], which can be expressed as

(19) ${\left.\frac{\delta \mathrm{\Omega }}{\delta \rho (\mathbf{r})}\right|}_{\rho ={\rho }_0}=0.$(19)

Once the equilibrium density profile is known, the grand potential of the system follows immediately

(20) $\mathrm{\Omega }=\mathrm{\Omega }[{\rho }_0(\mathbf{t})].$(20)

In the context of ion channels, C-DFT can be viewed as complementary to molecular dynamics, since it provides configurations that are averaged over the grand-ensemble rather than single frames. This can be important when system states are characterized by rare events or large fluctuations. C-DFT allows the exploration of larger time scales than those typically accessible to MD, even if the implementation of force fields is less obvious.

There are several studies of ion channels within the framework of C-DFT so far, related to gating, selectivity, and ion transport. In gating, the relationship between the gate configuration and the probability of a bubble formation was studied in Refs. [44,128]. The gate configuration was accounted for by the external potential $V_{ext}(\mathbf{r})$ and the bubble formation followed from the density profile $\rho (\mathbf{r})$ of a water-like simple fluid. In order to illustrate the application of C-DFT to gating in an ion channel, we show in Figure 3(a,c) the simplified geometry of a channel and its change from a closed configuration in (a) to an open one in (c). This geometrical change of the channel geometry is employed as an external potential $V_{ext}(\mathbf{r})$ in the DFT calculation. A cut through the resulting density profiles of a closed and an open channel is also shown in (a) and (c) [128]. While in the open configuration the water density in the channel remains liquid-like (c), a clear bubble is formed in the narrow region of the closed configuration (a). In addition to the structure, one obtains information about the energetics of the channel. In this application, we used $R_2$, the radius of the gate, as a control parameter to describe the geometrical configuration. For each value of $R_2$, we can ask, on the one hand, what the structure of the equilibrium configuration of water in the channel or the gate is. On the other hand, we can ask what is the structure of an open and of a closed configuration for a given value of $R_2$, that is, a density distribution with a liquid-like density, say ${\rho }_{op}(\mathbf{r};R_2)$, and one with a bubble in the gate, say ${\rho }_{cl}(\mathbf{r};R_2)$. Once these density profiles are known, it is straightforward to compute the corresponding grand potentials ${\mathrm{\Omega }}_{op}=\mathrm{\Omega }[{\rho }_{op}(\mathbf{r};R_2)]$ and ${\mathrm{\Omega }}_{cl}=\mathrm{\Omega }[{\rho }_{cl}(\mathbf{r};R_2)]$ and their difference $\mathrm{\Delta }\mathrm{\Omega }={\mathrm{\Omega }}_{cl}-{\mathrm{\Omega }}_{op}$. If $\mathrm{\Delta }\mathrm{\Omega }>0$ then the open state of the gate is the equilibrium configuration and the closed one is meta-stable and vice versa. Assuming a two-state system with probability $P_{open}(R_2)$ for being in the open state and $P_{closed}(R_2)=1-P_{open}(R_2)$ for being in the closed state, we can easily calculate from $P_{open}(R_2)/P_{closed}(R_2)=exp(-\beta \mathrm{\Delta }\mathrm{\Omega })$ the probability of finding the channel open:

(21) $P_{open}(R_2)=\frac{1}{1+e^{-\beta \mathrm{\Delta }\mathrm{\Omega }}}.$(21)

Figure 3. (a) The geometry of a simple model channel, that enters a C-DFT calculation as external potential and a cut through the corresponding equilibrium density profile of a closed configuration of the gate, that is blocked by a bubble. (b) The probability of finding the channel open as function of the gate radius $R_2$. (c) The geometry and the corresponding equilibrium density profile of an open channel that is filled with a liquid-like water density. Reprinted (abstract/excerpt/figure) with permission from [128]. Copyright 2017 by the American Physical Society.

The shape of this probability as function of $R_2$ is shown in Figure 3(b).

Here, the accurate knowledge of the state point of water in relation to the bulk-phase diagram has proven beneficial to study the formation of a bubble and the probability of opening and closing of the gate. The influence of hydrophobic gases and hydrostatic pressure was studied due to the relatively small computational costs [44]. Within this framework, a more complete model of an ion channel was formed that included voltage sensors and a hydrophobic gate [129].

Selectivity can be driven by specific (chemical) interactions or by physical mechanisms, such as a competition between the electrostatic interactions and entropy. In the latter case, C-DFT can help to understand the physics involved in selectivity by providing structure and energetics of the process [130–132]. Coupling classical dynamical DFT (C-DDFT)[125,127] and the electrochemical gradient across a membrane allows one to study the characteristics of ion flux through channels [105] adding microscopic details to the more macroscopic approaches described in the previous sections.

3. Continuum approaches: applications

Continuum methods, despite their limitations, lend themselves to a variety of applications owing to their flexibility. A complete survey of the possible applications to ion channels is impossible within the extent of this review, thus we will limit the discussion to a few selected examples without any claim to be exhaustive.

3.1. Capture models

The continuum model for capture (Equation (14a)(14a) $\frac{\mathrm{\partial }{\rho }_p(r,t)}{\mathrm{\partial }t}=-\mathrm{\nabla }\cdot \mathit{j}(r,t)$(14a) -(14b(14b) $\mathit{j}(r,t)=-D_p\mathrm{\nabla }{\rho }_p(r,t)+\mathit{v}{\rho }_p(r,t)+\beta D_p\mathit{F}{\rho }_p(r,t)$(14b) )) can be used to interpret some experiments in synthetic nanopores, for example, those of Larkin et al. [117] in which two globular proteins, ProK and RNase A, are captured by a solid-state hafnium dioxide (HfO₂) pore immersed in a 1 M KCl solution at pH 8.1 at which the HfO₂ pore is slightly negatively charged. The surface charge induces an EO flow, which superimposes to the EP of the external field. Equation (15)(15) $\nu =\frac{2\pi {\rho }_{p,\mathrm{\infty }}D}{\frac{D}{k{r}_e^2}{e}^{\varphi (r_e)}+{\int }_{r_e}^{\mathrm{\infty }}d\rho \frac{e^{\varphi (\rho )}}{\rho }}$(15) can be used to estimate the capture frequency of the two proteins in a wide range of voltages.

Figure 4, left panel shows the capture frequency when the entrance contribution was neglected, i.e. ${\nu }_e=0$ and $\nu ={\nu }_a$, for protein concentration ${\rho }_{p,\mathrm{\infty }}=1$ nM. The points refer to Equation (15)(15) $\nu =\frac{2\pi {\rho }_{p,\mathrm{\infty }}D}{\frac{D}{k{r}_e^2}{e}^{\varphi (r_e)}+{\int }_{r_e}^{\mathrm{\infty }}d\rho \frac{e^{\varphi (\rho )}}{\rho }}$(15) that includes EO, EP, and DEP, whereas lines indicate the case in which DEP is set to zero. The small displacement of symbols with respect to the solid lines suggests that, in the experiment, dielectrophoresis is negligible. Moreover, as the pore is negatively charged, for positive $\mathrm{\Delta }V$, EO is directed towards the pore so that it cooperates with EP. Theoretical capture frequencies are 4–5 times smaller than the ones observed in the experiments (see Figure 3 of Ref [117]), which is remarkable given that no fitting parameters were introduced.

Figure 4. Left) Capture frequency for the experiment of Larkin et al. [117]. Circles and squares refer to the theoretical capture frequency calculated via Equation (15)(15) $\nu =\frac{2\pi {\rho }_{p,\mathrm{\infty }}D}{\frac{D}{k{r}_e^2}{e}^{\varphi (r_e)}+{\int }_{r_e}^{\mathrm{\infty }}d\rho \frac{e^{\varphi (\rho )}}{\rho }}$(15) for RNase A and Protein K. Lines refer to the capture frequency when DEP is neglected. In both cases, protein bulk concentration is ${\rho }_{p,\mathrm{\infty }}=1$ nM and complete adsorption is assumed (${\tau }_e$ = 0 and $\nu ={\nu }_a$, transport limited regime). The pore is slightly negatively charged at the experimental pH = 8.1 while both proteins are positive, i.e. EO and EP cooperate at $\mathrm{\Delta }V>0$. Protein sketches are created using VMD [7]; blue and red colors identify to positively and negatively charged residues, respectively. Right) Experimentally observed capture frequency $f$, black points. The V-shaped plot indicates that capture occurs at both positive and negative voltages. The dashed line refers to the theoretical capture frequency estimated via Equation (15)(15) $\nu =\frac{2\pi {\rho }_{p,\mathrm{\infty }}D}{\frac{D}{k{r}_e^2}{e}^{\varphi (r_e)}+{\int }_{r_e}^{\mathrm{\infty }}d\rho \frac{e^{\varphi (\rho )}}{\rho }}$(15) with $\mathrm{\Delta }{\mathcal{F}}_0=22k_{B}T$ in the Arrhenius formula for the absorbing rate $k$ (see text), reproducing the typical V-shaped behaviour. Figures are adapted with permission from [115]. Copyright 2020 American Chemical Society.

The analytic capture model, based on the continuum equations of Section 2.3, developed by Chinappi et al. [115] also allows the interpretation of an experiment of DEP capture into an hemolysin channel (aHL) of a $\beta$-hairpin engineered to carry a permanent dipole, while remaining electrically neutral at pH 7. In practice, the peptide is a dumbbell of length $d=4$ nm, with zero global charge, and dipole $p=8e\cdot$nm. Figure 4, right panel shows the capture frequency into the aHL obtained from the ion current traces, highlighting that capture occurs independently of the voltage polarity. The figure shows that the experimental data are in agreement with the predictions of DEP capture model.

3.2. Channel permeability

Another application of the presented continuum theory, in particular, of Smoluchowski equation, deals with the computation of channel permeability when local inhomogeneity and chemical heterogeneity of the channel makes the assumptions of constant diffusion coefficient not viable. This is the typical case for ion channels. The passive transport of ions/molecules through a pore in the limit of dilute molecule concentration (as in the previous discussion) can be written in the simplest form as the Fick’s law (see Section 2.1 and Ref [94].):

(22) $j=-\mathcal{P}\mathrm{\Delta }\rho$(22)

where $\rho$ is either the ion or molecule density and $\mathcal{P}$ is the permeability coefficient containing all geometrical and chemical channel properties but the dependence from the concentration. An expression for the permeability coefficient can be derived starting again from the Smoluchowski diffusion-drift equation with respect to the axial transport coordinate, the position of the particle along the pore axis $x$, and the Kramers relations [119,133]

(23) $\mathcal{P}=1/{\left[{\int }_0^L\frac{exp\beta U\left(x\right)}{D\left(x\right)}dx\right]}^{-1}$(23)

where $L$ is the length of the pore, $D(x)$ the local diffusion constant, $U(x)$ the PMF characterizing the ion/molecule effective interactions with the pore. This expression is based on the assumption that the permanent ions/molecules crossing the membrane are uncorrelated and that there is no saturation of the occupancy of the pore. This assumption is reasonable at low concentrations, but may be violated if the pore is narrow and there are high-affinity binding sites along the pore.

Equation (23)(23) $\mathcal{P}=1/{\left[{\int }_0^L\frac{exp\beta U\left(x\right)}{D\left(x\right)}dx\right]}^{-1}$(23) shows two important features of transport inside a nanopore: its non-locality, due to the presence of an integral over the pore length $L$, and its statistical nature, due to the presence of the potential of mean force $U(x)$. This means that any local method, such as docking, would predict permeation only in a few limited cases, such as a single binding site. Furthermore, one needs to thermodynamically average all other variables for a correct estimation of the $U(x)$ profile. Reliable methods to evaluate $U(x)$ to be used in Equation (23)(23) $\mathcal{P}=1/{\left[{\int }_0^L\frac{exp\beta U\left(x\right)}{D\left(x\right)}dx\right]}^{-1}$(23) are based on advanced sampling techniques combined with classical molecular dynamics (see Section 4.2). Equation (23)(23) $\mathcal{P}=1/{\left[{\int }_0^L\frac{exp\beta U\left(x\right)}{D\left(x\right)}dx\right]}^{-1}$(23) is valid only at dilute concentrations or far from saturation, presumably with the pore occupied by a single independent particle. In case of saturation conditions, there are additional methods that can be applied to correctly estimate the flux, such as the Markov state model [134,135].

While for ions diffusing through specific channels the approximations of Equation (23)(23) $\mathcal{P}=1/{\left[{\int }_0^L\frac{exp\beta U\left(x\right)}{D\left(x\right)}dx\right]}^{-1}$(23) are quite severe due to strong correlations among ions and multiple occupancy, the method works well for small molecules diffusing through channels (weak molecule–pore interactions, single occupancy). Electrophysiology data provided a good test for Equation (23)(23) $\mathcal{P}=1/{\left[{\int }_0^L\frac{exp\beta U\left(x\right)}{D\left(x\right)}dx\right]}^{-1}$(23) in the case of a charged particle, the beta-lactam inhibitor avibactam. Its permeability through the general channel OmpF was measured with the reversal potential method and compared with that calculated via Equation (23)(23) $\mathcal{P}=1/{\left[{\int }_0^L\frac{exp\beta U\left(x\right)}{D\left(x\right)}dx\right]}^{-1}$(23) . Metadynamics was used to quantify $U(x)$, and the agreement was remarkable [136].

3.3. Hydration patterns

The analytical approach to characterise the hydration patterns in nanopores [28] and based on Equation (17)(17) $\frac{{\rho }_w(r)}{{\rho }_{w,0}}=g_{iw}(|\mathit{r}-{\mathit{r}}^i|)\prod _{j=1}^{N_p}{g}_{pw}(|\mathit{r}-{\mathit{r}}_j^p|).$(17) opens an avenue for a number of potential applications. For example, one can hypothesise that the ionic transport rate and the permeation pathway can be tailored by iteratively choosing the pore geometry and the level of externally applied strain. First, such engineering of the selective and conductive properties appears of high demand in water desalination technologies, where the maximisation of the output water quality and energy efficiency have to meet the environmental sustainability criteria [137]. Secondly, the maximised energy-yield during blue energy harvesting has to exceed 5 Watt/m² to be economically viable while minimising Joule heating [138,139]. Thirdly, next-generation fast DNA sequencing would benefit by applying the same logic of optimisation to slow down and guide the DNA during its translocation through a nanopore [140,141]. In that regard, the classical Coulter-principle-based nucleobase detection could be improved by maximising the current blockage optimising the geometrical features of the nanopore. Fourthly, the Materials Genome Initiative would benefit from cataloguing 2D and 3D nanopore isomers [142]. Finally, transverse (corrugated) waves, travelling along a carbon nanotube (CNT), would invoke the ionic motion against the applied electrochemical gradient, thus realising a nanoscale ionic pump [143]. Thus, the method supports the design and optimization of controllable nanoionic devices with on-demand selective and conductive properties, finding its applications in biophysics, industry, and nanotechnology.

4. Atomistic approaches

4.1. Overview of the methods

Intrinsic to the continuum approaches introduced in Section 2 is the mean field approximation, which neglects all fluctuations and correlations which might be relevant at the scale of ion channels. In particular, ion–ion correlations beyond the mean-field approach of Poisson-Boltzmann (see Section 2.1) are particularly important in narrow channels: a notable example is the selectivity filter of potassium channels discussed in the Introduction. Although, in principle, the geometrically and chemically complex molecular structure of proteins may be accounted for in a continuum model, it is generally preferred to directly resort to particle-based approaches, which more naturally encompass both fluctuations and the atomic structure. The modelling at the atomistic level of ion channels is the most refined available description, providing all the biochemical details necessary to understand complex phenomena such as gating [144], hydrophobic control [145], effect of mutations [146], protein–membrane interaction [147], and the pivotal role of water molecules through H-bonds network formation and local dielectric screening [145].

Atomistic methods explicitly including all atoms and all inter-atomic interactions allow a spatial resolution of a single hydrogen atom and the time-resolution of a single femto-second, so that Molecular Dynamics (MD) simulations [148] can be considered as an atomic resolution microscope [77,149,150]. Of particular interest is the capability of MD to bridge structural information on ion channels with their dynamics and, with some challenges, to their function, with spatial and temporal resolution which is yet unmatched in experimental approaches [151]. In the following, we describe the conceptual framework of MD, showing the equations of motion, the definition of the force field, which prescribes the interactions between atoms, and the challenges specific to ion channel simulations.

Three main ingredients are contributing to the success of MD as extremely refined and powerful tool to investigate conformational structure and changes and the relationship between structure and function in proteins. First, the development of highly parallelized and optimized hardware and software tools allows the simulation of systems of increasing dimensions and with increasing trajectory lengths. The second ingredient is the development and refinement of dedicated force fields, including not only amino acid groups but also lipids, solvent molecules, and specific substances as glycans, whose fundamental role has emerged clearly, for example, in the case of the glycosylate spike Covid-19 protein [152]. Last, but not least, the development and improvement of a wealth of methods and algorithms to improve the statistical sampling of the conformational space (see Section 4.2).

The mathematical framework of MD is conceptually simple: for each of the atoms (the same holds in the case of coarse grained models), MD trajectories are generated by integrating Newton’s equations of motion²:

(24) $m_i{\mathit{a}}_i(t)=-\frac{\mathrm{\partial }V(\mathit{x}(t))}{\mathrm{\partial }{\mathit{x}}_i(t)}$(24)

relating the mass and acceleration of the $i$-th particle to the force acting on it, with $\mathit{x}$ the vector of the $3N$ cartesian atomic coordinates and $V$ a classical empirical potential function. In this representation, interactions are established between atom nuclei, while the electronic degrees of freedom are averaged out. A typical analytical expression of the interaction potential $V$ reads

$V=\sum _i^{\mathrm{b}\mathrm{o}\mathrm{n}\mathrm{d}\mathrm{s}}\frac{k_{l,i}}{2}(l_i-l_{0,i}{)}^2+\sum _i^{\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{l}\mathrm{e}\mathrm{s}}\frac{k_{\alpha ,i}}{2}({\alpha }_i-{\alpha }_{0,i}{)}^2$

$+\sum _i^{\mathrm{t}\mathrm{o}\mathrm{r}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}}\left\{\sum _k^M\frac{V_{ik}}{2}[1+cos(n_{ik}{\mathit{\theta }}_{ik}-{\mathit{\theta }}_{0,ik})]\right\}$

(25) $+\sum _{i,j}^{\mathrm{p}\mathrm{a}\mathrm{i}\mathrm{r}\mathrm{s}}{\epsilon }_{ij}\left[{\left(\frac{r_{0,ij}}{r_{ij}}\right)}^{12}-2{\left(\frac{r_{0,ij}}{r_{ij}}\right)}^6\right]+\sum _{i,j}^{\mathrm{p}\mathrm{a}\mathrm{i}\mathrm{r}\mathrm{s}}\frac{q_i{q}_j}{4\pi {ϵ}_0{ϵ}_r{r}_{ij}}$(25)

where the first four terms represent bonded interactions due to covalent bonds, which enforce chain connectivity as well as the secondary structure, while the latter two, the Lennard-Jones and Coulomb potentials, quantify non-bonded interactions of van der Waals and electrostatic nature. In particular, the first two terms are harmonic potentials allowing only small oscillations around the specified equilibrium values ($l_{0,i}$, ${\alpha }_{0,i}$) of bond lengths and bond angles, respectively. The intensity of these restraining potentials is quantified by the force constants $k_{l,i}$ and $k_{\alpha ,i}$. In the third term, the torsional potential, $V_{ik}$ and $n_{ik}$ represent the barrier height and the multiplicity, respectively, while ${\mathit{\theta }}_{ik}$ and ${\mathit{\theta }}_{0,ik}$ are the current and equilibrium value of the torsional angles. In the Lennard-Jones potential, ${\epsilon }_{ij}$ represents the depth of the energy well whose minimum is located at $r_{ij}=r_{0,ij}$. Finally, in the Coulomb term, $q_i$ and $q_j$ are the atom charges, while ${ϵ}_r$ and ${ϵ}_0$ are the local and vacuum dielectric constants, respectively.

The sets of parameters of the empirical potential in Equation (25)(25) $+\sum _{i,j}^{\mathrm{p}\mathrm{a}\mathrm{i}\mathrm{r}\mathrm{s}}{\epsilon }_{ij}\left[{\left(\frac{r_{0,ij}}{r_{ij}}\right)}^{12}-2{\left(\frac{r_{0,ij}}{r_{ij}}\right)}^6\right]+\sum _{i,j}^{\mathrm{p}\mathrm{a}\mathrm{i}\mathrm{r}\mathrm{s}}\frac{q_i{q}_j}{4\pi {ϵ}_0{ϵ}_r{r}_{ij}}$(25) are referred to as force field. Some of them are computed through quantum simulations of small organic molecules exploiting the concept of transferability, while others are estimated so as to match thermodynamic measurements. The most popular force fields for biomolecular simulation include AMBER [153], CHARMM [154], OPLS [155], and GROMOS [156] and, along with amino acids, they also include parameters for nucleic acids, lipids, carbohydrates, and ions. The parameterization of small organic compounds, for example, molecules of pharmaceutical interest, is more challenging, but the task has been simplified thanks to generalized force fields (GAFF [157], CGenFF [158]) and toolkits like Antechamber [159], Paramfit [160], and GAAMP [161] with global optimization of parameters [162]. MD atomistic force fields are continuously revised and assessed [163]. Even water models are of pivotal relevance for a proper description of mechanisms based on hydrophobicity [164].

The integration of Equation (24)(24) $m_i{\mathit{a}}_i(t)=-\frac{\mathrm{\partial }V(\mathit{x}(t))}{\mathrm{\partial }{\mathit{x}}_i(t)}$(24) can be performed with a variety of schemes [148]. For instance, the velocity Verlet algorithm updates positions and velocities according to:

(26a) ${\mathit{x}}_i(t+\delta t)\approx {\mathit{x}}_i(t)+{\mathit{v}}_i(t)\delta t+\frac{1}{2}{\mathit{a}}_i(t)\delta t^2$(26a)

(26b) ${\mathit{v}}_i(t+\delta t)\approx {\mathit{v}}_i(t)+\frac{1}{2}[{\mathit{a}}_i(t)+{\mathit{a}}_i(t+\delta t)]\delta t.$(26b)

The timestep $\delta t$ plays a key role, since, to keep the integration stable, it has to be tuned to the period of the fastest oscillating modes in the system. Typically, in biological system, these are associated with the covalent bonds, with $\delta t$ of $1$ or $2$ femtoseconds. Fixing the most rapidly vibrating bonds [165], or reassigning the mass of hydrogens to the associated heavy atoms [166] can allow for $\delta t$ up to 5 femtoseconds. $\delta t$ thus limits the achievable length of the simulated trajectories of biological systems which, currently, is of the order of few microseconds on high-performance supercomputers.

In order to evolve Equation (24)(24) $m_i{\mathit{a}}_i(t)=-\frac{\mathrm{\partial }V(\mathit{x}(t))}{\mathrm{\partial }{\mathit{x}}_i(t)}$(24) , initial conditions for the atomic positions and velocities need to be specified. For the latter, a Maxwell-Boltzmann distribution is typically used. For the atomic coordinates and the additional information on the protein topology needed to setup the bonded potentials in Equation (25)(25) $+\sum _{i,j}^{\mathrm{p}\mathrm{a}\mathrm{i}\mathrm{r}\mathrm{s}}{\epsilon }_{ij}\left[{\left(\frac{r_{0,ij}}{r_{ij}}\right)}^{12}-2{\left(\frac{r_{0,ij}}{r_{ij}}\right)}^6\right]+\sum _{i,j}^{\mathrm{p}\mathrm{a}\mathrm{i}\mathrm{r}\mathrm{s}}\frac{q_i{q}_j}{4\pi {ϵ}_0{ϵ}_r{r}_{ij}}$(25) , one typically resorts to structural data deposited in databases such as the Protein Data Bank [167]. The provided structures also include information on atom types, which is required to setup the non-bonded interactions. Typically, even when the entire sequence of the ion channel is resolved, structural data need to be complemented with appropriate information for (a portion of) the cellular membrane and the surrounding aqueous solution. Initial ‘equilibration’ runs are needed in order to relax the structure in simulated conditions as close to the physiological ones as possible. During these runs, hydration and relaxation of the structures to the locally stable structure are achieved. In principle, equilibrium simulations for timescales comparable to experimentally relevant ones (milliseconds) should yield information on the dynamics of the protein and on its function, for example, the gating process of an ion channel, starting solely from structural information. Such brute force approach, however, is computationally inapproachable, except using ad-hoc supercomputer architectures discussed in Section 8. In addition, atomistic trajectories spend most of the time in well-defined regions of the phase space, only rarely jumping to another basin, which makes brute-force sampling extremely inefficient; the problem of rare events and enhanced sampling techniques are discussed in Section 4.2.

A typical issue in ion channel simulation is the calculation of ion currents and conductances. The importance of performing all-atom simulations is related to the pivotal role that the details of the chemical structure play in determining the behaviour of ion channels, gated by voltage, ions concentration, or ligand binding. Even a point-like mutation, or a residue chain rotation can be able to significantly change the channel conductance, by tuning or even impeding ion passage. In simulations, ions can be driven through the channel with a constant electric field $E$ perpendicular to the membrane plane that generates a membrane potential across the channel driving ion movement (as in Figure 5, left panel). In this method, the potential applied to the membrane is $V=E{L}_z$, where $E$ is the constant field and $L_z$ is the length of the simulation box in the direction of the membrane normal [168,169]. The constant electric field acts on all charged particles throughout the simulated system and affects their spatial distribution–especially in the bulk solution. This distribution results in a nontrivial transmembrane potential in the region of the pore that differs from $V$. The usage of periodic boundary conditions means that there is a single compartment because ions passing through the channel can be recycled avoiding charge build-up on either side of the membrane. However, it is also possible to impose appropriate conditions at the edge of the periodic box to maintain asymmetric ionic concentrations, allowing the simulation of the reversal potential yielding zero ionic current [170].

Figure 5. Two typical simulation setups in computational electrophysiology: applying an electrical field (left) and using two reservoirs with different concentrations. Adapted from [173], Copyright 2020, with permission from Elsevier.

An alternative method consists in simulating a system with two membranes, dividing the system into two compartments (as in Figure 5, right panel). By adjusting the ion concentrations in each compartment a net charge imbalance can be created directly to model the desired electrochemical gradient [171–173]. In this method, the potential applied to the membrane is $V=Q/C$, where $Q$ is the charge imbalance and $C$ is the membrane capacitance. In this case, ions must be swapped between compartments after passing through the channel to maintain the desired potential through the charge imbalance.

The membrane potential $V$ and the charge imbalance $Q$ are conjugated thermodynamic variables [168]. In the constant field method, $V$ is under control and $Q$ can fluctuate, whereas in the compartment method, $Q$ is under control and $V$ can fluctuate. Provided they are appropriately applied, using an external electric field or a charge imbalance can yield very similar results [174]. Although early simulations examined only one or two permeation events [12], increases in computational power now allow for many permeation events to be captured in multi-ms simulations, see for example [13,16,175]. Agreement with experimentally measured currents cannot be guaranteed, due to many reasons, including the state and condition in which the channel structure was obtained as well as limitations in the MD force fields themselves.

One limit of the atomic description is that its high level of detail currently precludes, for computational limitations, the simulation of mesoscopic biological systems like whole organelles, vesicles, or viral particles. For instance, the simulation of an action potential would require the explicit simulation of a large stretch of membrane hosting a population of many different types of ion channels. This system is way too large for classical atomistic force fields. One possible approach to the problem is Brownian Dynamics, which is, in a way, a hybrid between mean field and atomistic approaches, see Section 4.3. A different approach more akin to the atomistic one, is that of coarse-grained force fields where the number of interaction centers (and thus the number of interactions) is reduced by lumping several atoms in a single Coarse-Grained (CG) particle (see for example, [176–178]). A particularly successful CG tool for the study of protein/lipid systems is the MARTINI model [179] which combines the speed-up benefits of simplified models with the resolution obtained for atomically detailed models. Indeed, the MARTINI model allows the simulation of very large membrane systems, for example, an entire viral envelope with 29 millions of CG particles [180], allowing to access much larger time scales than atomistic simulations. Despite this high speed-up the MARTINI modeling retains a significant amount of molecular detail. In the MARTINI model on average four heavy atoms plus associated hydrogens are represented by a single interaction center. Mapping of water is consistent with this choice, with four real water molecules mapped to a CG water particle, while ions are modelled as a single particle, which represents both the ion and its first hydration shell. The MARTINI model has four main types of particles: polar, non-polar, apolar, and charged. Within each type, there are sub-types based on hydrogen-bonding capabilities or the degree of polarity, giving a total of 18 particle types or building blocks. These CG particles interact with one another through classical bonded potentials, Lennard-Jones interactions, and Coulomb interactions. The force field parameters have been tuned to reproduce a number of experimental thermodynamic data, in particular, the water/oil partitioning behaviour of a variety of compounds. The MARTINI model has been successfully used in a wide range of lipid-centered applications including the characterization of lipid membrane properties, lipid polymorphism, protein-lipid interplay and membrane protein oligomerization (see [179,181] and references therein), and lipid nanoemulsions [182]. The MARTINI model assumes rigid protein secondary structures. However, recent refinements extends its applicability to cases where large conformational changes are involved, for example, ion channel gating [183], where the Martini 3 force field was used [184] in combination with a G$\bar{\mathrm{o}}$-like model for proteins [185].

It should be noted that both all-atom and coarse-grained approaches encompass some degree of empiricity; for both methods this is true for the force fields which need to be tuned with some ‘first principle’ data. However, moving away from the atomic description inherently requires additional assumptions on the location, interactions, and dynamics of the CG particles. In other words, improvements in the accessible time and length scales come at the price of an additional level of empiricity, which should be assessed on a case-by-case basis.

4.2. Rare events in ion channels

As introduced in Section 4.1, the requirement of a small timestep in the integration of the equations of motion limits the length of the simulated trajectories to a few microseconds. This limitation is particularly problematic in the case of biomolecular systems, such as ion channels, where many important functional processes (such as gating or ligand binding) take place on much longer timescales.

If we associate the different functional states of a channel to distinct minima in an energy landscape, in a standard MD simulation, the trajectory spends most of the time close to the local minimum in which it was initialised, and, only seldom, jumps to a different local minimum owing to sufficiently large thermal fluctuations (see Figure 6). If the thermal energy of the system is lower than the energy barrier separating the minima, the activated process happens very infrequently. However, when it occurs, it is quite fast. Sampling such rare events is thus a major challenge for MD simulations, which motivated the development of several specific methods to speed up such events and enhance conformational sampling [80]. The crux is that, according to Eyring Equation (16)(16) $k=\kappa \frac{k_{B}T}{h}exp\left(-\frac{\mathrm{\Delta }\mathcal{F}}{k_{B}T}\right),$(16) , the waiting time $k^{-1}$ between rare events scales exponentially with the height of the free-energy barriers $\mathrm{\Delta }\mathcal{F}$ measured in thermal units $k_{B}T$. As a result, the trajectories spend most of the time close to minima, while a comparatively short one in the transition state (on top of the barrier). Because of the temporal limitations of MD, normally very few barrier-crossing events can be observed and the sampling of the transition state is very poor.

Figure 6. Left) Illustration of the problem of rare events in biophysical simulations via a cartoon of the protein free-energy landscape: during typical simulations only local barriers can be overcome, while a number of configurational changes or conduction events lies beyond the current computer capabilities. Reproduced from [186], under creative commons licence CC BY 4.0. Right) Conformational changes giving rise to opening and closing of the CRAC channel, see Ref [187].

In the specific case of ion channels, there are three main classes of thermally activated processes of paramount interest: i) structural changes, which switch between the functional states of the channel: open, closed, or inactivated (gating), ii) ion conduction, in which each ion has to overcome a free-energy barrier in order to translocate through the channel, and iii) binding of ligands or drugs to the channel. In the following, we will describe few enhanced sampling techniques developed to overcome the rare event problem and cite some examples on i) and ii). Problem iii), instead, is common to different kinds of proteins and is reviewed elsewhere [188].

Enhanced sampling techniques can be broadly classified into three categories [189]: (i) Techniques that do not focus on specific reaction coordinates; (ii) Techniques that aim to determine the path between known end points; and (iii) Techniques that enhance sampling or reconstruct the free energy along one or more selected reaction coordinates. While not truly an enhanced sampling technique, a final approach that has extensive use in ion channel simulations is that of alchemical perturbation discussed below. The actual implementation of enhanced sampling techniques is currently managed either by built-in modules to common MD codes or by specific packages, such as Plumed [190], which works together with some of the most popular MD engines, as Amber [191], NAMD [192], and Gromacs [193].

4.2.1. Techniques that do not focus on specific reaction coordinates

The first group includes techniques that do not focus on specific reaction coordinates (also referred to as Collective Variables, CVs) but simultaneously heat all the degrees of freedom of the system to increase the frequency of barrier crossing. These approaches have the advantage that they do not require predefined knowledge of the system to define reaction coordinates of interest.

The archetypal techniques of this group are Temperature and Hamiltonian Replica Exchange Molecular Dynamics (T-REMD and H-REMD) [81,194]. In T-REMD [81] a number of replicas are run in parallel at different temperatures, and exchanges are accepted or rejected according to a Metropolis-like criterion ($P_{acc}=min[1,exp(-\mathrm{\Delta }\beta \mathrm{\Delta }E)]$). A replica can thus easily evade from a kinetic trap when it is warmed up due to an exchange of atomic coordinates with a neighbouring replica. This method thus addresses the multiple-minima problem arising from the ruggedness of the free-energy landscape that threats the ergodicity of MD simulations. The most important limit of T-REMD simulations is their poor scaling with the size of the system, as the number of replicas is proportional to the square root of the number of atoms. This occurs because energy is an extensive property so that for large systems, the energy difference between neighbouring replicas will be large imposing a small temperature difference in order to have a sufficiently high overlap in energies of each replica to obtain a high enough acceptance probability. This, however, increases the number of replicas and thus the computational cost of the simulation. Heating only some parts of a system, such as the protein but not the solvent as in replica exchange with solute tempering (REST), can help alleviate this issue but can reduce the enhancement of conformational sampling [195].

This kind of problems can be avoided in Hamiltonian Replica Exchange [194] where replicas are simulated at the same temperature but with different energy functions. The different energy functions can be attained by splitting the energy function in different terms weighted by different scaling coefficients. The terms with unitary weights cancel out in the expression of the acceptance probability that therefore will depend only on the energy of a subset of degrees of freedom. Despite this improvement, H-REMD is computationally demanding for the simulation of membrane proteins and not frequently used in this field. Indeed, this technique has been mainly used for the simulation of peptides and small organic molecules [196–198].

An alternative to running multiple simulations is to alter the potential energy function of the entire system so as to speed up transitions as in accelerated MD (aMD) [199] or Gaussian accelerated MD (GaMD) [200]. In these, the potential energy of wells below a certain threshold value are increased, while those above are not. This speeds up transitions between energy wells. In theory, re-weighting of the simulation trajectories can be achieved to reproduce the unbiased results. Re-weighting is prone to noise, something that is reduced in GaMD.

4.2.2. Techniques to determine the path between known end points

The most known of these approaches is transition path sampling (TPS) [84]. In TPS, a statistical description of reactive trajectories is introduced via the definition of the transition path ensemble as the set of trajectories connecting the end states of a transition. Based on this description, reactive trajectories are directly obtained via specifically designed sampling procedures, and analyzed together to obtain information on the reactive process.

A key quantity in the approach is the reactive path probability, that is, the probability that a dynamical trajectory $\mathit{x}(\tau )$ started in a metastable set (A) at time zero will end up in a different metastable set (B) at a successive time $\tau$, defined as

(27) ${\mathcal{P}}_{AB}[\mathit{x}(\tau )]=h_A({\mathit{x}}_0)\mathcal{P}[x(\tau )]h_B({\mathit{x}}_{\tau })/{\mathcal{Z}}_{AB}(\tau ).$(27)

In (27) $\mathcal{P}[\mathit{x}(\tau )$] is the probability of the dynamical path, $h_A$ and $h_B$ are the characteristic functions of the corresponding regions (i.e. $h_M$ = 1 if $\mathit{x}\in M$ and $h_M$ = 0 otherwise), ${\mathit{x}}_0$ and ${\mathit{x}}_{\tau }$ are the states of the trajectory at times 0 and $\tau$, respectively, and ${\mathcal{Z}}_{AB}$ is a normalization factor.

Various deterministic and stochastic methods have been proposed over the years to sample this probability and, thus, the reactive path ensemble. The most commonly used rely on a Monte Carlo importance sampling procedure in the space of paths: starting from a given transition path, a new one is generated and accepted or rejected according to a Metropolis rule based on the reactive path probability. In practice, it is of relevance how the new paths are generated, since the acceptance rate, and thus the efficiency of a transition path sampling, will depend on how this is done. The most widely used is the so-called shooting algorithm, whose basic idea is to randomly pick a point along the old path, perturb the momenta of the system and propagate it forward and backward in time until it reaches the reactant or product state.

Another set of approaches to directly generate reactive trajectories is that of the string method and its variants, a class of algorithms based on the Transition Path Theory developed by E and Vanden-Eijnden [201]. The zero-temperature string method [202,203] allows to identify the minimum energy path (MEP) in the potential landscape $V(\mathit{x})$ of a system, that is, the curve $\varphi$ that connects two minima of $V(\mathit{x})$ via a saddle point and whose tangent is everywhere parallel to the force, that is

(28) $(\mathrm{\nabla }V{)}^{\mathrm{\perp }}(\varphi )=0.$(28)

An efficient procedure to obtain the MEP is to discretize an initial path (the string) into a number of points (called images) and evolve them according to steepest descent dynamics, while enforcing a chosen parametrization of the whole curve. Examples of commonly used parametrization schemes are by equal or energy-weighted arc-length.

A different algorithm, named finite-temperature string method [204,205], was designed to calculate the so-called principal curve, that is, the curve joining the mean positions of the system on the hyperplanes perpendicular to the curve itself. Because of its definition, the location of this curve depends on specific features of the underlying energy landscapes, and is thus appropriate to be used when the landscape is rugged, or variations in the width of the reaction channel are relevant. The curve can be obtained via successive iterations of i) constrained sampling in the orthogonal hyperplanes, as in the original formulation [204], or sampling in the Voronoi tessellation of the space generated by the discretized images along the curve itself [205], ii) update of the images to the new centers of the distributions, and iii) a reparametrization step.

While the zero- and the finite-temperature string methods were first formulated in the space of the Cartesian coordinates of the system $\mathit{x})$, versions using CVs have also been developed [85,205]. These allow to calculate the minimum free-energy path (MFEP), or the principal curve, in free-energy space, and permit the use of a large number of CVs, being thus useful for applications involving complex systems with many degrees of freedom such as ion channels.

Over the years, variants of the string method in CVs have been developed, such as the on-the-fly [206], the swarms of trajectories [207], or the dynamical string method [208]. From a theoretical perspective, these algorithms calculate different paths [201,209,210], and as a consequence, they also differ on the technique used to evolve the images. Concerning ion channels, applications of the string methods include the conformational transition of a potassium channel voltage sensor domain [211] and the gating of the bacterial GLIC channel [42,212].

4.2.3. Techniques that employ selected reaction coordinates

Perhaps, the most commonly applied category of enhanced sampling techniques in ion channel simulations are those that aim to enhance sampling or reconstruct the free-energy (or Potential of Mean Force, PMF) landscape of the system with respect to a selected set of CVs, defined as

(29) $\mathcal{W}(\mathit{z})=-k_{B}Tln\frac{\int e^{-\beta U(\mathit{x})}{\mathrm{\Pi }}_{i=1}^n\delta [{\theta }_i(\mathit{x})-z_i]{\mathrm{d}}^N\mathit{x}}{\int e^{-\beta U(\mathit{x})}{\mathrm{d}}^N\mathit{x}},$(29)

where $\{{\theta }_1(\mathit{x}),..,{\theta }_n(\mathit{x})\}$ are $n$ CVs and $\mathit{z}=\{z_1,..,z_n\}$ is a set of their values.

These techniques can be further classified into equilibrium (like Umbrella Sampling [213] or the Blue Moon Ensemble technique [214]), out-of-equilibrium (like Steered Molecular Dynamics [215] and Targeted Molecular Dynamics [216]), and close-to-equilibrium techniques (where the bias is small and incremental), like Metadynamics [217] or Adaptive Biasing Force [218] to cite but a few. In the last two techniques, the biasing potential is history-dependent, and favours the exploration of previously unvisited regions of the conformational space.

Good CVs should be low-dimensional descriptors of the functional dynamics of the system (often captured by the slow modes) [86]. If a system is biased to explore the whole range of values of the CV, it will perform a uniform sampling of all relevant metastable states as well as the Transition State Ensembles separating them. Given the structural complexity of biomolecular systems, as membrane-embedded ion channels, the choice of the most appropriate CV is not trivial at all. A direct solution would be to bias a huge number of CVs to make sure to include the relevant ones. Unfortunately, the volume of the space of CVs that the simulation needs to explore grows exponentially with the number of CVs so that, for practical purposes no more than 2–3 CVs can be simultaneously biased. To properly identify the CVs, knowledge is required of the system degrees of freedom that are relevant for the free-energy barrier separating the minima. In complex biochemical systems, this choice is typically performed in a system-dependent manner. Yet, these approaches are widely used for the study of conformational transitions in proteins. For the specific case of ion channels, they are useful both to investigate conformational transitions and to evaluate the PMF for ion translocation.

Umbrella Sampling and Metadynamics have proven to be particularly useful in the field of ion channel research. In the first approach, the interval of CVs values to be explored is discretized, and several copies of the system are simulated, each restrained via biasing potentials to the different values. In this way, equally accurate sampling can be obtained across the whole CV range, including high-free energy regions. The CV distributions obtained from the restrained simulations are then unbiased and combined together to yield the proper full distribution, and hence the free energy as a function of the CV, for example, using the Weighted Histogram Analysis Method (WHAM) [219]. Umbrella sampling has been frequently used to study ion permeation and selectivity [10,25,220–222] and ligand binding [223,224], for example.

In metadynamics, the energy function is supplemented with an history-dependent biasing potential composed by a sum of Gaussian functions centered in the visited regions of the CV space. In this way, the biasing potential discourages the visit of regions already explored and pushes the simulation towards unexplored areas of the conformational space. In the limit of a very long simulation, the biasing potential fills all energy minima and the free-energy profile can be simply attained by reversing the sign of the biasing potential. Ideally, a metadynamics simulation should be stopped when all energy minima have been filled. Since it is difficult to identify this moment, there may be problems of overfilling that limit the accuracy of the free energy calculation. This kind of problems can be avoided by using well-tempered metadynamics [225] where the height of the energy Gaussians decreases exponentially during the simulation and it is possible to tune the fraction of the energy basins that will be filled. Metadynamics has found multiple applications in ion channel research, for example, in studying permeation [226,227], gating related conformational changes [228], and ligand binding [229,230]

Collective Variables (CVs) allow for a reduced and physically meaningful representation of the process under investigation. If few CVs satisfactorily describe a process it is possible to use them to reconstruct a low-dimensional free-energy landscape, for example, via umbrella sampling or metadynamics, or to drive an accelerated exploration of the phase-space [231,232]. For ion channels, this approach can be pursued only in some special cases, for example, ion permeation [233] or inactivation [234]. Unfortunately, in the general case of gating, it is not easy to find simple descriptors for the complex structural rearrangements undergone by ion channels. High-dimensional CVs, however, preclude the reconstruction of large portions of their free-energy space. In those cases, path methods, such as the string method in collective variables [85] and related ones mentioned above, are more efficient because they entail a local optimization of the path. This efficiency, on the other hand, comes at the cost of exploring only a region close to the initial guess path. Overall, determination of CVs in gating problems is still an open challenge: in Section 6 we make some examples, which explore several rare event problems in ion channel research and discuss the choice of suitable descriptors.

4.2.4. Techniques that employ alchemical perturbations

Alchemical perturbations, such as the addition or removal of atoms, or the mutation of one amino acid residue to another during a simulation can be considered as an example of techniques that exploit the knowledge of end states to compute the free-energy difference between them. The most commonly employed approaches for this are the related methods of Free-Energy Perturbation (FEP) [235] and Thermodynamic Integration (TI) [236]. The approach exploits the fact that free energy is a state function so that the free-energy difference between end states can be attained as the sum of the free energies of a number of reactions arranged so as to close a thermodynamic cycle. Moreover, since the path is of no interest here, some or all of these steps may be un-physical which justifies the name of alchemical methods. In practice, the transition is broken into multiple steps (characterised by the parameter $\lambda$ that varies from 0 to 1), each of which is run as a separate stage of the simulation or as separate simulations (i.e. atoms slowly appear or disappear across this transition). The total free-energy change can be summed from the individual steps.

A typical application of alchemical methods is drug design, in which a key parameter is either the absolute free energy of binding of a ligand to a protein, or the relative-binding free energy of two different ligands [224,237]. In the former case, the reaction of interest is the migration of the ligand from bulk to the protein-binding site, which can be calculated from summing the free energy of ligand annihilation in bulk and ligand creation in the site. For relative free energies, the annihilation of one ligand in the binding site can be done simultaneously with the creation of the other. This approach has also found significant value in understanding the basis of ion selectivity, where the free energy to swap one ion type to another is computed at different points in the permeation pathway or in model-binding sites [20,238].

For the sake of completeness, we will cite an alternative method for the computation of free energies of binding, MM-PBSA or MM-GBSA [239]. In these methods, molecular mechanics energies (representing binding energies in vacuo) are combined with solvationfree energies (of ligand, receptor, and complex) computed by solving the linearized Poisson-Boltzmann or generalized Born equation. Using these contributions it is possible to build a thermodynamic cycle to compute binding free energies in solution. Of note, the direct calculation of the binding free energy of solvated reactants would be highly inaccurate because energy would be dominated by solvent–solvent interactions and the fluctuations of the total energy would be an order of magnitude larger than the binding energy. MM-PBSA and MM-GBSA can be considered intermediate in both accuracy and computational effort between empirical scoring and strict alchemical perturbation methods. In the last years, these approaches have been applied to a large number of systems with varying success [240–242].

4.3. Hybrid methods: Brownian Dynamics

Brownian Dynamics can be considered as a compromise between the computational efficiency of continuum methods and the molecular accuracy of atomistic ones. Finite size and ion–ion correlation effects are recovered in Brownian Dynamics (BD) [243], a technique where ions are explicitly modeled as particles while water, lipids, and proteins are still treated as continuous dielectric media. BD allows to observe statistically relevant numbers of permeation events in a time range from the microsecond to the millisecond and it thus represents the method of election for the computation of ion currents. Ion trajectories are computed integrating Langevin equation [244]:

(30) $m_i{\ddot{\mathit{x}}}_i(t)=-\frac{\mathrm{\partial }\mathcal{W}(\mathit{x}(t))}{\mathrm{\partial }{\mathit{x}}_i(t)}-{\gamma }_i{\stackrel{.}{\mathit{x}}}_i(t)+{\mathit{f}}_i(t).$(30)

The right-hand side of Equation (30)(30) $m_i{\ddot{\mathit{x}}}_i(t)=-\frac{\mathrm{\partial }\mathcal{W}(\mathit{x}(t))}{\mathrm{\partial }{\mathit{x}}_i(t)}-{\gamma }_i{\stackrel{.}{\mathit{x}}}_i(t)+{\mathit{f}}_i(t).$(30) represents the forces acting on particle $i$. The term $\mathcal{W}$ is an effective potential, whose negative gradient with respect to ${\mathit{x}}_i$ is the force stemming from the interaction with other particles, with the fixed charge distribution of the protein and membrane, and with the reaction field. The role of the solvent is implicitly accounted for by the friction term $-{\gamma }_i{\stackrel{.}{\mathit{x}}}_i(t)$ and the gaussian fluctuating force ${\mathit{f}}_i(t)$, which models random collisions.

The dynamics of ions in water is well approximated by the over-damped regime where the inertial term can be neglected and the Langevin equation reduces to

(31) ${\stackrel{.}{\mathit{x}}}_i(t)=-\frac{D_i}{k_{B}T}\frac{\mathrm{\partial }\mathcal{W}(\mathit{x}(t))}{\mathrm{\partial }{\mathit{x}}_i(t)}+{\zeta }_i(t)$(31)

where $D_i$ is the ion diffusion coefficient and ${\zeta }_i(t)$ is a Gaussian random noise. In principle, the diffusion coefficient should be position-dependent but most simulation works employ a single position-independent diffusion coefficient for all atoms of the same species. The effective potential $\mathcal{W}$ is split in a number of contributions that account for ion–ion interactions as well as the interactions of the ion with the pore, solution, and other biomolecules [97,245]:

(32) $\begin{array}{cc}\mathcal{W}\left(x\right)=& \sum _i\sum _{j\ne i}{u}_{ij}\left(\left|x_i-x_j\right|\right)+\sum _i{U}_{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}}\left(x_i\right)\\ & +\sum _i{q}_i\left[{\varphi }_{\mathrm{s}\mathrm{f}}\left(x_i\right)+{\varphi }_{\mathrm{r}\mathrm{f}}\left(x_i\right)\right]\end{array}$(32)

where the first contribution on the RHS models ion–ion interactions and includes a Lennard-Jones and a Coulomb term, the second contribution is a confinement potential that models the geometry of the pore representing the ion-accessible region, ${\varphi }_{\mathrm{s}\mathrm{f}}$ represents the static field, that is, the electrostatic potential generated by the fixed charge distribution of proteins and lipids, and ${\varphi }_{\mathrm{r}\mathrm{f}}$ is the reaction field, that is, the electrostatic potential generated by the polarization of the dielectric boundaries [246]. Notable extensions of the basic BD algorithm include that of Roux and coworkers [97,247], who designed a Grand Canonical Monte Carlo/Brownian Dynamics scheme that allows to model a variable number of ions and asymmetric ion concentrations, and that of Chung and colleagues [248], who included mobile protein elements.

5. Brownian Dynamics: applications

The accuracy of Brownian Dynamics can be enhanced if the effective potential is computed from fully atomistic Molecular Dynamics simulations. Using this approach Berneche and Roux [249] analyzed the conduction mechanism in the KcsA potassium channel. In this work, the effective potential was considered as a function of the axial coordinates of the three ions occupying the selectivity filter and it was decomposed as

(33) $W_t(x_1,x_2,x_3)=W_{eq}(x_1,x_2,x_3)+\sum _{i=1}^{3}q{\varphi }_{mp}(x_i)$(33)

where $W_{eq}$ is the equilibrium potential that was derived from all-atom umbrella sampling simulations and ${\varphi }_{mp}(x)$ is the transmembrane potential profile that was computed with Poisson-Boltzmann equation. In particular, the Brownian simulation was implemented as a continuous-time Markov chain with discrete states corresponding to the ions positions. Forward and backward transition rates were computed as

$\begin{array}{cc}{k}_{(z_1,z_2,z_3)\to (z_1\pm \delta z,z_2,z_3)}& =\\ \hspace{0.17em}& \left(\frac{D(z_1)+D(z_1\pm \delta z)}{2\delta z^2}\right)exp\left(-\frac{W_t(z_1\pm \delta z,z_2,z_3)-W_t(z_1,z_2,z_3)}{2k_{B}T}\right)\end{array}$

Using this approach the authors measured a maximal conductance of the channel in excellent agreement with experimental data [17] (550 and 360 pS for outward and inward fluxes, respectively). More importantly, the large statistics of crossing events, accrued during the long BD simulations, enabled the authors to clarify the mechanism of outward permeation which proceeds along the sequence of states $[S_3,S_1]\to [S_4,S_3,S_1]\to [S_4,S_2,S_0]\to [S_4,S_2]\to [S_3,S_1]$, where $S_i$ with $i=0,\cdots ,4$ identifies the binding sites inside the selectivity filter of KcsA. Inward permeation was found to occur in the reverse order. When this information was integrated with the free energy maps it became clear that the key step of the process was the $[S_4,S_3,S_1]\to [S_4,S_2,S_0]$ transition because configuration $[S_4,S_3,S_1]$ was less stable than $[S_4,S_2,S_0]$. This finding thus justified the difference between the outward and inward conductance, challenging the traditional view of isoenergetic states and barrierless conduction.

The use of ion channels and nanopores for DNA sequencing is another application where the use of MD simulations to tune the parameters of Brownian Dynamics can lead to invaluable results. The idea to sequence DNA exploiting the signature currents of the four nucleotides as the chain is electrophoretically pulled through a nanopore, dates from the beginning of the 2000s [250,251]. However, it soon became clear that the current cannot depend on a single base but more likely a stretch of a few nucleotides whose sequence varies in time in a sequence-dependent manner. Since no experimental device can currently image the conformation of a DNA chain as it is pulled through a nanopore, Molecular Dynamics algorithms have become the tool of election. Unfortunately, the computational measurement of the small currents flowing through ion channels requires prohibitively long simulations. For instance, a typical 100 pA ion current corresponds to the observation of approximately 10⁴ ion crossing events which requires a 16 $\mu$s MD simulation. Simulations of this length can be routinely run using Brownian Dynamics but the accuracy of the simulation requires a careful tuning of the effective potential using information deriving from atomistic MD. In the protocol designed by Comer and Aksimentiev [252] (named Atomic Resolution Brownian Dynamics – ARBD) the effective potential is decomposed in an ion–ion interaction term and a term of interaction of the ions with nanopore and DNA:

(34) $W({\mathit{x}}_1,{\mathit{x}}_2,\cdots )=\sum _{ij}{W}_{ij}^{ion-ion}(|{\mathbf{x}}_i-{\mathbf{x}}_j|)+\sum _i{W}_i^{sys-ion}({\mathbf{x}}_i)$(34)

The $W^{ion-ion}$ term was computed through three sets of Umbrella Sampling simulations that determined the PMF of the K⁺-K⁺, K⁺-Cl ^– and Cl^–-Cl ^– pairs as a function of the distance between the two ions. The $W_i^{sys-ion}$ were also computed from Umbrella sampling simulations where a single potassium or chloride ion was restrained to the center of the cells of a hexagonal lattice surrounding a single nucleotide restrained to the typical B-form of DNA and immersed in a water box. Comer and Aksimentiev showed that the 3D-PMFs generated from isolated nucleotides can be combined together to generate the PMF of large systems like oligonucleotides and base pairs. As an example, the reconstructed PMF of the A$\cdot$T base-pair, $w_{AT}$, can be computed from the PMF of the isolated A and B nucleotides ($W_A$ and $W_B$ respectively) as:

(35) $w_{AT}(\mathbf{X})=W_A(T_A^{-1}\mathbf{X})+W_B(T_B^{-1}\mathbf{X})$(35)

where $T_A$ and $T_B$ are the rigid-body transformations that map the coordinates ${\mathit{x}}_i^A$ and ${\mathit{x}}_i^B$ of the isolated A and B nucleotides to the coordinates ${\mathit{X}}_i$ of the same nucleotides in the A$\cdot$T base-pair. In a similar fashion, Comer and Aksimentiev showed that the PMF of a single-stranded DNA chain (ssDNA) can be attained by summing the PMF of the individual component nucleotides (with appropriate geometric transformations). Moreover, if a nucleotide is replaced with another one, the PMF of the ssDNA can be updated by appropriately replacing the PMF of the mutated base. The limit of this approach is that it cannot be used to build the PMF of double-stranded DNA (dsDNA) because the interaction of potassium ions for the minor groove is significantly under-estimated. The PMF of dsDNA must therefore be computed through atomistic umbrella sampling simulations. Mutations, however, can then be applied just as with ssDNA. It is important to stress that the computational efficiency of the method by Comer and Aksimentiev is critically reliant on the ability to combine elementary PMFs to generate the PMF of larger systems. Indeed, despite the power of Brownian Dynamics, there will be no computational saving if for each new system an expensive umbrella sampling simulation must be run to produce the effective potential. The accuracy test showed that ARBD generates ion distributions consistent with those predicted by atomistic simulations. Moreover, at low salt concentrations, ARBD computes currents in quantitative and qualitative agreement with those of atomistic MD. Finally, at higher concentrations, the agreement between ARBD and atomistic simulations is only qualitative. The authors, however, observed that the currents computed by ARBD in bulk solution are in better agreement with experimental data than the currents predicted by atomistic simulations. It remains to be tested whether this better agreement is also preserved in the confined environment of nanopores partially occluded by nucleotide chains.

6. Atomistic approaches: applications

The most notable properties of ion channels are ion selectivity, high permeation rates, and gating which allow them to orchestrate complex biological functions, such as neuron firing and muscle contraction. Because of the problem of rare events explained in Section 4.2, the characterisation of each of these properties requires enhanced sampling techniques. In the following we review few selected examples in which different techniques were used to assess permeability, selectivity, and gating in ion channels or models thereof. These examples are by no means exhaustive and reflect the authors’ expertise.

6.1. Permeability and selectivity of NaChBac

Guardiani et al. [253] used metadynamics to study the permeability of a bacterial voltage-gated sodium channel, NaChBac. Interestingly, while a PMF derived from a long equilibrium simulation exhibited four minima (in agreement with a study by Catterall and coworkers on NavAb [254]), the PMF derived from a metadynamics simulation biasing a single Na⁺ ion to explore the Selectivity Filter (SF) featured a single minimum so deep that the ion would remain trapped inside. The mismatch turned out to be due to an energy barrier in the middle of the SF that an ion could overcome only after receiving an electrostatic kick from a second incoming ion. Indeed, this knock-on mechanism was clearly identified when the metadynamics simulation was repeated biasing two ions instead of one. These results clearly highlight a requirement common to metadynamics and all collective-variable driven methods, including umbrella sampling [213], i.e., these approaches need some preliminary knowledge of the number and kind of collective variables relevant to describe the process at hand; in this case, this meant that some rough idea of the permeation mechanism had to be known a priori.

Metadynamics is a powerful tool also for the study of selectivity. Indeed, equilibrium and metadynamics simulations [255] revealed that a single calcium ion could access the SF of NaChBac but then remained stuck inside. This is due to the larger binding free energy of Ca²⁺ as compared to Na⁺ and to the unlikelihood of a Ca²⁺/Ca²⁺ knock-on mechanism. Indeed, in metadynamics calculations with a fixed and a mobile ion, the free-energy profile experienced by the incoming ion featured a very high barrier due to the electrostatic repulsion of the resident ion.

6.2. Bias exchange metadynamics simulations of NavAb

The NaChBac simulation study reported in Section 6.1 highlights a limit of metadynamics that is shared by many other enhanced sampling techniques: the necessity of a preliminary knowledge of the type and number of ions involved in the permeation mechanism. As also mentioned in Section 7, in Bias Exchange (BE) metadynamics several replicas of the system are run in parallel. Each replica runs a metadynamics simulation with bias on different collective variables, and structural exchanges are attempted at regular time intervals. The exchanges are accepted with probability $p_{acc}=min[1,exp\beta \mathrm{\Delta }V]$ where $\mathrm{\Delta }V=[V_1({\mathbf{x}}_1)+V_2({\mathbf{x}}_2)]-[V_1({\mathbf{x}}_2)+V_2({\mathbf{x}}_1)]$. This technique naturally lends itself to the study of multi-ion permeation processes. Indeed, if $N_{max}$ is the maximal number of ions that are expected to be involved in permeation, it is sufficient to use $N_{max}$ replicas each biasing the axial position of a single ion (plus a neutral replica where no biasing potential is applied). Since ions are free to leave or remain inside the pore, the number of ions taking part in conduction will also vary. Moreover, while each replica enhances the dynamics along the single biased CV, the exchanges between replicas enhance the sampling along the whole set of CVs biased in all the replicas. Domene et al [256] applied the method to the study of a simplified model of the NavAb sodium channel. Indeed BE-metadynamics yielded 1D- and 2D-PMF in excellent agreement with those generated by umbrella sampling at a fraction of the computational cost of the latter technique: while fullfree energy convergence took 300 ns in BE-metadynamics, more than 1 $\mu$s was required by umbrella sampling. Furthermore, in BE-metadynamics, the number and type of ions involved in conduction events did not need to be predefined, and conduction events where different numbers of ions are involved could be analyzed and compared using a single simulation.

6.3. Structure of claudin-based paracellular channels

Recently, computational studies on claudins, transmembrane proteins that are part of tight-junctions between adjacent epithelial or endothelial cells, have been reported. Using structural modelling and all-atoms molecular dynamics simulations Alberini et al [257] refined the channel architecture proposed by Suzuki et al [258,259]. To validate the model, the authors verified its experimentally known selectivity properties for cations [260] by employing the Voronoi-tessellated Markovian Milestoning method for free energy and kinetics calculations [205,261]. Milestoning allows reconstructing the long-time dynamics of a system from the crossing statistics of its trajectory through a set of hypersurfaces (the milestones) placed along the reaction coordinate [262,263]. It can be employed to obtain mean-first-passage-times and in turn permeability coefficients (compare with Section 3.2) across membranes or channels, even without relying on the so-called inhomogeneous solubility-diffusion model [264]. Voronoi-tessellated Markovian Milestoning is a version of the method that leverages independent simulations confined within a set of cells spanning the reaction coordinate. Transition path theory [265] is then used to obtain the kinetic properties of the full reaction from hitting statistics at the cell boundaries, identified as milestones.

6.4. Structures and permeability of nAChR

Enhanced sampling techniques are also very important for the functional annotation of Ligand Gated Ion channels (LGIC) that requires the investigation of the ion translocation through the pore crossing high free-energy barriers. Within this framework Chiodo et al [266–269] investigated the $\alpha$7 nicotinic acetylcholine receptor (nAChR), a widely expressed LGIC of the brain related to schizofrenia and Alzheimer’s disease [270,271]. Performing analysis of the profile and hydration in extensive molecular dynamics simulations of a homology model, the authors were able to associate different conformations to four different functional states: the open active state [266], a desensitized state [267], a closed-locked, non-conductive state [268] and the apo-resting state [268]. A comparison with the recent experimental structures of the channels in different states [272,273] confirms the accuracy of the models.

To provide a refined description of the process, the ion permeation across the all-atoms full-length (TMD+LBD) human $\alpha$7 channel model, was studied both in native and in the E-1’A mutant [274]. The single-ion PMF (i.e. the PMF of one ion while no other ions are present in the pore, see $U(x)$ in Equation (23)(23) $\mathcal{P}=1/{\left[{\int }_0^L\frac{exp\beta U\left(x\right)}{D\left(x\right)}dx\right]}^{-1}$(23) ) and the ion translocation kinetics were reconstructed for sodium and chloride by using the milestoning method with Voronoi tessellation [205,261].

The 1D PMF profiles shown in Figure 7 allowed to identify the structural determinants of ion translocation, i.e. the key residues responsible for the formation of energy barriers and kinetic traps, along the protein structure represented in background in Figure 7. The mean first passage time to traverse the full channel is smaller for sodium than chloride consistent with the experimentally determined cationic nature of wild type $\alpha$7. Moreover, results indicate that the free-energy barriers in the transmembrane domain play the major role in ion permeation, in agreement with results from simulations on GLIC [275].

Figure 7. PMF for the permeation of sodium (black) and chloride (red) ions through the full $\alpha$7 native channel. The curves are shifted along the y axis so that their values match at the intracellular end (on the left). On the background: two protein subunits are shown, represented in cartoon, for the sake of clarity. Lipids are in line representation; water represented with points. Dots indicate the centers of the Voronoi cells. Figure reproduced with permission from [274]. Copyright 2020 American Chemical Society.

6.5. Gating of hERG

So far examples focused on permeation and selectivity but ion channels exhibit another key property, gating, i.e. the ability to open and close in response to specific stimuli. A particularly interesting case is represented by hERG, a potassium channel involved in the repolarization of cardiomyocytes [276], whose mutations and nonspecific interactions with a wide range of drugs leads to the Long QT Syndrome [277,278]. In order to limit the incidence of this disease, the FDA has promoted the Comprehensive in vitro proarrhytmia assay [279], a screening program to exclude from the pharmaceutical development pipeline drugs exhibiting unintended interactions with hERG.

From the structural point of view hERG is formed by four identical subunits each comprising a Voltage Sensor Domain (VSD), a Pore Domain (PD), and a carboxy-terminal domain [280], see Figure 8(a,b). hERG responds to changes in the membrane potential with a movement of helix S4 of the VSD. This motion is somehow transduced to the PD where the gate is located. In order to clarify the mechanism of electro-mechanical coupling Costa et al [52] used a network theoretical approach already applied to the study of allosteric enzymes [49,50] and ion channels [51,281,282]. The channel is modeled as a network where each node corresponds to a residue, arcs correspond to interactions between residues and edge weights quantify the efficacy of motion propagation. Once the network is built, the minimal paths connecting all residues of a source and a sink region are computed using Dijkstra’s algorithm. Using this approach two families of VSD-PD pathways were identified, see Figure 8(c). Interestingly, several mutations known to affect hERG gating could be mapped on the discovered paths suggesting that mutations act by interrupting the electro-mechanical coupling between VSD and PD.

Figure 8. Structure of the hERG channel (b) with definition of the VSD, PD, and CTD (a). Pathways connecting the VSD and the PD (c) identified for a representative subunit by a network theoretical approach, which considers the residues as the nodes and the distances related to their correlations. Figures reproduced with permission from [52] under CC BY 4.0 license.

6.6. Hydrophobic gating in model channels

While hERG exhibits the typical voltage-gated steric occlusion of the channel, in other channels ion currents are blocked by hydrophobic gating [38,44,283,284]. Hydrophobic gating involves the formation of a vapour phase – a bubble – nucleating from the liquid phase inside the hydrophobic portion of the cavity of the biological channel [268]. Bubble formation in confinement has been shown to be deeply dependent on several characteristics of the nanopore, such as hydrophobicity, geometry, and size [45].

Simplified models of ion channels are a useful tool to disentangle the competing factors affecting the energetics of bubble formation in ion channels [44,283] and to provide microscopic interpretations of experimental results on nanoporous materials [285]. An example of this physical approach to hydrophobic gating is the study of the drying mechanism of model nanopores in the presence of an hydrophobic solute dissolved in water [47]. A cylindrical nanopore with the typical hydrophobicity of nonpolar aminoacids was studied by means of a rare event method, Restrained Molecular Dynamics [231], in order to characterize the kinetics of drying. Two collective variables were used: one controlling the water filling of the pore and the other the axial position of the gas. The bubble gating process was shown to be crucially modified by the physical presence of the hydrophobic solute inside the channel, see Figure 9. In particular, a single atom of argon was enough to abate the drying free-energy barrier between the wet and the dry state of the hydrophobic nanopore, causing the instantaneous drying of the whole nanopore for sizes comparable to the inner section of biological channels. Moreover, it was shown that the gas particle stabilizes the dry state of the nanopore and that the most probable mechanism of drying for this system corresponds to the infiltration of the gas particle inside the wet nanopore followed by the drying event. Such mechanism becomes relevant in connection to the context of general anesthetics, and in particular of volatile anesthetics, such as noble gases [286] and to the hypothesis that bubble formation may be involved in general anesthetics [44].

Figure 9. Illustration of the effect of a single hydrophobic particle on the dewetting of a model nanopore: the dewetting barrier is abated when an Argon atom enters the water-filled pore, thus accelerating the process of hydrophobic gating. Reproduced with permission from [47]. Copyright 2020 American Chemical Society.

6.7. Identification of functional states of K_v and Na_v channels

Voltage-gated ion channels undergo conformational changes in response to membrane potential, cycling between closed and open states. While X-ray crystallographic structures of K_v1.2, K_v1.2/K_v2.1 chimera and Na_vAb channels [254,287,288] in the active state became available already in the early 2000s, the resting state structure of the VSD resisted for very long to experimental determination. As of 2022 the situation has not particularly improved. In 2014, the VSD of Ciona intestinalis was crystallized in both the active and resting states [289], but the resting states of most other K_v and Na_v channels are still unknown.

This experimental bottleneck prompted massive computational efforts to predict the resting state. Early models were obtained using the Rosetta program [290] and refined through atomistic simulations [291]. Other models were generated through restrained MD simulations enforcing a wide array of constraints deriving from experiments of mutagenesis, cross-linking, fluorescence, and resonance-energy transfer. These techniques indeed identify residue–residue interactions that can be readily translated into distance constraints. A notable work in this field was performed by Vargas et al [292] who computationally engineered metal bridges between pairs of residues experimentally predicted to be in close proximity. An alternative approach is based on simulations incorporating the effect of the membrane potential. In this case, the ambition is not only to predict the resting state of the VSD but also to reconstruct the complete sequence of events of the active $\to$ resting conformational transition induced by membrane hyperpolarization. A particularly important contribution was given by Jensen et al [14] who ran hundreds of microseconds of MD simulations using a non-physiological potential (−750 mV) to accelerate the transition.

Remarkably, despite the wide variety of computational methods used by these investigators, the research results were largely in agreement with each other, yielding a ‘low-resolution’ consensus model [293] (see Figure 10) that can be used as a benchmark to discuss the idealized gating mechanisms proposed in the literature [294–297]. The consensus model [293] of the resting state of the VSD lends strong support to the helical-screw/sliding-helix model [294,295]. Indeed, all the resting state models show a significant roto-translation of helix S4 that undergoes a helical displacement of $\sim$ 10 Å. This displacement is significantly smaller than that predicted by the paddle model (20–25 Å) [296], but larger than that implied by the transporter-like model [297]. Even if extreme caution should always be exercised when dealing with computational predictions, the consensus model clearly shows the power of computational modeling and MD simulations in the solution of key problems of membrane physiology.

Figure 10. Main elements of secondary structure of the VSD in the active and resting state. (a) The VSD in the active conformation (taken from the x-ray structure; PDB accession no. 3LUT). (b) A superposition of different models of the resting- state configuration of the VSD obtained by independent research teams using different constraints and methodologies. The four helices, S1 (gray), S2 (yellow), S3 (red), and S4 (blue), are displayed. The spheres correspond to the C$\alpha$ atoms of E1 (Glu226) and E2 (Glu236) along S2, and R1 (Arg294) along S4. ©(2012) Vargas et al. Originally published in Journal of General Physiology [293].

7. Outlook and perspectives

The possibility to perform biomolecular simulations critically relies on the availability of good-quality structures. This necessity can be summarized in the statement: no structure, no simulation. This problem is particularly serious in the case of membrane proteins that are too large and hydrophobic to be resolved through NMR spectroscopy and are notoriously difficult to crystallize [298]. Indeed, membrane proteins represent between 20 and 30% of the proteome of most organisms, yet, to date, less than 1% of the entries of the Protein Data Bank are membrane proteins, and only a small fraction of these are eukaryotic. Eukaryotic membrane proteins cannot be expressed in prokaryotic systems because they require specific eukaryotic machinery for post-translational modifications and to be trafficked to the membrane. This is why eukaryotic membrane proteins cannot be produced recombinantly through standard genetic engineering techniques but they need to be purified from native sources. This, however, leads to another problem, because membrane proteins are unstable [299] even when treated with the mildest detergents needed for their extraction. The necessity to use detergents leads to high-volume crystals containing large amounts of solvent. These crystals are fragile, diffract at low resolution and suffer from irradiation damage. Over the last few years, cryo-electron microscopy (cryo-EM) has emerged as a powerful alternative to X-ray diffraction for the structural resolution of membrane proteins [71]. Indeed, cryo-EM requires small volumes and lower concentrations so that it can be applied to samples purified from natural sources. Moreover, cryo-EM does not require removal of flexible loops or the addition of stabilizing substrates to induce crystallization. Finally, while X-ray structures may suffer from artifacts due to crystal packing forces, the rapid quenching to the temperature of liquid nitrogen required by cryo-EM freezes the protein in a set of conformations more likely to represent functional states [300].

Although the overall number of available ion channel structures is still underrepresented as compared to other protein families, the improvement in the experimental techniques provided in the past five years has led to a remarkable increase in their availability. This is, for example, the case of human nicotinic receptors and voltage gated sodium channels for which a wealth of eukaryotic structures have been very recently deposited [272,273,301–309].

The scarcity of experimental ion channel structures prompted in the years the development of a series of computational tools to predict ion channels directly from protein sequences. Homology modeling combined with atomistic simulations provides a valuable approach to investigate the structure-dynamics-function relationship in several ion channels [310,311]. Methods based on machine learning algorithms have been recently developed in ion channel research, from sequence-based predictions to the topology of the transmembrane region, including structure-activity relationship [312–314]. Moreover, machine learning shows further promise for structure prediction, see the AlphaFold [315] successes at CASP13 and CASP14.

Nowadays the most popular technique for the computational study of membrane proteins is represented by Molecular Dynamics simulations. Despite impressive improvements over the last few years, MD simulations still suffer from a number of limitations that in 1990 were highlighted by Karplus and Petsko [316] in a sentence that still holds true today: ‘Two limitations in existing simulations are the approximations in the potential energy functions and the lengths of the simulations. The first introduces systematic errors and the second statistical errors’.

In order to overcome the approximations of the potential energy functions, massive efforts are being performed in the development of more accurate force fields. A major limitation of most currently used force fields is that they employ fixed atomic charges corresponding to the electrostatic environment of bulk solvent. This assumption is not always correct. For instance, non-polarizable force fields under-estimate by a factor of 2 the dielectric constant of the lipid bilayer [317] doubling the energy barrier of an ion crossing the membrane. The lack of lipid polarization is reflected in the free-energy profile of an ion crossing the gramicidin A channel [318–320]. Non-polarizable force fields also lack accuracy in the modelling of divalent cations like Ca²⁺ and Mg²⁺ [321] that, owing to their high charge density, tend to polarize the surrounding water shell. At the moment, the most commonly used polarizable force fields are the Drude oscillator [322,323] and the AMOEBA force field [324]. Even if promising results have been reported, the development of these force fields is still a work in progress hampered by the limited availability of model parameters and the increased computational cost [325]. For instance, using NAMD the computational cost of a simulation with a polarizable force field is twice that of the non-polarizable counterpart [326,327]. The increasing use of GPU cards for biomolecular simulations will likely provide new momentum to the development of polarizable force fields [328]. Besides the issue of computational cost, the modeling of electrostatic interactions also needs further refinement. Polarizable force fields are more accurate than classical force fields in the simulation of heterogeneous environments featuring regions with different dielectric constants or different hydrophobicity.

The problem of the short length of MD trajectories, already highlighted by Karplus and Pesko [316], lies at the heart of any simulation method. Indeed, MD simulations have the ambition to connect structure with function through dynamics but, while current trajectories routinely cover hundreds of nanoseconds, relevant biological phenomena occurs in milliseconds to seconds. For instance, this is the case of electrophysiological recordings that are measured over several milliseconds. In order to rationalize at the molecular level the experimental patterns, MD simulations need to reach longer time-scales. A brute force approach to overcome this limitation is to increase computing power. While in the 1980s supercomputers contained just a few powerful Central Processing Units (CPUs) in the 1990s the trend emerged to combine several CPUs into a cluster [72]. Nowadays supercomputers of top high-performance computing centres may contain hundreds of thousands to millions of CPU cores. Another trend involves the use of Graphical Processing Units (GPUs) originally designed for gaming [78]. This hardware is particularly well suited for highly repetitive operations like the force calculations that represent the most demanding part of MD simulations. The bottleneck here is the necessity to translate MD codes from Fortran or C to a more appropriate programming language like CUDA. However, given the small cost of GPUs it can be foreseen that in the near future ordinary researchers will afford to buy multi-GPU workstations providing the opportunity to do high-performance computing at their desk. An interesting alternative to CPU- and GPU-based general purpose computers is the design of a special purpose machine with an architecture tailored for MD simulations like the Anton supercomputer [79]. Anton was named after the early microscopist Anton van Leeuwenhoek with the declared ambition to develop a computational microscope with a resolution unattainable by any experimental method. Anton, that owes part of its high performance to a parallel architecture and optimized inter-node communications, has led to several important discoveries. For instance, Anton has helped to clarify the mechanism of transition between open and closed states in voltage gated potassium channels [13] and has allowed simulations of drug-binding in G-protein coupled receptors [329]. Moreover, the long simulations allowed by Anton revealed the flaws of a force field that would remain latent in shorter simulations. The third generation of Anton supercomputers has been announced in 2021 [330], which claims to perform on 512 nodes $100\mu$s per day for a million-atoms system.

In many cases, even the use of high-performance hardware does not enable sufficient sampling – anyway not efficient – calling for advanced enhanced sampling techniques. In section 4.2 we described several enhanced sampling approaches. Two of them, Metadynamics and Replica Exchange Molecular Dynamics have proven to be particularly effective, but they still present a number of limitations. The most critical problem of metadynamics is represented by hidden degrees of freedom, which may not be accurately described by the chosen CVs and can frustrate the exploration of the phase space, sometimes limiting the extent of convergence and accuracy of results, see, for example, the rightmost panel in Figure 11 in which biasing along a single CV, say S₁, does not allow the correct reconstruction of the free-energy landscape. REMD simulations, on the other hand, do not need any preliminary choice of CVs and they enhance the sampling of all degrees of freedom through high temperatures. Unfortunately, the necessity to use a number of replicas proportional to the square root of the number of atoms, makes the technique too computationally demanding for the simulation of large systems as a membrane-embedded ion channel. It is thus clear that the strengths and weaknesses of CV-based methods and REMD are complementary and their combination generates a new class of promising methods, see also [331,332]. A first combination of the two techniques is represented by Parallel Tempering Metadynamics [82]. In this algorithm, a number of metadynamics simulations biasing the same CV are run in parallel at different temperatures and exchanges between neighbouring replicas are performed with the Metropolis rule as in REMD. In this scheme, metadynamics ensures an effective sampling along the chosen CV while the exchanges with high-temperature replicas allow sampling along orthogonal CVs. Another hybrid algorithm combining metadynamics and REMD is Bias Exchange Metadynamics [83] where a number of replicas are run in parallel with bias on different CVs and exchanges are performed at regular time interval with the Metropolis criterium.

Figure 11. Illustration of machine learning applications to the analysis and enhancement of MD simulations. High dimensional data, e.g. the protein trajectory, is used as the input of an artificial neural network (ANN). The ANN is trained to project the input to a low dimensional space minimising a loss function. Depending on the structure of the ANN and on the loss function, the low dimensional representation captures different features that are considered to be important, such as slow modes. In some methods, the features learnt are used to further enhance the sampling of MD simulation as shown by the arrow in the bottom. Reproduced from [340], Copyright 2020, with permission from Elsevier.

Powerful hybrid techniques arise not only from the combination of different enhanced sampling techniques but also from the combination of simulation methods with different resolutions. These methods are particularly important when a part of the system needs to be treated with atomistic detail while the rest of the system is very large and just acts as the environment of the relevant part, so that it is more convenient to treat it at a coarse-grained level. This situation is particularly relevant for the study of ion channels. For instance, the simulation of a whole vesicle is computationally very expensive so that membrane lipids could be described using a CG model while the ion channel or the few membrane proteins of interest could be described at atomistic level. This would also account for the fact that the MARTINI force field currently models lipids more accurately than proteins. Such a framework is provided, e.g., by the adaptive resolution scheme (AdResS) in its different flavours [333–335] where the atomistic and coarse grained domains are connected by an intermediate region that allows for a continuous transition between the two domains. The passage between regions modeled at different resolutions is described by a sigmoid transition function $\lambda (x)$ that grows monotonically from zero (CG region) to one (atomistic region) taking values $0<\lambda <1$ in the hybrid region. Assuming a mono-dimensional case, the transition function allows the calculation of the force acting on residue $\alpha$:

(36) ${\mathit{F}}_{\alpha }(x_{\alpha })=\sum _{\beta \ne \alpha }[\lambda (x_{\alpha })\lambda (x_{\beta }){\mathit{F}}_{\alpha \beta }^{AA}+(1-\lambda (x_{\alpha })\lambda (x_{\beta })){\mathit{F}}_{\alpha \beta }^{CG}]+{\mathit{F}}_{th}(x_{\alpha })$(36)

where the atomistic ${\mathit{F}}_{\alpha \beta }^{AA}$ and coarse grained force ${\mathit{F}}_{\alpha \beta }^{CG}$ are just the negative gradients of the corresponding potentials

(37a) $V_{\alpha }^{AA}=\frac{1}{2}\sum _{\beta \ne \alpha }^N\sum _{ij}^n{V}^{AA}(|{\mathit{x}}_{\alpha i}-{\mathit{x}}_{\beta j}|)$(37a)

(37b) $V_{\alpha }^{CG}=\frac{1}{2}\sum _{\beta \ne \alpha }{V}^{CG}(|{\mathit{X}}_{\alpha }-{\mathit{X}}_{\beta }|)$(37b)

with ${\mathit{x}}_{\alpha i}$ and ${\mathit{x}}_{\beta j}$ the position vectors of atom $i$ of residue $\alpha$ and atom $j$ of residue $\beta$ while ${\mathit{X}}_{\alpha }$ and ${\mathit{X}}_{\beta }$ are the position vectors of the center of mass of residues $\alpha$ and $\beta$. It is important to note that the equations of state of the atomistic and coarse grained systems are very different so that a given value of the density corresponds to different values of the pressure. Therefore, in order to simulate the multiscale system at uniform density without dissipating the pressure gradient between the atomistic and coarse grained domains, it is necessary to add a thermodynamic force (the term ${\mathit{F}}_{th}(x_{\alpha })$ in Equation (36)(36) ${\mathit{F}}_{\alpha }(x_{\alpha })=\sum _{\beta \ne \alpha }[\lambda (x_{\alpha })\lambda (x_{\beta }){\mathit{F}}_{\alpha \beta }^{AA}+(1-\lambda (x_{\alpha })\lambda (x_{\beta })){\mathit{F}}_{\alpha \beta }^{CG}]+{\mathit{F}}_{th}(x_{\alpha })$(36) ). This force is determined iteratively starting from the initial choice ${\mathit{F}}_{th}^0=M/\rho \mathrm{\nabla }p(x)$. The method has been initially used [336] for the calculation of chemical potentials and solvation-free energies of molecular systems. For instance, a characterization of the solvation shell of fullerenes revealed that surface-induced structuring of water does not go beyond the second hydration layer. The AdResS method, however, is extremely flexible and recent applications include quantum-classical coupling [337], non-equilibrium systems [338], and advanced free-energy calculations [339].

A class of techniques that is finding increasing application in the arena of biomolecular simulation is represented by Machine Learning (ML). ML is leading a revolution in all fields of science and technology with applications ranging from image and speech recognition to the prediction of drug action. In particular, we will discuss the application of deep feedforward neural networks [341]. A neural network is just a set of activating functions called neurons and organized in layers (Figure 11). A neuron in a layer receives as input the weighted sum of the outputs of the previous layer, transforms this total input through a non-linear function, and passes the output to the next layer. The most important feature of deep neural networks is that, through a training stage, they can automatically ‘learn’, that is, adjust their weights to fit any function of interest. Indeed, the error function, quantifying the distance between the solution proposed by the machine and the correct one, can be seen as a multidimensional function of the weights and can be minimized through an algorithm similar to steepest descent called stochastic gradient descent [342]. In order to apply machine learning to biomolecular calculations the neural network must be fed with inputs respecting the physical constraints [343]. For instance, if we consider an oxygen molecule in vacuo with a six-dimensional coordinate vector $\mathbf{x}$, the energy will be invariant to roto-translation, while the force will be equivariant:

(38a) $U(\mathit{R}\mathit{x}+\mathit{T})=U(\mathit{x})$(38a)

(38b) $-\mathrm{\nabla }U(\mathit{R}\mathit{x}+\mathit{T})=-\mathit{R}\mathrm{\nabla }U(\mathit{x})$(38b)

where $\mathit{T}$ is a translation vector and $\mathit{R}$ is a rotation matrix. In order to retain the invariances and equivariances, the network should not be fed with a vector of cartesian coordinates but with a set of distances and angles which significantly reduces the dimensionality of the problem. Other restraints that should be accounted for are: permutational invariance, for example, of solvent molecules, probability conservation, and detailed balance. ML can be applied to a wide range of problems of interest for biomolecular simulation. The list we provide below has no claim to be exhaustive.

1. Calculation of potential energy surfaces. A number of classical methods already exists [344,345] to fit energy values derived from quantum mechanical calculations, but these often suffer from poor transferability. A viable alternative is to train a deep neural network to minimize an energy or a force loss function [346,347]:

(39a) $L_{ene}=\sum _i[\hat{U}({\mathit{x}}_{\mathit{i}},\mathit{\theta })-U_i{]}^2$(39a)

(39b) $L_{force}=\sum _i|\mathrm{\nabla }\hat{U}({\mathit{x}}_i,\mathit{\theta })+{\mathit{f}}_i{|}^2$(39b)

where $\hat{U}({\mathit{x}}_{\mathit{i}},\mathit{\theta })$ is the energy predicted by the network for configuration ${\mathit{x}}_{\mathit{i}}$ using parameter vector $\mathit{\theta }$, while $U_i$ is the energy computed ab initio for the same conformation. In a similar way, $-\mathrm{\nabla }\hat{U}({\mathit{x}}_{\mathit{i}})$ is the force predicted by the neural network and ${\mathit{f}}_{\mathit{i}}$ is the force computed through quantum mechanics. This approach can be used to optimize the parameters of a force field.

• (2) Calculation of free-energy surfaces. This problem is formally similar to the estimate of Potential Energy Surfaces, but the training is done on classical MD trajectories of the rare event of interest. Indeed, the neural network [348] must be trained to minimize the free-energy loss

(40) $L_{free-ene}=\sum _i|{\mathcal{F}}_i-\hat{\mathcal{F}}({\mathit{y}}_i,\mathit{\theta }){|}^2$(40)

where $\hat{\mathcal{F}}({\mathit{y}}_{\mathit{i}},\mathit{\theta })$ is the free energy predicted by the network using parameter vector $\mathit{\theta }$ at point $y_i$ in the space of collective variables while ${\mathcal{F}}_i$ is the free energy of the same configuration computed with a method able to sample the equilibrium distribution $\mu (\mathit{y})$: $\mathcal{F}(\mathit{y})=-log\mu (\mathit{y})+const$.

• (3) Coarse graining. ML can be used to learn the effective energy of a coarse grained molecular model [349,350]. Let us consider a transformation $\mathit{y}=\mathit{E}\mathit{x}$ that associates a vector of cartesian coordinates $\mathit{x}\in R^{3N}$ to a vector of coordinates of the coarse grained particles $\mathit{y}\in R^{3n}$ with $n<N$. The transformation is assumed to be linear because normally the CG coordinates are a subset of the atomistic ones. For instance the CG particle representing a residue is normally placed in the position of the $C_{\alpha }$ atom of that residue. The effective energy of the CG model is given by the free energy as a function of variables $\mathit{y}$: $\mathcal{F}(\mathit{y})=-log\mu (\mathit{y})+const$. If the coarse grained model is to be thermodynamically consistent with the atomistic one, the force predicted by the neural network $-{\mathrm{\nabla }}_y\hat{\mathcal{F}}(\mathit{E}\mathit{x},\mathit{\theta })$ must match the force on CG particles ${\mathit{f}}^{\mathit{y}}(\mathit{x})$ computed from the atomistic forces $\mathit{f}(\mathit{x})$. The loss function that the network must learn to minimize is thus

(41) $\sum _t|{\mathrm{\nabla }}_{\mathit{y}}\hat{\mathcal{F}}(\mathit{E}{\mathit{x}}_{\mathit{t}},\mathit{\theta })+{\mathit{f}}^{\mathit{y}}({\mathit{x}}_t){|}^2$(41)

where the sum runs over all conformations $t$ of a training set.

• (4) Automatic discovery of Collective Variables. As discussed, many enhanced sampling methods critically rely on the appropriate choice of CVs. It is thus important to clarify what is meant by ‘good’ collective variables [86]. Good CVs should correspond to the important dynamical motions of the system so that a simulation biased along these CVs should overcome the free-energy barriers (rarely sampled in unbiased simulations) separating the relevant functional states. Given the complexity of biomolecular systems, identifying the good CVs through chemical intuition is often a daunting task. This is why a number of numerical techniques have been designed for the automatic identification of good CVs. A first group of techniques seeks to project a MD trajectory in a lower dimensional space spanned by base vectors aligned along the directions of maximal variance. This is the case, for example, of the popular Principal Component Analysis (PCA) technique [351]. A second class of methods looks for a low-dimensional CV space spanned by vectors aligned along the directions of maximal auto-correlation. These directions correspond to the slowest molecular motions most likely capturing the functional dynamics of the protein. A typical representative of this class of techniques is the Time-lagged Independent Component Analysis (TICA) [352]. Both families of techniques can benefit from the fitting power of deep neural networks.

An example of application of neural networks to the discovery of high-variance CVs is the Molecular Enhanced Sampling with Autoencoders (MESA) [353]. The MESA method requires a neural network with a special architecture including an encoder that maps the 3N-dimensional vector of cartesian coordinates $\xi$ to a lower dimensional vector of CVs, $\mathit{x}$, and a decoder that approximates the reverse mapping from $\xi$ to $\hat{\mathit{x}}$. The network is trained to reconstruct its input, that is, the weights are tuned so that $\hat{\mathit{x}}\cong \mathit{x}$. The algorithm performs cycles of CV identification and biased simulations using the discovered CVs until the value of the CVs remains stable. The method was successfully applied to the alanine dipeptide and Trp-cage systems [353,354].

We now provide an example of discovery of maximally autocorrelated CVs. This approach was originally developed by Ferguson and coworkers [87] who named it State-free Reversible Variational Approach to Markov Processes Networks, and further improved by Parrinello’s group (deepTICA) [88]. The approach is based on the observation that a Molecular Dynamics simulation can be considered as a dynamical process that evolves a density distribution $p_t(\mathit{x})$ at time $t$ towards the Boltzmann distribution $\mu (\mathit{x})$. Mathematically, the evolution of the density distribution is performed by the transfer operator ${\mathcal{T}}_{\tau }$. The eigenfunctions associated to the largest eigenvalues of the transfer operator are characterized by the slowest relaxation to equilibrium and are thus good CV candidates. In complex biomolecular systems, the analytical diagonalization of the transfer operator is not possible. However, exploiting a variational approach, it is possible to compute lower bound approximations ${\tilde{\lambda }}_i$ of the true eigenvalues ${\lambda }_i$ and the corresponding approximate eigenfunctions ${\tilde{\psi }}_i$. To further simplify the calculation, the variational eigenfunctions ${\tilde{\psi }}_i$ are expressed as a linear combination of trial descriptors $d_j(\mathit{x})$

(42) ${\tilde{\psi }}_i=\sum _j{\alpha }_{ij}{d}_j(\mathit{x})$(42)

and the set of parameters $\{{\alpha }_{ij}\}$ is sought that maximizes their autocorrelation (thus also maximizing ${\tilde{\lambda }}_i$). This leads to a generalized eigenvalue problem

(43) $C(\tau ){\alpha }_i={\tilde{\lambda }}_{i}C(0){\alpha }_i$(43)

where

(44a) $C_{ij}(\tau )=⟨d_i({\mathit{x}}_t)d_j({\mathit{x}}_{t+\tau })⟩$(44a)

(44b) $C_{ij}(0)=⟨d_i({\mathit{x}}_t)d_j({\mathit{x}}_t)⟩$(44b)

The neural network is fed with pairs of time-lagged descriptors $\mathit{d}({\mathit{x}}_t)$, $\mathit{d}({\mathit{x}}_{t+\tau })$ that are mapped to their latent space counterparts ${\mathit{h}}^{\mathit{\theta }}(\mathit{d}({\mathit{x}}_t))$ and ${\mathit{h}}^{\mathit{\theta }}(\mathit{d}({\mathit{x}}_{t+\tau }))$, where $\mathit{\theta }$ is a set of parameters. The solution of the generalized eigenvalue problem using the latent space variables leads to an estimate of the ${\tilde{\lambda }}_i$ parameters. The network is then trained to minimize the loss function $\mathcal{L}=-{\sum }_i{\tilde{\lambda }}_i^2(\mathit{\theta })$. Once the loss function is optimized, the deep-TICA collective variables can be computed as $s_i={\sum }_j{\alpha }_{ij}{h}_j^{\mathit{\theta }}$ and an enhanced sampling simulation can be run biasing these CVs.

8. Conclusions

The discussion in this review shows that simulation methods are powerful tools for the analysis of biomolecular systems and, in particular, for ion channels. Specifically, molecular dynamics simulations just need three ingredients: an initial state, a reliable force field, and suitable algorithms for the numerical integration of Newton’s equations of motion. In principle, these few elements should suffice to fully reproduce the dynamical behaviour of biomolecules and their function, in agreement with the 1812 statement by Pierre Simon de Laplace [355]:

… an intelligence which could comprehend all the forces by which nature is animated and the respective positions of the beings which compose it, if moreover this intelligence were vast enough to submit these data to analysis … to it nothing would be uncertain, and the future as the past would be present to its eyes.

Unfortunately, as discussed above, both continuum and atomistic methods suffer from a number of limitations and Laplace’s optimism ought to give way to a more cautious attitude along the lines of the 1929 statement by Dirac [356]:

The underlying physical laws necessary for the mathematical theory of a large part of physics and the whole of chemistry are thus completely known, and the difficulty is only that the exact application of these laws leads to equations much too complicated to be soluble. It therefore becomes desirable that approximate practical methods of applying quantum mechanics should be developed, which can lead to an explanation of the main features of complex atomic systems without too much computation.

What Dirac said about quantum mechanics also applies to classical molecular simulations and still resists the test of time, with the notable exception that much computer power seems always in high demand. This review, highlighting strengths and weaknesses of the different methods in computational ion channels research, can guide the reader in the choice of the most appropriate approach to the problem of interest. Moreover, we selected a few rapidly developing research lines ranging from hybrid enhanced sampling methods to multi-scale algorithms and machine learning that will enable to by-pass the problem of ‘too much computation’. While a perfect knowledge of ‘the future as the past’ is still precluded to mankind even in the restricted area of quantitative sciences, a clever application of the approaches discussed here will certainly push forward the boundary of our quantitative understanding of biology.

Acknowledgments

We thank G. Ciccotti for inspiring discussions at FICN2021 and for encouragement. MLB is grateful to W.A.T. Gibby, D.G. Luchinsky, P.V.E. McClintock, and A. Smolyanitsky for valuable discussions. LC, GC, and LM thank T. Malliavin for the long and fruitful collaboration on the $\alpha$7 nicotinic receptor. CG and AG thank F. Costa for useful discussions on the hERG channel.

Disclosure statement

No potential conflict of interest was reported by the author(s). Nundefined

Additional information

Funding

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme [grant agreement No 803213]. MLB has been partially supported by a Ph.D. Scholarship from the Faculty of Science and Technology of Lancaster University, UK, by the Engineering and Physical Sciences Research Council UK (grants EP/M016889/1 and EP/M015831/1), and by a Leverhulme Trust Research Project Grant RPG-2017-134. BR is supported by the National Institutes of Health (NIH) through grant R0- GM062342

Notes

1. We recall that the Reynolds number $\mathrm{R}\mathrm{e}=\rho vL/\eta$ quantifies the relative magnitude of inertial to viscous forces.

2. These equations refer to a system with constant number of particles, volume, and energy. Interactions of the (small) simulated system with external thermal baths (constant temperature), pistons (constant pressure), mass exchange, etc. can be accounted for by suitably modifying the equations of motion such that the correct probability distribution is sampled [357].