Insights into channel dysfunction from modelling and molecular dynamics simulations

ABSTRACT Developments in structural biology mean that the number of different ion channel structures has increased significantly in recent years. Structures of ion channels enable us to rationalize how mutations may lead to channelopathies. However, determining the structures of ion channels is still not trivial, especially as they necessarily exist in many distinct functional states. Therefore, the use of computational modelling can provide complementary information that can refine working hypotheses of both wild type and mutant ion channels. The simplest but still powerful tool is homology modelling. Many structures are available now that can provide suitable templates for many different types of ion channels, allowing a full three‐dimensional interpretation of mutational effects. These structural models, and indeed the structures themselves obtained by X‐ray crystallography, and more recently cryo‐electron microscopy, can be subjected to molecular dynamics simulations, either as a tool to help explore the conformational dynamics in detail or simply as a means to refine the models further. Here we review how these approaches have been used to improve our understanding of how diseases might be linked to specific mutations in ion channel proteins. This article is part of the Special Issue entitled ‘Channelopathies.’ HighlightsStructural templates are now available for many ion channel proteins.Homology modelling can be used to generate models of human ion channel proteins.The effects of disease‐causing mutations can be understood in a structural context.Molecular dynamics simulations can provide deeper insight into such mutations.


Introduction
Channelopathies are defined as disease states where the primary cause is a dysfunctional ion channel and include conditions such as epilepsy, pain syndromes, cystic fibrosis, cardiac arrhythmias, myotonia and many others (Ashcroft, 2006).
Many of these conditions can be attributed to mutations that lead to alteration of ion channel function in some way (Cannon, 2007).If we are to be able to offer some kind of drug-based treatment for these conditions in the future, a complete molecular understanding of the effects of the mutations will provide the most rational footing from which to proceed.Many single nucleotide polymorphisms (SNPs) associated with diseases of ion channel function are fatal, for example where the mutation means that the channel or auxiliary protein will simply not even get expressed or correctly trafficked to the appropriate membrane.We will not discuss those mutations.Rather, we will focus on cases where molecular modelling has led to the most insight.
There has been spectacular progress in the structural biology of ion channels in recent years (Du et al., 2015;Li et al., 2017a;Li et al., 2017b;Long et al., 2005; Morales- Perez et al., 2016;Park et al., 2017;Payandeh et al., 2011;Sobolevsky et al., 2009;Whicher and MacKinnon, 2016;Wu et al., 2015;Zhang and Chen, 2016) and this can provide a step change in our understanding of function, but obtaining such structures is not trivial.Even with the structure in hand there is still the problem that ion channels by definition are very dynamic entities.Thus, understanding their transitions between states also becomes necessary, especially in the context of drug-design -one surely has a better chance of designing a compound when you know which state to target.Computational approaches can provide very useful complementary tools to the structural and electrophysiological approaches that are currently applied to the study of ion channels and their mutations.Two approaches in particular have given, and will continue to give, clear insight into the effects of ion Homology modelling is the primary route when there is no structure of the channel.
The single most important factor that dictates the quality of the model is the alignment to a suitable template (a structure that is homologous and has reasonable sequence identity).Sometimes the alignment can be very obvious, for example where the sequence is highly conserved as is the case for the selectivity filter region, but if one is using a template which only has low sequence identity (less than 20%), then researchers often try to utilize additional data that may exist in literature -for example knowledge of whether residues are exposed to a water-accessible pore or are lipid exposed.Another common issue is that although the structure has been "solved", some parts, typically the more mobile loop regions, are simply not visible from the electron densities.Related to this, quite often crystallographers use a construct that has artificially short linker regions wherefore one has to decide whether to use the artificial linker or whether to try modelling the native, but usually longer and more flexible loop region.Sometimes these loop regions can be functionally important so modelling them correctly may be a high priority in some cases.
In many cases, just building the homology model itself can be remarkably insightful and lead to new hypotheses concerning function.If the model is of high enough quality and there is good confidence in it, then in some cases it can be taken forward to explore conformation stability and/or dynamics as indeed might be done with a crystal structure.The models typically used in MD simulations try to reproduce the native environment (or that of an experiment) as closely as possible.Thus, simulations typically incorporate the surrounding lipid molecules, mimicking the cell There is huge modelling literature that makes reference to mutations that are associated with channelopathies and there is not enough space to review them all here.Thus, in this review, we focus on those studies where homology modelling and MD simulations were used to obtain a direct understanding of specific mutations associated with channelopathies.We have tried to highlight how the predictive and the explanatory powers of these methods pair very well with the more traditional experimental methods used to characterize ion channel function.

Voltage-Gated Potassium Channels
Potassium channels share a common topological architecture with sodium and calcium channels (Fig. 1), but as the first ion channel structure to be solved by X-ray crystallography was of a potassium channel (Doyle et al., 1998), it is perhaps not surprising that modelling in this area is the most mature.Much of the computational literature surrounding potassium channels is focused on understanding the properties of selectivity, conduction and activation, but more recently efforts have focused on mutations associated with channelopathies.
One potassium channel family that has received a lot of modelling attention is the voltage-gated K v 11.1 (from the KCNH2 gene), better known as the hERG1 channel (Perez et al., 2016).Mutations in this gene lead to Long QT Syndrome (LQTS), a disorder that gives individuals a predisposition to life-threatening arrhythmias.An additional problem is that promiscuous binding of drugs to a pocket formed by the S6 helix (Fig. 1B) blocks conductance of the channel and leads to an acquired form of LQTS (Sanguinetti and Tristani-Firouzi, 2006), also causing potentially fatal arrhythmias.The latter is a particular problem in pre-clinical drug safety trials and is a 6 common reason why many drugs fail at that stage.Thus, there is considerable interest in this particular potassium channel, but there is currently no crystal structure.Many groups have built models of the hERG1 channels based on available crystal structure templates, but most of these studies focused on channel block, which is not the focus of this review (the interested reader is referred to recent reviews (Dempsey et al., 2014;Stary et al., 2010) for insight into this specific area).indeed that is where the channel blocking effects of promiscuous drugs were thought to interact.However, in terms of genetic predisposition there may be positions away from the central pore that could still lead to disease states.One of the first works to attempt to build a full model in this regard was that reported by Subbotina et al. (Subbotina et al., 2010).These authors used a template-driven approach in conjunction with de novo design for regions missing from the template and were able to rationalize apparently conflicting experimental data particularly concerning the interaction of the voltage sensor and the pore region.Furthermore, it has become apparent that an emerging therapeutic strategy might be to develop small molecules that act as channel openers (Grunnet et al., 2008), an aspect that has also been explored via molecular modelling and MD simulations (Guo et al., 2015;Guo et al., 2014;Perissinotti et al., 2015).
In vivo, many potassium channels function as part of larger complexes.For example, KCNE subunits are known to modify the properties of several K-channels (including KCNQ1, ERF, HCN3, K v 1.1, K v 1.3, K v 2.1, K v 3.1 and K v 4.3) in ways that often reproduce currents that have native characteristics as reviewed by Strutz-Seebohm in 2011(Strutz-Seebohm et al., 2011).Interestingly, mutations in either KCNE or the KCNQ1 (K v 7.1) protein can also lead to LQTS (Abbott et al., 1999).It has been proposed that disease-associated mutations in KCNE1 or KCNQ1 might prevent the usual and necessary interactions between the two, leading to LQTS (Seebohm et al., 2001).Smith et al. (Smith et al., 2007) showed how mapping over 85 disease-linked mutation sites to models of both open and closed states of the KCNQ1 channel can be used to generate structure-based hypotheses for disease phenotypes and were able to highlight the role of certain interfaces within the structure that may serve as hot-spots for mutation.Strutz-Seebohm et al. (Strutz-Seebohm et al., 2011)  containing four α-subunits (based on the K v 1.2 crystal structure (Long et al., 2005)).
The KCNE1 transmembrane (TM) domain of the β-subunit was modelled de novo as an α-helix and several conformers were docked to the homology model of the channel.The conformation that gave the best agreement with the double mutant cycle data was simulated with MD to examine the conformational stability and indeed was found on this timescale (50 ns) to be a stable conformation.The model suggests that the closed state would be further stabilized, as indeed observed by slower activation rates, in the heteromeric complexes.The model's stoichiometry of four KCNQ1 α-subunits and two KCNE1 β-subunits is compatible with reports of varying stoichiometry (more or fewer β-subunits could exert similar effects).Although no direct data for the mutant was presented in this study, it certainly provides a working model with which to rationalize the effects of SNPs in the S4-S5-S6 region.
The KCNA (K v 1.1) gene has been linked to the autosomal dominant potassium channelopathy episodic ataxia type 1 (EA1).The C185W SNP in the conserved region of the S1 voltage sensor helix has been examined (D'Adamo et al., 2015) and interpreted with structural modelling (a homology model based on the K v 1.2/2.1 chimera (Long et al., 2007)) as a template followed by MD simulations (Abraham et al., 2015)).The authors suggest that the substantial kinetic changes arising from the mutation are due to steric clashes created by the presence of the bulkier tryptophan at this position.
Mutations in the voltage sensor region of K v 7.2 are associated with neonatal epilepsies.Miceli et al. (Miceli et al., 2013) (Yang et al., 2013) in their study of K v 10.2 showed that the R327H mutation (also linked to epilepsy) stabilizes the channel in the closed state.Thus, changing the stability of the closed, open or inactivated states appears to be an emerging theme of many mutations in K v channels that have been linked to epilepsy.

Voltage-Gated Calcium Channels
The point mutation I745T in the S6 segment of the L-type calcium channel Ca v 1.4 causes severe visual impairment by causing a large negative shift (-34 mV) in the voltage dependence of the channel activation (Hemara-Wahanui et al., 2005).
Moreover, this mutant channel exhibits slower deactivation kinetics than the wild type (WT), perhaps by reducing the energy barrier necessary to open the channel (Hemara-Wahanui et al., 2005).Interestingly, several amino acid substitutions of the equivalent residue in another L-type calcium channel, Ca v 1.2, also resulted in negative shifts in the voltage dependence and slower deactivation kinetics (Hohaus et al., 2005).Altogether, this suggests that a common gating mechanism involving the channel-lining S6 segment of these proteins may exist.To study that hypothesis, Stary et al. (Stary et al., 2008) performed MD simulations of the WT and mutants of a homology model of Ca v 1.2 embedded in a lipid bilayer.The comparison between the simulations of the WT and systematic mutations of residues of the S6 segment in the four domains of Ca v 1.2 revealed the dependency of this protein on the hydrophobic interactions between its S6 segments for pore gating.
This hypothesis was further explored in the study of the de novo mutations G402S and G406R, located in the S6 segment of Ca v 1.2 (Depil et al., 2011), which are known to lead to the Timothy syndrome, a channelopathy consisting of physical,

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
10 neurological and developmental defects (Splawski et al., 2004).These missense mutations have been found to cause a reduction of the voltage-dependent inactivation, resulting in prolonged calcium entry.To gain some structural insights into the functional effects of these two mutations, a different homology model of an open Ca v 1.2 channel was used (Kudrnac et al., 2009), and compared to an equivalent model of the same channel in its closed state (Depil et al., 2011).The visualization of the open and closed states of the model of Ca v 1.2 suggested that the G402S and G406R mutations altered the hydrophobicity of a conserved motif of small hydrophobic residues located in the lower section of the four S6 segments of the channel.As this segment constitutes the activation gate, the authors hypothesised that mutations altering this tight "hydrophobic seal" located in the S6 segments would lead to higher channel opening rates.Indeed, such structural observations could perhaps be extended to the other L-type calcium channels due to their high sequence similarity.

Voltage-Gated Sodium Channels
Mutations in voltage-gated sodium channels (Na v s) have been linked to various disorders including epilepsy and heart arrhythmias (Imbrici et al., 2016b).Na v s are also implicated in different forms of chronic pain.Na v 1.7 (SCN9A) mutations in particular have been linked to various different pain disorders, from congenital insensitivity to pain through to inherited erythromelalgia and paroxysmal extreme pain disorder (Lampert et al., 2010).
Despite the lack of structures for human Na v s, several studies employing Na v homology models based on the bacterial KcsA potassium channel (Doyle et al., 1998) or the bacterial Na v Ab channel (Payandeh et al., 2011) (Fig. 2A) have added to our understanding of Na v based channelopathies.Furthermore, chimeric structures formed by a combination of bacterial and human sequences have shown to be useful in the design of new treatments for pain (Ahuja et al., 2015) and may serve as useful templates in the future.(Lomize et al., 2006).(B+C) Model of the hydrophobic gate in Na v 1.7 with the WT F1449 (B) and the mutant F1449V (C), causing an increased channel diameter.Na v channels are formed from a single polypeptide chain (Fig 1C), and in the WT channel, the hydrophobic gate is formed by a tyrosine residue from DI and three phenylalanine residues from each of DII, DIII and DIV.F1449V is found in DIII.The models are based on the KcsA structure (Doyle et al., 1998), only the channel-lining S6 helix is included and the model is viewed from the intracellular side.
The F1449V mutation in Na v 1.7 causes inherited erythromelalgia and has been suggested to increase channel activity by increasing the channel diameter in the closed state (Lampert et al., 2008).A homology model of the Na v 1.7 pore domain in the closed state was generated using the KcsA potassium channel structure (Doyle 12 et al., 1998) as a template for both the WT and different F1449 mutants.F1449 is found on the S6 helix of the third domain and faces into the pore (Fig. 2B).
Therefore, mutations of the large and hydrophobic F1449 side chain were suggested to affect the pore diameter and alter the stability of the closed state (Lampert et al., 2008).Mutation to a smaller side chain would be expected to increase channel activation, whereas mutation to a larger side chain is expected to reduce channel activation.Models were constructed for WT, F1449V, F1449L, F1449Y and F1449W and the pore diameter was measured for each model.The WT model showed that F1449 forms part of a hydrophobic gate (Fig. 2B) with a WT channel diameter of 0.6 Å, compared to 1.1 Å for the F1449V mutant (models in Fig. 2B and 2C).For F1449L the diameter was 0.76 Å whereas for F1449Y and F1449W the pore diameter is even smaller than WT at 0.5 Å and 0.45 Å, respectively (Lampert et al., 2008).Thus, the side chain at this position can contribute to regulating the pore diameter and was predicted to alter the stability of the closed state (Lampert et al., 2008).These predictions were confirmed by experiment.The F1449V opened faster than WT, F1449L activated similarly to WT, and both F1449Y and F1449W took longer time to reach peak current than WT.The F1449V mutation (but not the others) also appeared to open at more negative potentials, which could contribute to its association with erythromelalgia.Thus, this mutant has a less stable closed state and opens faster that WT (Lampert et al., 2008).
Another Na v 1.7 mutant, A1632T, which also causes inherited erythromelalgia, did not alter the voltage dependence for channel opening.However, steady-state fast inactivation was observed at more depolarised potentials than for the WT channel (Eberhardt et al., 2014).A homology model based on the Na v Ab channel structure (Payandeh et al., 2011) illustrated that A1632 is located in the S4-S5 linker of the fourth domain, so the most likely explanation for the effects of the mutation is that it alters the dynamics of inactivation (Eberhardt et al., 2014).

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
13 Na v 1.7 mutants involved in pain syndromes are often resistant to pharmacotherapy.
However, a study has illustrated how combining structural modelling and thermodynamic analysis might aid in the identification of mutants with increased response to known drugs (Yang et al., 2012).Using the Na v Ab channel (Payandeh et al., 2011) as a template, it was clear from a homology model that the S241T mutation (in the S4-S5 linker of the first domain) was only 2.4 Å away from the V400M mutation (in S6 of the first domain), which is known to be susceptible to modulation by carbamazepine (Yang et al., 2012).Thermodynamic mutant cycle analysis was used to explore this further.The experiments supported a shared mechanism between the two mutations and suggested that S241T may also be affected by carbamazepine.Indeed, carbamazepine depolarised the voltage dependence of S241T activation and reduced the hyperexcitability of Na v 1.7 S241T neurons on a cellular level (Yang et al., 2012).
Mutations in Na v 1.4 and Na v 1.5, (SCN4A and SCN5A genes, respectively), have also been implicated in various channelopathies.The I141V mutation (in the S1 helix of the first domain of both channels) has been studied both by MD simulations and functional experiments (Amarouch et al., 2014).Functional experiments showed that I141V in both Na v 1.4 and Na v 1.5 caused a shift to more negative potentials in the voltage dependence of activation.Homology models of Na v 1.4 WT and I141V were constructed based on the Na v Ab structure and embedded into a membrane bilayer and studied by MD simulations for 100 ns (Amarouch et al., 2014).
The Na v Ab structure was solved with the voltage sensor in an activated conformation but the channel in a closed conformation (Payandeh et al., 2011).R225, one of the conserved positive charges in the S4 helix of the voltage sensor, is located very close to V141 in the models, believed to illustrate an activated voltage sensor.For

M A N U S C R I P T A C C E P T E D
ACCEPTED MANUSCRIPT 14 the mutated structure, a hydrogen bond between the Y168Hη atom and the backbone oxygen atom of R225 was seen to form during the simulations.This hydrogen bond was not observed in the WT simulations, however, the reduced steric hindrance from the valine side chain relative to the original isoleucine appears to be enough to allow the formation of this extra hydrogen bond.The active state could therefore be stabilised by this additional hydrogen bond in the I141V mutant, increasing the activation of the channel in agreement with experiments (Amarouch et al., 2014).However, mutations on Y168 designed to test the hypotheses further were difficult to completely reconcile with the simulations (Amarouch et al., 2014).
Mutations in Na v 1.5 are associated with phenotypes involving arrhythmias and dilated cardiomyopathy.Moreau et al. studied the effect of two Na v 1.5 mutations, R222Q and R225W, both found on S4 in the voltage-sensing domain of the first domain (Moreau et al., 2015).Again, homology models based on the Na v Ab structure (Payandeh et al., 2011) were generated but in three different functional states (Moreau et al., 2015).To generate models of the different functional states of the voltage sensor, the alignment for S4 was manually shifted by zero (state α), three (state β), and six (state γ) residues towards the N-terminal, with α being the most activated state and γ being the state closest to deactivation.The different models were embedded in a lipid bilayer, solvated, and their dynamics studied for 100 ns using MD simulations.
Electrophysiological experiments with the two mutants suggested that both might allow for so-called gating pore current, since the observed currents were also present in the presence of tetrodotoxin, a neurotoxin and potent channel blocker (Moreau et al., 2015).In each of the three states for the WT, studied by MD simulations, the water extended from the membrane surface and down to a constriction zone, approximately in the middle of the membrane, and water also protruded somewhat 15 into the membrane from the cytoplasmic side.In each state one of the conserved charged residues on S4 were located at this constriction zone; R228 in α, R225 in β and R222 in γ.The introduction of R222Q in the γ state and R225W in the β state disrupted the salt bridge stabilising these residues in the constriction zone.Following from this disruption, the four-helix bundle expanded and a solvated pathway was observed all the way through the voltage-sensing domain (Moreau et al., 2015).
Thus, the simulations support the existence of gating pore currents in the R222Q and R225W Na v 1.5 mutants, and these cation leaks, in principle, reflect gain of function and thus increased channel activity, and are consistent with the formation of a channel that cannot be blocked by tetrodotoxin.In the simulations, the pore was wider in the R225W mutant than in the R222Q mutant, potentially further explaining the experimentally measured larger current through the R225W relative to the R222Q mutant (Moreau et al., 2015).
Mutants of Na v 1.9 (SCN11A gene) that are implicated in pain have also been identified, e.g., the R222H mutant (Han et al., 2016).R222 is the first arginine on the voltage-sensing S4 helix of Na v 1.9, and fully conserved among human Na v s.To gain insight into the structural explanation for the increased channel activity, Han et al. performed multistate modelling of the Na v 1.9 channel (Han et al., 2016).Homology models of Na v 1.9 in different states were constructed based on templates extracted from earlier MD studies of voltage-gated potassium channels.The modelled states included, in the order of appearance in the functional cycle, an activated state, early loss-of-activation state, late loss-of-activation state, resting/closed state, early activation state and late activation state, illustrating the channel opening and closing movements and the change in interactions of important residues between these different states (Han et al., 2016).
Arginine is positively charged at physiological pH, whereas histidine most likely exists in a predominantly neutral state at physiological pH.R222 was observed to be involved in charge-charge interactions between the voltage-sensing helix and other helices in the voltage-sensing domain, and therefore mutation to a predominantly neutral side chain may have large consequences for these interactions and thus affect gating of the channel.The multistate modelling procedure suggested that R222 was interacting with negatively charged residues on the S2 and S3 helices in all states except for the late activation step and the fully activated state (Han et al., 2016).Thus, the presence of an arginine most likely stabilises the states in which R222 interacts with a negatively charged side chain, whereas substituting the arginine for the more neutral histidine side chain is likely to favour the two open states in which residue 222 is not interacting with a negatively charged side chain, i.e. the late activation state and the activated state.Overall, R222H should therefore show an increased stability of the active state, which is in agreement with the experiments showing that the R222H mutant opens more readily than the WT channel (Han et al., 2016).

Potassium Inward-Rectifier (K ir ) Channels
K ir 1.1 channels are expressed in the kidney, and their function is to secrete excess K + ions into the urine in response to changes in intracellular pH at physiological levels.Mutations in the K ir 1.1 pH-sensing domain are known to cause type II Bartter syndrome, a channelopathy characterised by hypokalaemia (Hebert, 2003).
Rapedius and co-workers have investigated the role of the 'RKR' triad using MD simulations of a K ir 1.1 model, paying special attention to the pK a values of these residues (Rapedius et al., 2006).The resulting simulations showed that K80 (of the RKR triad) takes part in transmembrane inter-helix interactions with A177, controlling the ability of the channel to gate in response to pH change.Interestingly, the mutation A177T also causes Bartter disease in an undefined manner (Peters et al.,

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
17 2003), while the mutation of the equivalent residues in K ir 6.2 causes neonatal diabetes (Ashcroft, 2005).Therefore, such protein segments seem not only structurally important but also constitute a "hotspot" for channelopathies.Additionally, the triad arginine residues were also observed to maintain highly conserved interand intra-subunit interactions (Rapedius et al., 2006).Therefore, despite not constituting the pH-sensing domain itself, the RKR triad was found to be indirectly involved in pH-related channel gating.
The K ir 6.2 channel is an ATP-sensitive channel and mutations surrounding the ATPbinding site have been associated with neonatal diabetes and hyperinsulinemia (Ashcroft, 2005).Moran et al. (Moran et al., 2013) calculated the free energy of binding of ATP to WT and mutant K ir 6.2 channels.By combining these results with experiments designed to explore the influence of mutations on the free energy, the authors were able to propose the most likely mode of ATP binding.The results should help to understand the origin of ATP-sensitivity for this channel.
The K ir 6.2 channel is also regulated by phosphoinositol bisphosphate (PIP 2 ).Some mutations associated with neonatal diabetes lie close to the PIP 2 binding site and this aspect has been explored with homology modelling and MD simulations (Haider et al., 2007).The authors were able to rationalize how mutations affecting the PIP 2 binding might indeed give rise to impaired insulin secretion.Homology models have also been used to help explain the E23K polymorphism in K ir 6.2 that has been associated with congestive heart failure (Reyes et al., 2009) and also a heterozygous mutation, W68R, in the K ir 6.2 subunit, found in a patient with transient neonatal diabetes (Männikkö et al., 2011).

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
18 Mutations in the CNGA3 and CNGB3 subunits of the cone photoreceptor CNG ion channels have been linked to a reduction in cone function (achromatopsia).For example, D277N in CNGB3 results in day blindness in dogs and this mutation is located on the same locus where human achromatopsia mutations have been identified.Sequence analysis has revealed that this aspartate residue is located in the second transmembrane domain S2, and is surrounded by a further two conserved aspartate residues (Tanaka et al., 2014).As mutagenesis studies showed that the mutation of these residues results in incorrect localization of the channels, it was hypothesized that these D/N mutations might affect inter-subunit assembly.To further investigate the structural consequences of mutation of the Tri-Asp motif, Tanaka and colleagues (Tanaka et al., 2014) constructed a homology model of the WT transmembrane domain using the K v 1.2/2.1 channel (Long et al., 2007) as a template.The model was then embedded in a lipid bilayer and simulated for 65 ns at physiological conditions.Throughout the simulation the authors observed a number of salt bridges between these aspartate residues and positively charged residues located in S3/S4 and they concluded that the substitution with asparagine disrupts such interactions, which appeared to be key for correct folding.The recent appearance of a 3.5 Å resolution single-particle electron cryo-microscopy structure of a CNG channel from Caenorhabditis elegans (Li et al., 2017a) will undoubtedly lead to further modelling and simulation efforts.

Transient Receptor Potential Ankyrin (TRPA) Receptor
The gain-of-function N855S mutation located in the S4-S5 linker of the TRPA-1 receptor has recently been related to a familial human pain disorder (Kremeyer et al., 2010).To understand the effect of this mutation on the TRPA-1 channel structure and dynamics, Zíma et al. (Zíma et al., 2015) produced a homology model of the protein and performed MD simulations.These authors proposed that it is a change 19 in the electrostatic environment experienced by nearby charged residues that causes a depolarized shift in the voltage-dependence of the channel and that this is what governs the behavior of the N855S mutation.To that end they performed various in silico and electrophysiological studies to examine this underlying hypothesis.For example, after a 5-ns MD equilibration, the mutation N855R was introduced manually and simulated for another 15-30 ns within a POPC bilayer.During both WT and N855R mutant simulations a stable inter-subunit salt bridge was formed between the highly conserved and neighbouring residues E854 and K868, leading the authors in turn to proceed with the electrophysiological characterisation of the charge reversal E854R and K868E mutants.They observed a reduced channel response upon chemical and voltage stimuli.However, although this data indicates that these neighbouring electrostatic residues might be involved in sensing voltage changes in the membrane, the observed increase in Ca 2+ -induced currents observed upon mutation of the N855 residue led the authors to suggest that the N855S mutation might increase Ca 2+ activation efficacy rather than increasing ion channel function as such.It will be interesting to see if similar gain-of-function mutations emerge that are located in this region.

Cys-Loop Receptors
Even though the Cys-loop family of receptors is quite large and has been implicated in many channelopathies (Ashcroft, 2006), the absence, until recently, of high resolution structural data has precluded in-depth analysis of SNPs that might affect the gating properties.Indeed, one of the largest problems with Cys-loop receptors for structural biology is actually assigning the corresponding physiological state (i.e. whether the structure represents a resting or desensitized state).Modelling of channelopathies has therefore mainly remained at the level of building homology

M A N U S C R I P T A C C E P T E D
ACCEPTED MANUSCRIPT 20 models to gain some kind of structurally informed hypothesis that can be tested experimentally.
For example, several studies have used homology modelling (using GluCl as a template) to map GlyR mutations related to hyperekplexia in order to investigate the structural impact of mutations (James et al., 2013) and to find inspiration for biochemical experiments probing the surrounding residues (Lynagh et al., 2013).
Furthermore, the Torpedo muscle-type nicotinic acetylcholine receptor (nAChR) was used as a template (Chung et al., 2013;Chung et al., 2010) to create models of WT GlyR and a range of hyperekplexia mutants identified from systematic DNA sequencing of the GLRA1 gene in 88 unrelated human hyperekplexia patients.The models were used to assess the structural impact of the mutations and to correlate this with findings from functional assays.Modelling the positions of such mutations can provide useful working hypotheses (Bode and Lynch, 2014).With the recent structural efforts (Du et al., 2015;Huang et al., 2015;Huang et al., 2016), we now have very good starting points for more elaborate investigations of disease-related mutants.For example, (Safar et al., 2017) used the zebrafish α 1 GlyR structure as a template for the human α 1 GlyR.A previous model based on GluCl had revealed that the side chain of E103, which when mutated to a lysine results in hyperekplexia, is located close to the side chain of R131, and that these might form a salt bridge at the back of the binding site.The new model was used to perform MD simulations that confirmed the presence and conformational stability of the proposed salt bridge.
Arias et al. tried to explore how some treatments for slow-channel congenital myasthenic syndrome (a disease of the nAChR receptor), such as the open channel blocker fluoxetine, can be rationalized via docking to a model of the muscle-type nAChR (Arias et al., 2010).Using the α 1 subunit from the same receptor as a template, Changeux and colleagues (Taly et al., 2006)  or at the interface between rigid blocks within a subunit.These interfaces are significantly modified during a quaternary twist motion believed to be important for the channel gating, which might explain how these mutations exert their effect.The recent structure of the α4β2 nicotinic acetylcholine receptor (Morales- Perez et al., 2016) should provide an improved template on which to base future models on.
Genome sequencing has revealed that there are a substantial number of mutations in GABA A receptor genes that are associated with a range of diseases such as epilepsy, autism, schizophrenia and addiction (Yuan et al., 2015).A recent structure of the human β 3 GABA A receptor allowed the authors to map mutations related to epilepsy, febrile seizures and insomnia to several distinct structural elements throughout the structure (Miller and Aricescu, 2014).With the recent structure at hand, we can expect to see more modelling studies in this area in the near future.

The Cystic Fibrosis Transmembrane Regulator (CFTR)
The cystic fibrosis transmembrane conductance regulator, or CFTR, belongs to the ATP Binding Cassette (ABC) superfamily of proteins.Although the vast majority of ABC proteins are transporters, CFTR is an ATP-gated chloride channel (Linsdell, 2006;Muallem and Vergani, 2009), which makes it unique in this large superfamily.
Mutations in CFTR can cause cystic fibrosis, which is one of the most common lethal autosomal recessive disorders in the Caucasian population.The most common mutation is a deletion of F508 located in the first nucleotide-binding domain (NBD1) (Fig. 3).The F508del mutant protein is generally misfolded and targeted to the endoplasmic reticulum-associated degradation pathway.However, the few molecules   (Mornon et al., 2015).The locations of some of the disease-causing mutations in the NBD1 of CFTR (I507, F508, G551) are shown in red surface representation.(C) A closer look at the location of the disease-causing mutations (I507, F508, G551) mentioned in this review.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
23 Until recently (Zhang and Chen, 2016), only individual NBDs have been solved to atomistic resolution (Atwell et al., 2010;Lewis et al., 2003;Lewis et al., 2010;Lewis et al., 2005;Thibodeau et al., 2005).Nonetheless, these have been used in several theoretical studies to compare the stability and conformational changes of WT and mutant NBDs by molecular modelling and MD simulations.
In a study that used 2 ns long MD simulations of WT and F508del CFTR, the authors observed noticeable changes in the intra-domain dynamics between the two systems (Warner et al., 2007).However, docking of the small ligand CPX (8-cyclopentyl-1,3dipropylxanthine) into the mutant NBD1 changed the profile of dynamic movements to be closer to the WT dynamics.Another study of longer (20 ns) MD simulations of WT and F508del NBD1 revealed that the mutant shows more conformational flexibility, which led the authors to suggest that the mutant protein is targeted by the degradation pathway due to a larger amount of exposed hydrophobic residues (Wieczorek and Zielenkiewicz, 2008).In yet another study, Bisignano and Moran did not observe any major differences between WT and F508del NBD1 in terms of changes in entropy and enthalpy, but found that the mutant had a less favourable solvation energy and a significantly lower (30-fold) binding energy with ATP (Bisignano and Moran, 2010).
A graph theory approach was used on data from discrete MD simulations of WT, F508del and I507del NBD1 to determine which residues form dynamic coupling between structurally distinct areas of CFTR.S492 was identified as a key residue in allosteric interactions between the Regulatory Insertion (RI) and the F508del and I507del mutations (Proctor et al., 2015).Replica Exchange MD simulations of WT and F508del NBD1 showed that the mutant exhibits higher overall fluctuations, causing exposure of hydrophobic residues in the vicinity of F508, reduced ATP binding capacity and impairment of interactions with other domains (Zhenin et (Aleksandrov et al., 2012).Since F508del causes misfolding, its folding pathways have been explored using MD simulations.One such study observed that F508del affects folding dynamics of CFTR since the metastable states that appeared along the folding pathways were differently populated for the WT and mutant NBD1 even in the presence of correcting mutations (Serohijos et al., 2008b).Similar conclusions were reached in a study that used high-temperature MD simulations to unfold the WT and mutant NBD1 (Estácio et al., 2016).Estácio et al. observed the formation of unique stable intermediate states that were misfolded.The second-site rescue mutations did not completely block the formation of the misfolded states, suggesting only limited potential for correction.
The availability of the full-length ABC exporter structures (Dawson and Locher, 2007;Ward et al., 2007) allowed the construction of models of the transmembrane domain (TMD):NBD assembly.The first models suggested that atoms from the main and side chain of F508 interacted with intracellular loop (ICL) 4 (Mornon et al., 2008;Serohijos et al., 2008a).Both models were based on the Sav1866 crystal structure (Dawson and Locher, 2007) in the outward-open state, which is thought to represent the open channel state for CFTR.Later, homology models of both human CFTR in open and closed channel states were created using MsbA crystal structures (Ward et al., 2007) as templates (Mornon et al., 2009).No changes between the two states were observed in the TMD:NBD coupling interfaces that include F508.

25
The dynamics of the MSD:NBD assembly were assessed in MD simulations for the WT and two mutants (F508del and G551D) (Belmonte and Moran, 2015).While no large conformational changes were observed for the MSDs and NBDs, the interactions between the NBDs were affected in both mutants.While F508del showed an increase in the amount of ubiquitously dispersed contacts between NBDs, G551D had less internal contacts between the domains along with altered interactions with the ICLs, potentially explaining the gating defect caused by this mutation.In another study, rapid conformational changes lead to the first full-open channel state, which was reached after the F508-bearing alpha-subdomain of NBD1 moved to an alternative stable position (Mornon et al., 2015).This observation potentially highlights the role of F508 in the protein, since opening and closing of the channel likely require only localised changes.Analysis of the full-open model also revealed that most well-characterised missense mutations are clustered in hotspots such as pore-lining residues, the ICL bundle, the NBD1:ICL4 interface and the canonical ATP binding site.

ClC Chloride Channels
Dietary sodium chloride intake sensitivity is present in approximately half of hypertensive patients (Weinberger, 1996), with renal Cl -and Na + homeostasis playing a major role in hypertension pathogenesis (McCallum et al., 2015).In particular, two kidney ClC proteins, namely ClC-Ka and ClC-Kb, together with the accessory Bartter protein, govern chloride absorption and urine concentration.Accordingly, mutations in these proteins can lead to reabsorption problems or Bartter syndrome.
Until recently (Park et al., 2017), there was no mammalian structure of a chloride channel and most modelling efforts used the 3.5 Å resolution X-ray structure of a bacterial ClC (Dutzler et al., 2002) chloride channel as a template.For example

M A N U S C R I P T A C C E P T E D
Gradogna et al. (Gradogna et al., 2012) generated a model of ClC-Ka based on ClC and subsequently used it to perform a docking campaign over a new series of highaffinity benzofuran-derivative compounds (Liantonio et al., 2008).An iterative docking protocol allowed the authors to hypothesize a putative binding site, which was further tested using MD simulations of the best-scoring poses.Importantly, the resulting binding site did not overlap with two common polymorphisms associated with hypertension, rationalizing why the drug inhibition level was well maintained in ClC-Ka proteins that present such genetic variants.Furthermore, MD simulations were also employed to rationalise the results of multidisciplinary studies of mutations affecting the voltage-gated CLC-1 chloride channel (Imbrici et al., 2016a;Imbrici et al., 2015).In particular, the effects of dominant mutations F484L (Imbrici et al., 2015) and T335N (Imbrici et al., 2016a) have been characterised through medical diagnosis, patch-clamp electrophysiology, Western blot and quantitative PCR.In these studies, both homomeric and heteromeric channels were shown to maintain the same expression levels as WT, while showing a reduction of the chloride currents with a positive increase on the dependency on the voltage applied, thus reducing the probability of channel opening at physiological voltage levels.To gain insights into the structural basis of this reduced channel function, each point mutation was introduced in a homology model.
Both WT and mutants were embedded in a POPC bilayer and simulated for 25-50 ns, finding that both mutations induced changes in the H-bond network: while F484L seemed to indirectly facilitate a closed configuration of the pore by creating a stable inter-dimer H-bond between the pore residues E232 and Y578 (Imbrici et al., 2015), the T335N mutation favoured intra-dimeric interactions instead, highlighting the complex conformational role of this protein region on channel gating (Imbrici et al., 2016a).

On the interpretation of modelling data
In this review we have tried to illustrate how a variety of modelling and MD approaches have been used to provide insight into channelopathies.In some cases, a simple homology model can, as exemplified above, be very helpful in generating hypotheses that can be explored experimentally.In other cases, dynamical modelling might be needed to gain further insight.The depth of detail in the approaches can range from quite superficial through to something that could only be applied with expert knowledge.However, the validity of these approaches, and how useful they are, depends to a large extent on the question being asked and is often highly specific to a particular study.These days it is very easy to generate homology models.The single biggest factor that dictates the quality of homology models is the sequence alignment, but assessing the quality of the model afterwards is critical in deciding how much reliance can be placed upon it.The simplest measures include analyzing the Ramanchandran distribution to check there are no "violations" of known ϕ or Ψ angles or analyzing the potential energy of the model.Whilst the latter is definitely sensible, it is only really useful to provide a ranking of models.It does not provide an absolute measure of model quality.The latter can however be assessed by programs such as QMEAN (Benkert et al., 2008), although this is currently focused on soluble proteins.
In the context discussed here, modelling and MD are more often than not used in what might be described as a qualitative "hypothesis generator" rather than a truly quantitative predictor of functional changes.In that case, a simple homology model that helps a researcher visualize the location of important residues in threedimensional space can be just as useful as a more complete study that used MD simulation to refine the model and examine some of the underlying dynamic properties.The dynamic processes of proteins exist over a large timescale range (femto-second to hours and beyond), thus it is not particularly useful to provide ball-M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 28 park guidance for how long a simulation needs to be in order for it to be useful.
Typically the use of MD in the contexts above has been to investigate conformational stability, but it can provide much more quantitative insight if desired.However, the reader should ask themselves to what extent have the authors shown that the property under consideration has converged or at least if an error estimate has been provided?

Future Directions
One reason why modelling and MD approaches have been used in the way described above has been that matching experimental timescales has been difficult.This, however, is changing and simulations are starting to appear that, for some processes at least, are of comparable timescales.In addition to this, there is a growing expectation that simulations will be repeated multiple times, just as "wet" experiments are.In the past this has been difficult because the resource to do that as a matter of course was simply not available.As compute resource becomes more readily available, and multiple repeats are performed routinely, it will become easier for non-experts to judge how much weight they can assign to the modelling results.As computer power increases we can expect to see not only better assessment of errors but also more quantitative predictions per se.
Another aspect that is likely to see expansion as computer power increases, is the tighter integration of experimental data as input to drive computational methods.For example, MD can be combined with cryo-EM maps to help in the refinement process (Singharoy et al., 2016) and this is likely to become increasingly important as the use of cryo-EM also becomes more widespread.Metadynamics and similar methods allow a route to readily incorporate different types of experimental data into an MD simulation ranging from small angle X-ray scattering (Kimanius et al., 2015) to NMR data (Granata et al., 2013).

Conclusions
We have focused on simulations of proteins where the paradigm is "one mutation causes one phenotype" as this is the most straightforward place to begin modelling studies.However, as more structural information starts to appear, it is becoming increasingly apparent that an understanding of how mutations influence the stability of the different conformational states will become crucial.Modelling and simulation approaches can provide a unique tool to help connect what is observed structurally with the dynamics reported via electrophysiological experiments, especially single channel recordings.
As the field advances, we can look forward to structures of complexes that more closely reflect the in vivo situation and in particular the role of auxiliary subunits.Indeed, in some cases the role of these subunits in channelopathies has already been appreciated (O'Malley and Isom, 2015).In the long run, being able to move between different time and length scales, from molecules to tissue level, will become necessary and indeed there is already progress in this area (Kharche et al., 2012;Miceli et al., 2013).
The development of sophisticated multiscale models is likely to be especially important as the "one mutant, one phenotype" paradigm does not hold true in many cases.For example, large exome sequencing (Klassen et al., 2011) revealed that rare missense variation was as complex in patients associated with sporadic idiopathic epilepsy as unaffected patients.The authors conclude that computational modelling of genetic variation in cell models will become even more important in future years in order to assess risk to excitability disorders.A further study appears to support this where no single rare variant can be linked to idiopathic generalized epilepsy (Heinzen et al., 2012).

M A N U S C R I P T A C C E P T E D
) homology modelling and ii) molecular dynamics (MD) simulations.
molecules and ions and can thus provide very detailed atomicresolution insight into the problem.

Figure 1 .
Figure 1.Schematic of voltage gated ion channel topologies.(A) Topology of K v and CNG channels.(B) Topology of the K ir channel, which can sometimes associate with an SUR1 protein -mutations in SUR1 are known to influence K ir 6.2 (K ATP ) channel behaviour.(C) General membrane topology of Na v and Ca v channels, which are both expressed as a single protein.The voltage-sensing domain is shown in purple and used an alanine scan with double mutant cycle analysis to identify putative interaction sites KCNE1.The results from these experiments suggest that KCNE1 binds to a region formed by helices S4, S5 and S6 of KCNQ1.These experimental observations were interpreted with a homology model of KCNQ1 combined electrophysiology and different modelling approaches to rationalize two mutations at the R213 position (R213W and mutations were found to destabilize the open state.The homology modelling revealed that R213 likely plays a critical role in stabilizing the activated voltage sensor configuration.Conversely, Yang et al.

Figure 2 .
Figure 2. (A) Snapshot of the Na v Ab channel structure(Payandeh et al., 2011) embedded in a lipid bilayer with a solution containing sodium and chloride ions in the and mapped 27 spontaneous mutations known to cause congenital myasthenia and autosomal dominant nocturnal frontal lobe epilepsy onto the structure.Most of the mutations were located at the interfaces between subunits, fate suffer from lower activity and membrane stability (for a review see(Lukacs and Verkman, 2012)).

Figure 3 .
Figure 3. CFTR structure and disease-causing mutations.(A) Topology of the CFTR mutations, such as V510D and the Teem (G550E, R553Q, and R555K) mutations, reduce fluctuations in the NBD1, thereby dampening the flexibility of this domain.Insertion of proline residues in the key dynamic regions can also lead to a decreased flexibility of the NBD1